Solving Large/Infinite Results w/Yee 1D FDTD Code

[angelos_physik]

Hello to all!
I have written this 1D FDTD code but the results I take are very large and after some timestep they get to infinite/NaN... I really don't know what's wrong. The expressions I have used are the exact given solution for the Maxwell equations.

Thank you for your help.

program yee
    implicit none    real :: Delta, tau,omega,n,mu,time,Eprev,Hprev
    real, dimension(:), allocatable :: sigma,eps,A,B,C,D,J
    real, dimension(:), allocatable :: E,H
    real, parameter :: pi=3.1415926
    integer :: i,m=10000, Length=5000,ls,k    Allocate(E(Length),H(Length),sigma(Length),eps(Length))
    Allocate(A(Length),B(Length),C(Length),D(Length))
    Allocate(J(Length))

!Initialization
    open(unit=1,file="Elektrfeld.txt")
    Delta=0.02d0
    tau=0.9d0*Delta
    omega=2.0d0*pi
    n=1.46d0
    mu=1
    ls=1000
    E(:)=0.0d0
    H(:)=0.0d0
    A(:)=0.0d0
    B(:)=0.0d0
    C(:)=0.0d0
    D(:)=0.0d0
    J(:)=0.0d0
    Eprev=0.0d0
    Hprev=0.0d0    do i=1,Length
        if (i.le.6/Delta) then
            sigma(i)=1.0d0
        else if (i.gt.6/Delta .and. i.lt.((Length*Delta-6.0d0)/Delta))then
            sigma(i)=0.0d0
        else if (i.ge.(Length*Delta-6.0d0)/Delta .and. i.le.Length) then
            sigma(i)=1.0d0
        end if
    end do    do i=1,Length
        if (i.lt.Length*0.5d0) then
            eps(i)=1.
        else if (i.ge.1100 .and. i.lt.(1100+2)) then
            eps(i)=n**2.
        else if (i.ge.(1100+2)) then
            eps(i)=1.
        end if
    end do    do i=1,Length
        A(i)=(1.-sigma(i)*tau/(2.*mu))/(1.+sigma(i)*tau/(2.*mu))
        B(i)=(tau/mu)/(1.+(sigma(i)*tau)/(2.*mu))
        C(i)=(1.-sigma(i)*tau/(2.*eps(i)))/(1.+sigma(i)*tau/(2.*eps(i)))
        D(i)=(tau/eps(i))/(1.+sigma(i)*tau/(2.*eps(i)))
    end do    write(*,*) "time"
    read(*,*) time

!----------------Main calculation----------------
    do k=1,m        !time
        E(1)=0.0d0
        E(Length)=0.0d0
        J(ls)=sin(2.*pi*k*tau)*exp(-((k*tau-30.)/10.)**2.)
        do i=2,Length
            Eprev=E(i)
            E(i)=D(i)*(H(i)-H(i-1))/Delta+C(i)*E(i)-D(i)*J(i)
        end do
        do i=1,Length
            Hprev=H(i)
            H(i)=B(i)*(E(i)-E(i-1))/Delta+A(i)*Hprev
        end do
        if (k.eq.time ) then
            write(1,*) "-----",k
            do i=800,1200
                write(1,*) E(i)
            end do
        end if
    end do

    close(1)
end program yee

 

[phinds]

we have "code" tags that get your code properly formatted. It would be a good idea if you modified your post to use them.

 

[gsal]

Well, at quick glance, the second i-loop goes from 1 to Length and it contains the term E(i-1), this means that for i=1, you are calling on E(0) which is not defined...it is outside your program altogether given the way you allocated the arrays, their indexes start at 1.