（计算数学专业论文）dtn重构算法的有效性分析.pdf

摘要 本文主要研究了以声波导为背景的h e l m h o l t z 方程，用d t n 重构算法步进计 算波的传播性态的有效性问题。对于含有弯曲内部界面的声波导，首先需要进 行坐标变换，将波导中弯曲的界面拉直，把h e l m h o l t z 方程转化成为一个可以用 步进方法求解的形式。再用d i r i c h l e t - t o - n e u m a n n 映射，将边值问题转化为初值 问题。 在数值实施d t n 重构算法时，涉及到微分、积分、求逆、特征值等问题，因 此会产生计算误差，需要对其有效性进行分析。本文在数值模拟中发现，对于 缓慢变化的内部界面，得到的结果较稳定。但用内部界面模拟海底对称陡坡时， 用r a y l e i g h 迭代方法求解特征问题时会发生不稳定现象，使得到的波形振幅过 小。使用重置方法可以解决这个问题，虽然计算时间较长，但计算过程中产生 的误差不大，所得声波形状变化较火是由于内界面函数变化过快引起的。另外， 当在界面函数的极点附近进行坐标变换时，迭代过程中某被积函数可能出现分 母为零的情况。我们通过研究被积函数的性质，利用其对称性，可以将其转化 为分段积分的形式，其中包含极点牙的区域恰好积分为零。这样，在保证计算正 确性的基础上，解决了问题。 关键词：h e l m h o l t z 方程，坐标变换，算子离散，d t n 重构算法，内部界面， 误差分析 a b s t r a c t t h i sa r t i c l ed i s c u s s e st h ev a l i d i t yo ft h en u m e r i c a lm 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 n ( d t n ) m a p f o ra c o u s t i cw a v e g u i d e sw i t hac u r v e di n t e r f a c ei nt w o l a y e rm e d i a ，w ep r o p o s eac o o r d i n a t et r a n s f o r m a t i o nt of l a t t e nt h ei n t e r f a c ea n d c h a n g et h eh e l m h o l t ze q u a t i o ni n t oas o l v a b l ef o r m t h e nw eu s eal a r g er a n g e s t e pm e t h o dt od i s c r e t i z et h er a n g ev a r i a b l ea n dr e d u c et h eb o u n d a r yv a l u ep r o b - l e mt oa ni n i t i a lv a l u ep r o b l e mb yd t nm a p i nn u m e r i c a ls i m u l a t i o n ，t h ep r o c e s so fd i f f e r e n t i a l ，i n t e g r a l ，i n v e r s ea n d e i g e n v a l u ep r o b l e m sm a yi n d u c ec a l c u l a t i o ne r r o r ，s oi ti sn e c e s s a r yt oe v a l u a t e t h ev a l i d i t y w ef o u n dt h a tf o rt h ew a v e g u i d ew i t hi n t e r n a li n t e r f a c e sw h i c h c h a n g es l o w l y , t h er e s u l t sa r es t a b l e b u ti ft h ei n t e r f a c es i m u l a t e sas y m m e t r i cs t e e ph i l li nt h es e a b e d ，t h ea m p l i t u d ei st o os m a l ld u et ot h ei n s t a b i l i t yo f r a y l e i g hi t e r a t i o n w bc a nu s et h er e s e tm e t h o dt os o l v et h i sp r o b l e m t h o u g h i tm a yt a k em o r et i m e o u rr e s e a r c hd e m o n s t r a t e st h a tt h ee r r o ri nn u m e r i c a l s i m u l a t i o ni sn o ts e r i o u s ，t h er e a s o no fv i b r a t i o no fw a v e f o r mi st h a tt h ei n t e r - f a c ec h a n g e st o oq u i c k l y i na d d i t i o n ，w h e nw et r a n s f o r mt h ec o o r d i n a t en e a r t h ep o l e ，t h ed e n o m i n a t o ro ft h ei n t e g r a n di nt h ei t e r a t i v ep r o c e s so ft h ec o o t - d i n a t et r a n s f o r m a t i o nm a yb e c o m ez e r o a f t e rr e s e a r c h i n gt h ec h a r a c t e ro ft h e i n t e r g r a n d ，w em a k eu s eo fs y m m e t r yo ft h ei n t e r f a c ef u n c t i o nt os e p a r a t et h e i n t e g r a lr e g i o ni n t os e v e r a lp a r t s ，a n dt h ei n t e g r a lo ft h ep a r tw h i c hc o n t m u st h e p o l ei se x a c tz e r o s ot h ep r o b l e mi ss o l v e do nt h eb a s i so fa c c u r a c y k e y w o r d s ：h e l m h o l t ze q u a t i o n ，c o o r d i n a t et r a n s f o r m ，o p e r a t o rd i s c r e t i z a t i o n ， d t nr e - f o r m u l a t i o n ，i n t e r f a c e ，e r r o ra n a l y s i s 浙江大学研究生学位论文独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作及取得的 研究成果。除了文中特别加以标注和致谢的地方外，论文中不包含其他人已经发 表或撰写过的研究成果，也不包含为获得迸婆盘堂或其他教育机构的学位或 证书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文 中作了明确的说明并表示谢意。 学位论文作者签名： 签字日期：年月日 学位论文版权使用授权书 本学位论文作者完全了解迸鎏盘堂有权保留并向国家有关部门或机 构送交本论文的复印件和磁盘，允许论文被查阅和借阅。本人授权澎婆盘堂 可以将学位论文的全部或部分内容编入有关数据库进行检索和传播，可以采用影 印、缩印或扫描等复制手段保存、汇编学位论文。 ( 保密的学位论文在解密后适用本授权书) 学位论文作者签名：导师签名： 签字口期：年 月e l 签字日期： 年月 日 致谢 值此论文完成之际，谨在此向多年来给予我关心和帮助的老师、同学、朋 友和家人表示衷心的感谢! 首先，我要对导师朱建新教授致以最崇高的敬意。感谢朱老师这几年在学 业上的谆谆教导，以及在生活上的关心照顾。朱老师对科学的钻研精神和严谨 的治学态度将是我一生的榜样。 其次，我要对师姐陈芝花表示诚挚的谢意，感谢她在论文撰写中给我的帮 助。也感谢和我一起参加讨论班的宋仁成，张学仓，沈浙奇，以及已经毕业的师 兄师姐，还有新来的师弟师妹们，在和你们的交流中让我获益匪浅。还有感谢 潘一力同学，和我一起度过两年生活，共同学习，共同成长。 最后，感谢我的父母和家人这么多年对我的关心，并支持我完成学业。 谨把本文献给所有那些曾经帮助过我的人们! 第一章引言 1 1 应用背景 在过去的几十年间，世界形势的快速变革开创了建模与仿真工作国际合作 的新途径。电子通信的发展极大地促进了信息之间的交流与新技术的诞生。其 中，水声学在描述水下概貌，借助海洋波导传递信息，测量海洋特性等方面起 着突出作用。历史上，海军曾支持声呐技术人员开创和发展水声建模工作来改 进声呐系统设计与评估效果。由于军事安全上的限制，很多研究成果鲜为人知， 但仍有大量文献公开发表。这些文献大部分涉及到以海洋建模为背景的数值程 序的开发与优化，并逐渐成为一个被称为计算海洋声学的新兴学科。随着建模 技术的成熟和向公众领域的移植，大量成果应用于民营工业，对地球物理探索 和声呐技术等领域有所促进。近来，声学海洋学者已经使用水声模型作为反演 技术( 从海洋物理特性的直接测量中获取海洋的有用信息) 的辅助工具，用于 获得大面积海洋的概貌，以及i i , 在n 海洋的长周期变化 2 7 】。 1 2基础理论 波动方程是从较基本的状态方程、连续方程和运动方程导出的，这在很多 物理基础教科书 1 l 上都有严格的推导，本文不再重复，直接从波动方程出发进 行数学研究。 波动方程的具体形式随基本假设条件和具体应用场合的不同而不同。对于 大多数的应用，通常都采用简化的与时间有关的双曲型二阶线性偏微分方程： v 2 垂= 石1 等 ( 1 1 ) 其中v 2 = ( a 2 l a x 2 ) + ( 0 2 l o y 2 ) + ( 0 2 1 0 2 2 ) 是拉普拉斯算子，圣是势函数，c 是声 速，t 是时间。 为了获得与时间无关的h e l m h o l t z 方程，又引入了谐和解作进一步简化。即 假定势函数圣的谐和解是 垂= c e 一础( 1 2 ) 第一章引言 2 其中是与时间无关的势函数，u 是声源频率( 2 7 r f ) 。于是波动方程( 1 1 ) 就简化成 v 2 - t - 舻咖= 0 其中托= u c = 2 r a 是波数，入是波长。 本文考虑声波在水下传播的情形，h e l m h o l t z 方程可以写成这种形式： ( 1 3 ) t 正- 4 - u ：4 - k 2 ( ，z ) u = 0 ，0 z + o 。，一o 。 z + o 。( 1 4 ) 其中z 方向是垂直向下的，z = 0 表示海平面，z 方向为横向传播方向，u 为声压经 过傅立叶变换之后的形式。 对上述方程( 1 4 ) ，由于k 在水平x 方向上是弱衰减的，故可以把z 方向分成三 部分，即z l 。当z l 时，假设仡为常数，可由分离 变量方法得到方程的精确解。本文主要考虑0 z l 区域中，h e l m h o l t z 方程 的求解问题。 因为在0 z l 区域，k 并不是常数，而与水平距离z 和深度z 等因素有关， 故无法直接给出解析解，需要数值计算进行求解。 1 3发展现状 在数值求解h e l m h o l t z 方程的过程中，首先遇到的是边界问题。声波在无限 深度是无法进行数值计算的，所以必须引入适当的吸收边界条件把计算空间截 断，并且应保证在截断边界处只有向外传播的波而没有反射波。从吸收边界条 件的研究历史看，大致分为两个阶段【3 8 】。 第一阶段是2 0 世纪7 0 8 0 年代，共提出了四大类吸收边界条件：基 于s o m m e r f i e l d 辐射条件的b a y l i s s - t u r k e l 吸收边界条件；基于单向波动方程 的e n g q u i s t m a j d a 吸收边界条件；利用插值技术的廖氏吸收边界条件；以及梅 一方超吸收边界条件。这些吸收边界条件的反射系数大约在0 5 一5 之间。 第二阶段是2 0 世纪9 0 年代，由b e r e n g e r 提出了完全匹配层( p e r f e c t l ym a t c h e d l a y e r ，简写为p m l ) 的理论模型f 2 】它在求解区域边界上的反射系数只有上述 各种吸收边界条件的1 3 0 0 0 ，可以吸收各个频率，各个方向的波，被认为是臼前 最好的吸收边界条件之一。 其次，关于波的传播计算方法，主要有以下几种方法3 8 1 ： 第一章引言 3 一是利用解析法。根据分离变量法来求解偏微分方程，将解表示为己知函 数的形式，可以得到精确的结果，并可作为近似解和数值解的检验标准。但是 这种方法使用范围非常有限，仅能求解很少一部分方程，不具有普遍性。 二是利用常规数值方法，如有限差分法和有限元法。由于限制少，这是目 前应用最j “泛的数值方法。但这些数值方法都属于全局方法，将导致高阶系统， 会占用大量的c p u 内存及计算时间。若遇上非稀疏矩阵，更在存储空间上碰上 了很大困难。目前最新发展动向是研究高效的并行数值算法。 三是利用模式匹配方法【3 ，4 】：主要有正规模的方法和耦合模的方法。当折 射率与纵向无关的时候，使用正规模方法：当折射率随纵向变化的时候，使用耦 合模方法，这是迄今为止比较成熟的理论方法。它利用一组正交完备组的模式 展开，将问题的解由相应的一组耦合系数来表示。但是这种方法公式推导繁琐， 且随着模式的增加，计算量比较大。 四是基于光传播的单向光束法和双向光束法【5 ，6 ，7 ，8 ，9 ，i 0 ，i i 】。这种方法 来源于近似h e l m h o l t z 方程的o n e - w a y 方程也u = i 、理+ 仡3 礼2 “。o n e - w a y 近似是 使用最广泛的一种抛物方程近似法。若将z 看作时间变量，则抛物方程可以当作 初值问题来求解。波导出射端的初始算子可以通过逼近平方根算子取得，一般 采用有理p a j c l d 逼近。 再者，如果存在弯曲边界或弯曲的内部边界f 3 9 1 ，在进行波的传播计算时， 传统的做法就是使用“阶梯”逼近f 1 2 ，1 3 ，1 4 1 。这需要用很小的步长进行模拟， 影响运算速度。因此，海洋声学的研究者们提出了一些数值方法来避免粗糙 的”阶梯”逼近。文献【1 5 ，1 6 使用了一种全局坐标变换的方法，它是一种保角变 换，需要很大的计算量。另外一种方法，局部坐标变换，也出现在了各种各样的 实际应用之中f 1 7 ，9 ，1 8 ，1 9 ，2 0 】。其中，局部正交变换在拉直弯曲边界或者界面 的同时，变换后的h e l m h o l t z 方程中不存在关于传播方向和纵方向变量交叉导数 项，这使得可以对改进的h e l m h o l t z 方程使用大步长方法求解【2 0 】。 1 4 算法概述 本文主要研究的是基于狄利克莱一纽曼映射( d i r i c h l e t t o - n e u m a n nm a p ) 的 d t n 重构算法在含有弯曲内部界面的声波导中的应用 2 0 】。这种算法是1 9 9 3 年由 f i s h m a n 首先提出的 7 】。1 9 9 6 年l uy a y a n 和m c l a u g h l i n 提出了基t - d i r i c h l e t - t o - n e u m a n n 映射的m c c a t i 方法，将边值问题转化为初值问题【1 1 】。这种方法所需内 第一章引言 4 存空间非常小，是一种非常实用的方法。 针对如下问题： u 2 z + u 2 名+ k 2 ( z ，z ) u = 0 考虑声波在0 z l 区域传播的情形，深度方向用z = d 1 截断( 注：更好的 方法是使用完美匹配层( p m l ) 进行截断，本文为了方便计算，只使用这种简单 的截断) 。函数z = 危( z ) 把这个区域分为两层，当0 z 危( z ) 时，密度为p 1 ， 当h ( x ) 名 1 ) 时，密度为p 2 。 对于含有弯曲内部界面的声波导，这里使用局部正交变换将弯曲界 面z = ( z ) 拉直 2 0 】。令 童= f ( x ，z ) ，三= g ( x ，z ) ， 其中( z ，z ) 为原来波导所在的坐标系，称为物理坐标系，反映了波导的实际物 理形状，( 岔，2 ) 为变换之后得到的新坐标系，称为计算坐标系。这样，波导的 弯曲界面被平滑化，相当于把区域【( 。，名) i 一 z o 。，0 z 危 ) ) 映射 到 ( 圣，三) l 一。 金 o o ，0 三1 ) ，【( z ，z ) l o 。 z ，h ( x ) z d i 映 射到 ( 圣，三) i 一 2 o 。，1 三d i 。简便起见，这里选取g 为线性插值函 数，可以通过正交性条件厶如+ 五g z = 0 确定。 令札( z ，z ) = w ( z ，z ) y ( z ，z ) ，方程( 1 4 ) 化为 其中选取 + a ( 仝，三) k 童+ 卢( 岔，2 ) k + 7 ( 圣，乏) y = 0( 1 5 ) w ( x ，z ) = w ( x ，名) = w ( x ，名) = 0 乏 1 器生h ( x 掣萨d ，l 三 。 7 ( 盒) + ) 【 一九( z ) 】2 一、一 九( 2 ) d h ( x + ) 7 ( z ) 系数口( 圣，2 ) ，( 金，2 ) ，y ( ，三) 在附录a 中具体给出。 d 0 ，i = 1 ，2 ，一1 选取对角矩阵d = d i a g ( d 1 ，d 2 ，d ) ，使得d a = s ，s 为对称矩阵。 可以推出a ：d 一 ( d 一 s d 一) d 其中j d 一= ( d ) 一1 ，d = d i a g ( x 丽l ， 瓦，、石) m a 。d 一；s d 一 知矩阵a 存在n 个加权线性无关的特征向量y = 【，k ，】， 满足 a = d i a g ( a 1 ，入2 ，入) a 巧= 巧， j i = 1 ，2 ， 将区间【o ，纠关于圣离散为 0 = 奎o 2 1 童8 _ 1 圣8 仝s + 1 圣m = l 定义算子q 和y ，满足 = q ( 圣) kv ( l ，) = y ( 圣) y ( 2 ，) ( 1 6 ) 代入方程( 1 5 ) ，可得 y 7 = 一y q ( 动，q = 一q 2 一( q 霹+ p 侥+ 一y ) 边界条件为 q ( l ) = i 、q ( l ，2 ) 秀+ 7 ( l ，乏) ，y ( l ) = i i 为单位矩阵。 动 z ，il7 + 以 动 z ，- 8+ 谚 ，- 、 a 子 算将 第一章引言 6 令算子q ，y 在筑分别为q 。，k ，通过在区间( 以一1 ，筑) 的中点圣s - 1 2 = ( 圣。一1 + 筑) 2 处满足的方程和界面条件可以推出饥，k 和仉一1 ，k l 之间的关系。这样， 就可以从q l ，y l 推出q o ，y o 。当给定在圣= 0 处的初始值v ( o ，) 时，由y ( l ，) = y ( o ) v ( o ，) 可得出波传播到圣= l 处的值。 具体步骤如下：在区间( 丸一1 ，豇) ，令7 = 以一岔，1 ，已知q 。，k ，在圣。一1 2 得 到对应的矩阵也。，并得到特征值矩阵人和特征向量矩阵y 。设： 嘭砌= ( 圣外1 2 ) ，曙仳伽= k ( 丸一1 2 ) ，j = 1 ，2 ， 根据特征向量之间的加权正交性，可以得到矩阵w = ( w t j ) n n ，z = ( 锄) x n 使 得 n 嘭删= w , j v , ( 1 s e l o ) 嘭一) = z , j y , 似田 i = 1l = 1 将算子仉，k 表示为对应矩阵国。，磁，可通过如下公式求得国。一l ，e 一1 国= w q 。z p = ( q + i v a ) 一1 ( - q + i a ) _ 一 广- r ：e 打瓶p e i r 瓶 国。一l = i c x ( i 一只) ( ，+ r ) 一1 览一1 =e z ( ，+ p ) e 打e x ( i + r ) 一1 在上述运算过程中，可以用nxn 矩阵来逼近算子q ，y ，也可以用截断的 特征函数展开方法，即取出a 矩阵最大的n 个特征值所对应的特征向量，得到一 个较小的nxn 矩阵来代替，佗的下限通常为传播模的个数。由于n n ，所以这 种方法更加简便有效。 在数值计算中，从数学问题转化为数值问题会产生截断误差，由于计算机 的字长有限，运算过程中的舍入会导致舍入误差。在数值实施d t n 重构算法时， 涉及到积分，求逆，特征问题等计算，因此产生计算误差。在实际操作中，我们 发现对于不同的内部界面，横向步长，坐标变换，计算精度等诸多因素都会影响 算法的有效性。本文将对此进行分析，根据具体情况的不同，对一些细节方面 的算法进行改进，使d t n 重构算法得到进一步完善。 第二章有效性分析 在求解实际问题中，从建立数学模型，将数学问题转化为数值问题，到设 计算法，在计算机上运行得到结果，每一步都可能产生误差。数学模型是通过 对被描述的实际问题进行抽象，忽略次要因素简化而得到的。一般说来它是原 始问题的近似，其误差称为模型误差。在数学模型中还会涉及到一些观测得到 的物理量，比如温度，质量，密度，电压等等，这些参量也包含误差，称为观察 误差。上述两种误差不属于数值计算的研究范围。数值分析主要研究从数学问 题转化为数值问题的算法时所产生的误差，称为截断误差或者方法误差。另外， 由于计算机的字长有限，原始数据的保存和运算过程中的舍入也可能产生误差， 统称为舍入误差。有些误差可能会互相抵消，有些误差可能无法抵消，甚至会 增长。我们关心计算结果的可靠性，但对于大规模数值计算，在理论上对误差 进行定量估计有很大难度f 2 8 ，2 9 1 。 在实际操作中，数值计算的有效性应结合具体数值方法来进行分析。对于 不同的算法也有各种各样分析方法。通常情况下，一个算法是否可靠与误差传 播密切相关。如果输入数据有扰动( 即误差) ，而计算过程中误差不增长，就称 此算法是数值稳定的，否则是数值不稳定的。对于任何输入数据都是稳定的算 法，称为无条件稳定。而对某些数据稳定，某些数据不稳定的算法，称为条件稳 定。通常对条件稳定算法可通过重新设计算法或增加计算机的精度，使算法的 可靠性得到改善。本文针对步进算法在模型本身，数值近似，计算方法等方面 进行有效性分析。 2 1 不同的界面函数 考虑声波在水下传播的情形，内部界面函数 ( z ) 相当于海水与海底之间 的分界面。与大陆板块一样，海底地形地貌也具有多样性的特点。本文选 取了几个典型的函数来模拟海底的形状。在计算过程中，均选取金= 0 处算 子口( 圣，乏) 谚+ p ( 圣，三) 如+ ，y ( 圣，三) 离散成的矩阵的最大特征值所对应的特征函 数v ( o ，乏) 乘以w 为波的初始值，取截断特征函数的个数为r t ，纵向第一层0 乏1 等分为l 份，第二层1 2 d 作2 等分，第三层d 三d 1 作3 等分， 第二章 有效性分析 8 横向步进时的步长为r 。首先考虑水平的海底界面，这相当于取函数危( z ) 等于常 值函数。 例1 h ( x ) = 1 ，l = 1 0 ，d 1 = 4 0 ， 1 6 ，乞2 = 1 1 2 ，g l = 1 0 0 ， n 2 = 5 0 ， 算传播到圣= l 时的结果如图2 1 所示： ( a ) u ( 厶名) 的实部 d = 1 5 ， p 1 = 1 ，, 0 2 = 1 7 ， k 1 = n 3 = 2 5 0 ，佗= 4 0 ，r = 1 2 5 6 ，计 ( b ) u ( l ，：) 的虚部 图2 1 ：例1 当h ( x ) = l 串 l u ( l ，z ) 的波形 为了验证结果的准确性，本文对步长7 - 分别取1 4 ，1 8 ，1 1 6 ，1 3 2 ，1 6 4 ，1 1 2 8 ， 并以7 - = 1 2 5 6 为标准，比较经上述步长进行计算时，所得结果u ( l ，z ) 与丁= 1 2 5 6 时所得结果的相对误差。简记步长为7 - 时计算结果为乱r ，定义相对误差函 数： 。 e ( 7 i ) ：坚；掣 一 l iu a 2 5 6l | 通过图2 2 表明，随着步长的减小，相对误差逐渐缩小，u ，是收敛的。在后面的例 子中均可采用这种方法来证明计算结果的收敛性。 下面考虑海底为楔形的情况。我们把海底形状抽象为线性函数h ( x ) = k x + b ，根据水声建模与仿真【2 7 】中对海底状况的描述，取k = 一o 0 5 ，b = 1 ( 图2 3 ) ， 然后针对下面的例子进行计算： 第二章有效性分析 9 图2 2 ：乱，( l ，z ) 的相对误差 r a d g e z 图2 3 ：含内界面h ( x ) = 1 0 0 5 x 的声波导 一一。粤“mo n 茸苞o q 第二章有效性分析 1 0 例2 h ( x ) = - 0 0 5 x + 1 ，l = 1 0 ，d 1 = 4 0 ，d = 1 5 ，p l = 1 ，j d 2 = 1 7 ，k 1 = 1 6 ，k 2 = 1 1 2 ，1 = 1 0 0 ，n 2 = 5 0 ，n 3 = 2 5 0 ，佗= 4 0 ，r = 1 2 5 6 ，计算传播到圣= l 时的结果如图2 4 所示： ( a ) u ( l ，：) 的实部( b ) u ( l ，：) 的虚部 图2 4 ：例2 当危扛) = 1 0 0 5 x 时“( 厶名) 的波形 可以看出，对于楔形海底界面d t n 重构算法也很稳定。再考虑海底变化时 非线性的情况，比如h ( x ) = 1 - a r c t a n ( x ) 8 ，相对于传播距离来说， ( z ) 变换非 常缓慢( 图2 5 a ) ，但局部放大之后，会发现这种函数在前半部分变化较大，后 半部分几乎是平的( 图2 5 b ) 。通过例3 可以知道，这种海底界面d t n 重构算法也 是很有效的。 第章有效性分析 1 1 ( a ) 完整的声波导( b ) 局部放大的声波导 图2 5 ：含内界面h ( x ) = 1 - a r c t a n ( x ) 8 的声波导 侈9 3 h ( x ) = 1 一a r c t a n ( x ) 8 ，l = 1 0 ，d 1 = 4 0 ，d = 1 5 ，p l = 1 ，以= 1 7 ，k 1 = 1 6 ，仡2 = 1 1 2 ，1 = 1 0 0 ，n 2 = 5 0 ，n 3 = 2 5 0 ，扎= 4 0 ，7 = 1 2 5 6 ，计算传播到童= l 时的结果如图2 6 所示： ( a ) u ( 三，z ) 的实部( b ) 缸( 厶2 ) 的虚部 图2 6 ：例3 当忽0 ) = 1 - a r c t a n ( x ) 8 时u ( l ，名) 的波形 可以发现，对于非对称缓慢变化的界面，步进算法计算效果很好。下面考 第一章有效性分析 1 2 虑对于对称的情形，比如海底呈山状，用函数危 ) = 1 6 e 一口( 芒一) 2 来模拟对称 的山坡，图2 7 a 中= 0 2 ，盯= 2 0 ，表示平缓的山坡，图2 7 b 中= 0 2 ，口= 6 0 0 ， 表示较陡的山坡。 o2468 o r w ， ( 8 ) = 0 2 ，盯= 2 0 、 ( b ) = 0 2 ，口= 6 0 0 图2 7 ：含内界面九( z ) = l e e 一仃( 量一 ) 2 的声波导 例4 对于图2 仡中的声波导，h ( x ) = 1 一e e 一矿( 詈一丢) ， e = 0 2 ，仃= 2 0 ，l = 1 0 ，d 1 = 4 0 ，d = 1 5 ，p l = 1 ，p 2 = 1 7 ，k 1 1 6 ，k 2 = 1 1 2 ，1 = 1 0 0 ， n 2 = 5 0 ，n 3 = 2 5 0 ，佗= 4 0 ，下= 1 2 5 6 ，计算传播到金= l 时的结果 如图2 8 所示： 第二章 有效性分析1 3 ( a ) u ( l ，z ) 的实部 ( b ) t ( 厶z ) 的虚部 图2 8 ：例4 当危( z ) = 1 一e 一玎( 量一 ) 2 ，= 0 2 ， 盯= 2 0 时u ( l ，z ) 的波r d 例5 对于图2 7 6 中的声波导，h ( x ) = l e e 一矿( 暑一 ) ，= 0 2 ，口= 6 0 0 ，l ： 1 0 ，d a = 4 0 ，d = 1 5 ，p l = 1 ，晚= 1 7 ，k l = 1 6 ，物= 1 1 2 ，n 1 = 1 0 0 ， n 2 = 5 0 ，n 3 = 2 5 0 ，礼= 4 0 ，r = 1 2 5 6 ，计算传播到2 = l 时的结果 如图2 9 所示： ( a ) u ( 厶z ) 的实部 ( b ) t ( 厶z ) 的虚部 图2 9 ：例5 当允( z ) = 1 一g e 一口( 量一言) n ，= 0 2 ， 盯= 6 0 0 时u ( l ，z ) 的波形 第章有效性分析 1 4 从波形图可以看出，例4 的计算结果较稳定，而例5 的计算结果与我们之前 算 h 的波形明显不同，振幅非常小。再计算例5 在声波传播之前，童= 0 时的波 形： ( a ) u ( o ，名) 的实部 图2 1 0 ：例5 当 0 ) = 1 一e e a ( 考一 尸， = 0 2 ， 仃= 6 0 0 时u ( o ，z ) 的波形 由图2 1 0 和图2 9 对比可知，例5 中声波在传播之后，波形和振幅变化较大， 这一方面与内部界面较陡有关，另一方面也可能在模拟计算中存在某些问题， 导致结果不稳定。 在步进过程中，每一步都需要涉及到特征问题。因为在传播方向声波导性 质变化缓慢，对应算子离散成的矩阵变化不大，所以通常使用r a y l e i g h 迭代方法 来求解特征值和特征向量，以节约时间。r a y l e i g h 迭代的过程主要为 4 0 】： 1 取上一步得到的特征向量u ，规范化后作为初始值u o ，计算初始特征 值a o = u $ s u o ，s 为需求解特征问题的矩阵； 2 求解( 入知，一s ) u k + l = 乱七， k = 0 ，1 ，并规范化u 七+ l = 横孙； 3 a 知+ 1 = u ：+ 1 s u k + l ； 4 若0 丸+ 1 一a 七l l i i 入七l i ，为控制精度，或者达到一定循环次数，则退 出循环，否则令k = k + l ，跳到步骤2 继续迭代。 第_ 章有效性分析 1 5 当盯= 6 0 0 时，经逐步调试程序发现，在使用r a y l e i g h 迭代时发生了不稳定 现象，特征值的求解出现异常，如图2 1 1 。 ( a ) 在圣= 6 1 2 5 的前4 0 个较大特征值( b ) 在圣= 3 8 7 5 的前4 0 个较大特征值 图2 1 1 ：例5 在岔= 6 1 2 5 和圣= 3 8 7 5 处特征值的对比 正常情况下，在步进中矩阵相应的特征值应呈稳定变化状态，而且对于对 称的声波导，在关于畲= 寺对称的横向坐标的特征值也应一样( 后文第三节中 有相应描述，如图2 1 7 ) 。但从图2 1 1 中可以看出，在金= 6 1 2 5 和金= 3 8 7 5 处特 征值相差很大。从圣= 6 1 2 5 步进计算到圣= 3 8 7 5 的过程中可以发现，这是由于 使用r a y l e i g h 迭代时产生计算误差累积而引起的。为了改善这种状况，本文采 用m a t l a b 软件自带的函数，对每一步的算子离散后形成的矩阵分别求解特征值 和特征向量，这称为重置方法。虽然在速度上有所影响，但可以得出正确结果， 如图2 1 2 。 第_ 章有效性分析 1 6 ( a ) u ( 厶。) 的实部( b ) u ( 厶z ) 的虚部 图2 1 2 ：例5 使用重置方法求解特征问题，当步长7 - = 1 2 5 6 时u ( l ，名) 的波形 取步长7 分别为1 4 ，1 8 ，1 1 6 ，1 3 2 ，1 6 4 ，1 1 2 8 ，1 2 5 6 ，其它参数不变。假 定r = 1 2 5 6 时为精确解，记步长为7 时计算结果为u r ，相对误差为e ( 7 - ) = 骂导业掣，则相对误差和计算时间如下表： i i u l l 2 5 6 l i 步长7 相对误差e ( 7 )计算时间( 秒) 1 4 3 8 9 5 1 9 0 1 0 17 1 4 8 8 0 8 1 0 1 1 8 1 6 4 7 7 1 9 1 0 11 4 1 1 0 0 0 1 0 2 1 1 6 3 2 1 8 3 4 9 1 0 22 7 9 2 8 8 5 1 0 2 1 3 2 4 4 0 1 7 8 4 1 0 35 5 6 4 5 3 3 1 0 2 1 6 4 3 8 1 7 3 2 7 1 0 31 1 1 2 2 1 7 1 0 3 1 1 2 8 4 1 1 3 7 2 7x1 0 32 2 2 1 7 9 9 1 0 3 1 2 5 6 0 0 0 0 0 0 04 4 3 5 2 9 9 1 0 3 表2 1 ：例5 中步长变化时相对误差和计算时间对比 从表2 1 可以看出，给定步长7 - 时，所需计算时间与步进次数呈线性关系。将 区f , - j 0 ，驯关于圣离散为o = 圣o 童1 主：s - 1 丸 丸+ 1 盒m = l ， 从丸步进计算至以一1 大约需要1 7 秒，随着剖分的加细，需要的时问越来越长，而 第二章有效性分析1 7 精度也越来越高。 当口= 2 0 的时候，使用r a y l e i g h 迭代方法和重置方法求解问题都可以进行 计算，并且结果基本一样，相对误差只有1 0 一。 ( & ) t ，z ) 的实部( b ) , 4 l ，。) 的虚部 图2 1 3 ：例4 分别使用r a y l e i g h j , 塞代方法和重置方法( 7 - = 1 2 5 6 ) u ( l ，z ) 的波形 对比 再对比两种方法在不同的步长所需要的时间，其中时间单位为秒。 步长丁 r a y l e i g h 迭代方法计算时间 重置方法计算时间 1 4 4 2 4 7 9 6 2 1 0 17 1 4 8 8 0 8x1 0 1 1 8 8 0 0 6 0 8 2 1 0 1 1 4 1 1 0 0 0x1 0 2 1 1 6 1 5 2 7 2 2 8 1 0 22 7 9 2 8 8 5 1 0 2 1 3 2 3 0 1 5 9 3 9 1 0 25 5 6 4 5 3 3 1 0 2 1 6 4 6 0 5 2 7 6 5 1 0 21 1 1 2 2 1 7 1 0 3 1 1 2 8 1 1 9 9 9 5 0 1 0 3 2 2 2 1 7 9 9 1 0 3 1 2 5 6 2 3 0 9 0 0 4 1 0 3 4 4 3 5 2 9 9 1 0 3 表2 2 ：例4 中步长变化时使用r a y l e i g h 迭代方法和重置方法计算时间对比 从表2 2 中可以发现，对于相同的步长，使用r a y l e i g h l 塞代方法时间较短，比 第- 章有效性分析 1 8 重置方法节省近一半的时间。但这种方法对于较陡的界面则不适用，比如例5 ， 在使用范围上有一定的局限性。另外，我们注意到，对于例5 这种较陡的界面， 使用重置方法求出的波形震荡也较大，为了验证数值模拟的正确性，下面本文 对数值计算中两个重要过程进行分析。 2 2 线性方程组分析 通过分析步进算法的具体步骤，我们发现在求解线性方程组 p = ( 国+ i v y ) 一1 ( 一国+ t 天) 国。一1 = i 五( ，一只) ( ，+ r ) 一1 e l= e z ( j + p ) e r 瓶( ，+ r ) 一1 ( 2 1 ) ( 2 2 ) ( 2 3 ) 的过程中，系数矩阵的条件数很大( 图2 1 4 ) ，而且变化不规律，有可能出现病 态现象，下面尝试对这种情况进行改善。( 这里取矩阵的2 条件数进行分析。) c o n d ( a ) = i ia1 1 2 l ia - 1 | 1 2 = 对原方程作某些预处理，可以降低系数矩阵的条件数。因为c o n d ( q a ) = c o n d ( a ) ， 显然不能通过每一个方程都同乘以一个常数q 的方法来处理。但是可以采用将 矩阵的每一行( 或每一列) 分别乘以适当常数的方法。即找到可逆的对角阵d 1 和d 2 ，使得方程组a x = 6 化为 d 1 a d 2 y = d l b z = d 2 y 这称为矩阵的平衡问题。理论上，最好选择合适的对角阵石l ，百2 满足 c o n d ( d 1 a d 2 ) = m i nc o n d ( d 1 a d 2 ) 其中的m i n 是对所有属于可逆对角阵的d 1 ，d 2 取的，这在实际工作中很难实现。 一般矩阵的平衡问题要针对具体问题具体分析，而不是寻求般适用的策略。 d 1 是用来平衡系数矩阵的，现是用来平衡未知数的。一种简单的方法是 令d 2 = i ，再选择d 1 使得d l a 的每行的o o 一范数大致相等，这可避免消元过 程中小数与大数的相 j 1 1 1 2 9 。 第二章有效性分析 1 9 本文对方程( 2 1 ) ，采用行平衡法，对于系数矩阵a = 【a i j 】n n = q + i 、a ， 取d i = m a ，xi a i j l ，i = 1 ，2 ，n ，令d t = d i a g ( 1 d l ，1 d 2 ，1 厶) ，则原方程 由a x = b 的形式转化为d 】a x = d 1 b ，这时c o n d ( d 】a ) c o n d ( a ) 。对于方 程( 2 2 ) 和( 2 3 ) ，也可以使用类似的列平衡法。 针对例5 ，在计算过程中，增加对矩阵国+ 狐和，+ 刷挣平衡预处理，然 后再求解方程组。图2 1 4 为例5 在预处理之前，相关矩阵的条件数。其中横轴表 示声波的传播方向，把( 0 ，l ) 进行2 5 6 0 份等距剖分，0 = 2 0 仝l 圣。一l 筑 圣曩+ 1 圣2 5 6 0 = l ，纵轴表示在筑处( s = 1 ，2 ，2 5 6 0 ) ，求解方程 组( 2 1 ) 一( 2 3 ) 时，系数矩阵的条件数。图2 1 5 为将例4 进行预处理之后，得到的条 件数。从图2 1 4 与图2 1 5 对比中可以发现，预处理后的矩阵条件数明显降低了， 但总体来看，依然分布杂乱，没有规律。 ( a ) 口+ i 瓜的条件数( b ) i + r 的条件数 图2 1 4 ：例5 进行预处理之前相关矩阵的条件数 第二章有效性分析 ( a ) 审+ t 瓜的条件数( b ) i + r 的条件数 图2 1 5 ：例5 进行预处理之后相关矩阵的条件数 我们选取在进行预处理之前，国+ z 瓶的条件数较大的几个点，和周围的几 个点，来对比改善情况。 步数预处理之前的条件数预处理之后的条件数 1 3 5 09 8 1 3 7 86 0 7 2 1 1 3 5 12 1 4 8 7 5 11 2 5 4 1 1 1 3 5 25 5 1 7 1 8 13 3 0 5 1 0 1 3 5 31 9 0 4 3 0 81 2 2 3 8 0 1 3 5 41 0 9 1 2 5 3 7 7 0 2 7 1 3 5 57 7 5 0 5 06 0 7 9 9 1 3 5 66 0 7 0 2 84 8 2 4 4 表2 3 ：例5 计算过程中预处理前后国+ i 以的条件数对比 从表2 3 中可以看出，经过预处理之后条件数有所降低，但原先国+ z 瓶的 条件数较大的几个点，与其它点相比仍然较大。再比较预处理前后，方程( 2 1 ) 的 残差i i ( 国+ i v 伍) p 一( 一国+ i v x ) i i 。 第一章有效性分析2 1 步数预处