摘要：在用一维时域有限差分方法计算波的传播时，为了防止出现发散的结果，要求时间步长Δt<■，但具体取多大的值，没有定论。通过编程试验，时间步长取Δt<0.75■，而时间步数Nt大于波从波源传到观察点所需要的时间，又小于反射波到达观察点的时间。而且计算的时间格点和位置格点全部都可以取整数。
　　关键词：FDTD；有限时域差分方法；参数
　　中图分类号：G642.4？摇 文献标志码：A？摇 文章编号：1674-9324（2013）03-0155-02
　　时域有限差分方法（FDTD）1966年由Yee提出后，经过许多年的发展，在电磁波散射方面得到了广泛的应用。现在，在光子晶体能隙的计算方面也得到了很多应用。笔者在一维、二维的FDTD编程计算中，发现时间步长的选择不同，会导致结果的不同，下面进行讨论。
　　一、一维FDTD差分方程表达形式
　　一维波动方程■+■■=0，在一维无限长介质中，其解为u=f（x-vt），是以速度v向x轴正向传播的行波。用差分形式，把波动方程进行改写：
　　由u（x+Δx）=u（x）+■Δx+■■（Δx）2，u（x-Δx）=u（x）-■Δx+■■（Δx）2。两式相减，u（x+Δx）-u（x-Δx）=2■Δx，于是，■=■。同理：■=■。于是，一维波动方程■+■■=0写为：■+■■=0，整理得到：un+1（i）=un-1（i）-■[un（i+1）-un（i-1）] ①
　　二、初始、边界条件的确定以及编程计算的结果
　　边界条件：u（t，1）=0，u（t，Nx）=0，对t=1到Nt。
　　初始条件：由于在差分形式的递推公式①中，要计算t=3的值，需要知道t=1，2时的值。本文规定初始条件为u（t，i）=0，对t=1，2.
　　下面是两组不同的参数，用Matlab编程计算的结果：
　　第1组：Lx=1×10-4，Δx=1×10-6，Nx=100，Nt=100，v=c=3×108，Δt=■
　　信号源u（t，10）=5sin（■t），观察点在i=80，也就是x=80Δx处。
　　第2组：Lx=1×10-4，Δx=1×10-6，Nx=100，Nt=1000，v=c=3×108，Δt=■
　　信号源u（t，10）=5sin（■t），观察点在i=80，也就是x=80Δx处。参数中时间长度Nt与第1组不同，其他参数相同。
　　第3组：Lx=1×10-4，Δx=1×10-7，Nx=1000，Nt=3000，v=c=3×108，Δt=■■
　　信号源u（t，10）=5sin（■t），观察点在i=500，也就是x=500Δx处。
　　从图中可以看出：
　　t=1到600，在观察点无信号，信号还没有到达；
　　t=600到900，是过渡阶段，与数值化的计算有关，在观察点信号尚不稳定；
　　t=900到2400，观察点信号与用连续方程得到的结果一致；
　　t=2400以后，观察点的信号的幅度比源的信号还大，是因为文中取一维两端边界为零，反射波到达观察点，与正向波形成干涉。
　　用快速傅里叶变换得出的频谱透射率，由于已经包含了开始的过渡阶段和以后的反射波的干涉的影响，在0.8左右，与波的无衰减传播结果接近。
　　三、结论
　　在一维的FDTD计算中，直接使用整数时刻和整数空间坐标，而不需要使用半整数时间和半整数空间坐标，简化了坐标的标记形式；时间步长取Δt=■■，时间总长Nt的取值，要求在这个时间信号源的信号可以到达观察点，但又没有界面反射波到达观察点，这样的参数比较理想，可以得到与连续方程相同的结果。
　　参考文献：
　　[1]葛德彪，阎玉波.电磁波时域有限差分方法[M].第2版.西安：西安电子科技大学出版社，2005.
　　[2]曾辉，杨亚培.FDTD法与平面波展开法在光子禁带计算中的差异分析[J].电子科技大学学报，2005，34（6）：901-904.
　　[3]张国华，袁乃昌，付云起.FDTD方法分析光子带隙能带结构[J].微波学报，2001，17（4）：14-17.
　　[4]张亮，寇晓艳.FDTD的Matlab语言的实现[J].延安大学学报：自然科学版，2009，28（2）：57-59.
　　[5]朱章虎，卢万铮，冯奎胜.FDTD计算中PML的简化应用及编程实现[J].空军工程大学学报（自然科学版），2006，7（2）：55-57.
　　作者简介：陈义万，湖北工业大学理学院物理学副教授，研究方向：光子晶体及性质。