∃ Λ ∃ 2 2 2 第 17 卷　第 4 期 2006 年 12 月 　JOU RNAL O F GUAN GX IUN IV ER S IT Y O F T ECHNOLO GY　　　 　　　　　　　　　　广 西 工 学 院 学 报 　　　　　　　　　　　V o l 17　N o 2006 D ec 4 文章编号　1004 6410 (2006) 04 0043 04 基于M ATLAB 的 FD TD 算法编程 赵　嘉 (广西工学院 计算机工程系, 广西 柳州　545006) 摘　要: 介绍了时域有限差分 (FD TD ) 法的基本原理, 推导了二维 TM 模 Yee 算法的 FD FD 表达式, 并结合算例阐述 了基于M A TLAB 编程的基本方法。 关　键　词:M A TLAB; FD TD; Yee 算法; 编程 中图分类号: T P321　　　　文献标识码: A 0　引言 S 时域有限差分 (FD TD ) 法是六十年代由 K Yee 提出并首先用于求解电磁散射问题, 其主要思路是在 三维空间和时间轴上对场量离散, 并且用中心差分代替偏微分, 将麦克斯韦方程组转化为差分方程, 通过在 时间轴和空间轴上采取蛙跳法 (leapfrog) 逐步推进地求解, 最终求出一定边值与初值条件下的空间场解。 随 着计算机技术的发展, 近年来 FD TD 计算技术也得到了越来越多的应用。对于 FD TD 算法的编程求解, 最常 用的有V C 和 FOR TRU N , 而M A TLAB 作为一种可视化效果很好的科学计算软件, 在 FD TD 计算中能充分 发挥编程简单、可视化程度高、能显示动态场 效果的特点。 本文首先推导二维 TM 模 Yee 算 法 的 FD FD 表 达 式, 并 讨 论 如 何 结 合 FD TD 算法边界条件的特点, 用 M A TLAB 语言进行编程的方法和应注意的问题。 1　二维 TM 模 Y ee 算法 1 1　算式推导 在自由空间中, 对于二维问题, = 0, 对 于 TM 波, H z = E x = E y = 0,M A XW ELL 的 两个旋度方程可分解为 1 : E z y H x t = - , 　　 z 0 图 1　二维 TM 波Y EE cell E z x = H y t 0 , 　　 0 E z t = H y x - H x y 　　构造二维 TM 波 YEE cell 如图 1 所示: 按 YEE 元胞对上式偏导用中心差商代替, 可得: 　　　　　 E z ( i, j + 1 2 , n) - E z ( i, j - 1 2 , n) y = - 0 H x ( i, j + 1 2 , n+ 1 2 ) - H x ( i, j + 1 2 , n- 1 2 ) t (1) 收稿日期: 2006 作者简介: 赵　嘉 (1965 08 20 ) , 女, 广西柳州市人, 广西工学院计算机工程系讲师。 

44 广西工学院学报　　　　　　　　　　　　　　　　　　　　第 17 卷　 　　　　　 E z ( i+ 1 2 , j , n) - E z ( i- x 1 2 , j , n) H y ( i+ = 0 1 2 , j , n+ 1 2 ) - H y ( i+ 1 2 , j , n- 1 2 ) t 　　　 0 E z ( i, j , n+ 1 2 ) - E z ( i, j , n- 1 2 ) t 1 2 1 2 ) H y ( i+ 1, j , n+ ) - H y ( i, j , n+ H x ( i, j + 1, n+ - H y ( i, j , n+ 　　　　= 　　方程中出现了半个网格和半个时间步, 为便于编程, 可将上面差分式改写为如下 FD TD 算式: 　　　　　　　　　H x (i, j , k+ 1) = H x ( i, j , k) - [E z (i, j + 1, n) - E z ( i, j , n) ] y x - t 1 2 1 2 ) (2) (3) (4) (5) (6) 0 y t 0 x 　　　　　　　　　H y (i, j , n+ 1) = H y ( i, j , n) + [ E z ( i+ 1, j , n) - E z ( i, j , n) ] 　　　E z ( i, j , n+ 1) = E z (i, j , n) + 1 2　数值稳定性条件 t 0 H y ( i+ 1, j , n) - H y ( i, j , n) x - H x ( i, j + 1, n) - H y ( i, j , n) y y、 x、 t 与空间步长 FD TD 方法是以一组有限差分方程来代替M A XW ELL 旋度方程进行数值求解的方法, 在执行形如上 式 FD TD 算法时, 随着时间步的增长, 保证算法的稳定性是一个很重要的问题。数值解是否稳定主要取决于 时间步长 z 间的关系。按照 Cournan t 稳定性条件, FD TD 算法中空间和时间间隔之 间应满足的关系为: 1 　　　　　　　　　　　　　　　　　C y 2+ 其中 C 为真空中的光速。 对于 TM 模的二维情形, 令式中 C 是可变时, 即 应满足下式: (7) y , 此时, y 不等时, 则应取二者中的最小值。当沿两个轴向的网格元 y 分别是 i, j 的函数, 则应取每个轴向上的最小值, 再选二者中之最小者。总之, 时空步长 1 z 2 z = ∞, 网格单元尺寸通常可取 2 , 一般可选 C t≤ 1 2。当 x 2+ t≤ x = x , x , t= - 1 2 x x 　　　　　　　　　　　　　　　　　　　 t= 才能确保算法在较长时间步长上运行的稳定性。 m in ( x m in· y m in) 2C (8) 2　基于M ATLAB 的编程方法 1　算例 2 　　考察区域为正方形的二维自由空间, 区域的边 界为理想电导体, 为了激励外向辐射波, 设在区域的 中心有电场分量 E z 随时间按高斯变化, 求解该区 域内电场随时间的变化, 场域如图 2 所示。 2 2　M ATLAB 编程方法及结果 图 2　算例计算区域 基于M A TLAB 对本算例 FD TD 算法进行编程 求 解可归纳为以下几个步骤: (1) 正确建立二维 TM 模 Yee 算法的数学模型。 这个工作包括正确地 划分网格、推导出正确的 FD TD 算式等。 本算例的算式已在上节给出, 下面介绍网格的划分: 根据算例的尺 寸, 可将区域划分为 200×200 个网格, 其中空间步长取 2C , 其中, C 为 自由空间的光速。 t 的选取必须使它们的关系满足 Cournan t 稳定性条件。尽管时间和空间步长越 小, FD TD 运算的精度会越高, 但太细的网格会使计算机运算时间太长, 同时, 在网格分辨率 (即波长与空间 步长之比) 小到一定程度时, 数值误差基本就不变了。对于采用M A TLAB 编程, 需要特别注意时间步长和空 间步长的比值, 如果比值选取不当, 有可能在伪彩图中没有激励的波形, 一般可选取该比值小于 10 倍光速。 005, 时间步长取 y 和 y = 0 x、 x = t= x 

　第 4 期　　　　　　　　　　　　　　赵　嘉: 基于M A TLAB 的 FD TD 算法编程 54 　　 (2) 确定边界条件和设置激励源。 FD TD 方法可求解电磁学中的散 射及辐射等问题, 为了在有限的空间内模拟无限大空间的问题, 必须在计 算时考虑吸收边界条件。 常用的吸收边界M ur 和 PM L 吸收边界。 一般 情况下, 对角区要求不高的场合, 可用M ur 边界条件, 而对计算精度要求 高且计算机性能较好时, 可采用 PM L 吸收条件。边界条件的 FD TD 算式 可参阅相关的资料, 在编程中可将边界条件的 FD TD 算式直接写在程序 里。 对于本算例, 求解的是二维硬边界条件, 可在程序中直接在四个边界 上将电场值赋 0, 同时将 E x、E y 和 H z 赋 0 (TM 波) 即可。此外, 激励源的 设置在 FD TD 算法中至关重大, 常用的激励源主要有正弦激励和高斯激 励, 本算例中采用微分高斯脉冲, 其表达式为: t0) 2 2 　　　　　　　E i ( t) = ( t- (9) 　　微分高斯脉冲的特点是不含零频率分量, 为使M A TLAB 显示的二 维伪彩色效果明显, 需选择合适的 t0 和 值, 本算例中, 经实验取 t0= 50 nS, t0) exp (- = 5 S. (3) 按要求绘制M A TLAB 程序流程图, 编写程序。M A TLAB 具有 强大的可视化功能, 在本算例中, 可以通过编程的方式使运算结果以二维 或三维伪彩图的形式输出, 程序流图如图 3 所示。 2 3　算例运行结果 ( t- ) 图 3　程序流图 图 4 是运算时间步分别为N = 137、N = 159 和 N = 174 时场 E Z 的伪彩色显示图。 其中, N = 137 步时, 波未到达边界, N = 159 步时, 电磁波刚好到达硬边界, N = 174 步时, 显示电磁波在硬边界的反射情况。 显 然,M A TLAB 运行能很好地显示 E Z 随时间步变化的动态效果图, 能很好地反映电磁波在区域内传播和在 边界反射的情况。 在 N = 120 步时, 给出的是M A TLAB 显示的 E Z 三维效果图。 

2 64 3　结束语 广西工学院学报　　　　　　　　　　　　　　　　　　　　第 17 卷　 本文结合算例介绍了使用M A TLAB 对 FD TD 算法编程求解电磁场问题的思路和方法。 FD TD 在电磁 场数值分析中, 具有概念清晰、不用对大型矩阵求逆, 因而编程相对简单的特点, 而M A TLAB 强大的数据处 理功能和图形处理功能, 二者结合起来, 能使 FD TD 算法的编程更加简明扼要, 编程效率更高, 结果显示更 直观。 参　考　文　献: 电磁波时域有限差分法[M ] 西安: 西安电子科技大学出版社, 2002 时域有限差分法的M A TLAB 仿真[J ] 现代电子技术, 2005, 8 (199) : 45 46 1 葛德彪, 闫玉波 2 郭春波 3 A 4 石澄贤, 郑明芳 T aflove Com putational E lectrodynam ics D ifference T im e M A TLAB 语言及其在电子信息工程中的应用[M ] T he F in ite Dom ain M ethod[M . A rtech Hourse, 1995 北京: 清华大学出版社, 2004 FD TD ar ithmetic programm ing ba sed on M ATLAB language ZHAO J ia (Com puter Engineering D epartm en t, Guangx iU n iversity of T echno logy, L iuzhou 545006, Ch ina) Abstract: T he basic p rincip le of FD TD w as in troduced and the 2 derived. T he p rogramm ing m ethod based on M A TLAB w as illustrated w ith exam p le. Key words: M A TLAB; FD TD; Yee A rithm etic; P rogramm ing D TM m ode FD TD arithm etic exp ression w as (责任编辑　赖君荣) (上接第 38 页) 参　考　文　献: 1 张立明. 人工神经网络的模型及其应用[M . 上海: 复旦大学出版社, 2001. 163 2 胡守仁. 神经网络导论[M . 湖南: 国防科技大学出版社, 1993. 3 3 丛　爽. 面向M A TLAB 工具箱的神经网络理论与应用[M . 合肥: 中国科学技术大学出版社, 1998. 11. 4 广西统计局. 广西统计年鉴[M . 北京: 中国统计出版社, 2005. 5 陈桂明, 戚红雨, 潘　伟. M atlab 数理统计 (6. X) [M . 北京: 科学出版社, 2002. 3. 164. 10. The appl ica tion of neura l network in foreca sting the eff ic iency of util iz ing energy FAN Song hai, L ID E w ei (D epartm en t of F inance and Econom ics, Guangx i U n iversity and T echno logy, L iuzhou 545006, Ch ina) Abstract: In view of the non linear characteristic that the develop ing trend of the efficiency of utilizing energy in Guangx i has p resen ted, it w as sim ulated and fo recasted based on BP neural netw o rk. T he resultsw ere con trast ed w ith the trend of the efficiency of utilizing energy in our coun try. A nd the p ractical sign ificances of the fo re casted o r sim ulated results w ere analyzed. Key words: neural netw o rk; efficiency of utilizing energy; fo recast (责任编辑　李　捷)