（计算数学专业论文）步进算法的稳定性分析.pdf

摘要 本文主要研究了基于狄利克莱一纽曼( d i f i c h l e t - t o - n e u m a n n ) 映射的步进算 法在计算波的传播性态时的稳定性分析众所周知，好的微分方程数值解法要 具有：好的稳定性和收敛性基于狄利克莱一纽曼映射的步进算法在计算波的 传播形态时是收敛的，但目前缺乏对其稳定性的研究 实际上，在运用该法计算波的传播时涉及到矩阵t o + i 瓜) 的求逆，发现在 某些情况下该矩阵会产生病态，从而造成数值的不稳定。通过进一步分析本 文认为造成不稳定的主要原因是步进过程的中的误差积累特别的，当步进的 步数取值较大或者深度方向的离散不恰当时，误差积累更快。 本文对于矩阵( 亘+ i 人) 出现病态的情况，给出了这样的改进方法：改善 深度方向离散；将正则化方法和均衡法有机的结合起来，选择适当的正则化因 子对病态矩阵进行处理数值模拟表明，上述改进方法能使步进算法保持快速 稳定计算，并得到符合实际的波的传播性态 关键词：步进算法，亥姆霍兹方程，正则化，均衡法，数值试验 n a b s t r a c t s t a b i l i t ya n dc o n v e r g e n c ea r et h em a i n l ys t a n d a r d st ot h en u m e r i c a ls o l u t i o n so f d i f f e r e n t i a l e q u a t i o n m a r c h i n gm e t h o db a s e do nd i r i c h l e t - t o - n e u m a n nm a pi s w i d e l yu s e dt oc o m p u t et h ew a v ep r o p a g r a t i o ni na c o u s t i c w a v e g u i d e i t s c o n v e r g e n c eh a v eb e e np r o v e d b u tt h e r ei sl i t t l es t u d yo ni t ss t a b i l i t y t h i sp a p e ri s m a i n l ys t u d yi t ss t a b i l i t y d u r i n gt h ep r o c e s s , w eh a v et oc o m p u t ei n v e r s eo fm a t r i x ( 垂+ i x ) b u tt h e m a t r i xm a yb e c o m ei l l - c o n d i t i o n e d w h e ns t e pn u m b e r sa l o n gm a r c h i n gd i r e c t i o ni s t o ol a r g eo rd i s c r e t e n e s sa l o n gt h ed e p t hd i r e c t i o nu n r e a s o n a b l y , t h em a r c h i n gm e t h o d w i l lb en o tn u m e r i c a l l ys t a b l e t w om e t h o d sa r eg i v e nt oi m p r o v et h ei 1 1 c o n d i t i o n e dm a t r i x ( q + i f a ) t h e f i r s tm e t h o di sd i s c r e t i z et h ed e p t hd i r e c t i o na v e r a g e l y i ft h i sm e t h o d d i s a b l e d t h ee q u i l i b r i u mm e t h o di si n t r o d u c e dt oi m p r o v et h ec o n d i t i o nn u m b e r su n d e r t h e a v e r a g ed i s c r e t e n e s s t h er e s u l ti n d i c a t e dt h a t ，t h ep r o p o s e dm e t h o di sh i g h l ye f f e c t i v e ，s t a b l e ，a n dw e c a ng e tt h er e a lw a v ep r o p a g a t i o n k e y w o r d s ： m a r c h i n gm e t h o d ，h e l m h o l t z eq u a t i o n ，r e g u l a r i z a t i o n ， e q u i l i b r i u mm e t h o d ，n u m e r i c a ls i m u l a t e 1 1 1 浙江大学研究生学位论文独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作及取得的 研究成果。除了文中特别加以标注和致谢的地方外，论文中不包含其他人已经发 表或撰写过的研究成果，也不包含为获得迸江盘鲎或其他教育机构的学位或 证书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文 中作了明确的说明并表示谢意。 学位论文作者繇卢王、两 签字嗍勰年6 月7 日 学位论文版权使用授权书 本学位论文作者完全了解逝姿盘堂有权保留并向国家有关部门或机 构送交本论文的复印件和磁盘，允许论文被查阅和借阅。本人授权逝姿盘堂 可以将学位论文的全部或部分内容编入有关数据库进行检索和传播，可以采用影 印、缩印或扫描等复制手段保存、汇编学位论文。 ( 保密的学位论文在解密后适用本授权书) 学位论文作者签名：卢王、二为 签字日期：跏扩年f 月7e l 导师签名： 毕 签字嗍腓石月7 日 致谢 在论文完成之际，谨向年多年来给与我关心和帮助的老师，同学，朋友 和家人表示衷心的感谢。 首先，感谢我的导师朱建新教授! 在浙江大学学习的这两年里，朱老师 在学 - - j 上诲人不倦，让我学到了许多知识；科研上，悉心指导，并用他锐意 进取，严谨踏实的精神熏陶我；生活上，无微不至的关心和支持我在此， 对朱老师给与的教导，鼓励和支持表示衷心的感谢! 其次，感谢陈芝花，李鹏，宋仁成，张学仓，朱敏，邓小毛，沈浙奇， 沈晶雪，潘一力等在课题研究和论文撰写中给与我的帮助。特别感谢陈芝花 师姐在程序设计方面给与的指导! 感谢父母，弟弟，妹妹以及其他亲人多年来对我始终如一地关心、鼓励 和支持没有他们的支持，我也很难取得今天的成果。 谨以此文献给所有关心、支持和帮助过我的朋友们 浙江大学硕士学位论文第一章引言 第一章引言 在声学、电磁学、地震学及其它许多应用领域中的大规模的波传播问题的 研究中，通常需要在一个长度比波长还要大许多的区域中求波场的分布。以海 洋中的声波的传播为例，在探索开发海洋的过程中，海洋声学起着特别突出的 作用，因为只有声波才能在水中进行较长距离的传播，从而进行水下测距、定 位，通讯和遥感等操作。 海洋是一种极其复杂的声学介质海洋介质最具特征性的现象是非均匀性 这种非均匀性强烈地影响着海洋中的声场同时海洋波导的特性随着水平距离 也有变化，但变化非常缓慢。声速随深度的有规律变化会形成水下声道。利用 水下声道可以进行数百，甚至数千公里的远程声传播，但受海水声吸收的限制， 要求用于传播的声波频率c o 很低，即相应的波长，要很大才行。 若波源是单一频率的时间调和的，则波场的决定方程是以下h e l m h o l t z 方程， + “+ r 2 ( z ，z ) u = 0 ，d x 4 - o o ，0 d l 形= = ， 即波长远小于波导的深度，且水平距离也远大于波导的深度。 对上述的h e l m h o l t z 方程，由于七在水平z 方向上是弱衰减的，故可以把z 方 向分成三部分x o ，0 工l ，x l ，其中后边这两部分由水平的o u t g o i n g 条件连 接。当x o ，x 上时，可由分离变量的方法得到方程的精确解，这里主要考虑 0 x s l 这部分的h e l m h o l t z 方程的求解问题 当r 是常数，并且定解区域是矩形区域时，以上h e l m h o l t z 方程有精确平面 浙江大学硕士学位论文第一章引言 波解和球面波解但实际的大部分问题中，r 并不是常数，而是与水平距离x 和 深度z 等有关的，故无法直接给出h e l m h o l t z 方程的精确解 求解h e l m h o l t z 方程有许多直接的数值方法，如有限元和有限差分法等 2 5 8 o 但方程定解区域的水平距离与波长相比非常的大，而且波场的内部振荡也很大， 每一个波长范围内部需要用很密的格点( 或基函数) 来表征波场，对于这个水 平距离比波长还要大得多的区域来说，若用有限差分或者有限元方法来处理这 样一个水平区域很大的h e l m h o l t z 方程时，产生的线性系统的阶数将非常的大， 导致相当大的存贮空间，计算的代价也很高昂同时这些系统常常也是不定的， 或非对称的，这就使得方程的求解更加困难 就海洋中声波导问题来说，一方面，在波传播方向上求解区域的长度，要 比纵向( 垂直于传播方向) 的要大得多( 都远比波长大) ；另一方面，波导介质 沿着传播方向有变化，即和区域是相关的，但是很细微。显然一个好的数值方 法应该充分利用声波导的这两方面的特性。利用波导对区域水平距离依赖很弱 的特性，有很多更合适的方法 实际应用中的声场区域的底边，即水下声道的中间界面，通常是弯曲的 针对这种弯曲的底边，有两大类的处理方法，一是用“阶梯状”的折线段( 实 际上是不连续的平行线段) 来近似代替原弯曲的底边n 引。在每一段小区间内( 水 平方向的) ，视作是平坦的底边，这样原来整个大区域( 弯曲底边的) 求解问题， 就转化为一系列连续的平坦的小区间上的求解问题。这种近似方法很容易实现， 但这也同时带来了一对矛盾。首先相邻阶梯底边的不连续性会带来误差，特另l j 是水平距离很大的场合，这误差会积累得很大；另一方面，若要很好的近似弯 曲底边，减小误差，要求水平方向的分划很细，印步长要取得很小才行，但这 就失去了计算的大步长的优点 另一类方法是通过坐标变换，将带有弯曲的中间界面的区域变成易于求解 的平坦的矩形带状区域有三种坐标交换，1 ) 全局变换；2 ) 局部变换( 非正 交) ；3 ) 局部正交变换 2 浙江大学硕士学位论文第一章引言 全局变换是用一种共变映射( c o n f o r m a lm a p ) 拉平原弯曲底边n 们，同时也 保持了原方程的简单形式，但对于这种全局变换来说，当波导范围很大，且边 界很复杂的时候，计算起来非常困难广泛使用的局部坐标变换相比之下计算 起来要容易得多n9 2 ，但对于含有法向导数的底面边界条件来说，变换后得到的 边界条件是水平方向和深度方向的偏导数的组合形式。水平方向的偏导数将给 数值计算上带来很大难度 局部正交坐标变换不但能保持原方程的简单的形式，同时原的法向导数边 界条件，变换后也不合有关于水平方向的偏导数项，适合进一步的m a r c h i n g 求 解新旧坐标的的转换也可以方便的用n e w t o n 迭代来实现。 在下一章的叙述中我们可以看到，在声波导中用步进算法计算含有一层， 或者多层弯曲界面或者底边的波的传播，已经得到了很好的发展，但是在步进 算法的稳定性方面，还没有很好的进行研究基于此种目的本文主要对单测重 构的步进算法在计算波的传播形态时的稳定性做了研究。 发现在实际的解h e l m h o l t z 方程的计算中，并不是每一种离散都可以得到一 个很好的解在经过大量的数值试验之后得出，当波的传播方向和深度方向的 离散不恰当时，波的传播会变得很差经分析得，产生这种现象的原因是在计 算的过程中某个矩阵的病态导致的 众所周知，离散化技术的关键在于如何剖分曲线所在的区域为若干个小区 域，再将方程中的各阶偏导数，用相应的离散点处未知函数的差商取代后，将 微分方程转化为代数方程组，通过解代数方程组获得未知函数在离散处的离散 值，作为原微分方程的数值解，此解能否取代原方程的解取决于剖分的方法和 稠密程度。分割的技巧在于如何处理下述矛盾：分割太细，计算量大，计算误 差的积累和传递影响数值解；分割粗了，获得的数值解又不能完全反映原问题 的解因此本文也指出了针对文中条的解h e l m h o l t z 方程时如何进行剖分。 3 浙江大学硕士学位论文第二章步进方法基础 第二章步进方法基础 本章将回顾总结正向基本解算子步进方法( f f o m m ) 的发展和主要的细节 f f o m m 是y y l u 等在19 9 6 年【4 】首次提出的，之后这一方法得到了很大程度的 发展，从最初的只能求解水平分层的波导中的波的传播，到后来的可以处理具 有弯曲底部边界或者一层内部界面的波导 l 3 3 s 1 1 9 9 9 年，y y l u 的o n e - w a yl a r g er a n g es t e pm e t h o d ( o l r s m ) 方法进一步 解决了大步长计算的问题b ，们o l r s m 方法利用d t n 映射将原h e l m h o l t z 方程的 边值问题，转化为一个算子的初值问题( 将水平方向的坐标变量看作是时间变 量) ，然后再进一步用水平方向上的m a r c h i n gm e t h o d ( m m ) 方法来求解数值 实现时，若深度z 仍用刀个点来离散，q 可以用n xn 的矩阵来近似，总的内存的 需求就是o ( n 2 ) 更进一步可以用局部的特征函数系展开的方法，只需保存一个 特征值中的前个最大的特征值，通常可以取到1 0 n 刀这样原本用n x ，l 阶 的矩阵近似算子q ，可用阶数要小得多的n xn 阶矩阵来代替，所需的内存空间 将更小而且运算所需的存贮空间只和深度方向的分划大小有关，和水平距离 的大小及分划无关。因此，在海洋这种水平距离远比深度及波长要大得多的声 场环境中，o l r s m 方法与传统的直接差分法和有限元方法相比，是一种非常实 用的方法 本章的以下各节将讲述用步进方法求解带有弯曲中间界面的h e l m h o l t z 方程 的实现过程。首先给出基本方程，用局部正交变换把中间界面拉直，使得原来 求解h e l m h o l t z 方程得问题转化为一个可以用步进方法求解得形式。其次，根据 边界条件将转化后的方程进行离散，再用步进算法对其进行求解最后本章 给出了步进算法进行求解的具体过程。 4 浙江大学硕士学位论文第二章步进方法基础 2 1 基本方程 我们从二维h e l m h o l t z 方程开始 “。+ ”。+ r 2 ( 工，z ) u = 0 ( 2 1 ) 其中砌sxs 4 = o o ，0 z q，处于0 z j i l ( x ) 的一层密度为岛，处于 j j l ( x ) z 1 以下给出上端和底端的边界条件：材脚= 0 ，；n = 0 连接条件为： l i mu ( x ，z ) = l i mu ( x ，z ) z h ( z ) -z _ o ) + 上l i m o u ( x , z ) ：上l i m o u ( x , z ) , o l 。+ 6 ( 。) 踟 岛z 6 ( 。) + 锄 其中以是z = h ( x ) 的法向量，假设在x 0 和x l 时，波导是水平距离无关的 设在x 0 时，x ( x ，z ) = x o ( z ) ，h ( x ) = ；在x l 时，x ( x ，z ) = k ( z ) ，h ( x ) = 九 如果没有波从佃传过来，则在x = 三处的边界条件( 散射条件) 是 以= f 柝事雨，其中江行，在x = 0 处最简单的边界条件就是在x = o 取 “= u o ，其中u o 是一个关于z 的给定的函数假设( 2 1 ) 的解是存在且唯一的。 2 2 局部正交坐标变换 在0 z h ( x ) 印第一层进行局部坐标变换，令量= f ( x ，z ) ，三= g ( x ，z ) ， 且满足眠z ) 卜 x 佃，o z 厅( x ) o q 眠三) 卜 殳 佃，o s 三l ，其 中o ，z ) y o1 日坐标，( 量，三) 为新坐标。因为我们假设是由f ，g 产生的正交坐标变换， 浙江大学硕士学位论文 第二章步进方法基础 所以令： 六。+ 正2 0 ，g ( x ，o ) = o ，g ( x ， ( x ) ) = 1 取g ( 石，z ) 2 南 ( i ) ：当i l ( 戈) 0 时， r 静+ 三( 狮) ) 2 一o ( 2 2 ) 对于给定一新坐标( j ，三) ，可以通过( 2 2 ) 式用牛顿叠代来求解工，从而 z = 三厅o ) ； 如果给定一1 日坐标( 五z ) ，则三= 熹，同时曼也可以由( 2 2 ) 式得 h 石l 到。 ( i i ) ：如果x 满足厅( 工) = o ，令舅= 石，三= ( x ) 附注一：( x ，z ) 卜z ：q ( 量，1 ) 且z = h ( x ) 在办( z ) z q 即第二层进行局邵坐标燹挟，其中1 d 4 令： 童= 7 ( x ，z ) ，主= 季( x ，z ) ，满足： ( x ，z ) i 硼 x 佃，厅( x ) zs d + z = 王专 ( 舅，三) l 咖 o ( i = l ，2 n - i ) 现在的问题归结为求解矩阵a 的特征值，求解的这种非对称的三对角矩阵 “ 力r 0 咖 职 o m 瞰躲 姆m i t 圪 = 一亡、加生c 似卜一上 咖 咖 o 阶 阶 g 耵 g i q姆躲 力 x 联k 磐躲 l i = 力畦 m 磐 g 讲 n p 、 浙江大学硕士学位论文 第二章步进方法基础 的特征值，常用的方法是直接调用l a p a c k ，为此我们需要将a 修改为对称矩阵， 以下为修改的方法： 因为a 为三对角矩阵，故存在一个对角矩阵d = d i a g ( d 。，畋，式) ，使得 l l l l d a = s ，a d2 s d2 ，矩阵( d ：s d ：) 仍是实对称矩阵，因此矩阵 有刀个线 形无关的特征向量，且可以对角化a = 1 7 a 1 7 一， e - - z , ，吒】，a = d i a g ( , 毛，五，厶) ，且有4 巧= 乃巧 2 4 步进算法 设有两个算子q ( 曼) 和】，( 曼) ，定义为： 攻= q ( y c ) v 矿( 厶) = l ，( 主) 矿( 量，) ，于是，若能求得r ( o ) ，就能从给定的初 值v ( o ，) ，得到矿( 三，) = 】，( o ) 矿( o ，) 。 算子q 和】，满足： y ，= 一y q ( 量) q = 一q 2 一( 动i 2 + 8 + 厂) 并所需要满足的边值条件： q ( 曼) k = f 扛面獗厕， r ( o = i 其中： i 为单位矩阵 我们通过以下方法来求得】，( o ) ： 同样，我们也是在传播方向即x 方向进行离散，x o ：0 ，：三，i i l ：旦， m 置= s h ，s = 0 ，l ，m 浙江大学硕士学位论文第二章步进方法基础 以( t 巾t ) 为例，令y i ；，工皇 y ，圪，q ， n q 一，巧劬”1 巧， i = l 皇( 矿d ) m 其中令： 哆枷) 全巧( t ) ， 巧删) ：兰巧一) 矿= ( ) 删 并假定 缸) 皇( 劬p 1 ) m ， 垒巧( t ) 匕”- - zw , y , 枷) 2 = ( 乃) m 注：我们将算子转换成在工= 三处基矿( ) 下的表示，这样可以减少数值计算的误 差，否则会导致曲线在界面处的断裂 由此我们就可以通过龟，求得龟一。，e 一。，以下为具体的过程： q = 吧z p = ( 亘+ i 仄) 。1 ( 一圣+ i , g ) r ：p i r 4 x p e i r 压 龟一。= i 天( ，一尺) ( ，一r ) - 1 e 一。= e 2 ( ，+ 尸) 7 瓜( ，+ 尺) - 1 其中，x a = d i a g ( x 夏- ，压，历) ，f = 毫一o 1 3 形 卜 一 = 形 卜 吩 l 一 亿 浙江大学硕士学位论文第三章稳定性分析 第三章稳定性分析 3 1 稳定性分析基础 评判一种算法的优劣的一个重要的标准是稳定性对于一个数值计算方法， 如果对于计算中产生的舍入误差能得到有效的控制，则这种算法是数值稳定的 否则，如果这种误差得不到控制，在计算中误差得到传播甚至放大，则这种算 法就是数值不稳定的。很多方法在数学上是等价的，但用来计算，结果却不同， 也就是说它们对舍入误差的敏感程度是不同的 条件数是线性方程组d x = 6 的解x 的稳定性的有效度量。若条件数很大，则 称厶= 6 是病态方程组病态问题的一个后果就是计算解的严重失真衡量方程 组出= 6 的计算解的另一个重要尺度是计算残余向量，= 6 一a x a ，其中x a 为计算 解。若，在某种意义下是小的，则称算法是稳定的。但是，“小”并不能保证计 算解能很好的接近精确解x 】 在本文中我们计算矩阵的条件数是用如下条件数计算： c o n d ( a ) ：= 1 1 4 ：la 1 1 2 = 其中k ( 么) 和k ( 彳) 分别是半正定矩阵a r 么的最大和最小特征值 对于改进某个矩阵彳的病态性，我们可以通过改进以么为系数的线性方程 组a x = 6 来得到。 3 2 病态线性方程组的解法 对病态线性方程组求解的方法有很多，比较好的有：求解病态线性方程组 的正则化方法1 3 4 】；迭代法m ；预处理方法【3 8 l 等 1 4 浙江大学硕士学位论文第三章稳定性分析 在实际试验中，上述方法都采用过，发现前两者的效果并不好如：用迭 代法进行改善时，残余向量，= 6 一以虽然很小，但是改善后的矩阵的条件数几 乎没有什么改变，波的传播的形态也几乎没有得到改进。因此本文采用的是第 三种方法：预处理方法 预处理的基本思想如下： 对于”阶线性方程组： a x = b ( 3 2 ) 其中系数矩阵彳的条件数比较大，也就是( 3 2 ) 为病态线性方程组，在求解之 前，现对系数矩阵彳的条件数进行改善人们提出将方程组( 3 2 ) 化为易于求 解的等价方程组 匆= 6 ( 3 3 ) 其中：j ，= c x 彳仍然保持原来的特性，c 为预处理矩阵，且容易从方程组c x = y 解的x 。若j 的条件数c o n d ( j ) ：比么的条件数c 伽i d ( 彳) ：小，则用相应的算法解方程组( 3 3 ) 的收敛速度会比方程组( 3 2 ) 的要快【删 下面介绍具体的预处理方法： 因为需要改进的矩阵的元素的数量级相差比较大，所以我们采用的是行平 衡法设彳= ( ) r ，计算s = m 。；，a s ，x a , ， ( i = 1 ，2 ，刀) ，a d = d i a g ( s 。，s ：，s n ) ， 于是求出= 6 等价于求( e d ) a x = ( 占d ) 6 ，或者五= 舌( 其中占也称作正则化因 子) 这时五= ( e d ) a 的条件数可得到改善对于s 须根据具体情况取定，不能 太大也不能太小，若取的太大，会改变方程的解，若取得太小，对方程的解不 会有改善 浙江大学硕士学位论文第三章稳定性分析 3 3 步进计算的数值稳定性 对于步进计算，是否产生病态以及病态的程度，与步进的步数和深度方向 如何离散有关。当深度方向的离散确定时，若步进的步数超过一定的数值，则 会使得亘+ i 天的产生病态，且超过的步数越多，垂+ i x 的病态程度越严重 当步进的步数一定时，若深度方向的离散取的不恰当，亘+ i x 也会产生病态 但后面这种情况产生的亘+ i 人的病态程度没有前面那种情况严重。由于会出现 后面的这种情况，因此，在深度方向进行离散时，最好采用等分的方式。 经过大量的数值计算发现：对于产生病态的亘+ i 天的条件数，并不是没有 规律的。通常在步进计算的前半部分是没有病态的，囝+ i 石的条件数与不产生 病态时几乎是相同的。产生病态也仅仅是步进到一半步数以后，亘+ i 天的条件 数开始变大的，最后变得没有规律产生这种现象得原因是在步进计算进行到 一半以后时，由于误差积累的原因，使得矩阵委+ i 天受到扰动所引起的 因此，用步进算法时要避免上述的问题，对深度方向离散时，采用等分取 节点，若发现在等分时仍不能改善病态的时候，则在深度方向等分基础上在步 进计算过程种采用行均衡的预处理方法 1 6 浙江大学硕士学位论文 第四章数值模拟 第四章数值模拟 下面我们通过在f o r t r a n 中计算的数值例子验证文中所提到的问题，及对出 现的问题给出修正的方法 本章主要按照矩阵( 0 + i 4 - 天) 的病态程度，采用不同的方法对其进行改善 即：当( q + i 4 a ) 的病态不是很严重的时候，我们通过改变深度方向的离散， 改善其病态性之后再求解。当病态性加重至不能通过前面的方法改善时，我们 通过在步进计算中对病态矩阵进行均衡法的处理，从而得到好的结果，获得稳 定的数值解本章中所取的精确解图形与参考文献【3 】，【3 5 】中的相同 ( 1 ) 当矩阵( q + i x a ) 的病态性不是很严重的时候 例1 ： 小一= 【o 1 6 7 , 啪，；冀箔 三= 1 0 ，d = 1 5 ，q = 4 0 ，, o l = 1 ，p z = 1 7 ，i = 2 0 0 ，2 = 1 0 0 ，3 = 1 0 0 ， 办( x ) ：l 一一褂，占：0 2 ，盯：2 。，f ：而1 传播所得的结果如下： 1 7 浙江大学硕士学位论文第四章数值模拟 d e p t h d e p t h f i l 4 1 例1 中l = 2 0 0 ，2 = 1 0 0 ，3 = l o o ，f = 上1 7 0 时的”( 厶z ) 我们取f = 上1 2 8 时的解作为精确解，对f = 去步进过程深度方向的离散进行 浙江大学硕士学位论文第四章数值模拟 改进，将l = 2 0 0 ，2 = 1 0 0 ，3 = 1 0 0 ，分别改进为：i = 1 0 0 ，2 = 5 0 ，3 = 2 5 0 ， 得到较好的结果，如下： d e 0 t h 图4 2 n i = 1 0 0 ，2 = 5 0 ，m = 2 5 0 时f = 击和f = 去时z ，( 厶z ) 的比较 浙江大学硕士学位论文第四章数值模拟 改变深度方向离散后矩阵( + i 忑) 的条件数如下： 籁 生 蠕 图4 3 改变深度方向的离散后( 囝+ i 天) 的条件数 产生这种现象的原因是f = 而1 的步进过程中的( 亘+ i 天) 的病态引起的， 为了更直观的表示矩阵( 亘+ i 天) 的条件数是如何变化的，从而更清楚的知道 病态是从什么地方开始的，我们将下面图4 4 中矩阵( 亘+ i 天) 的条件数小于 1 2 的部分进行放大，见图4 5 浙江大学硕士学位论文 第四章数值模拟 籁 芒 熊 藉 圣 媒 剖分数 图4 4 改变深度方向的离散起矩阵( 委+ i 天) 的条件数 剖分数 图4 5 图4 4 中纵轴小于2 2 的部分 浙江大学硕士学位论文 第四章数值模拟 经过大量的数值计算，我们发现即使是在深度采用等分离散，当步进的步 数过多时，步进算法还会变得不稳定，并且步进的步数越多，矩阵( q + i 石) 的条件数到后面就会变得越大但是当步进步数不是很大时，在某些深度方向 非等分离散的情况下，矩阵( 亘+ i 万) 病态也不会产生病态对于本例，在其 他条件不变的前提下，仅仅取相对小一点的步进步数f 2 丽1 ，矩阵( 囝+ i 灭) 不会产生病态，波的传播性态也不会改变。如下为l = 2 0 0 ，2 = 1 0 0 ，3 = 1 0 0 ， f 2 而1 时的条件数 糕 盘 熊 图4 6 例1 中l = 2 0 0 ，2 = 1 0 0 ，m = l o o ，f = 上1 4 0 时的条件数 这说明，步进算法的不稳定性是由步进方向和深度方向的误差积累造成 的，并且深度方向的不同离散方式，对误差积累的速度的影响是不同的 浙江大学硕士学位论文 第四章数值模拟 ( 2 ) 矩阵( 亘+ i x ) 的病态性加重时 当矩阵( q + i x a ) 的病态性加重时，改变深度方向的离散已经不能改进其 病态性时，我们在等分深度方向的基础上采用均衡法对矩阵( q + i x a ) 进行处 理，如： 例2 ： 小一书1 啪6 , ，z 渊x ) l = 1 0 ，d = 1 5 ，q = 4 0 ，肛= 1 ，岛= 1 7 ，l = 1 0 0 ，2 = 5 0 ，3 = 2 5 0 ， 办( x ) 小韶一删一0 2 产而1 对矩阵( 0 + i 4 灭) 的病态程度进行改进前，波的传播如下： 图4 7 例2 中矩阵( 耍+ i 瓜) 的病态程度改进前的u ( l ，z ) 的实部 浙江大学硕士学位论文 第四章数值模拟 兮 e 图4 8 例2 中矩阵( q + i x a ) 病态程度改进前的u ( l ，z ) 的虚部 经分析发现，原因还是由于矩阵( 西+ i x a ) 的病态引起的，我们给出矩 阵( 0 + i x - a ) 的病态程度进行改进前的条件数为了更直观的看出矩阵 ( 亘+ i 万) 的病态性是从什么地方开始的，我们将图4 9 中( 亘+ i x ) 的条 件数小于1 2 时的部分进行放大，如下图： 浙江大学硕士学位论文第四章数值模拟 籁 盎 礞 藏 跫 媾 图4 9 例2 中矩阵( 亘+ i x ) 的病态程度改进前的条件数 剖分数 4 1 0 图4 9 中纵轴小于1 2 时矩阵( 豆+ i 天) 的条件数 浙江大学硕士学位论文第四章数值模拟 由此我们可以看到：( q + i , f a ) 开始产生病态是从步进步数为2 0 0 0 到2 2 0 0 开始的，而非一开始就是病态的 在步进过程中采用均衡法对( 0 + i f a ) 进行处理，当正则化因子占= o 0 0 5 时得到较好的结果，如下为改进后( q + i 4 a ) 的条件数： 籁 圣 礞 图4 1 1 例2 中矩阵( o + i , z ) 病态程度进行改进后的条件数 糊酬郴) _ 们1 6 , 啪：三锚 三= 1 0 ，d = 1 5 ，d l = 4 0 ，岛= 1 ，p z = 1 7 ，n = 4 0 0 l = 1 0 0 ，2 = 5 0 ，3 = 2 5 0 ， 厅( x ) ：l 一钟- d 生- q ， ，占：。2 ，f ：函1 为精确解，则对例2 所给出的情况，将进行改善后的波的传播与精确解作如下 的比较： 浙江大学硕士学位论文第四章数值模拟 d e p t h d e p t h 图4 1 2 例2 中( q + i 一f a ) 病态程度改进后的u ( l ，z ) 与精确解的比较 2 7 浙江大学硕士学位论文第五章总结 第五章总结 本文以带有弯曲的中间界面的声波导为背景，对步进算法的稳定性进行了 分析本文指出影响步进算法稳定性的主要因素有：步进的步数及深度方向的 离散当前者都比较小时，步进算法是数值稳定的，但是，步进的步数变大， 或者深度方向的离散不恰当时，就会破坏这种稳定性。产生这种现象的原因是 误差积累造成的 本文分另0 从步进的步数对步进算法稳定性的影响，深度方向的离散步进算 法稳定性的影响，这些方面进行说明，指出造成这种现象的原因是，深度方向 的不同离散对误差积累的影响是不同的因此，我们应采用等分的方式进行离 散，以使得深度方向的离散对误差的积累影响最小对等分离散时仍然出现不 稳定的传播给出了正则化方法种的一种：行均衡的改进方法并通过数值模拟， 说明步进算法的不稳定性是存在的，所给出的原因是正确的，对不稳定的情况 进行改善的方法是有效的，即改进后能使得步进算法保持快速稳定的计算，并 能得到真实的波的传播。 附注 附注 当h ( r ) - - o 时，我们定义口( 瓦三) = ，l i m ，口( 曼，三) ，( - ，三) = ，l i m 膏( 叠，三) ， y ( i ，三) = 磐y ( 曼，三) ，则第一层中： 口( 五三) = l e 一，器 即一一器p 。锪 啊铲 器i 锱卅+ 一乍鬻， + 手曙仁中，需 第二层中： 俐= 丽( d - i ) 2 p 跏牡m 。 俐= 筹产p 龉郴弘批。 俐= p 辫郴脚们。 器r 一一黑+ 3 尚( 叫咿2 d ：+ ：2 ) 卜 d h ( i l 。 三f 里：篓堕：p 器邶2 d - z - h ( y t ) ) + 4 d 一向( i ) 办( i ) 。 器。t 纠+ 糯卜固州啦d ：+ ，】) 浙江大学硕士学位论文附注 第三层中： i 口( 瓦芝) = e - 川叫 l ( i ，三) = o 卜铲；器十c 纠+ 器卜审) + 捣 l 华一抖叫咿+ 丢帮p 一- ) 更进一步，如果办( i ) = o 且有矿( i ) = o ，则： 雠晦箭 卜铲面( d - 丽1 ) 2 ( 五三) = o 旧一2 ( 纠+ 型等业龋 口( i ，2 ) = 1 ( 墨三) = o 厂( i ，2 ) = 茁2 ( i ，三) + 1 4 d 硝4 ( i ) 浙江大学硕士学位论文 参考文献 参考文献 【1 】l uy y ；h u a n gj a n dm c l a u g h l i nj r ，l o c a lo r t h o g o n a lt r a n s f o r m a t i o na n d o n e w a y m e t h o d sf o ra c o u s t i c s w a v e g u i d e s w a v em o t i o n , 2 0 0 1 ， 3 4 ( 2 ) ：19 3 2 0 7 【2 】l uyy ，o n e w a yl a r g er a n g es t e pm e t h o d sf o rh e l m h o l t zw a v e g u i d e s j c o m p u tp h y s ，19 9 9 ，15 2 ：2 31 2 5 0 【3 】z h uj x ；l uy y = ，l a r g er a n g es t e pm e t h o df o ra c o u s t i cw a v e g u i d ew i t ht w o l a y e rm e d i a p r o g r e s si nn a t u r a ls c i e n c e ，2 0 0 2 ，12 ( 11 ) ：13 18 【4 】l uy y ，m c l a u g h l i nj r ，t h er i c c a t im e t h o df o rt h eh e l m h o l t ze q u a t i o n j a c o u s ts o ca m ，19 9 6 ，10 0 ( 3 ) ：14 3 2 - 14 4 6 【5 】5 l uy y ，l a r g er a n g es t e pm e t h o df o rh e l m h o l t zw a v e g u i d e ，i n ：j a d c s a n t o ( e d ) ， m a t h e m a t i c a la n dn u m e r i c a la s p e c t so fw a v ep r o p a g a t i o n ，s i a m ，p h i l a d e l p h i a ， p a ，19 9 8 ，p p 6 2 6 - 6 2 8 【6 】l uy = y ，e x a c to n e w a ym e t h o d sf o ra c o u s t i cw a v e g u i d e s m a t h e m a t i c sa n d c o m p u t e r si ns i m u l a t i o n ，19 9 9 ，5 0 ：3 7 7 - 3 9 1 【7 】李荣华，冯果忱，偏微分方程数值解法( 第三版) ，高等教育出版社，1 9 9 6 年1 2 月 【8 】董光昌，陈仲慈，汤国桢，数学物理方程( 第一版) ，浙江大学出版社 【9 】潘祖梁，非线性问题的数学方法及其应用( 第一版) ，浙江大学出版社，1 9 9 8 年 【l o 】徐树方，矩阵计算的理论与方法( 第一版) ，北京大学出版社，1 9 9 5 年 【l l 】高等代数( 第二版) ，高等教育出版社，1 9 8 8 年 12 1g e o r g ev f r i s k , o c e a na n ds e a b e da c o u s t i c s 【1 3 】g e n eh g o l u b & c h a r l e sf v a nl o a n ，m a t r i xc o m p u t a t i o n s ( t h i r de d i t i o n ) ，t h e 3 l 浙江大学硕士学位论文参考文献 j o h nh o p k i n su n i v e r s i t yp r e s s ，19 9 6 【1 4 】l a r a h a m s s o n ，h 一o k r e i s s ，n u m e r i c a ls o l u t i o no ft h ec o u p l e dm o d ee q u a t i o n s i nd u c ta c o u s t i c s ，j c o m p u tp h y s ，1 9 9 4 ，1 1 1 ：l 一1 4 【15 】c a b o y l e s ，c o u p l e dm o d es o l u t i o nf o rac y l i n d r i c a l l ys y m m e t r i co c e a n i c w a v e g u i d ew i t har a n g ea n dd e p t hd e p e n d e n tr e f r a c t i v ei n d e xa n dat i m ev a r y i n g r o u g hs e as u r f a c e j a c o u s ts o ca m ，1 9 8 3 ，7 3 ：8 0 0 8 0 5 【l6 】m d c o l l i n s ，a p p l i c a t i o n sa n dt i m ed o m a i ns o l u t i o n so fh i g h e r - o r d e rp a r a b o l i c e q u a t i o n si nu n d e r w a t e ra c o u s t i c s j a c o u s rs o ca m ，19 8 9 ，8 6 ( 4 ) ：10 9 7 110 2 【17 】t w d a w s o n ，j a f a w c e t t , ab o u n d a r yi n t e g r a le q u a t i o nm e t h o df o ra c o u s t i c s c a t t e r i n gi naw a v e g u i d ew i t hn o n - p l a n a rs u r f a c e s j a c o u s ts o ca m ，19 9 0 ， 8 7 ( 3 ) ：1l1 0 - 11 2 5 【l8 】l d i e c i ，n u m e r i c a li n t e g r a t i o no ft h ed i f f e r e n t i a lr i c c a t ie q u a t i o na n ds o m e r e l a t e di & s u e s s i a mj n u m e r a n a l ，1 9 9 2 ，2 9 ( 3 ) ：7 8 1 - 8 1 5 【19 l f i s h m a n ，a k g a u t e s e n ，z s u n ，u n i f o r mh i g h f r e q u e n c ya p p r o x i m a t i o n so f t h es q u a r er o