（计算数学专业论文）辛紧致格式fdtd方法的研究.pdf

辛紧致格式f d t d 方法的研究中文摘要 辛紧致格式f d t d 方法的研究 中文摘要 传统的时域有限差分法 f i n i t ed i f f e r e n c et i m e d o m a i n f d t d 即y e e 格式 在 时间空间导数离散上都采用二阶中心差分格式 格式精度较低 色散耗散误差较大 对电大问题作电磁波传播长期响应分析时 由于误差的积累 往往造成波形的严重失 真 这是传统的f d t d 方法的固有缺陷 针对这些问题 本文采用时间辛算法与空间紧致格式相结合的方法 构造出时间 空间均达到四阶精度的f d t d 格式 新格式理论上不产生耗散误差 相比其他一些典 型高精度格式有更低的色散误差 和更好的稳定性 在具体格式分析过程中 讨论了 紧致格式的色散误差及其各向异性特性 介绍了h a m i l t o n 系统及其辛性质 辛格式 的主要构造方法 介绍了辛算法在m a x w e l l 方程中的应用 和辛f d t d 格式的构造 给出了辛紧致f d t d 格式结合的行为分析 包括稳定性 色散分析 以及同 误差限 制下的存储量与计算量估计 数值实例中证实了本文提出的辛紧致f d t d 格式对长期 响应分析的有效性 另外还将辛格式应用到电磁目标的散射计算中 得n j 与解析解 相吻合的结果 关键词 时域有限差分法 辛算法 紧致格式 色散分析 稳定性 作者 傅平 指导老师 凌国平 教授 f h es t u d yo nt h es y m p l e c t i cc o m p a c tf d t dm e t h o d a b s t r a o t t h es t u d yo nt h es y m p l e c t i cc o m p a c tf d t dm e t h o d a b s t r a c t t h es t a n d a r df i n i t ed i f f e r e n c et i m e d o m a i nm e t h o do fm a x w e i le q u a t i o n i e t h ey e es c h e m e w h ic hu s e st h es e c o n d o r d e rc e n t e rd if f e r e n c es c h e m ein b o t hl i m ea n ds p a c e i sal o w o r d e re x p l i c i ts c h e m e w i t hg r e a t e rd i s p e r s i o n a n dd i s s i p a t i o ne r r o r s f o re l e c t r i c a l i yl a r g ed o m a i n sa n dl a t e t i m ea n a l y s i s i tw o u l dr e s u l ti nt h ew a v e f o r md i s t o r t i o n d u et ot h ea c c u m u l a t i o no fn u m e r i c a l e r r o r s i ti st h ei n t r i n s i c1 i m i t a t i o nf o ry e es c h e m e i nt h i sp a p e r at e m p o r a ls y m p l e c t i ci n t e g r a t o rp r o p a g a t o ra n das p a t i a l c o m p a c td i f f e r e n c es c h e m ea r ec o m b i n e di n t oan e wm e t h o d w h i c hi s f o u r t h o r d e ra c c u r a t ei nb o t hs p a c ea n dt i m e w i t hn od i s s i p a t i o ne r r o r l o w e r d i s p e r s i v ee r r o ra n dl a r g e rs t a b i l i t yc o m p a r e dw i t ho t h e rc l a s s i ch i g h o r d e r d i f f e r e n c es c h e m e s f i r s t d i s p e r s i v ee r r o ro fc o m p a c td i f f e r e n c es c h e m ei s d i s c u s s e d t h e n w ei n t r o d u c et h ea p p l i c a t i o no fs y m p l e c t i cm e t h o di n t o m a x w e i ie q u a t i o n a n dg i v et h ec o n s t r u c t i o no fs y m p l e c t i cf d t ds c h e m e l a t e r w es h o wt h eb e h a v i o r a la n a l y s i sf o r t h es y m p l e c t i cc o m p a c tf d t dm e t h o d a t l a s t n u m e r i c a le x a m p l e sd e m o n s t r a t et h es u p e r i o rp e r f o r m a n c ei nl a t e t i m e a n a l y s i so ft h em e t h o d a n dt h ee l e c t r o m a g n e t i cs c a t t i n gf i e l dc a l c u l a t i o n i sa g r e e dw e l l w i t ht h ea n a l y s i sr e s u l t i i w r i t t e nb yf u p i n g s u p e r v i s e db yp r o f l i n g g u o p i n g y 7 8 1 2 8 6 苏州大学学位论文独创性声明及使用授权声明 学位论文独创性声明 本人郑重声明 所提交的学位论文是本人在导师的指导下 独立进行研究工作所 取得的成果 除文中已经注明引用的内容外 本论文不含其他个人或集体已经发表或 撰写过的研究成果 也不含为获得苏州大学或其它教育机构的学位证书而使用过的材 料 对本文的研究作出重要贡献的个人和集体 均已在文中以明确方式标明 本人承 担本声明的法律责任 研究生签名 盘聋 日期 塑苎 竺 学位论文使用授权声明 苏州大学 中国科学技术信息研究所 国家图书馆 清华大学论文合作部 中国 社科院文献信息情报中心有权保留本人所送交学位论文的复印件和电子文档 可以采 用影印 缩印或其他复制手段保存论文 本人电子文档的内容和纸质论文的内容相一 致 除在保密期内的保密论文外 允许论文被查阅和借阅 可以公布 包括刊登 论 文的全部或部分内容 论文的公布 包括刊登 授权苏州大学学位办办理 研究生签名 导师签名 日期 玉哩 尘 日期 五 竺 辛爨数格式f d t d 方法的研究 第一章引言 第一章引言 1 1 传统y e e 格式的缺点及改进 传统的时域有限差分法 f d t d 即y e e 格式 1 1 1 2 1 作为一种典型的全波分析时域 算法 在电磁问题的众多领域中得到广泛的应用 展现了广阔的发展前景 然而 传 统的f d t d 方法也有其自身的不足之处 至少表现在以下几个方面 一 格式在空间 时间导数离散上都采用二阶中心差分 精度较低 二 格式的色散误差较大 同时由于时间项是以二阶精度显式推进 受到c o u r a n t 稳定性条件的限制较大 也就是时间步长的取值由空间步长的最小值决定 时间步长 的减小 使得迭代步数和运算时间大量增加 三 长期响应分析结果较差 对电大问题作电磁波传播长期响应分析碍寸 髓着时 间计算步的增加 格式离散所产生的色散和耗散误差会慢慢积累 最终导致计算波形 的严重失真 为提高精度 降低色散误差 一般采用高阶差分格式来离散m a x w e l l 方程 空 间离散的传统高精度格式有 f a n g 3 1 提出的典型的四阶中心差分格式 然而f a n g 4 t h 格式由于其使用网格点数较多 对吸收边界条件的处理适应性差 6 j 另一种高精度离 散格式是有限紧致差分格式 与传统的高精度差分格式相比 在相同糟度下 紧致差 分格式有着网格基架点少和分辨率高的优点 4 并使得格式在边界处理要比传统的高 精度差分格式更加容易 在时间推进上 传统的高精度格式有f a n 9 1 3 1 四阶格式 它利 用m a x w e l l 方程 对二阶中心差分格式的截断误差瑷透过空闯的高阶导数来近似 使 时间推进精度达到四阶 另外还有典型的四阶r u n g e k u t t a v 以t 简称r k 方法 7 然 而 f a n g 4 n d 格式对非均匀媒质 需要求解介电系数或磁导系数的二阶导数近似 处理复杂 并由于其需要求解高阶的空间导数导致对边界问题和表面处理的计算复杂 1 6 1 r k 方法不仅产生色散误差 更有耗散误差的引入 相对其他高精度时间格式 还需要更多的内存空间来存储中间变量 而且上述两种格式对长期响应都存在问题 计算后期容易导致波形的失真 近年来由基于h a m i l t o n 系统提出的辛算法在m a x w e l l 兰茎堡垡丛 堡 里查鳖璺磷究 第一章引言 时域差分法中得到应用 该方法主要特征是保持系统的相空间体积不变和总能量守恒 j 这使得它能克服上述长期响应的波形失真问题 并且显式的辛格式有着和y e e 格 式一样的简单计算形式 仅增加了计算步数 对边界处理容易 相比r k 方法而言 辛格式则不需要存储中间变量 降低了存储量 在本文中 将采用时州辛算法与空间紧致格式相结合的方法 并应用于m a x w e l l 方程的求解 提出一种时间 空间均达到四阶精度的辛紧致f d t d 方法 新格式理论 上不产生耗散误差 相比其他典型高精度格式有更小的色散误差 和更好的稳定性 而且保持了能量的守恒 能克服长期响应波形失真问题 这方匠的工作 尤其是辛算法在m a x w e l l 方程时域有限差分法中的应用 目前在 国外已有少量文献出现 国内在这方面的研究还刚刚起步 2 0 0 4 年才见有这方面的 文章 该课题的研究是比较新的 并且有很好的实际应用价值 1 2 本文内容安排 第二章介绍了传统时域有限差分法y e e 格式的基本知识 分析了它的稳定性和色 散误差 第三章介绍了空间导数离散的紧致格式 给出格式的色散分析和各向异性特性 并与其他一些常用格式进行比较 第四章介绍了h a m i l t o n 系统及其辛算法的一般知识 重点介绍了一种显式辛算 法的构造方法 p r k 法 并给出简单算例说明辛算法的保能性 第五章把h a m i l t o n 系统辛算法引入到m a x w e l l 方程中来 给出无耗和有耗媒质 情况的四阶p r k 计算格式 并进行了辛格式的稳定性及色散分析 第六章对辛紧致时域有限差分格式的联合行为 包括色散性 稳定性等进行了分 析 并与其他 些高精度格式进行了比较 还给出了同一误差限制下不同格式的计算 量与存储量估计值的比较 第七章给出数值算例 计算了一维的高斯脉冲波的长期传播波形 证明了辛紧致 格式的有效性 再将辛紧致格式应用到电磁目标的散射计算中 给出二维金属圆柱散 射电场幅值波形 说明辛算法的适用性 2 辛蜷致格式f d t d 方法的聊f 究第二章时域有限差分法 f d t d 第二章时域有限差分法 f d t d 2 1 传统时域有限差分法 y e e 格式 描述电磁场的m a x w e l l 方程为 v 雷 丝 了 v 屉一票一无 口 u v 雪 p v 西 p 各向同性线性媒质中的本构关系为 d e e 一 b 一 l 冠 j g e 一 j m a m 叠 其中豆 疗 西 豆分别是电场强度 磁场强度 电位移矢量和磁感应强度 1 7 1 7 分别是 电流密度和磁流密度 a p 分别是电荷密度和磁荷密度 它们分别是空间位置和时 间的函数 弘占分别是媒质的磁导系数和介电系数 d 盯分别是导磁率和电导率 对于各向同性媒质 它们是标量 而对于各向异性媒质 它们是张量 对于均匀媒质 它们是常量 而对于非均匀媒质 则它们随空间位置而变化 在直角坐标系下 可以将m a x w e l l 方程写为六个标量方程 等一警 鸣 占拿 1 t 5 仃e 二十占 一 种 d z 警一警 0 1 哆5 s 鲁 b z 曲 土一 三 s ll 厶z a l o z蕊 o z 警一等 鸣 占睾u t似卯 3 辛紧致格式f d t d 方法的研究 第 章时域有限差分法 f d t d 拿o z 一熹 吒以 掣c t口v 等一警 髟 孥o t b z n 2 一 爿y 2 2 b 出出 熹一警 1 t 1 牟孕o t 亍一 2 牟 鲫 y e e 1 1 首先将空间按立方体分割 电磁场的六个分量在空间的取样点分别放在立 方体的棱边和表面中心点上 如图2 1 所示 图2 1 基本y e e 阏格以及其中的电磁场分量配置 从图中可以看出 电场分量位于网格棱边中心并且平行于棱边 每个电场分量环 绕着四个磁场分量 磁场分量位于网格表面中心并且垂直于表面 每个磁场分量环绕 着四个电场分量 在空间取样上 电场和磁场分量在任何方向上始终相差半个网格步 长 这种场量的配置不仅允许旋度方程作中心差分近似 也符合法拉弟 f a r a d a y 感应定律和安培 a m p e r e 环路定律的自然几何结构 在时间取样上 磁场分量与 电场分量相互错开半个时间步 这样 增量的时间导数也容易用中心差分近似 以 2 2 a 中的第 个方程离散为例 1 1 e f 寺 七 c 爿 删 e n t l 十i l 工j c b m 竺 l 二圭 圭 竺二竺 l 二圭 三 竺一兰 圭 垫二竺 i 二皇 a y a z 式中嘶 器器卷 c b 埘 面万2 丽a t 4 芝苫 z 辛紧数格式f d t d 方法的研究第二章时域有限差分法 f d t d 上式括号中m 1 2 j 对 2 2 a 中的冥他方程有类似的离散表达式 2 2 b 中第一个方程离散为 彬 2 f 吉 女 了1 c p m 晖 2 f 吉舢百1 一c q m 霹 1 t 百1 一e i i 1 髟 露 1 一蟛 十 七 r 土一 j 二 自 一 a y 丛 其中凹 m 器糌 c 鼬 丽厩2 a 面 上式括号中聊 f 1 2 k 1 2 对 2 2 b 中的其他方程有类似的离散表示式 由于采用了中 心差分近似 y e e 格式在空间和时间h 均具有二阶精度 2 2y e e 格式的色散性 稳定性 假定媒质是均匀无耗p 0 盯 0 从m a x w e l l 方程可导出 在直角坐标系中 电磁场任意分量均满足下列齐次波动方程 c 导 嘉 导一吉 y 丽 萨 萨一7 列y 0 式中c 是均匀媒质中光的速度 上式相应的有限差分近似为 1 矿 1 f 意 一2 9 t 工七 矿 1 f j k 7 一 y f l k 一2 f j k y f 一1 七 一 i 丁 一一 妒 f j l k 一2 y f k 妒 f 一l k 石广一一一 吵 f j k 1 一2 y 七 f j k 一1 一一一 z r 一 一 在空间网格点 f k 上考虑平面波解可表示为 妒 f j k e x p j k i a x k y j a y t j z 一g o n a t 式中丸 k y 砭分别是波数矢量的三个坐标分量 国为角频率 代入波动方程离散式可 得 兰墨塾登三u 里 里堕堡盟 型签 茎三童堕燮查堡墨坌垩 望 堡 击 2 s i n 2 竽 古s 聍 竽 古s 掰 竽 古咖z 竽 但 3 2 3 给出了波动方程离散后平面波波矢量惫 k x k y 砭 与频率国之间应满足的关 系式 即数值色散关系 平面波在连续均匀无耗媒质中的解析色散关系式为 f 竺 砖 蟛 砖 2 4 可以证明 当a t g x 缈 缸趋于零时 式 2 3 变为 2 4 因此 式 2 4 是数值色散 关系式 2 3 的极限式 两者的差异即为色散误差 当网格尺寸足够小时 这种数值色 散可以减小到所要求的程度 一般要求 应用此格式 即得到 其中向 曰 目 q i 鳊 7 令 则把上面的r k 格式拆分如下 z n l t z f 妒 砖 j 乙 f 勺f t l t 夏d q k 翼o h 测心 黝 警 一豢 鲰吼 1 烈力o 鼻 p f a u f p q j q f q r 口筘 弓 g j l p 1 p f 6 f 厂 只 q i l q 1 q r 蟾 只 q 江l s 为了使变换具有典则变换的特性 要求有如 f 限帝d t l l f 1 2 b a j q 一6 f 屯 0 i j l j 一般情况下 由上面的限制条件构造的辛差分格式都是隐式的 因此在计算过程 主鲨璧壁堑垒羔旦旦l 查鲨堕型苎l 一 釜型童 型 翌i 塑里蔓堡生主蔓些 中不可避免的需要矩阵求逆 这将消耗大量的时间 下丽对可分的h a r n i l t o n 系统来建 立 般的显式差分格式 i i 1 若h a m i l t o n 函数能表示成下列的形式 n z h p q r p v q 此时称h a m i l t o n 系统是一个可分的h a m i l t o n 系统 这是正则方程组 4 3 可写成 f亟 a r g 仇1dt 孰 一豢 尥 l 讲 对这样可分的h 锄j t 彻系统 在应用r k 方法时可以对p 和g 采用不同的推进步 称为分块龙格一库塔法 p a n i t i o n e dr u n g e k u t t a 简称p r k t 1 如下所示两个b u t c h e r 列阵 第 个是对应p 的推进系数 第二个是对应q 的推进系数 a 1 a i 每 蟊 i 临 坠l 堡l 矗 6 l b si 云 葭 可以得n 般的表达式如下 以h 表示时间步长 只 少 a 口f q v i q f q 毛g 0 i p 1 p 6 l 厂 q l l q 1 q 五g 霉 l l 同样要求有如下限制 格式才满足典则性 辛变换 b i 磊q b n i b 1 i 0 i j l s 下面具体介绍如何从上述的一般表达式得到我们所需要的显式的辛差分格式 对 pn q 的推进系数取如下两b u t c h e r 列阵 可以证明系数的设定是满足辛变换的要求 i6 1 s 兰鉴塾整垫 望 望立鳖塑墅塑 c 1 0 00 c lc 2 0 c c 2 i j q l 0 c 1c 2 c j 一1c 3 i i i 蔓些至婪 苎 塑 墨壅兰兰曼堕 o0 00 d l 0 0 d 吒 i 0 0 l 刍生 生 旦 d ld 2 d s lt 现以一个二步二阶精度的辛p r k 为例说明 c 0 00 p 的推进瓢i 羞cc 2g 的推进瓶i 糟6 2 l矗l 按照p r k 方法一般表达式展开 得到 q 1 q 4 5 a e p h b f q 1 4 5 b 幺 q 4 娲g p o 4 5 c b p h b q 1 6 2 q 2 只 h b 2 厂 q 2 4 5 d p 1 p h b q 1 b 2 f 0 2 4 5 e q q h b j g p 1 岛g p d 易 8 2 9 县 4 5 f 把 4 5 b 式中的q l 直接代之以q 并注意到 4 5 d 4 5 e 相同 最 p 并略去q i 改写q 2 为翁 简写上面六式得到下面的四式 月 p h b l f q q q 氓g 日 p e l 十 6 2 鸟 q 1 q i h b 2 9 p 1 这样就构造了一个二步的显式辛格式 选取合适的q 吐值 可以使它达到二阶精度 c 1 c 2 l d l d 2 l c 2 d l 1 2 类似有四步的显格式如下 只 p h c f 碍 最 日 h c 2 f q 1 b h e 3 f q 2 1 9 q 1 矿 h d g e 1 q 2 q l m g 只 q 3 蜴 h d 3 9 p 3 辛紧敛格j cf d t d 方法的研究第四章h a m i l t o n 系统与辛算法 p 1 b h c 4 q 3 g q 3 h d 4 9 p 1 当系数取值满足以下关系式 格式的精度可达到四阶 2 c l c 3 1 d l d 2 l 2 一1 2 彰 1 4 8 0 从迭代格式可以看出 上述迭代可以只保存当前的一步结果 而无需像传统的 r k 方法那样保存所有的中j 目l 临时变量 所以辛格式减少了存储量 同时也提高了计 算精度 4 4 数值算例 为说明辛算法的保能性质 以简单的振动方程 1 7 s i m p l eh a r m o n i co s c i l l a t o r s h o 为例 f 粤 口 西 罄一p 它的一个准确解为 p t s i n t q t c o s t 其h a m i l t o n 能量函数为 h p 2 r q 2 t 1 作为比较 该方程同时用辛算法 四阶辛p r k 方法 y e e 格式和非辛算法 四 阶r k 法 计算 计算时间步f 0 8 执行时间r 0 至f 2 5 0 图4 1 给出了它们随 迭代时间步的数值能量h 的比较 2 0 辛紧毁格式f d t d 方法的研究 第四章h a m i lj o n 系统与辛算法 工 姗i 盘 j 型 裁 能量误差比较 圈4 1r k 格式 辛格式 y e e 格式计算所得的能量误差比较 从图4 1 中可以看出 辛格式的能量值始终在1 附近振荡 四阶辛格式要比y e e 格式的效果好 震荡的幅度较小 而非辛格式 能量的衰减是很显然的 同时 给出 相空间中p q 值随计算时间变化的轨迹图 表示为图4 2 辛算法的运动轨迹是比较 清晰的 而r k 方法则没有可分辨的运动轨迹 有很严重的耗散效应 r 潜式 4 r i d 辛格式 4 n q 图4 2 r k 格式 辛格式计算 p q 随时间运动的轨迹比较 2 l 叽 们 蚰一 的 雌 们 们o 辛紧致格式f d t d 方法的研究 第五章辛算法在m a w e l l 方程中的戍用 第五章辛算法在m a x w e i 方程中的应用 5 1m a x w e ii 方程组的h a m ii t o n 表述 对一般的m a x w e l l 方程 v x 雷 堡 j 研 v 肛一署一无 定义h a m i l t o n 函数 1 8 1 9 l 蜃 西 昙4 雪 v 后 喜西 v 1 3 一了 豆 五 西 上 s 则m a x w e l l 方程能表示为下面h a m i l t o n 典则方程的形式 a do h a t8 b a 豆o h o t 0 1 这样就可以将适用于h a m i l t o n 系统的辛算法应用到求解m a x w e l l 方程的f d t d 方法 中来 5 2 辛f d t d 的计算格式 5 2 1 无耗媒质的情况 首先考虑磁导率盯 0 电导率盯 0 的简单情况 这时上面定义的h a m i l t o n 函 数为可分 可以利4 3 节的p r k 的构造方法来建立显式的辛差分格式 为了和传统的f d t d 格式取得一致的表现形式 即用豆 疗来取代西 后 则上述 m a x w e l l 方程的典则形式可改写如下 n 饵 两 当 疗 v 心 豆 v 锺 二 2 2 辛紧数 各af d t d 方法的研究第五章辛算法在m a x w e l l 方程中的应用 盟 一鸳 一l v 豆 o to nu 丝 罂 v 厅 8 ta h 则应用4 3 2 节的二阶显式p r k 方法的结果 则有 以 豆 疗 作为开始迭代量 丘 h 一一坐c 1 v 豆一 磊 豆n 6 t d v 霞 占 雷 露一a tc 2 v 磊 雷 磊 a t d v 厅 例如取系数 c l c 2 1 2 4 1 吐 0 则 蜃 h 一 a t v 旁 5 1 a 1 2 豆 豆一 a t v 露 5 1 b 占 驴1 再一尝 v 槽 1 5 1 c 考虑这样的迭代过程 可以合并 5 1 a 5 1 c 式 消去曰项 从而得到 丘 露一a t v 豆一 豆一 豆一 坐 v 厅 丘 丘一a t v 雷川 毋 譬 v 丘 若令昂 曰n l 2 则得到上面的迭代格式如下 辛紧致格式f d t d 方法的研究第五章辛算法在m a x w e i i 方程中的应用 疗 2 疗 2 一a t v 豆一 x 豆 豆一 a t v 疗 2 从形式上这就是传统的f d t d 方法 即y e e 格式二阶时间交错格式 可以看出 传统 的f d t d 方法是最简单的一类辛算法 又如取四步显辛p r k 法如下 丘 疗 n 一坐q v 豆n 磊 霹 丝吐 v 异 恿 p 竺c 2 v 磊 磊 磊 a t d v 丘 占 丘 丘一生岛 v 龟 霄 n 丘一坐c 4 v 磊 磊 磊 垒以 v 丘 占 豆 磊 坐以 v 豆n 占 例如取系数为 i c如4 1 6 2 训a 咖c2 c一3 1 6 1 a 1 3 2 1 3 12 a d 4 0 舯口 扼 瓶碣 鸸 口 哎 一 实际计算过程中合并q c 4 的两步迭代 精度达到了四阶 5 2 2 有耗媒质的情况 对m a x w e l l 方程磁导率 0 电导率盯 0 的情况 由于其h a m i l t o n 函数是不 可分的 原则上不可以用4 3 节的方法来建立显式的辛格式 但仍然可以把上述构造 辛格式的过程形式地应用到一般的m a x w e l l 方程中来 以 2 2 a 中第一个方程为例 堡一堡 盯e 占堡 印 z 4 a f 在一个辛时间步d t a t 内离散上式 则有 扭 砂 o h y a z 叫珊 坐盘掣 壁赵掣 辛紧毁格式f d t d 方法的研究第五章辛算法在m a x w e l l 方程中的应用 层 j 珊 扣栌酬 嗍 c b 坝等一警 k中ca m 垫然 cb 咖丽而2djat2e m r m d at面 2 占 胧 盯 m f 上式括号中m i 1 2 j 七 对 2 2 a 中的其他方程有类似的离散表达式 左端项 没有用类似传统的f d t d 展开 因为这里说明的仅是时间显式辛格式建立的过程 并 且在本文中 空间导数离散应用的是紧致差分格式 而非二阶中心差分格式 同样对 2 2 b 中第一个方程可离散为 2 扣尹1 c p 脚孵m 扣争c 鼬噜一争 其中凹 咖篇镂 c q m 丽而2 c f a t 砸 匕式括号中m f 1 2 k 1 2 1 对 2 2 b 中的其他方程有类似的离散表示式 在本文下面的计算分析中 总以取这样一组四阶精度的辛系数作为讨论对象 t 9 1 t 1 0 1 1 2 l fc i 0 1 7 3 9 9 6 8 乞 0 1 2 0 3 8 5 0 4 岛 o 8 9 2 7 7 6 2 9 9 c 4 c 2 岛 c i i 吐 o 6 2 3 3 7 9 3 吐 0 1 2 3 3 7 9 3 2 吃 吐 吐 吐 以 0 实际计算时 由于以 0 可以合并q c 5 的计算过程 5 3 辛f d t d 格式的稳定性分析 由上面的分析得知 辛格式的构造可建立在r k 方法上 所以对辛格式的稳定性 分析亦可以由r k 方法的稳定性分析推广而来 下面首先介绍对单一变量的r k 的稳 定性分析方法 i 由此再继续讨论p r k 方法的稳定性分析 5 3 1 单一变量的眯稳定性分析 设某数值方法在节点以处的数值解应为n 但实际计算时因为有舍入误差 算得 辛紧致格式f d t d 方法的研究 第五章辛算法在m a x w e l l 方程中的应用 的值为或 记皖 见 y 称为第n 步数值解的扰动 若某一种数值方法在节点矗处的数值解为y 伴随有大小为 的扰动 由瓯的传 j 瓯j i 皖i m n l z 2 1 咒 以 五 酏 f 毛 虬 h j t z 岛 令 a 呈 七 主 e 爿 三 芝 儿 y n h a b 7 j j 5 2 由式 5 3 求得k i m 爿 t 圪代入 5 2 式得到 令z h 2 0 于是 称r z 为稳定性函数 可求得 h 1 i h 2 b 7 i h a a p 帆 r z z b 7 一z a 一e 儿 r 0 当陋 z isl 时 格式为绝对稳定 r z d e t 蕊l z 百a z e b r 辛紧敛格式f d t d 方法的研究第五章辛算法在m a x w e l l 方程中的应用 5 3 2 辛格式 p r k 的稳定性分析 对p r k 格式的稳定性分析 因为l a r k 格式含两个变量 且取不同的推进步 b u t c h e r 列阵不一致 上述的对单变量的i l k 方法的稳定性分析已不适用 同时注意 到上述r k 方法稳定性分析 一般给出z 与实轴的交集 这样对试验方程 可以简化 为 y m y r 国 0 对试验方程 埘 r 甜 0 或写为 鲁 孑 j j 为单位辛阵 用p r k 方法得 曰 p 吩g q q 一矗国 a u p j p 1 p h c o 岛q q g 一砌 配 改写上面的迭代为矩阵形式 孙砌冰荆纠 阱阶国繇批 令只 多 一 三 以 苫i o j 马 兰 方程 5 4 5 5 写为 5 4 5 5 p 蕙 埘 i j i 咖 盛由一出 辛紧致格式f d t d 方法的研究第五章辛算法在m a x w e l l 方程中的应用 z y 而珊 呜出 5 6 虬 i 虬 向珊 e 如 5 7 醐廿e h 东 以 i 兑 h o b 7 露 5 9 y n l s y s 扭7 i z 4 一1 e 三 h c o r z t l o z l o 类似稳定性分析 辛格式 p r k 有迭代形式 y s y 其递增矩阵的特征值为 扎 圭f 跏 f 1 1 t r s v 2 v h c o 2 f s v f 容易知道i i 1 即格式是无耗散误差的 格式的色散误差可由上面的定义表示 v v a r g 2 或以归一化相速表示为妒 v 里丛型 则可以利用稳定性分析求得的递增矩 阵的迹s v 来给出相速分析 作为比较同时给出四阶r k 方法递增因子为 删一罢 知心一争 3 0 辛紧致格式f d t d 方法的研究第五章辛算 法在m a x w e l l 方程中的应用 时间色散误差 图5 1 时间离散格式的相速分析比较 图5 2 在基本稳定区间 0 2 内以归一化相速误差给出了辛格式的时间离散色散 分析 同时给出了其他一些格式的相速误差作为比较 从图中看出 辛格式 4 n d 比其他一些格式更接近准确解 色散误差更小 rll剃婆挈 驭 辛紧致格式f d t d 方法的研究 第六章辛紧致f d t d 格式的行为分析 第六章辛紧致f d t d 格式的行为分析 6 1 辛紧致f d t d 格式的色散分析 本节将分析辛紧致f d t d 格式的时间 空间离散联合的色散特性 以一维 m a x w e l l 为例 写成矩阵形式 a e 西 8 h 疣 a h 8 e o t a x 8 e a h o t占0 x o 三旦 占盘 一三旦 o 咖 6 1 首先用f o u r i e r 方法分析空间导数离散化的数值色散误差特性 对平面波电磁场 的各分量可以表示为列矢量 f f t e 一淑 为波数代入 6 1 对空间离散化并 以差分代替微分后可以将它化为以下形式的关于时间的矩阵微分方程 掣 q 弛 6 2 d t 该方程具有e 山形式的解 五是系数矩阵q 的特征值 以标准的y e e 格式计算为例 得 g o 12 s i n k t x 2 1 2 s i n k a x 2 占 x 0 其中c o c o 为空间离散的数值波数 警 一占 警 一 主璧塞堡丛 望查些盟堕壅 篁盔童圭鉴墼婴堡堑垄盟堑塑坌堑 然后按照上面对时间的色散分析 可求得递增矩阵 1 2 m 鸩 嘲 鹕 三a t 3 m 鸩2l 1 a t 2 m 1 m 其中 一i ii c o c o m 2 掣u x s c 特征值为 丢 r 2 m 如土r i j 丽 以m 鸠代入得 u 1 i 1 面a t 2 国研士f 一吉y 国k 国 2 土 i 石丽 其中y 当 坐为c o u r 趾t 稳定性值 8 厶x 2 奴 定义归一化相移为 v 竺笔釜堂 实际上分析比较格式的色散 相速 误差是比较买际相移c o a t 与数值稠移 a r c s i n 1 一丢 2 一伊脚 2 2 的逼近程度 则色散方程定义为 础 a r c s i n l 一三 2 叫2 国 2 2 为了使数值相移有简单的形式 则重新写为 s i n 2 础 l 一 2 v 2 c o 缈 2 2 础 1 一掣 s i n 2 芋 三嘲甜 2 v 2s i n 2 竽 6 3 6 3 式即为所熟知的一维y e e 格式的色散方程 对其他维数t 其他格式可进行类似的 色散分析 得到如下色散方程 以三维情况列出 y e e 格式 2 n d 2 t h s i n 2 警 却蚴2f 专s i n 2 竽 矿1s i n 2 竽 十万1s i n 2 竽 其中t k y 屯分别为x y z 方向的波数 f a n g 格式 2 n d 4 t h s i n 2 譬 邛 f 2 古s i n 2 竽 c o s 2 竽 c o s 2 竽 寿s i n 2 竽 c o s 2 竽 c o s 2 竽 虿1s i n 2 睁 耐 竽卜 竽 f a n g 格式 4 竹d 4 t h 咖2 竽 却埘b n 2 竽n 1 2 一s i n 2 竽 扣2 竽 一l p 2 s i n 2 扣2 竽 灿n 2 竽圳 其中 p 出y f 壶甜 竽 古耐降 吉血2 降 辛格式 4 n d 4 n 1 咖2 竽h 1 丫1 回 其中 护 回 2 喜昂 c 2 2 旌 刃 记 r g c d j r c i 叱 c d j d d b o d i c l s j 5 i t 2 j j j p e m 1 q h j 2 s5 0 p j m 玎 儿 f 竺笋 占 缸 缈 位 为空间导数离散f o u r i e r 分析所得到的各方 向的数值波数 例如对空间四阶f a n g 格式 t 塑螋瓮掣 2 7s i n k y a x 2 一s i n 3 k y a x 2 巩刮 i 面 一 珑 j 2 7 s i n k a x 瓦2 石 s i n 3 k a x 2 对空间四阶紧致格式 盯f 丝型坐盟塑堂警型坚竺幽 x 2 a s i n k y a y 2 2 6 s 缸3 一缈 2 缸c o s 缈 1 叫 忑f 一 仉 f 2 a s i n k a z 2 2 b s i n l 3 k f a z 2 2 a c o s 一k a z i 一 作为比较同时给出r k 法与四阶空间离散格式的联合色散关系 表达为递增因子形 式 仃 f 一塑萼世 丝掣 t c a t r l f l y 1 1 一塑警型1 根据上述方程 图6 l 给出了一维情况下各釉格式的色散误差比较 设 血 a y a z 以归一化相速误差表示 并以每波长分割数作为横坐标 其中c o u r a n t 稳审幛值靴取0 5 擎紧致格式f d t d 方法的研究第六章辛紧致f d t d 格式的行为分析 蚓 罂 挈 1 鱼 不同格式的联合色散分析 图6 1 一维情况下联合格式的相速误差比较 由图6 1 可以看出 辛紧致格式比传统的y e e 格式和其他的高精度格式的色散 误差都小 有更好的数值分辨率 6 2c o u r a n t 稳定性 通过上面的分析 得到了各种格式的色散方程 根据这些方程 我们很容易给 出这些格式的c o u r a n t 稳定性值 例如对y e e 格式色散方程 s i n 2 警 邓 2 古s t n 2 竽 寿s n 2 竽 古s n 2 竽 要求右端项小于1 则有 若令缸 缈 厶z 万 警为c u r 8 n t 稳定性值 对y e e 格式则有 墨 忑1 其中门为空间维数 其他格式可用类似的方法得到相应的c o u r a n t 稳定性值a 表6 1 给出了上述一些计算格式的c o u r a n t 稳定性值 可见四阶辛紧致格式的稳定性值是最 辛紧毁格式f d t d 方法的 阱究 第六爹辛紧致f d t d 括式的行为分析 大的 一般情况下都取0 5 作为实际计算时的c o u r a n t 稳定性值 表6 1 不同格式的c o u r a n t 稳定性值 s m a x y e e f a n gf a n gr k f a n gr k c o m ps y f a n g s y c o m p 2 n d 2 t h 2 n d 4 t h 4 n d 4 t h 4 n d 4 t h 4 n d 4 m 4 n d 4 t h 4 n d 4 t h 一维 10 8 5 7 1l i 1 5 0 3 l 1 5 8 31 2 9 0 61 2 9 9 5 二缝0 7 0 70 6 0 60 7 0 7o 盘1 3 5o 8 1 9 009 1 2 60 9 i s 9 三维05 7 70 4 9 4 60 5 7 70 6 6 4 10 6 6 8 70 7 4 5 1 0 7 5 0 3 6 3 辛紧致格式解的效率 为了确定不同格式解的效率 这里给出不同格式在同一色散误差限制下的存储 量与计算量 6 3 1 最小网格分割数 按照s h l a g c r se ta l 2 6 的定义 23 6 0 0 砉一1 为每波长色散误差 以度为 单位 其中一称为数值波长 与数值波数的关系为 孕c o s 庐 s i n 8 a 为自由 阿 空阃波长 若色敷误差限制在1 度的范围内 表6 2 绘出不同格式的每波长最小分割 数 表6 2 误差1 度限制下的每波长最小分割数 占m y 毫e f a n gf a n g r k f a n g r k c o m p s y f a n gs y c o m p 2 n d 4 m 2 r i d 4 t 1 1 o n d 4 t 1 1 4 n d 4 t h 4 r i d 4 t h 4 n d 4 t h n d 4 t h f 二维1 8 1 66888 7 三维2 01 2798 87 由此结果 下面将给出在此色散误差限制下各种格式的存储量与计算量比较 6 3 2 存储量与计算量 给出不同格式的存储量 首先要知道每个计算单元网格上的存储量 这里设 辛紧致格式f d t d 方法的研究 第六章辛贤致f d t d 格式的行为分析 s t o r a g e c e l l 为单元网格存储量 2 8 为每波长网格分割数 s t o r a g e 2 2 为每波长平 方的存储量 则二维 三维存储量表达为 二维 s t o r a g e 2 2 s t o r a g e c e l l 2 1 6 2 三维 s t o r a g e 3 3 1 一 s t o r a g e c e l l 五 8 3 表6 3 6 4 分别给出了上述格式的存储量比较 可以看出辛紧致格式的存储量在众 多格式中是有明显优势的 表6 3 二维情况格式每平方波长的存储量比较 y e e f a n gf a n gr k f a n gr k c o m ps y f a n gs y c o m p 2 n d 2 t h 2 n d 4 t h 4 n d 4 t h 4 n d 4 1 1 i 4 r i d 4 t h 4 n d 4 t h 4 n d 4 t h s t o r a g e c e l 3 3 3 1 5 1 533 棚1 8 1 6 6 8 887 f s t o r a g e 2 9 7 2 7 6 81 0 89 6 09 6 0 1 9 2 1 4 7 y e e f a n gf a n gr k f a n gr k c o m ps y f a n gs y c o m p 2 r i d 2 t h 2 n d 4 t h 4 n d 4 t h 4 n d 4 t h 4 r i d 4 t h n d 4 t h 4 n d 4 t h s t o r a g e c e l l 6663 03 066 竹坩 2 0 1 279887 s t o r a g e 3 4 8 0 0 0 1 0 3 6 8 2 0 5 82 1 8 7 01 5 3 6 03 0 7 2 2 0 5 8 对计算量的讨论 需要知道格式的计算步数 当然计算机处理对实数的加减 乘 除的时间消耗是不一样的 为讨论方便 这里不再加以区别 都作为一个算时考虑 这里设 o p e r a t i o n 丑2 a t 为每平方波长 每单位时间的算时数 o p e r a t i o n c e l l f 为每网格 每单位时间的算时数 则二维 三维计算量表达为 二维 o p e r a t i o n 2 2 a t o p e r a t i o n c e l l a t x 1 s t 8 3 三维 o p e r a t i o n 2 3 a t o p e r a t i o n c e l l f 1 s m a 2 8 4 表6 5 表6 6 分别给出格式的计算量 从中可以看出 辛紧致格式的计算所消耗时 间是比较少的 在二维情况下 和f a n g 4 n d 4 t h 格式几乎一致 在三维情况时间消 耗是最少的 辛紧致格式f d t d 方法的研究第六章辛紧致f d t d 格式的行为分析 表6 5 二维情况格式每单位时间每平方波长的计算量比较 y e e f a n gf a n gr k f a n gr k c o m ps y f a n g s y c o m p 2 n d 2 t h 2 n d 4 t h 4 n d 4 t h 4 n d 4 t h 4 n d 4 t 1 1 4 r