（计算数学专业论文）高频波在声波导中传播的计算.pdf

摘要 众所周知，高频h e l m h o l t z 方程的波传播问题很难求解，是热门的研究领域。 产生困难的原因是高频时，解是高振荡的，为达到一定精度，给定计算区域上的 离散问题计算量很大。各种不同类型的方法被应用于求解高频h e l m h o l t z 方程， 如有限元方法、边界元配置方法，渐近分解方法及离散奇异卷积方法，但都存 在缺陷。鉴于改进步进方法在大范围传播的缓变波导中的有效性，本文将此方 法应用于高频波传播计算问题。理论分析及数值实例显示，如果离散格式能够 精确逼近二阶偏导算子的特征值，步进方法对高频波传播计算是有效的。传统 的应用二阶及四阶离散格式的步进方法对高频波失效，原因就在此。因此，本 文将对二阶偏导算子的特征值逼近度高的切比雪夫拟谱方法应用于步进方法， 数值实例显示此时步进方法对计算高频波传播问题有效。 关键词：h e l m h o l t z 方程，高频波，改进步进方法，切比雪夫拟谱方法 a b s t r a c t a s w ea l lk n o w , t h ew a v ep r o p a g a t i o np r o b l e m sf o rh e l m h o l t z e q u a t i o ni n v o l v i n gl a r g ew a v e n u m b e r sa r ed i f f i c u l tt os o l v en u m e r i c a l l y i ti sav e r yh o tr e s e a r c h f i e l d t h eo b s t a c l eo ft h i sp r o b l e mi st h a tw h e nt h ew a v e n u m b e rki sl a r g e t h es o l u t i o ni sh i g h l yo s c i l l a t o r y i no r d e rt og e tac e r t a i nd e g r e eo fa c c u r a c y , m u c hm o r e u n k n o w n sa r en e e d e d a sar e s u l t , t h ec o m p u t a t i o n a lc o m p l e x i t yo ft h i sp r o b l e mi s q u i t el a r g e av a r i e t yo fm e t h o d sa r ea p p l i e di ns o l v i n gt h i sp r o b l e m ，s u c ha sf e m ， b e m ，a da n dd s c ，b u tt h e ya l lh a v et h e i rd i s a d v a n t a g e s s i n c et h er e - f o r m u l a t i o n m a r c h i n gm e t h o di sv a l i di nl o n gr a n g ew a v ep r o p a g a t i o np r o b l e m sw i t hl o wa n d g r a d u a lv a r i a t i o nw a v e n u m b e r s ，i ti sa p p l i e di nh i g hf r e q u e n c ya c o u s t i c a lw a v e g u i d e s p r o p a g a t i o np r o b l e m s t h e o r e t i c a la n a l y s i sa n dn u m e r i c a le x a m p l e ss h o wt h a tt h e m a r c h i n gm e t h o di sv a l i di ft h ee i g e n v a l u e so ft h et w o - o r d e rp a r t i a ld i f f e r e n t i a t i o no p e r a t o ri se x a c t l ya p p r o x i m a t e d n u m e r i c a le x a m p l e ss h o wt h a tw h e nt h ec h e b y s h e v p s e o d o s p e c t r a lm e t h o di su s e di na p p r o x i m a t i n gt h et w o - o r d e rp a r t i a ld i f f e r e n t i a t i o n o p e r a t o ht h em a r c h i n gm e t h o di se f f i c i e n ta n dv a l i df o rc o m p u t i n gt h eh i g hf r e q u e n c y w a v ep r o p a g a t i n gp r o b l e m k e y w o r d s ：h e l m h o l t ze q u a t i o n ，h i g hf r e q u e n c yw a v e ，r e f o r m u l a t i o nm a r c h i n gm e t h o d ， c h e b y s h e vp s e o d o s p e c t r a lm e t h o d 浙江大学研究生学位论文独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作及取得的 研究成果。除了文中特别加以标注和致谢的地方外，论文中不包含其他人已经发 表或撰写过的研究成果，也不包含为获得迸婆盘堂或其他教育机构的学位或 证书而使用过的材料。与我，一同工作的同志对本研究所做的任何贡献均己在论文 中作了明确的说明并表示谢意。 学位论文作者签名： 签字日期：年月日 学位论文版权使用授权书 本学位论文作者完全了解迸姿态茔有权保留并向国家有关部门或机 构送交本论文的复印件和磁盘，允许论文被查阅和借阅。本人授权逝鎏盘堂 可以将学位论文的全部或部分内容编入有关数据库进行检索和传播，可以采用影 印、缩印或扫描等复制手段保存、汇编学位论文。 ( 保密的学位论文在解密后适用本授权书) 学位论文作者签名：导师签名： 签字日期： 年 月 日签字日期：年月 日 致谢 短短两年的研究生生活即将划上句号，回忆这两年来的学习和生活，我心 中充满感激。相比本科的我来说，研究生的生活让我迅速地成长起来。在此衷 心感谢我的导师朱建新教授。在读硕期间，特别是论文写作期间，朱老师始终 给予我细心的指导，他多次询问研究进程，并为我指点迷津，帮助我开拓研究思 路，还不厌其烦的多次修改我的论文，并提出意见和建议，让我能顺利的完成论 文和学业。从他身上，我学到了无论做人还是做学问都要严谨求是的态度。 同时，向我的同门、师兄师姐、师弟师妹表示诚挚的感谢。感谢师兄宋仁 成、张学仓，师姐陈芝花，同门沈浙奇，师弟黄文星、娄琪，师妹陈赠思，和你们 一起度过了快乐的研究牛牛活。特别是沈岛雪，和非常怀念和你一起自习，探 讨学术的时光，你开朗幽默的个性让人记忆深刻。最后，感谢我的父母对我无 私的爱和奉献，感谢他们给予我的支持和包容。 第一章背景介绍 考虑如下边值问题：- a u k 2 u = f ，为时谐波源，k 是波数。该问题有 很晕要的物理应用，特别是在海洋和电磁波传播问题中。当波长和计算区域 的长度相比非常小时，即当k 较大时，就变成高频问题。h e l m h o l t z 方程数值解 的精确度极大地依赖于波数圪。一般来讲，当k 很大时，h e l m h o l t z 方程的解是高 振荡的。高频h e l m h o l t z 方程由于求解困难，各种不同的方法被应用于求解此 类方程，传统的如有限元法( f e m ) 以及改进【l ，2 ，3 】，较新的如边界元配置方 法( b e m ) 【4 】，渐近分解法( a d ) 【5 】，离散奇异卷秘法( d s c ) 【6 1 ，是当前热门 的研究领域。这些方法存在的问题是，为了让数值解达到一定的精度，单位波 长的离散点个数要很大，实际模拟中的计算量就会很大，不仅如此，当保证单位 波长离散点数不变，即分辨率不变，随着k 增大，数值格式的误差也会大幅增加， 产生所谓的污染效应。因此，寻找一种合适的计算量小且稳定的数值方法是关 键。 1 1 传统方法 b a b u g k a 及他的合作者，对求解高频h e l m h o l t z 方程作了一定的研究，特别 是对于所采用的数值方法关于波数的稳定性作了探讨。他采用了有限元方 法( f e m ) 来求解高频h e l m h o l t z 问题，并提出了污染效应这一概念，即数值方 法关于波数的稳定性。虽然高频问题的计算精度取决于单位波长离散点的个 数，i i l j a h 越小，精确度越高。但即使a h 保持不变，即分辨率不变，随着k 增大， 数值格式的误差也会大幅增加，产生污染效应。所谓污染效应【l 】，即若存在常 数c ，7 ，s r ，t 0 ，使得数值解的误差e 满足烈c ，也就是说，当h 矿为 常数时，随着k 的增大，数值解的误差会越来越大，则说数值格式具有污染效 应。针对传统的g a l e r k i nf e m 方法，b a b u s k a i 及他的合作者提出一般化的 改进的f e m 方法，可以用来求解高频问题，但存在缺陷，在一维问题中，可避 免污染效应，但在二维或更高维问题中污染效应不可避免，但可设计达到最 小。特别地，在文献【7 ，8 】中，给出了在日1 半范数下的危矿有限元解的误差，满 足e 1 a p + q k p 2 ，其中p = ( t c h 2 p ) p ，c 1 ，q 为与k ，h 无关的常数，上式右端 第一章背景介绍 2 的第二项表示的就是污染项。对有限元方法而言，污染效应无法避免。 除了污染效应外，有限元方法的计算复杂度也很高。以三维问题为例，有 限元方法( f e m ) 涉及的未知元个数将达k o 的三次方，其中，c 是波数，a 是问题 的维数，事实上有可能更差【7 】。因此，b a b u g k a 又发展出对有限元方法的改进算 法。这些算法【2 ，3 】的共同特质是在设计数值格式的过程中，用到解高振荡的特 性，从而减少了问题的计算复杂度。 1 2 最新进展 s e o n g j a ik i m 5 ，采用渐近分解( a d ) 方法求解高频h e l m h o l t z 方程。对于 相对光滑的介质，达到相对精度单位波长只需采用4 5 个离散点。而二阶差分格 式理论上单位波长要采用1 0 - 2 0 个离散点，实际模拟中，单位波长要采用2 0 3 0 个 离散点。该方法大大地改进了标准二阶差分方法，对于二维及三维问题，计算 速度提高1 0 1 0 0 倍。 e i d a rg i l a d i 【4 】发展出一种针对高频h e l m h o i t z 方程积分形式的边界元配置 方法( b e m ) ，其中基函数是渐近得到的。每个基函数是光滑的振幅和高振荡的 相位因子的乘积，就像方程的渐近解。相位因子由几何光学( g o ) 和散射几何 理论( g t d ) 得到，而光滑的振幅是高阶样条函数。这是一种高阶的数值方法， 其中未知元的个数几乎与k 无关。离散问题可有具有d ( ) 个元的稀疏矩阵a 近 似，其中是基函数的个数。 g a n gb a o 6 提出了一种新的求解高频h e l m h o l t z 方程的方法，离散奇异卷积 方法( d s c ) ，即局部的谱方法。其优势为：其一，改进有限元方法，选取合适的 参数将完全避免污染效应，离散奇异卷积方法( d s c ) 方法本质上无污染，且污 染效应并不会随维数的增加而增加；其二，单位波长仅取两个离散点即可得到 高精度的数值解，计算量小。存在问题是，数值离散方程过程中将会用到求解区 域以外点的信息，对于不同的方程采用不同的边界条件【9 】。如f o k k e r p l a n c k 方 程，对于自然边界条件，用0 替代求解区域以外点处的函数值；对d i r i c h l e t 边界 条件，用边界处函数值替代；对于纽曼边界条件，用边界处的函数值及其导数替 代；对于特殊的具有周期性解的问题，可以采取平移的方法解决。但对于其它 的问题，区域外点处的函数值的替换将大大影响数值解的精度。 第一章背景介绍3 1 3 本文采用的方法 y y l u 1 0 ，针对缓变的大范围传播的波导问题作了研究，提出了单向步进 方法及改进的步进方法。传统的有限差分法和有限元法涉及大量的未知数不适 用于大范围传播问题。而改进步进方法方法对于求解缓变波导的低频波传播问 题非常有效，且在传播方向上可以采用大步长。 1 4 本文结构 本文第一章介绍高频波实际物理背景和现有的求解方法，第二章给出本文 采用的大步长步进方法的介绍及具体实现过程，第三章探讨大步长步进方法在 求解高频波时遇到的障碍及原因，采用三种不同离散格式，比较了它们在求解 高频波中的优劣，并得出了相应的结论，附以数值实例加以验证。 第二章大步长步进方法介绍 2 1 单向步进方法介绍 单向步进方法在缓变波导的波传播问题中有广泛使用。最原始的单向 模式由t a p p e r t 1 l 】在7 0 年代提出，应用在水下声波中。差不多同一时期，又 由f e i t 和f l e c k 1 2 弓i 入光波导巾。更精确的单向模式由之后的c l a e r b o u t 1 3 和 g r e e n e 1 4 提出，分别被应用在几何物理学和水下声波中。由于历史原因，单向 模式方法在声波和光波中分别被称为抛物方程法( p e ) 和波束传输法( b p m ) 。y yl u 在文献【1 5 】中，给出了单向步进方法应用于无限长缓变波导波传播问题的 介绍及主要思想。同时jxz h u 在文献【1 6 】中验证了单向模式在缓变波导中的有 效性，并作了进一步的推广，即考虑更为复杂的分层且交界面弯曲的缓变波导 的例子。下面简要介绍单向步进方法【1 5 】，考虑如下的t e 极化h e l m h o l t z 方程，其 中k 表示声波的波数： 磋岛+ 畿岛+ 圪2 ( z ，2 ) 岛= 0 ( 2 1 ) 如果波导在z 方向上不变，则可将波分解为乱= u + + u 一，其巾札= 岛，前面的 分解元让+ 精确满足如矿= i m u + ，其中m = 硬了虿，是根号算子。对缓变波 导，k = 代( z ，z ) 随z 方向缓慢变化。如果我们只对向z 轴正方向传播的波感兴趣， 则原始h e l m h o l t z 方程( 2 1 ) 可以用下面的单向模式近似： 免”i m u( 2 2 ) 方程( 2 2 ) 就称为单向h e l m h o i t z 方程，与初始h e l m h o l t z 方程( 2 1 ) 相比，此完美 的单向方程( 2 2 ) 更易于求解，但存在同有误差无法消除，因为它只是对原 始h e l r n h o l t z 方程的近似。求解该单向方程的核心在如何逼近根号算子m 。 在实际数值求解过程中，如果对m 算子用谱分解的方法计算，计算量非常 的大，所以我们通常应如一个参考波数b ，例如 1p 2 畿+ k 2 ( z ，z ) = 仡：( 1 + z ) ， z = 壶硬+ 告一1 ( 2 3 ) 第章大步长步进方法介绍 5 因此，m = k 。丽。如果我们令乱= v e i ，c ，则 塞= i k 。( v q + z _ 1 ) 钉 ( 2 4 ) 最简单的近似是 厕一1 昙 ( 2 5 ) 更高阶的近似为用咖纠阶p a d 6 有理近似： 厕一喜尚 ( 2 6 ) 其中p 是正整数，a l ，b l 有下列式子精确给出： 口f = 丽2s i n 2 岛，玩= 甜仇，仇= 鼎 ( 2 7 ) 对于k 是缓变的，根号算子m 可用上述有理逼近得到。对于k 是常数的情况( 也 是本文所考虑的例子) ，单向模式( 2 2 ) 本身就是原始h e l m h o l t z 方程( 2 1 ) 的精确 转化，方程的解直接可表示为： u ( x ，z ) = u ( o ，z ) e 1 讹( 2 8 ) 若逼近根号算子m 的离散矩阵厨可对角化，则假设a 为逼近算子理+ ，c 2 的离散 矩阵，且a = v a v 一1 ，于是厨= y 以v 一1 ，其中a 是矩阵a 特征值组成的对角 阵。于是( 2 8 ) 转化为 u ( z ，z ) = e y 仍一l u ( o ，z ) = v e l 厄y u ( o ，名) ( 2 9 ) 若给定初始波场u ( o ，名) ，则可由公式( 2 9 ) 得到波传播到z = l 处的波场u ( 厶z ) = v e l v x l v 一1 缸( 0 ，z ) 。 2 2 改进步进方法介绍 改进步进方法有别于单向步进方法，避开对于根号算子有理逼近的讨论， 转而弓 2 k d t n 算子q ，和基本解算子y 1 5 ，1 7 1 。以二维h e l m h o l t z 方程为例，假定 波导沿x 轴正方向传播，横截方向z 单位化为( 0 ，1 ) 区间，在边界z = 0 ，z = 1 处添 加均匀的和与传播方向无关的边界条件。对如下方程： 扣描“2 x , z ) u = 0 0q 碚 0 ，考虑，c 2 有一个小的扰动， k 2 一k 2 ( 1 + e ) ，e = o ( 1 ) 则a a + e r 2 ，分析范数比世暑。 事实上，本文下面考虑的用二阶精度和四阶精度离散格式得到的矩阵a 部 是实对称阵，一定可以对角化，即a = v a v 一，其中y 是实正交阵。假定对角化 后得到的是a 的精确特征值和特征向量，于是 a + e 咒2 i = v ( a + e 圪2 ) y 一1 号凸 e i - g l l l 2 = 蟛鬻= 嵴鬻 因为谱范数是酉不变的，故世书雾署= 世爵。我们可以转而分析 范数比世暑。 锝= 百i l e i c x l 网l l 2 ( 3 3 ) = ：一= - - = _ 一 ij j ， j i e l a l i l 2 0 e 1 、a l 0 2 、7 若假定i i e k 2 a 一1 i l l = d ( 1 ) ，则可对订忑丽作泰勒展开，即 v z l + e k 2 a - i - - - - - 1 + 互il e 仡2 a 一1 ) 一百1 l e 仡2 人一1 ) 2 + 。( i i ( k 2 a 一1 ) 2 1 1 1 ) 现只保留前两项，代入( 3 3 ) 后得 也龄= 螋挚 也雄挚= i l e 丽新1 1 1 2 = l i e i 扣2 脚1 1 2 ( 3 4 ) 第_ 三章大步长步进方法在高频波中的应用 g = e x p ( i 考e n 2 a 一1 2 ) 则 e x p ( i l 蔫) e x p ( i 考岳) e x p ( - i l 孺e k 2 ) e x p ( i l 矢) e x p ( 考罴) 1 0 i l a l l 2 5 压厕2 乒忑磊i 面鬲万丽 = m a x 1 ，唧( 鲁忐) ，e x p ( 考矗) ( 3 5 ) n ( 3 5 ) 代入( 3 4 ) 得 丙l l e i 厕l i l 2 1 2 = m a x 7 2 焉，2 禹) 当矩阵a 无负特征值时，范数比上界即为l ；当a 有负特征值时，若扰动 量e o ，范数比上界为e l 2 高，其 中a 七+ 1 是矩阵a 模最小的负特征值由于上式作泰勒展开的前提是i i e k 2 a 一1 1 1 1 = d ( 1 ) ，且注意到入七+ 1 是矩阵a 模最小的负特征值，于是 l 2 悬= 考蔫厂石铷) 厂丽。( 1 ) e l 2 1 k 一 h即 ， 1 z 叽 设 排 驰 崞 同 。 船、l 由 久 且 一h 一，f-。il、 孙 = o a 旭 记 为 ， 值 k 征 特 一 的 钏 一 阵 h 矩 k 阶 八i t ，、- ， 魄 。i 删 应 不 孙 清 第三章大步长步进方法在高频波中的应用 因此，当矩阵a 有负特征值时，范数比仍以l 为上界。 综上所述，当砖交大时，若假定扰动量e 满足| i e k 2 a 1 1 1 1 = o ( 1 ) ，无论矩阵a 是 否有负特征值，都有 螋i l e i 瓶翌l i l 2 1 号下l i e i 厕可l 2 1 e p , j , 扰动量并不会影响e i 、舭的运算，故e l x j l 稳定。 下面以数值实例验证上述理论分析。其中对- 阶偏导理以阶精度格式和 四阶精度格式离散，具体参见公式( 2 1 3 ) 。 假定l = 1 0 ，矩阵阶数佗= 1 0 0 ，扰动量e = 4 - 1 0 6 ，对不同的k = 5 0 ：3 0 0 ， 记p 1 = 惭2 a 一1 1 1 1 r l = | | e i 厕1 1 2 i i e i 仉i | 2 r 2 = i l e , 厕l i l 2 l l e i , 压l i l 2 ，二 阶精度格式和四阶精度格式的数值结果列于表3 1 、表3 2 。从表中可以看出， 当p 1 d ( 1 ) 时，r 1 您1 ，符合上述扰动分析结果。 表3 1 ：二阶精度格式：三= 1 0 ，佗= 1 0 0 ，a = d 2 + k 2 ， 表3 2 ：四阶精度格式：l = 1 0 ，佗= 1 0 0 ，a = d 4 + k 2 j 第j 章大步长步迸方法在高频波中的应用 1 2 3 2 1 2 改进步进算法的扰动分析 在改进步进算法中，参见步进递推公式( 2 1 2 ) ，同样涉及e i h * 讥的有效性分 析，同上不赘述。下面仅讨论( i 、a + s ) - 1 的有效性。同样假定a 可对角化，且对 角化精确。当k 2 有小扰动，即a _ a + e j ，a a + e k 2 ，。 我们考虑k 是常数的情形。对于k 是变化的，可以用分段常数逼近，然 后再在每一段上分析。根据步进算法，并注意到仡是常数，此时步进算法中 的a d d2a 仳t l ， s = v 一1 q 喇k 如d = i v - a - 磊a d ，a 眦仰= v a v 一1 令s = v 一1 i 佤v = i v 一1 y 抠y 一1 v = i 瓶 号( i 瓶+ s ) 一1 = ( 2 i v y ) 一1 因此只需讨论( i 厕一1 的有效性，即分析范数比必黯 瞬= 世耸铲烂铡黔产 = i i ( 1 + e k 2 a 。1 ) 一1 2 1 1 2 1 1 1 一e ，c 2 人。1 1 1 2 当e 0 时，若矩阵a 无负特征值，则范数比上界恒小于l ；若矩阵a 有负特征值， 则范数比上界为1 一；鲁，其中a k + l 为模最小的负特征值，上述范数比仍为有 界量，接近于1 。当e 5 0 0 时，四阶精度格 式误差逐渐变大，当n = 1 0 0 0 时，可达1 0 。总体说来，只要矩阵阶数n 5 0 0 ， 二阶、四阶精度格式的离散矩阵对角化稳定。 第一章大步长步进方法在高频波中的应用 1 6 手 三 亨 参 至 ( a ) 二阶精度格式 荨 参 圣 ( b ) 四阶精度格式 图3 1 ：矩阵a 阶数n = 1 0 0 ，横轴为仡，纵轴为i i v a v a l l 2 ( a ) 二阶精度格式 亨 参 至 b ) 四阶精度格式 图3 2 ：分辨率c = k 危= 2 ，横轴为矩阵阶数n ，纵轴为0 y 1 a v a l l 2 第三章大步长步进方法在高频波中的应用 1 7 手 三 章 ( a ) 二阶精度格式，纵轴为0 y a 2 v a l l 2 ( b ) 四阶精度格式，纵轴为l gi i v - 1 山y a l l 2 图3 3 ：分辨率c = k 危= 1 ，横轴为矩阵阶数佗 综合上面分析，在用二阶和四阶精度格式逼近算子硬+ 尤2 为相应的矩 阵a 后，只要矩阵阶数不大( n 5 0 0 ) ，矩阵对角化是精确且稳定的，从而步进 过程中矩阵求逆和开根号也稳定。因此理论上，用改进步进方法计算高频波传 播问题应该是有效的。但事实上，在后续的数值计算巾，我们发现用采用二阶 精度及网阶精度格式的步进方法对计算高频波传播问题失效。 3 2 3 数值实例 例如我们考虑方程( 2 1 0 ) 的波传播问题，令艽= 1 0 0 ，2 0 0 ，取不同的正整 数m ，即对应不同的初值 u ( o ，z ) = s i n 芸( 2 m 一1 纠 初始波的频率为，：2 m f - - 1 。我们将步进方法得到的传播到z ：l 处f g = u l 与精 确解u 管耐 钍警耐= e i 、k 2 一( 吾( 2 m 一1 ) ) ”s i n 孑- ( 2 m 一1 ) 2 比较，记解相对误差为e r ( u ) = l 乱 积一乱l i l 啦? 础| 。同时，我们还将算子精确 的特征值和特征向量应用于步进方法中，用以验证若特征值完全精确，改进步 进算法计算高频波是否有效，将结果列于表3 5 、3 5 中。其中我们分别取，特征 第一章大步长步进方法在高频波中的应用 1 8 值最大的传播模，特征值最小的传播模为特例，这里我们不考虑迅衰模，即对应 特征值为负数的情形，因为此时波衰减得很快，已经传播不到z = l 处。 n = 5 0 二阶四阶精确特征值 m 口 e r ( u )e r ( u )e r ( u ) l0 2 59 9 9 8 7 7 9 9 4 7 7 5 5 e 一0 0 64 8 0 9 9 7 5 e o l01 1 2 7 9 3 2 e 一0 1 3 3 21 5 7 51 4 3 8 3 51 2 5 1 8 2 6 e + 0 0 02 3 5 7 8 6 6 e 一0 0l1 0 4 7 7 4 1 e 一0 1 3 n = 1 0 0二阶四阶 精确特征值 ”z | 口 e r ( u )e r ( u ) e r ( u ) l0 2 59 9 9 8 7 72 5 11 8 0 7 e 一0 0 62 9 7 2 5 7 2 e 0ll4 0 9 4 0 8 0 e 一0 1 3 3 21 5 7 51 4 3 8 3 51 3 9 7 5 4 6 e 舢6 2 6 4 4 7 e 一0 0l4 9 1 7 3 4 3 e 0 1 3 礼= 1 0 0 二阶网阶精确特征值 仃l 口 e r ( u )e r ( u )e r ( u ) l0 2 5l9 9 9 9 3 81 2 5 5 7 8 6 e 0 0 61 6 8 9 1 5 0 e 0 ll1 3 9 0 3 5 9 e 0 1 3 6 43 1 7 51 4 2 5 7 91 3 6 0 7 7 9 e + 0 0 07 0 4 0 9 2 1 e o o l 8 7 1 5 7 8 9 e 一0 1 4 n = 2 0 0 二阶四阶精确特征值 m 8 e r ( u )e r ( u )e r ( u ) lo 2 51 9 9 9 9 3 83 1 5 5 1 5 9 e 0 0 73 9 0 8 4 9 3 e ol21 2 9 7 4 ll e 一0 1 2 6 43 1 7 51 4 2 5 7 99 8 4 1 4 0 4 e 一0 0 l5 13 5 7 3 4 e 一0 0 21 1 9 1 7 1 7 e 0 1 2 从表3 5 、3 6 中我们可以看出，其一，如果在步进过程中采用精确的特征值 和特征向量，步进方法对计算传播模式的高频波是有效，相对误差可达1 0 _ 1 3 。 其二，在计算高频波时，四阶精度格式优于二阶精度格式，且随着离散矩阵阶数 的增加，误差减小。其三，四阶和二阶格式对特征值大的传播模有效( 即较大的 传播常数) ，对特征值小的传播模失效( 即较小的传播常数) ，无法计算。 更直观地，我们取k = 4 0 0 ，n = 4 0 0 ，m = 1 0 ，2 5 ，用二阶精度计算的波和精 确波放在同一张图中，如下图3 4 所示。图中实线和点划线表示精确解的实部和 虚部，和+ 表示数值解得实部和虚部。从图中可以看出，m = 1 0 ，m = 2 5 分别 为对应较大的正特征值和较小正特征值的传播波模，采用二阶精度格式的步进 方法对前者有效，对后者失效。 第三章大步长步进方法在高频波中的应用 1 9 ( a ) 仇= 1 0 ，p = 3 9 9 9 9 6 9( b ) m = 2 5 ，卢= 3 9 2 5 2 4 9 图3 4 ：二阶精度，横轴是z 从上述的分析可以看出，步进方法在计算高频波传播问题时失效的原因应 该是，逼近算子雳+ k 2 的离散矩阵a 的特征值对原问题的特征值逼近度不够。 文献【2 0 】指出，虽然二阶中心差分对算子硬+ 圪2 逼近的误差阶为o ( 硬) ，但 对求算子谚+ k 2 的特征值的精度并非是二阶。以如下简单的s t o r m l i o u v i l l e 问题 为例： ， 一i f , + q y 2a y ( 3 6 ) l 爹( o ) = 爹( 7 r ) = 0 比如用m 阶对称三对角阵h - 2 a o + q 近似求解问题( 3 6 ) 的前m 个特征值，其中 a o = ( a i j ；o ) ，a i i ；o = 2 ，a i ，i + 1 ；0 = - 1 ，h = 丌( m + 1 ) q = d i a g ( q ( h ) ，q ( 2 h ) ，口何一 ) ) 即用二阶中心差分近似。将特征值以一定顺序排列p 1 p 2 ，则( 3 6 ) 的 第m 个特征值为= m 2 ，而h q a o + q 的第i 个特征值为p m = m 2 + q + o l m ( g ) ， 两者之间的误差为o ( m 4 h 2 ) ，随着特征值模的增大，误差越来越大。 因此，由离散矩阵a 求得的特征值已并非原来的算子鲤+ 托2 的特征值，这也 是步进算法无法算高频波传播问题的原因所在。 鉴于上述原因，本文在改进步进方法中，采用切比雪夫拟谱方法近似算 子雳+ k 2 。下面这一节将详细介绍切比雪夫拟谱方法，并将其与二阶、四阶格 式逼近算子砖+ k 2 特征值的精度进行比较分析，分析结果显示切比雪夫拟谱方 法对算子的特征值逼近度最高。 第一章大步长步进方法在高频波中的应用 3 3 大步长步进方法的改进 3 3 1 切比雪夫拟谱方法 切比雪夫拟谱方法即是将微分算子离散成切比雪夫微分矩阵。下面先简要 介绍切比雪夫微分矩阵是如何形成的 1 9 1 。先引入切比雪夫点 z j = c o s ( j = l n ) ，j = 0 ，1 ，n 对给定的定义在切比雪夫点上的网格函数v ，我们分两步得到离散导数叫： 锄是次数的多项式，唯一确定的且满足p ( 巧) = ，0 j n 令哟= p ( x j ) 由于运算是线性的，故可以用( + 1 ) ( n + 1 ) 阶矩阵表示，我们记为d ， 即叫= d v 。对于任意的，切比雪夫微分矩阵d 由如下定理给出： 对任意的1 ，a 4 - ( n + 1 ) ( n + 1 ) 阶切比雪夫微分矩阵巩的元的指标 从0 变到，则该矩阵的每个元可表示为： ( d n ) o o = 牮，( d ) = 一半， ( d ) 面。症务， j = 1 ，一l ， ( d ) 巧= 毒茬墨券，i 歹，t ，j = o ，一1 盹= r 岩扎 从上述定理可知，切比雪夫微分矩阵d 有如下性质 ( d ) 巍= 一路硝( 巩) 巧 ( d ) 巧= 一( d r ) 一t ，n j 利用切比雪夫微分矩阵，可以用来求解如下的二阶微分方程问题： 三-：u：ui+，三au玄o+b五u5=：fi 9 一一1 z 1 第t 章大步长步进方法在高频波中的应用 2 l 将离散点代入后得如下离散方程组，记( d ) 巧= d i j ，d 。( ，k ，k = 1 ，2 分别表示一 阶和二阶切比雪夫微分矩阵： 一u ( 鼢) + 口让知( 祝) + b u n ( x d = f ( x i ) i = 1 ，n 一1 口一i i n ( x n ) + 雎( z ) = g o c + 心n ( x o ) + 风乱；v ) = 9 + 尝o ( 一矿露+ o 露u n ( x i ) + b u n ( x ) = i ( x 1 ) i = 1 ，n 一1 口一u n ( x n ) + 阻甚ou d o ) ，“( 巧) = g o + 缸( 跏) + 辟岛戳缸s ( x a = 外 下面简要叙述用切比雪夫拟谱方法逼近算子谚+ k 2 ，形成矩阵a 舭b 的过程， 假定矩阵的元的指标从0 变到。针对本文所考虑问题( 2 1 0 ) ，考虑如下特征值 问题 ， j 栅“”( 3 7 ) 【u r ( 1 ) = 0 ，u ( o ) 20 记d 譬，歹= 1 ，2 分别为逼近一阶和二阶偏导算子的切比雪夫微分矩阵，其 中鲻= 【础】2 。由于切比雪大拟谱方法中的离散点取为白= c o s ( j ，r n ) ，j = 0 ，1 ，n ，因此易【- 1 ，1 1 ，而本文考虑的是勺【0 ，1 】，采用坐标变换易= 2 ( 乃一 ) ，于是u ：= 缸z 鬈= 2 让t ，u ! = 4 让；。因此对于本文中的特征值问题，记逼 近一阶和二阶偏导算子的矩阵分别为d ，j = 1 ，2 。于是 臻； ：i 磊，2 牡三三笋爿a 姚：d 妒+ 脚 d 嚣= 4 d 嚣= 【d 品】2 一“。船6 一“川。1 结合特征值问题的边界条件u ( 1 ) = 0 ，u ( o ) = 0 ，推出 u n ( z n ) = o ，岛础“( 巧) = 0 兮乱( 匈) = 一各1 呓u ( ) 删 将乱( 匈) 代同到a 窿6 中得 a 曲曲= a c o 船d 。( i ：n - i , i ：n - 1 ) 一a 黧。( 1 ：n - 1 , 0 ) 望紫 因为在步进方法中，为加快步进计算的速度和效率，一般会将矩阵对角化， 同二阶和四阶精度格式，我们在下面分析切比雪夫拟谱方法中矩阵a 抛6 的对 第三章大步长步进方法在高频波中的应用 角化问题。从切比雪夫微分矩阵d 的性质可以看出，矩阵a c 鼢是非对称的满 矩阵，故没有二阶和四阶精度格式下的良好性质。但考虑到当代增大到一定量 时，a 如6 可以是严格对角占优的，凶而可对角化。在实际计算过程中发现，切比 雪夫拟谱方法得到的离散矩阵确实是可对角化的。下面用数值实例观察切比雪 夫拟谱方法计算特征值是否稳定。 令矩阵a 如6 阶数为1 0 0 ，k 从5 0 变化到5 0 0 ，l i y o a 如b v a l l 2 如图3 5 ( a ) 所示。 令c = k = 2 ，矩阵阶数扎从5 0 变化到5 0 0 ，相应的k = 2 n ，i i v - 1 a 如6 y a l l 2 如 图3 5 ( b ) 所示。令c = k = 1 ，矩阵阶数n 从1 0 0 变化到1 0 0 0 ，相应的k = 扎+ 1 ，l gi i v _ 1 a 妇b v a l l 2 如图3 5 ( c ) 所示。从图中可以看出，当矩阵阶数大于5 0 0 后， 误差在1 0 _ 4 左右，切比雪夫拟谱方法矩阵的对角化精度已不高，所以在实际计 算过程中矩阵阶数不可取得过大，将会影响步进方法的有效性。 宇 至 芎 至 tn ( a ) 矩阵a 阶数n = 1 0 0 ，纵轴为0 y a 。h 。b v a i l 2c o ) 分辨牢c = 2 ，纵轴为8 y a 。h 。b v a l l 2 亨 ； n ( c ) 分辨率e = l ，纵轴为l gi i y a 。 。b v a l l 2 图3 5 ：切比雪夫拟谱方法 第j 章大步长步进方法在高频波中的应用 3 3 2 特征值精度分析 下面给出二阶、四阶、切比雪夫拟谱方法三种离散格式计算特征值的 精度的比较。对于特征值问题( 3 7 ) ，经过计算可知，该问题的第仇个特征值 为入m = 尤2 一( 吾( 2 m 一1 ) ) 2 。 令k = 1 0 0 ，2 0 0 ，3 0 0 ，佗= 2 0 0 ，将a 2 ，a ，a 砒b 的特征值与a 仇= 尤2 一( i ( 2 m 一 1 ) ) 2 的相对误差列于图3 6 ，其中特征值由大到小排列。将误差记为e r ( a ) = 世i a 型e i a e t t 。图中，实线表示二阶精度格式，虚线表示四阶精度格式，点表示切比 雪夫拟谱方法。 图3 6 ：特征值精度比较，横轴为第m 个特征值，纵轴为l ge r ( a ) 从图中可看出，四阶优于二阶，较大的正特征值算得好( 即传播模计算较 准) ，相对误差在1 0 书左右。这也是四阶及二阶格式对初始波低频，传播波方 向高频的高频波传播问题有效的原因。因为初始波低频，传播波方向高频的波 对应的就是起主导作用的传播模。二阶及四阶精度格式的特征值相对误差小 于1 0 - 4 的特征值个数为1 0 左右。而切比雪夫拟谱方法相对误差阶维持在1 0 q 4 ， 第三章大步长步进方法在高频波中的应用 多于一半的特征值可以算准。三种方法中切比雪夫拟谱方法的精确度最高。 3 4 高频波传播问题数值模拟 由上述的分析知，步进方法在计算高频波传播问题时失效的原因是对二阶 偏导算子特征值的逼近度不够。而切比雪夫拟谱方法对算子特征值的逼近度很 高，因此，若将此方法应用于步进方法，理论上计算高频波传播问题有效。下面 将针对高频声波在波导中的传播问题，比较三种离散格式计算高频波的优劣， 用数值实例验证上述结果。 3 4 1 k 为常数 首先，我们考虑的是k 为常数的情形。令k = 1 0 0 ，2 0 0 ，取不同正整数m ， 即对应不同的初值，我们将步进方法得到的传播到z = l 处解仳l 与精确 解够疵= e i 、舻一( 詈( 2 m - 1 ) ) 2 l s i n 【昙( 2 m 1 ) z 】l t 较，记解相对误差为e r ( u ) = i 乱哥耐一u , , i l 乱笋积i 。同时为了验证离散格式逼近算子特征值的精度和计算结果 之间的联系，引入特征值相对误差e r ( a ) = l a 疵一a i i 入疵f 其中给定m ，对 应算子爱+ k 2 的精确特征值记为a “以= k 2 一( 芸( 2 m 一1 ) ) 2 ，用离散格式得到的 对应特征值记为a 。传播常数为p = 入一删，初始波频率为f = ( 2 m 一1 ) 4 。将 解误差和特征值误差列于表3 7 、3 8 。 表3 7 ：高频波传播计算，k = 1 0 0 ，扎为矩阵阶数 仃= 5 0r 二阶四阶切比雪夫拟谱法 m | 口 e r ( u )e r ( a )e r ( u )e r ( a )e r ( u ) e r ( x ) io 2 59 9 9 8 7 79 9 4 7 7 5 5 e - 0 0 69 9 4 8 9 8 2 e - 0 0 94 8 0 9 9 7 5 c - 0 1 04 8 1 0 9 6 8 e m l 31 8 9 ( 1 5 3 c i e 加1 21 4 2 1 2 6 l e - 0 1 6 3 21 5 7 51 4 3 8 3 5 1 2 5 1 8 2 6 e + 0 0 07 3 5 8 9 1 8 e - 0 0 12 3 5 7 8 6 6 c _ 0 0 l 5 1 1 4 2 1 9 c a 0 l1 9 3 钙1 4 i e + 0 0 02 3 7 0 8 8 l e 0 0 2 n = 1 0 0- 二阶四阶切比劣夫拟谱法 m | 口 e r ( u )e r ( a )e r ( u )e r ( a )e r ( u )e r ( a ) i o 2 5 9 9 9 8 7 7 2 5 1 1 8 0 7 e 0 0 62 5 1 2 1 1 6 e - 0 0 92 9 7 2 5 7 2 e _ o l l2 9 8 4 6 4 8 e - 0 1 47 2 8 4 6 8 9 e - o l28 5 2 7 5 6 5 e - 0 1 6 3 21 5 7 51 4 3 8 3 51 3 9 7 5 4 6 c + 0 0 05 3 8 8 9 2 3 e 4 ) 0 i6 2 6 6 禾4 7 e - 0 0 l8 4 0 7 9 5 7 e - 0 0 29 8 9 6 6 4 l e - 0 1 21 2 1 0 2 9 6 = - 0 1 4 第蔓章大步长步进方法在高频波中的应用 表3 8 ：高频波传播计算，k = 2 0 0 ，n 为矩阵阶数 n = 1 0 ( i二：阶朋阶 切比霉天拟谱法 ，n | 口 e r f u )e r ( a )e