（计算数学专业论文）在多层无界区域中helmholtz方程的数值解法及其应用.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 方程的方法是设立假的边界条件将无界的求解区域有界 化，如设这些边界条件为第一、第二或第三类边界条件，这种做法只是将原方程 求解粗糙地近似化，所以无论是对此近似问题进行精确求解还是离散方程后进行 数值计算，解的偏差会变得很大，导致解的逼近效果变差。产生偏差的主要原因 在于这样的边界条件实际上并不符合无界区域中波传播的特性( 产生了较强的反 射波) 。基于此缺陷，j p b e r e n g e r 于1 9 9 4 年提出了完美匹配层( p e r f e c t l y m a t c h e d l a y e r ，简称p m l ) 的概念，1 9 9 7 年f r a n c i sc o l l i n o 针对上述添加的p m l ，在 数学上相当于对坐标做了一个复的伸展变换量= x + i xr r ( r ) d r ，把方程转换为一 圳 个新的复方程，我们称之为改进的复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 方程时，产生的线性系统的阶数将非常的大，导致相当大的存贮空 间，计算的代价也很高昂。同时这些系统常常也是不定的，或非对称的，这就使 得方程的求解更加困难。根据波导对区域水平距离依赖很弱的特性，本文采用在 波的传播计算中有效常用的步进方法，如o n e w a y 方法计算声波的传播。 最后，用步进方法计算波的传播问题将涉及方程的特征问题。曾有人做过界 面是平坦的光波导的特征值问题，而本文研究的是在带有弯曲界面的声波导中的 波传播性态。我们先选取适当的非线性坐标正交变换及方程变换将界面拉直，然 后从数值角度出发，离散上述问题的特征方程相应的算子，得到一复矩阵：为求 特征问题，利用r a y l e i g h 迭代方法的具有局部收敛和快速收敛的特点，构造多 重r a y l e i g h 迭代算法，成功地求出此复矩阵的特征值与特征向量，为波的传播 计算提供了基础。 结果表明，上述算法不仅可以较好地求出h e l m h o l t z 方程在无界区域中的模 ( p r o p a g a t i o n m o d e ，l e a k y m o d e ，b e r e n g e r m o d e ) 的特征分布，而且用此分布 可以较精确地求解无界区域上h e l m h o l t z 方程。我们所设计的方法具有易于数值 计算的优点，如保持三对角矩阵运算，极大地减少了存储空间和计算量。可在 光波、声波在无界区域上的传播计算中得到应用。 a b s t r a c t i nt h i sp a p e r ，w ep r o v i d eam e t h o dt os o l v eh e l m h o l t ze q u a t i o ni nu n b o u n d e d r e g i o n t h ef i r s t , w et r a n s f o r mu n b o u n d e dr e g i o nt ob o u n d e dr e g i o n m a n yk i n d so f b o u n d a r yc o n d i t i o n ss u c ha st h ef i r s tc o n d i t i o n ，t h es e c o n dc o n d i t i o n ，a n dt h et h i r d c o n d i t i o nh a v eb e e nu s e dt os o l v et h i s p r o b l e m t h e s ec o n d i t i o n s c a r l o n l y a p p r o x i m a t et h ee q u a t i o nr o u g h l y , t h i sr e s u l t i sn o tv e r yg o o d ，b e c a u s et h e s e b o u n d a r yc o n d i t i o n s a r en o tf i tf o rt h er e a lc a s e s a c c o r d i n gt ot h el a c ko ft h i s a p p r o x i m a t i o n ，j 。p tb e r e n g e rp u tf o r w a r dt h ep e r f e c t i ym a t c h e dl a y e r ( p m l ) i n1 9 9 4 f r a n c i sc o l l i n ou s e dac 。m p l e xt r a n 8 f o r m 量= x + ira ( f w rt 。c h a n g eh e l m h 。i t z e q u a t i o nt oa n o t h e re q u a t i o ni 1 11 9 9 7 ，a n dt h i se q u a t i o ni sc a l l e da sa ni m p r o v e d h e h n h o l t ze q u a t i o n 。 t h es e c o n d ，t h e r ea r em a n ym e t h o d st os o l v eh e l m h o l t ze q u a t i o ni nt h eb o u n d e d r e g i o nw i t hav e r yl a r g el e n g t h ，a sc o m p a r e dw i t ht h et y p i c a lw a v e l e n g t h s t a n d a r d m u m e r i c a lm e t h o d s ，s u c ha sf i n i t ee l e m e n ta n df i n i t ed i f f e r e n c em e t h o d ，g i v er i s et o v e r yl a r g el i n e a rs y s t e m s t h e s el i n e a rs y s t e m sa r ed i f f i c u l tt os o l v e ，s i n c et h e ya r e n o n s y m m e t r i ca n di n d e f i n i t e f u r t h e r m o r e ，w h e nt h e s eb o u n d a r yv a l u ep r o b l e m so f t h eh e l m h o l t ze q u a t i o na r ed i s c r e t i z e d ，t h el a r g en u m b e ro fu n k n o w n sa r es o l v e d t o g e t h e rg i v i n gr i s et oav e r yl a r g ed e m a n df o rc o m p u t e rm e m o r y f o rw a v e g u i d e p r o b l e m s ，t h el e n g t hs c a l ea l o n gt h em a i np r o p a g a t i o ni sm u c hl a r g e rt h a nt h e t r a n s v e r s el e n g t hs c a l e s ow eu s eo n e - w a ym e t h o dt oc a l c u l a t ea c o u s t i cw a v e p r o p a g a t i o n 。 t h el a s t ，w h e no n e - w a ym e t h o di su s e dt oc a l c u l a t ea c o u s t i cw a v ep r o p a g a t i o n ， t h ee i g e n v a l u e sa r en e e d e d ，t h er e s e a r c ho ns l a bo p t i c a lw a v e g u i d eh a sb e e nd o n e ，i n t h i sp a p e r ，w ei n v e s t i g a t et h ea c o u s t i cw a v ep r o p a g a t i o no nu n b o u n d e dr e g i o nw i t ha c u r v e di n t e r f a c e 。w eu s e al o c a lo r t h o g o n a lt r a n s f o r mt of l a tt h ec u r v e di n t e r f a c eo f t h ew a v e g u i d e t h e n ，w ed i s c r e t eo p e r a t o rc h a r a c t e r i s t i ce q u a t i o n ，a n dg e tac o m p l e x m a t r i x a n dw em a k eam u l i t - g e n e r a l i z e dr a y l e i g hq u o t i e n ti t e r a t i o nm e t h o dt o c a l c u l a t et h eo p e r a t o r se i g e n v a l u e s i tp r o v i d e st h eb a s eo fc a l c u l a t i n ga c o u s t i cw a v e p r o p a g a t i o n ， t h en u m e r i c a le x a m p l e ss h o wt h a to u rt r e a t m e n tc a nn o t o n l y c a l c u l a t e h e l m h o l t ze q u a t i o n se i g e n v a l u e s ( p r o p a g a t i o nm o d e ，l e a k ym o d e ，b e r e n g e rm o d e ) ， b u ta l s oc a l c u l a t ea c o u s t i cw a v ep r o p a g a t i o ni nu n b o u n d e dr e g i o n t h i sm e t h o dh a s s o m ea d v a n t a g es u c ha st r i a n g l e d i a g o n a l ，s a v i n gc o m p u t e rm e m o r ya n ds oo n i tc a l l b eu s e di nc a l c u l n i n gw a v ep r o p a g a t i o ni nu n b o u n d e d r e g i o nw i t hac u r v e di n t e r f a c e l 引言 在声学、电磁学、地震学及其它许多应用领域的大规模的波传播问题的研究中，通常需要 在一个长度比波长还要大许多的区域中求渡场的分布。以海洋的声波传播为例，在探索开发海 洋的过程中，海洋声学起着特别突出的作用，因为只有声波才能在水中进行较长距离的传播， 从而进行水下测距、定位、通讯和遥感等操作。 海洋是一种极其复杂的声学介质。海洋介质最具特征性的现象是非均匀性，这种非均匀性 强烈地影响着海洋中的声场同时海洋波导的特性随着水平距离也有变化，但变化非常缓慢 声速随深度有规律变化会形成水下声道利用水下声遭可以进行数百，甚至数千公里的远程声 传播，但受海水吸收的限制，要求用于传播的声波频率u 很低，即相应的波长 要很大才行 若波源是单一频率时间调和的，则波场的决定方程是以下h e l m h o l t z 方程 z 。+ u ；：+ k 2 ( 。，z ) u = 0 ，0 = + c o ，0sz + o 。 其中z 方向是垂直向下的，以海平面为z = 0 z 表示水平方向距离变量。 对上述h e l m h o l t z 方程，由于k 在水平z 方向上是弱衰减的，故可以把x 方向分成两部 分，即0 os l 和z l ，其中这两部分由水平方向o u t g o i n g 条件连接。当z 兰l 时，可由 分离变量的方法得到方程的近似解，本文主要考虑的是0 z 兰三这部分h e l m h o l t z 方程的求 解问题 当k 是常数，并且定解区域是矩形区域时，h e l m h o l t z 方程有精确平面波解和球面波解 但实际的大部分问题中，一并不是常数，而是与水平距离z 和深度z 等有关的，故无法直接给 出如下h e l m h o l t z 方程z + ：+ 萨( # ，。h = 0 ，0 。 + o o ，0 z l 的解析解。 首先，对于z 方向上的无限深度问题，以往的方法是简单地设立假边界条件将无界的求解 区域有界化，如设这些边界条件为第一、第二或第三类边界条件。这种做法只是将原方程粗糙 地近似化，所以无论是对此近似问题进行精确求解还是离散方程后进行数值计算，解的逼近效 果都较差产生误差的主要原因在于这样的边界条件实际上并不符合无界区域中声波传播的特 性 在这种情况下，必须引入适当的吸收边界层条件把计算空间截断，并且应保证在截断边界 处只有向外传播的波而没有向内传播的反射渡。吸收边界的效果直接关系到数值计算的正确性 和精确性，是影响计算品质的决定因素1 9 9 4 年，jp b e r e n g e r 首先提出了高效的二维理想 匹配层( p m l ) 吸收边界条件的概念【1 0 】。后来，又在理论上证明了该方法可以完全吸收来 自各个方向、各种频率的声波，而不发生任何反射。它被认为是目前最好的吸收边界层条件 1 9 9 7 年f r a n c i sc o l l m o 针对上述添加的p m l ，在数学上相当于对坐标做了一个复的伸展变换 女= z + i 学f ( r ) 打。把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 方程有许多直接的数值方法，如有限元和有限差分法 等但方程定解区域的水平距离与波长相比非常的大，而且渡场的内部振荡也很大，每一个波长 范围内部需要用很密的格点( 或基函数) 来表示波场，对于这个水平距离比波场还要大得多的 区域来说，若用有限差分或者有限元方法来处理这样一个水平区域很大的h e l m h o l t z 方程时， 产生的线性系统的阶数将非常的大，导致相当大的存贮空间，计算的代价也很高昂。同时这些 系统常常也是不定的，或非对称的，这就使得方程的求解更加困难 以有限差分法来求解以上的h e l m h o l t z 方程为例，离散时分别在z 方向和z 方向取个 点和m 个点，若步长相同，即有m ( l d ) n ，总的存贮空间的需求就是o ( l d ) n 2 ，离散 所得的系统是非对称的和不定的，并且它的带宽是，即使用带状矩阵的处理方法来求解，运 算所需空间的数量级将是o ( l d ) n 4 随着划分的加密，所需的存贮空间将是非常惊人的 就海洋中声波导问题来说，一方面，在波传播方向上求解区域的长度，要比纵向( 垂直于 传播方向) 的要大得多( 都远比波长大) ；另一方面，波导介质沿着传播方向有变化，即和区域 是相关的，但是很细微显然一个好的数值方法应该充分利用声波导的这两方面的特性利用 波导对区域水平距离依赖很弱的特性，有很多更合适的方法本文采用在波的传播计算中有效 常用的步进方法，如0 n e - w a y 方法计算声波的传播 最后，用步进方法计算波的传播问题将涉及方程的特征问题曾有人做过界面是平坦的光 波导的特征值问题【3 j ，而本文研究的是在带有弯曲界面的声波导中的波传播性态为求特征问 题，我们先选取适当的非线性坐标正交变换及方程变换将界面拉直，然后从数值角度出发，离 散上述问题的特征方程相应的算子，得到一复矩阵，利用p “y l e i g h 迭代方法的具有局部收敛和 快速收敛的特点构造多重r a y l e i g h 迭代算法，成功地求出此复矩阵的特征值与特征向量，为 波的传播计算提供了基础 以下第二章先介绍了声波导有关的二维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 方程的离散；第六章中介绍了利用多重广义l l a y l e i g h 商叠代法计算其对应 矩阵的特征值；第七章介绍了波的传播计算方法；数值例子将在第八章中给出。 我们从二维h e l m h o l t z 方程开始 2 基本方程 其中0 曼x l ，o z h b ，处于0 1 ，h b h d ，l d ，h ，h b 为常数，如图1 所示以下给出上端和底端的边界条件：u ：o = 0 ，i ；：日= 0 连接条件为 ：j 恕) _ “( z ，z ) = ：j 路) + “( z ，z ) ，击：跺) - 翌! 杀兰= 去；j 豫) + ! ! 舞型 ( 2 ) 其中n 是z = h ( x ) 的法向量。假设在zs0 和。l 时波导是与水平距离无关的设在z 0 时，k ( z ，z ) = k o ( z ) ， ( z ) = h o ；在z l 时，一( z ，z ) = k ( 兰) ：垒( 兰) 三垒妒如果没有波从+ o 。 传过来，则在z = l 处的边界条件( 散射条件) 是。= 以2 + k 。2 ( z ) u 。其中i = ，j 在 z = 0 处最简单的边界条件就是在z = 0 取u = u o ( z ) ，其中u o 是一个关于= 的给定的函数 假设解是存在且唯一的 2 h ( x ) 弋 l a y e t i i l a y e t i i 图1 ：基本方程示意图，其中= = ( 。) ：1 一 e 一。( 量一r ，# ：0 2 ， = 1 0 0 ，= 1 0 ，0 2 3 0 ，0 $ 曼1 0 3 局部正交坐标变换 在o ；s ( 。) 即第一层进行局部正交变换【2 令 ；：善爿，取9 ( z ，z ) 2 南，我 们还可以得到以下的公式： z 。器班+ 扣瑚2 - o 对于给定一新坐标( ，2 ) ，可以通过( 3 ) 式用牛顿叠代来求解。，从而。= oh ( z ) ；如果 给定一旧坐标慨= ) ，则2 = 志同时女也可以由( 3 ) 式获得 附注1 如果g 满足h i ( z ) = o ，则令= 2 = i 函 在 ( z ) 三= 兰d 即第二层进行局部正交变换【“令 ：21 j 。? ，我们可以取 z = 9 ，o j 如翻= 丽1 - d 蚪键糟 我们还可以得到以下的公式： e 帮班妒聊邓) _ d 用一o ( 5 ) 对于给定一新坐标( 圣，未) ，那么。= 弘一荆) 锚，z 可以由( 5 ) 式得到； 相反地，给出( 。，= ) ， 可以由( 4 ) 式得到，由( 5 ) 式用牛顿叠代可得到矿，女用牛顿叠代从 3 z 。端1 汗一。解得 附注2 如果z 满足 ( 。) = o ，则令岔= z ，童= 耥z + 与糟 在d 。h 即第三层进行局部拉伸变换1 2 1 。令 ：2 ，( 2 ，d ) ，对于给定一新坐标 ( ，1 ) ，那么# = 2 ，z 可以由( 5 ) 式得到，其中令0 = d ；相反地，给出( 。，z ) ，那么 = ；，由 ( 5 ) 式用牛顿叠代可得到矿，其中令z = 。，用牛顿叠代从z 。嘉譬出+ ；似矿) 】2 = 。解得 在日兰。蔓h b 即引入( p m l ) 的第四层进行局部拉伸变换令t ：3 ，( z ，d ) ，可知满 足 “z ，g ) l 一。 z + c o ，h z s h b 暂 ( 士，童) f 一。 圣 + 。，h 童h b 篙- 爨+ 甏塞= 吣( z ，硼= e 如，日印= 日b 对于给定一新坐标( ，1 ) ，那么= = ，# 可以由( 5 ) 式得到，其中令= d ；相反地， 给出( z ，z ) ，那么2 = z ，由( 5 ) 式用牛顿叠代可得到矿，其中令。= d ，2 用牛顿叠代从 z 。端弘1 ) 2 = 。解得 夕 ! z 图2 ：坐标变换示意图，其中叫表示( j ) ，t o ? 表示( x ，z ) ，且z = h ( z ) = 1 一e e 一。f 量一i r ， = o 2 ，e = 1 0 0 ，l = 1 0 ，o 曼 ： 3 0 ，0 # s 1 0 4 4 引人p m l 的方程变换 对于( 1 ) 式，由坐标的正交变换可转化为以下的形式： + a p ，鲁) k + 卢 ，辛) k + 7 ( 窖， ) y = o ， ( 6 ) 令( ，。) = w ( x ，# ) y ( z ，= ) 在第一层中忙】，选择( z ，z ) = 、黼，则( 1 ) 式可以简化为( 6 ) 式，其中( z ) 0 附注3 在第一层中，如果 7 ( i ) 2o ，我们定义o ( i ，j ) 21 觋a ( ， ) ，口( ，o ) 21 霉：p ( 。，。) 7 ( 2 t ) 。! i 。m ( 女，2 ) ，则 卜，沪南e 。帑洲础卜。z 器e “2 警， 卜沪髯o 。；爹卜e ，撕 更进一步，如果 7 忙) = o 且矿忙) = 。，则d ( 2 ，2 ) = 南，口( i ，) = 。，_ ( i ，1 ) = 一邓，一譬觜 在第二层舻】，选择肺，加1 端笺掣，其州舢舭) 5 “，j ) = 芦( 圣，动= 7 偿，= 附注4 在第。层中，如策 忙) 。o ，我们定义o ， ) = ! 骧“( 女，j ) ，p ( i ，童) = ! 警口( ，j ) 1 ( 2 ，) 21 噢7 ( ，e ) ，则 “( i ，2 ) = 厕( ) - 1 ) 2 e 篙岛( “( ) m 2 “， 滕，弘迎誉冬器燮s 龋妒q n 哪雄、) ， 化，护e 紫忙。2 。h n 吣) ) 等【罱蓊# 一 西警蒜十3 e 一箬黯妒_ 2 8 啪姆】一 ! 黑毯i e 拦裔( 一m ) ) ( 2 d m ) ) + 4 ( 劫f d 侮) 1 。 拦球涵计篆鬻卜e 龋妒。辞n 堆翻) 更进一步，如麓搿尊) 絮。豆酽穗) = 。，则a ( 是童) 裟踹， 芦( 聋，童) ：。；7 ( 墨 ) = 一2 ( 。，1 ) + ；- i ：二二! ；掣 t 4 ，( t ) 式 雀第三屡中【“，遮择w ( x ，# ) =其中掰和) 0 ，辩( 1 ) 蝣国= 簪端嬲， 俺国= 错卫群一1 器三鬻+ 丽而d 4 42丽x? 。，。j 3 焉舀莎一一+ 丽+ 一五i 孬可一丽+ 百i = 霹；巧y 疆i 可 ：f ( z + ) - x i ：；赫i ：；5 ；斋) i j ：j ；! ；i 2 + 一2 ( z ，s ) + 型嗡糍产燮+ 、拦裔+ 丛譬鼎逖 x 勰簧鞘器瓣j 2 t 附注5 在第三层中，如策睁) 2o ，我僻定义a ( e ， ) 。! 象8 ( ，声擘，2 ) 2l 臻p ( 2 ，o ) 8 磊巫州彗州巫 d ( i ，e ) = e - d h ( ”， 幛，牡j 3 酱掣慨卅i 3 器) e 。+ 丽掣丽 华一轫+ - i 觜小。 _ 1 ) ； 更进一步，如果 ( i ) = o 且( i ) = o ，则。( i ，= 1 ，p ( i ，。) = o ，7 ( i ， ) = k 2 ( 2 ，z ) + ： d ( 4 ) f 劫 在第四层中，由于p m l 的引入原方程变为： 丕( 坼) 笔) 柑( 叩) u = o 其中p ( 。) = 高，h z h b 同样令u ( 。，z ) = w ( 置! ! ：! 鱼生： 我们选择帅，加器。掣芒，其州舢则，式可以简化 r + k 2 ( z ，z ) + 性质1 在第四层中，如果和) = 0 ，我们定义血( 童，劫= ! 嗡口( 童， ) ，卢( 牙，2 ) = ! 粤；口( 套，o ) 7 ( j ，2 ) 。! 觋7 ( i ，j ) ，则 fa ( ，j ) _ p ( z ) 2 e - d h ( ”， ip ( ， ) = p ( z ) p ( z ) e - d h ( ”， 卜，牡i 酱十舻喊卅i 3 器) e 。+ 丽等厕 【 学一；卅i 1 觜 c 一1 ) ； 更进一步，如果h 国) = 0 且 ”( i ) = 0 ，则o ( i ， ) = p 0 ) 2 ，p 忙， ) = p ( z ) p ( z ) 7 ( j ，1 ) = 舻( g ，z ) + i i d h ( 4 ) ( i ) 7 鬓筘 pp他生。兰丛叫j 水圳弦崭笪黪露 5 特征值问题及其离散 我们考虑丑( t ， ) 的特征值问题，其中 r ( 者，j ) = 口( 士，害) k i + 卢 ，童) + 1 ( 未，o ) ef 8 ) 以下( g ) ，( 1 0 ) ，( 1 1 ) 式是对r ( 女，j ) 进行离散化的边界条件： 对于( 8 ) 式，它的两端边界条件如下：y j ：o = y b ：日日= o ，第一与第二层曲面交界条件 转化为如下的形式 2 1 ： ii 蟀w ( 舭) 。y ( 挪) = l i m + w ( z ，= ) y ( 叩) ， 击。l i r a w ( 叩) ( 渺( z ) 一2 毁擎 矿一掣) = 【去。骧w ( 。，z ) ( 。) + 2 告! 糟】矿一z ；赢f 1 + ，( 。) ) 。 k ) 对于( 8 ) 式，第二与第三层曲面交界条件转化为如下的形式 2 1 ； j 。l i r a 。一( 钔) y ( 舭) = l i m + ( 叩) 矿( 舭) ， i 。叶l i r a 。一阿扛，z ) 端砼之；璺磐o ，z ) - 对于( 8 ) 式，第三与第四层曲面交界条件转化为如下的形式： jj 骝y ( 舭) 1 i m 日+ v ( 叩) ，j _ h 、。 一日十 1 ” 【j + t i m 一畦2 。l i r a hh + k 、j - + 一 叶h + 1 根据( g ) ，( 1 0 ) ，( 1 1 ) 的边界条件，可把r ( 自，离散戒三对角的复矩阵，型如： 徉91 肚i ： i i如b b n 。- 1 麓1 我们把方程在2 方同进行离散，其中i = 0 0 = 1 进行n 1 阶划分，2 = 1 一 = d 进行n 2 阶划分，l = d 一 = 日进行n 3 阶划分， ：日。2 ：h b 进行n 4 阶划分 - = 矗砸，h ，= 弓帚， s = 与铲，k = 嬲，血为水平e 方向的第t 个离散点，白 为深度 方向的第个离散点；矩阵的具体数值如下： 当j = 1 时， f 山= 0 驴州血矧。2 掣； ( 1 2 ) 【。= 掣+ 掣； 8 卿 当1 j n 1 时 当j = 1 时 当= n 1 + 1 时 当n 1 十1 j n 2 时 9 ( 1 3 ) ( 1 4 ) f 1 5 】 ( 1 6 ) 掣拳 惹竿 铲耻铲 帝一岛h +旁 掣拳 蔫掣 铲驴耻 当j = n 2 时 当j = n 2 + l 时 ( 1 7 ) 。、m ：、2 击之胪 铲e 掣一掣，蔫雾； 确h 掣e 掣一酶，凳誊；( 1 8 ) 岛= 掣十掣； 当n 2 + 1 j n 3 时 当，= n 3 时 ( 1 9 ) 一【掣+ 掣 蹴； ( 2 0 ) 毒； f 山= 【掣一帮1 瓶； 卜懒诤掣+ 掣一掣，蹴； ， 【q = 掣+ 睁； 塑一 掣棼券 掣掣掣 掣一掣 掣攀 篱竿 铲驴耻 博攀 主! 惹二琴； 6 特征值i 目题的数值群痃 对于以上三对角的复矩阵，由于矩阵阶数大时，现有软件计算特征值结果不好，故我们用 以下几种叠代法进行求解比较： 对称r a y l e i g h 叠代法 1 0 h = 兰 1 2 d k x 垡k 。l l ，其中宴k 是矩阵a 的右特征向量的共扼转置，凰是矩阵j 4 的右特征向 量 2 。( a k l a 1 x ，k + l = x k ， 3 。凰+ 一横勘， 4 。a k + l = l 亟1 2 k + ! l 生x 亟k + 1 2 ， 5 。 | l 生i i 等址| | f ，其中是r a y l e i g h 叠代控制精度， 6 。 如果步骤5 。不满足，则 i + l = h ，x + l = x ，跳到步骤2 。，进行循环如果进行对称 ra _ y l e i g h 叠代到一定次数还不满足叠代控制精度，则跳出循环 非对称r a y l e i g h 叠代法 1 d h = 畿，其中x 是矩阵a 的左特征向量，其中k 是矩阵4 的右特征向量， 2 。 x 再- ( h 、- ，日a ) = x 宁，( 札i 直) k + l = 托， 3 。碍- 2 尚鲰一隅， a 舻毪等， 5 。 | i 地元者止| | 已其中f 是r a y l e i g h 叠代控制精度， 躲 站 驰 6 0 如果步骤5 。不满足，则a k + l k ，x 韪1 x ? ，k + l j k 跳到步骤2 0 ，进行循环如果 进行非对称r a y l e i g h 叠代到一定次数还不满足叠代控制精度，则跳出循环 多重非对称r a y l e i g h 商叠代法 2 】 令m 是一正整数，表示把c 分成m 等分给出较精确求且特征值的多重非对称r a y l e i g h 商 叠代算法具体如下： 第一步：令j = o ，a j ( z ) = 0 ，对应得到的矩阵山是实矩阵，求出其特征值 a ) ； 第二步：对j = 1 ，2 ，m 1 。计算q = 等，o a z ) = q r ，r = 看膏兰备，c r + ； 2 。 计算对应方程离散化后的对应矩阵a ，； 3 。 分别以a j 一1 得到的特征值 冲“) 为初值，进行非对称广义r a y l e i g h 商叠代，得到a ， 的特征值 ) ，；= 1 ，2 ， 注：矩阵a ，是由矩阵a ，一1 微小扰动得到的。根据有关定理，我们可得到，从山一处求 得的特征值可以作为在山处的非对称广义r a y l e i g h 商叠代的初值，从而保证它的局部收敛性 质。另外，由于山是一个实数范围内的三对角矩阵，通过把其相似对称矩阵，从而利用l a p a c k 中的程序包求出其特征值 7 波的传播计算方法 o n e - w a y 方法【1 1 一些基于e x a c to n e - w a yr e f o r m u l a t i o n s ( e o r ) 变换的近似方法有着比c p m 方法更好的优 点，其中广泛使用的抛物方程近似方法( p a r a b o l i ce q u a t i o n 堤多种对h e h n h o l t z 方程单向( o n e - w a y ) 近似中的一种，对于式( 6 ) ，其可转换为如下型式：矿( 奎州，三) 8 i h 、。明+ 肫+ y ( 莹。，2 ) 具体离散形式如下： 1 。v ( e 。，1 ) = v 0 ) ，2 ，k 骝2 ，一，k 骝2 陋1 ，n 2 ，口。 8 y ( 女。+ - ， ) = 【k ；i 。，k 搿。，k 骝，。 l a l , 岛，风】8 2 。 泸1 ，卢2 ，卢n 圩= e t i + h , l x 2 t je i - “i 1 五2 ，e i 认+ 1 1 2 ”j i 。十l ，。陋l ，口2 ，n n 】h 3 。 哦? i 脚k 翳2 ，。k 措。】归- ，胁，风 ”= k 2 ，噼j 伊，k 。 陋t ，a ”。n 。】” 4 。 根据3 。式求得【a 1 ，口2 ，o 。严 5 。 如果没计算到传播划分结束，则将其代入2 0 式，重复上述两步，直至传播划分结束 8 数值模拟结果 下面我们通过大量的数值模拟例子实现这篇文章中所提到的方法 第一个例子中，我们在p e k e r i s 情况下，即 ( z ) = 1 ，说明多重非对称l h y l e i g h 商叠代法 对此问题适用： 设0 z 1 时p t = 1 0 ，1 z h b 时, 0 2 = 1 7 ，d = 1 5 ，h = 2 0 ，h b = 3 0 ，水平距 离l = 1 00 ，处于0 2 1 时k = 1 6 0 ，处于1 z h b 时k = 1 4 4 ，且盯( 。) = c 南，r = 膏再! 岛，c r + ，以下1 一4 都代表离散的阶数! 首先，引出对于n e w t o n 叠代法中的，( ) ，并做个数值模拟： 1 2 f 西( z ) = e x p ( i 西f 五z ) 一e x p ( 一i 厢z ) = 2s i n ( 拆；口# ) ，0 = 1 ； 豢e l 篓餐拦型e x p ( ：篇未h 冀h b ， l ( z ) =- e 3 币0 k ；一 皇) + 。已一e 、碚一 毒) ， # 、吖 【l = 。+ t 戌口( t ) 出 及边界条件饵b ) = 0 和在z = 1 的界面条件 ( 2 6 ) 可得 ，( a ) = 击、石 = ic o s ( 、石i 。= 1 【1 一e x p ( 2 i 、，乍l ：万) ( h b 一1 + i 正芋8 盯( r ) c b ) 】) 击。厢s m ( 厢 1 + e x p ( 2 i 币) _ ( h b 一1 + i 詹8 口( r ) 如) ) = o 当不引入p m l ，由o u t g o i n gc o n d i t i o n 塞= t 污五， d n v 可得 9 ( ) = p 。历c 。s 历一p ，i 厢s i n 历= o r e ( 五) 图3 ：n e w t o n 叠代法n 1 = 2 0 0 ，n 2 = 1 0 0 ，3 = i 0 0 ，n 4 = 2 0 0 时的特征值 图( 3 ) 是取9 ( ) 为叠代函数得出的p e k e r i s 解的结果。此法对于p e k e r i s 情况下是可以实 现的，但由于其他的情况，如不能给出特征函数的解析解，不一定可以得到一个非线性特征方 程因此我们有必要构造多重非对称r a y l e i g h 商叠代求解特征问题 首先，我们验证多重非对称r a y l e i g h 商叠代的算法的正确性，以p e k e r i s 情况下，即分段 一的情况为例： 1 3 警上m 荆翔 + i | 馏掣似上m 味磐 r e ( 旯) 图4 ：7 一表示p e k e r i s 精确解1 , 1 表示多重非对称r a y l e 谵h 商叠代g ：1 0 0 0 0 的解 由图( 4 ) 比较可知，用多重非对称r a y l e i g h 商叠代算法所得的结果和p e k e r i s 问题的解吻 合，故多重非对称r a y l e i g h 商叠代算法对此问题也适用 其次，我们讨论有界的定解问题的改变和数值计算对p e k e r i s 问题解的逼近影响，本文从 下面三个方面考虑： ( 1 ) 剖分加密( 初值) 对特征值的影响 r e ( ) 图5 ：其中j t 表示n 1 = 2 0 0 n 2 = 1 0 0 n 3 = 1 0 0 ，_ 4 = 2 0 0 ，一表示n l = 1 2 0 ，n 2 = 6 0 ，n 3 = 6 0 ，n 4 = 1 2 0 c = 1 0 0 ，m = 1 0 0 ，i i f ( x ) l i = 1 0 2 1 4 三 善 r e ( 凡) 图6 ：其中v 表示n i = 2 0 0 ，n 2 = 1 0 0 ，n 3 = 1 0 0 n 4 = 2 0 0 一表示n 1 = 2 5 0 ，n 2 = 1 2 5 ，n 3 = 1 2 5 ，n 4 = 2 5 0 c = 1 0 0m = 1 0 0 1 i ，( x ) l l = 1 0 一。 通过图( 5 ) 图( 6 ) 比较可见，当划分到6 0 0 等份时，特征值基本上能稳定不变；对于p r o p a g a t i o nm o d e 和l e a k ym o d e ，划分到3 6 0 等份时，特征值基本吻合。 ( 2 ) 深度h b 增加对特征值的影响； 图7 ：其中川0 表示n i = 1 0 0 h b = 24 g = 1 1 5 5 0 3 = 5 0 n 4 r e ( a ) 图8 ：其中。表示n 1 = 2 0 0 ，n 2 = l o o ，2 1 3 = l o o ，n 4 = 2 0 0 ，h b = 3 0 ，表示n 1 = 4 0 0 ，n 2 = 2 0 0 ，n 3 = 2 0 0 ，n 4 = 6 0 0 ，h b = 5 , 0 ，g = 1 0 0 ，m = 1 0 0 ，i l i ( x ) 【i = 1 0 一2 通过图( 7 ) 图( 8 ) 比较可见，随着h b 的增加，越接近实际情况，特征值计算效果更好 不会漏解，但考虑到计算量及误差积累，h b 也不宜太大，取h b = 3 0 左右适当即可。 ( 3 ) a ( z ) 中的c 变大对b e r e n g e rm o d e 的影响 r e ( 兄) 图9 ：n 1 = 2 0 0 ，2 # 2 = 1 0 0 ，n 3 = 2 0 0 ，n 4 = 2 0 0 ，其中g = 1 0 0 ，m = 1 0 0 ，i l i ( x ) l i = l o 哪 r e ( 1 图1 0 ：n 1 = 4 0 0 ，n 2 = 2 0 0 ，n 3 = 2 0 0 ，n 4 = 4 0 0 ，其中g = 1 0 0 0 0 m = 1 0 0 0 0 ，i ( x h l = 1 0 3 由图( 9 ) 至图( 1 0 ) 可见，随着g 的增大，b e r e n g e rm o d e 会逐渐消失 第二个例子中，我们用多重非对称r a y l e i g h 商叠代来计算界面是弯的情况下的特征值： 设0 。 h ( x ) 时p l = 1 0 ，h ( x ) z h b 时m = 1 7 ，它们的界面曲线z = h ( x ) _ 1 一 曲一8 ( 量一 r ，d = 1 5 ，h = 2 0 ，h b = 2 5 ，水平距离l = 1 00 ，善= 0 2 ，子= 1 0 0 处于0 。( h ( z ) 的k = 1 6 0 ，处于 如) 。 日日的k = 1 4 4 ，且口( z ) = g i ，r = 音l i 乇，g r ，以下 n 1 ，n 2 ，n 3 ，n 4 都代表离散的阶数! 首先，我们取c = 1 0 0 ： 当= 1 0 时， r e ( k ) 图1 1 ：。表示多重非对称r a y l e i g h 商叠代的解，其中n 1 = 2 0 0 ，n 2 = 1 0 0 ，n 3 = 1 0 0 ，n 4 = 2 0 0 ，c = 1 0 0 1 7 当= 4 5 时 一 e r e ( 旯、 图1 2 ：7 一表示多重非对称1 h y l e i g h 商叠代的解，其中1 = 2 0 0 ，n 2 = 1 0 0 ，n 3 = 1 0 0 ，n 4 = 2 0 0 ，c =