【求助】看看我编的二维FDTD程序

lanzhiyu520

以下是我编写的二维TM波的FDTD 的FORTRAN程序
是自由空间中的电磁波的传播模拟，没有其他的限制条件，只是想实现模拟电磁波在二维真空中的传播情况，但这个程序的结果却是不稳定的，就是电场随时间无限增大，我看了下这个程序是符合时间稳定条件的，现在怎么也找不到毛病出在哪里了，请高手帮忙看看问题出在哪里：

program fdtd2D
real,dimension(:,:),allocatable::hx,hy,ez
real eps,mu,dt,dx,dy,c0,pi,f
pi=3.14159265
c0=3.0e8 !光速
f=3.e12  !激励源频率
dx=0.05  
dy=dx   
dt=dx/(2*c0)
allocate (hx(1:nn,1:nn))
allocate (hy(1:nn,1:nn))
allocate (ez(1:nn,1:nn))
hx=0.
hy=0.
ez=0.
eps=8.854e-12 !真空
mu=4.0e-7       !真空
nt=100      !时间步
nn=100     !网格数
  open(1,file='D:\MSDEV\ex.txt')
  !open(2,file='D:\MSDEV\cell.txt')
do n=1,nt  !外围循环是时间
   do j=2,nn
     do i=2,nn
       if (i.eq.2.and.j.eq.2) then
           ez(i,j)=sin(2*pi*f*n*dt) !激励源 !
       else
           ez(i,j)=ez(i,j)+(dt/eps)*((hy(i,j)-hy(i-1,j))/dx-(hx(i,j)-hx(i,j-1))/dy)
       end if
     end do
   end do

   do j=1,nn-1
     do i=1,nn-1
           hx(i,j)=hx(i,j)-(dt/mu)*(ez(i,j+1)-ez(i,j))/dy
           hy(i,j)=hy(i,j)+(dt/mu)*(ez(i+1,j)-ez(i,j))/dx
     end do
   end do
end do
  write(1,'(f20.10)') (ez(i,15),i=1,nn) !输出电场
  !write(2,'(i5)') (i,i=1,nn)
end

[ 本帖最后由 lanzhiyu520 于 2007-10-18 20:47 编辑 ]

lanzhiyu520

我把数组的初值重新定义了一下，就是赋初值0
do i=1,nn
  do j=1,nn
    hx(i,j)=0.
  end do
end do
do i=1,nn
  do j=1,nn
     hy(i,j)=0.
  end do
end do         
do i=1,nn
  do j=1,nn
     ez(i,j)=0.
  end do
end do
替换掉下面的代码：
hx=0.
hy=0.
ez=0.
结果还是会增大，只是增大的速度慢了

messenger

hx(i,j)=hx(i,j)-(dt/mu)*(ez(i,j+1)-ez(i,j))/dy
好象应为
hx(i,j)=hx(i,j)+(dt/mu)*(ez(i,j+1)-ez(i,j))/dy
-------------------------------------------------------
do j=1,nn-1
     do i=1,nn-1
          hx(i,j)=hx(i,j)-(dt/mu)*(ez(i,j+1)-ez(i,j))/dy
           hy(i,j)=hy(i,j)+(dt/mu)*(ez(i+1,j)-ez(i,j))/dx
     end do
   end do
end do

lanzhiyu520

回3楼的好心人：书上的是“-”号，不过我按你的建议试了下，结果还是会增大
公式应该没有错呀

lanzhiyu520

现在我把源改在中心了，网格也改成300nm了，并且把电场和磁场放在纵坐标内部循环了，结果还是变大，就是不收敛！！
        program fdtd2D
real*8,dimension(:,:),allocatable::hx,hy,ez
real eps0,mu0,dt,dx,dy,c0,pi,f
integer pmin,pmax
pi=3.14159265
c0=3.0e8  !光速
f=3.e12    !源频率
dx=3.e-7
dy=dx         
dt=dx/(2*c0)        
eps0=8.854e-12
mu0=4.0e-7
nn=100   !网格数
pmin=0
pmax=nn
  open(1,file='D:\MSDEV\ex.txt')
allocate (hx(pmin:pmax,pmin:pmax))
allocate (hy(pmin:pmax,pmin:pmax))
allocate (ez(pmin:pmax,pmin:pmax))
do i=pmin,pmax
  do j=pmin,pmax
    hx(i,j)=0.
  end do
end do
do i=pmin,pmax
  do j=pmin,pmax
     hy(i,j)=0.
  end do
end do         
do i=pmin,pmax
  do j=pmin,pmax
     ez(i,j)=0.
  end do
end do
nt=20 !时间步数
do n=1,nt
do j=pmin+1,pmax
   if(j>=pmin+1) then
     do i=pmin+1,pmax
           if (i.eq.50.and.j.eq.50) then
                 ez(i,j)=sin(2*pi*f*n*dt) !激励源        !
           else
                 ez(i,j)=ez(i,j)+(dt/eps0)*((hy(i,j)-hy(i-1,j))/dx-(hx(i,j)-hx(i,j-1))/dy)
           end if
     end do
   end if
   if(j<=pmax-1) then
      do i=pmin,pmax-1
              hx(i,j)=hx(i,j)-(dt/mu0)*(ez(i,j+1)-ez(i,j))/dy
              hy(i,j)=hy(i,j)+(dt/mu0)*(ez(i+1,j)-ez(i,j))/dx         
      end do
   end if
end do
end do
  write(1,'(f20.10)') (hx(i,52),i=pmin,pmax)
end

[ 本帖最后由 lanzhiyu520 于 2007-10-24 09:37 编辑 ]

sannxi

你不用花时间编二维程序，二维程序葛老师的那本书（西安电子科技大学出版社）后面已经编好了，你可以直接用了，

sannxi@163.com