（应用数学专业论文）波导计算中束传播法的一种改进.pdf

还 i l 北京化工大学位论文原创性声明 本人郑重声明：所呈交的学位论文，是本人在导师的指导下， 独立进行研究工作所取得的成果。除文中已经注明引用的内容外，本 论文不含任何其他个人或集体已经发表或撰写过的作品成果。对本文 的研究做出重要贡献的个人和集体，均己在文中以明确方式标明。本 人完全意识到本声明的法律结果由本人承担。 作者签名：塑垒日期： 关于论文使用授权的说明 州y 。；o 学位论文作者完全了解北京化工大学有关保留和使用学位论文 的规定，即：研究生在校攻读学位期间论文工作的知识产权单位属北 京化工大学。学校有权保留并向国家有关部门或机构送交论文的复印 件和磁盘，允许学位论文被查阅和借阅；学校可以公布学位论文的全 部或部分内容，可以允许采用影印、缩印或其它复制手段保存、汇编 学位论文。 保密论文注释：本学位论文属于保密范围，在上年解密后适用 本授权书。非保密论文注释：本学位论文不属于保密范围，适用本授 权书。 作者签名： 导师签名： l 粥诵 日期： 日期： 川，岁。 p 二 学位论文数据集 中图分类号 0 2 4 1 8 学科分类号 1 1 0 8 7 论文编号 1 0 0 1 0 2 0 0 7 0 6 6 1密级非保密 学位授予单位代码 1 0 0 1 0 学位授予单位名称北京化工大学 作者姓名周娟娟学号 2 0 0 4 0 0 0 6 6 1 获学位专业名称应用数学获学位专业代码 0 7 0 1 0 4 偏微分方程组解的性质 课题来源自选项目研究方向 及数值解法 论文题目波导计算中束传播法的一种改进 关键词 h e l m h o lt z 方程，r i c c a ti 方法，o n e - w a y 论文答辩日期 2 0 0 7 - 0 5 - 2 5论文类型 应用研究 学位论文评阅及答辩委员会情况 姓名职称工作单位学科专长 指导教师黄晋阳 教授 北京化工大学 评阅人l黄飞敏研究员中科院数学所 评阅人2许兰喜 教授 北京化工大学 评阅人3江新华 副教授 北京化工大学 评阅人4 评阅人5 答暇员会蛳黄飞敏研究员中科院数学所 答辩委员1黄晋阳教授北京化工大学 答辩委员2施小丁教授北京化工大学 答辩委员3许兰喜教授北京化工大学 答辩委员4江新华副教授北京化工大学 答辩委员5 注：一论文类型：1 基础研究2 应用研究3 开发研究4 其它 二中图分类号在中国图书资料分类法查询。 三学科分类号在中华人民共和国国家标准( g b t1 3 7 4 5 - 9 ) 学科分类与代码中 查询。 四论文编号由单位代码和年份及学号的后四位组成 tl。lj 摘要 波导计算中束传播法的一种改进 摘要 在声学、电磁学、地震学及其他许多应用领域的大规模的波传 播问题的研究中，通常需要在一个长度比波长还要大许多的区域中求 波场的分布，在激光的传输和光晶仪器的设计等问题研究中会遇到波 数对传输方向变量的依赖较强的情况，目前尚需改进和发展针对这些 问题的算法。本论文在前人工作的基础上，给出了对d t n 算子q 更好 的逼近，得到了更有效的算法。主要的工作有： 1 、讨论了描述波的传导的二维h e l m h o l t z 方程模型数值求解问 题，讨论了求解大范围波导问题的主要的现代方法。利用d t n 算子q 满足的r i c c a t i 方程，给出一种新的更精确的光束传播法( b p m ) 方 法的改进。我们新的算法对波数不用作限制，能满足散射边界条件。 2 、新算法直接采用了指数型算子函数的有理帕德逼近来得到近 似解，可以大大减少计算量和存储量，得到较为精确的结果，是一种 求解大范围波导问题有效的方法。通过引入理想匹配层( p m l ) 新算 法也可求解无限域的波导问题。 关键词：h e l m h o l t z 方程，r i c c a t i 方法，o n e - w a y l i ， 英文摘要 a b s t r a c t f o rl a r g es c a l ew a v ep r o p a g a t i o n p r o b l e m s i n a c o u s t i c s ，e l e c t r o - m a g n e t i s m ，s e i s m i cm i g r a t i o na n do t h e r f i e l d si ti so f t e nr e q u i r e dt os o l v e t h eg o v e r n i n ge q u a t i o n si nad o m a i nw i t hav e r yl a r g el e n g t hs c a l e ，a s c 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 f o rl a s e rt r a n s m i s s i o na n do p t i c a l c r y s t a le q u i p m e n td e s i g n ，u s u a l l yt h ew a v e n u m b e ri ss t r o n g l yd e p e n d e n t o nt h er a n g ev a r i a b l e i ti sr e q u i r e dt oi m p r o v ea n d d e v e l o pa l g o r i t h m sf o r s o l v i n gt h ew a v e g u i d ep r o b l e m s b a s e d o nt h ep r e v i o u sw o r k s ，w eg i v e ab e t t e ra p p r o x i m a t i o nt ot h ed t n o p e r a t o rq ，a n do b t a i nam o r ee f f e c t i v e a l g o r i t h m o u rm a j o r w o r ki sa sf o l l o w s ： 1 t h en u m e r i c a ls o l u t i o n so ft h et w od i m e n s i o n a lh e l m h o l t ze q u a t i o n f o rw a v e g u i d ep r o b l e m sa n dm o d e mm e t h o d st os o l v el a r g e - s c a l e w a v e g u i d ep r o b l e ma r ed i s c u s s e d u s i n gt h er i c c a t ie q u a t i o nw h i c h t h ed t no p e r a t o rq s a t i s f i e s ，an e w a n dm o lea c c u r a t ei m p r o v e m e n to f t h eb e a mp r o p a g a t i o nm e t h o d ( b p m ) i sg i v e n t h en e wa l g o r i t h mc a n s a t i s f yt h er a d i a t i o nb o u n d a r yc o n d i t i o nw i t h o u tr e s t r i c t i o no nt h e w a v e n u m b e r 2 t h en e wm e t h o du s e sr a t i o n a lp a d 6a p p r o x i m a t i o nt o e x p o n e n t i a l o p e r a t o rf u n c t i o nt og e tt h ea p p r o x i m a t es o l u t i o n t h i sm e t h o dc a n r e d u c et h ea m o u n to f c o m p u t a t i o na n ds t o r a g ea n dg i v em o r ea c c u r a t e r e s u l t s i ti sa ne f f e c t i v em e t h o df o rt h ew a v e g u i d ep r o b l e m sw i t h l a r g es c a l e t h en e wa l g o r i t h mc a na l s oa p p l yt o t h ew a v e g u i d e n -，：b-j- ， 英文摘要 p r o b l e m si na ni n f i n i t ed o m a i nb yu s i n gt h ep e r f e c t l ym a t c h e dl a y e r ( p m l ) m e t h o d k e y w o r d s ：h e l m h o l t ze q u a t i o n ，r i c c a t im e t h o d ，o n e - w a y i i i 气li 目录 目录 第一章绪论”1 1 1 研究背景1 1 2 研究的主要内容8 第二章h e l m h o l t z 方程8 2 1 基本方程：8 2 2o n e - w a yh e l m h o l t z 方程一1 0 第三章算子的计算及逼近1 1 3 1r i c c a t i 方法1 2 3 2 有理p a d 6 逼近1 4 第四章改进的算法及实现1 5 4 1 改进的算法1 5 4 2 算例1 8 第五章结论2 1 参考文献2 2 致谢2 4 研究成果及发表的学术论文2 5 作者及导师简介2 6 i v 气lr1 舟 ， 目录 c o n t e n t s c h a p t e r 1i n t r o d u c t i o n 1 1 1b a c k g r o u n do fi n v e s t i g a t i o n 1 1 2m a i ni n v e s t i g a t i o n 8 c h a p t e r 2h e l m h o l t ze q u a t i o n 。8 2 1 2 2 c h a p t e r 3o p e r a t o rc a l c u l a t i o na n da p p r o x i m a t i o n 11 3 1r i c c a t im e t h o d 12 3 1r a t i o n a lp a d 6a p p r o x i m a t i o n 1 4 c h a p t e r 4t h e i m p r o v e da l g o r i t h m a n dt h en u m e r i c a l s i m u l a t i o n 15 4 1i m p r o v e da l g o r i t h m 1 5 4 2c a l c u l a t i n ge x a m p l e - 18 c h a p t e r5c o n c l u s i o n s 2 1 r e f e r e n c e s 2 2 a c k n o w l e d g e 2 4 r e s e a r c ha c h i e v e m e n t 2 5 i n t r o d u c t i o no ft h ea u t h o ra n dt h em e n t o r 2 6 v 0 8 1 lr 0 1 j t 剐 一 叫 e 忱1 i n o 0 h 1 m t o l a e u h q e y a c w 1 上 一 s e a n b o 符号说明 符号说明 声波频率，h z 波长，m 传播方向( 设为水平方向) 距离变量 垂直于传播方向( 设为深度方向) 的变量 表示波数 z 的变化宽度 所要求的h e l m h o l t z 方程的解的传播方向的距离，m 平方根算子 参考波数 算子的特征值 对应于特征值允，的特征函数 z 方向离散节点数 v i 砂 力 国， x z讹d么三屯乃纵 ， 第一章绪论 1 1 研究背景 第一章绪论 在声学、电磁学、地震学及其他许多应用领域的大规模的波传播问题的研 究中，通常需要在一个长度比波长还要大许多的区域中求波场的分布。以海洋的 声波传播为例，在探索开发海洋的过程中，海洋声学起着特别突出的作用，因为 只有声波才能在水中进行较长距离的传播，从而进行水下测距、定位、通讯和遥 感等操作。海洋是一种极其复杂的声学介质【2 3 1 。海洋介质最具特征性的现象是非 均匀性，这种非均匀性强烈地影响着海洋的声场。同时海洋波导的特性随着水平 距离也有变化，但变化非常缓慢嗍。声速随深度有规律变化会形成水下声道。利 用水下声道可以进行数百、甚至数千里的远程声传播，但受海水吸收的限制( 2 5 1 ， 要求用于传播的声波频率国很低，即相应的波长z 要很大才行暖】。 在光通信中集成光学是一个重要的分支，光波导作为集成光学器件中一种 最基本的元件有着十分重要的研究价值。介质光波导是集成光路中普遍存在的。 随着光通信技术的不断发展和集成光路研究的不断深入，非常有必要对各种不同 波导结构的介质光波导的传播模式及电磁场分布进行详细的分析，建立起精确的 数值分析或近似计算方法，以期对光波导的结构进行优化，进一步提高器件的性 能。平板光波导是其中最常见的一种介质光波导。 介质光波导( 简称光波导) 是一种能够将光波限制在其内部或其表面附近， 引导光波沿确定方向传播的介质几何结构。在集成光学领域，介质光波导是在光 波导器件和集成光路中用以限制和传导光的元件。其结构有很多种，从截面来看， 有圆形，掩埋方形，脊形等等。不管什么类型的光波导，介质光波导都是由一个 折射率相对高的区域，与其周围折射率相对低的区域组成，将光波限制在其中进 行传播。光波导从折射率分布变化的大小可分为弱波导和强波导。强波导折射率 差可达5 8 ，而弱波导芯层的折射率只比包层大一点，比如标准的单模光纤芯层 折射率只比包层折射率大约0 3 。 目前对介质光波导的研究有很多解析和近似的方法。一是射线( 或光线) 北京化工大学硕士学位论文 分析法，它是建立在全内反射的原理之上，物理概念清晰，并能得出一些光在波 导中传输的基本特性，很直观。但缺点也很明显，它并不能得出波导中的场分布。 二是利用模式匹配方法( c o u p l e dm o d em e t h o d ) 。这是迄今比较成熟的理论 方法。它利用一正交完备组的模式展开，从而问题的解由相应的一组耦合系数来 表示。但是其公式推导繁琐，且随着模式的增加，计算量比较大。 三是解析法。根据分离变量法来求解偏微分方程，可将解表示为已知函数 的形式，从而计算出精确的数值结果。它还可以作为近似解和数值解的检验标准。 但是它仅能解决很少量的问题，不具有普遍性。平板光波导由于无限的边界条件， 还没有简洁的数学表达式。 四是数值求解波动方程。它是目前最常用的方法，在计算电磁学的各个领 域中得到了广泛应用。求解波动方程的方法可以有有限元方法、有限差分方法、 及积分方程方法等。这类数值方法对复杂问题有较强的适普性。这些数值方法都 属于全局方法，都将导致高阶系统，因此将会占用大量的c p u 内存及计算时间。 若遇上非稀疏矩阵，更在计算空间上碰到了很大困难。目前纯数值法正在发展当 中，大量研究工作都在探讨相应的计算技术及手段。最新的发展动向是研究高效 的并行数值算法。 一般来说，数值求解波动方程仍比其他方法来得更加容易，速度也更快， 而且与其他绝大多数波导器件不同的是，光波的波长很小而且光波导是开放式结 构，因此用有限网格划分空间并引进计算时，需要设置其边界条件。如果不对边 界条件进行处理，将会产生边界反射。在这种情况下，必须引入适当的吸收边界 层条件把计算空间截断，并且应保证在截断边界处只有向外传播的波而没有向内 传播的反射波。吸收边界的效果直接关系到数值计算的正确性和精确性，是影响 计算品质的决定因素。 光波导器件是靠芯层与包层折射率及结构上的微小差异来实现其各种功能 的【2 7 1 。随着工艺水平的不断提高，现在可以实现的器件截面结构以及折射率的分 布情况越来越复杂，这一方面满足了应用需求，但同时也给器件的设计带来很大 的困难，许多器件无法使用解析模型给予准确描述。使用计算机辅助手段，对光 波导进行数值分析，可以直观有效的实现器件的模拟，对光波导的设计起到重要 的作用。在光波导的仿真中，光束传播法以其简单有效，易于应用到复杂结构的 2 第一章绪论 器件而得到广泛的应用。 对光波导的分析方法可以分为两种：解析法和数值法。虽然解析法可以计 算出精确的数值结果，但由于光波导结构的复杂性，只有少量波导可以得到解析 表达式，因此更多的波导分析应用到数值方法。光束传播法( b e a mp r o p a g a t i o n m e t h o d ，简称b p m 方法) 【1 5 ，1 6 】就是一种被广泛应用到光波导仿真中的数值方法。 最初是用来模拟激光在大气中的传播，后来人们进一步将此方法运用到光纤研究 中，并广泛的应用于无源波导器件和光纤中电磁场模式特性的分析中，由此逐渐 发展出很多计算方法。最初的光束传播法都是基于快速傅立叶变换，称为快速傅 立叶变换光束传播法。它源于h e l m h o l t z 方程，因此只能处理一个极化分量，得 到标量场。由于其采用了快速傅立叶变换，也给它带来了先天的不足。首先它的 横向离散化格点宽度必须一致，而且不能取的太小，格点数还必须是2 的整数次 幂。另外快速傅立叶变换光束传播法受最大折射率变化和最大传播步长的影响， 计算效率低，而且在处理复杂波导时精度较差，其边界条件的设置也是相当麻烦 的，但快速傅立叶变换光束传播法在历史上有着重要的意义。 有限元光束传播法的原理是将波导横截面分成很多三角元，在每一个元内 的场用多项式来表达，然后加入不同元间场的连续条件，就可得到整个横截面的 场分布。有限元光束传播法采用的是不规则的节点元，这样在处理复杂几何曲面 ( 如弯曲波导、劈形波导) 时，它比快速傅立叶变换光束传播法或有限差分光束 传播法所采用的方格有效得多。尽管有限元光束传播法有很多优点，但由于很难 得到有限元光束传播法的矢量公式，很少使用矢量有限元光束传播法来计算三维 波导。另外要得到有限元光束传播法在两种介质界面处两个不同元间的介电常数 也存在着困难。 有限差分光束传播法将隐式差分格式应用于基于慢包近似的h e l m h o l t z 方 程，得到的方程是无条件稳定的。由于该算法在计算过程中涉及的矩阵是三对角 矩阵，可以通过优化算法大大降低算法所需的计算机存储量和计算时间，而且由 于该算法是在c r a n k - n i c h o l s o n 条件下离散化的，因此该算法无条件稳定，从而 可以选用较大的纵向的步长。 另外，光束传播方法根据光束传播方向与计算方向之间夹角的不同情况， 可以分为傍轴近似和广角光束传播法。在很多文献中，对于广角光束传播方法出 3 北京化工大学硕士学位论文 现了多种光束传播方法的变型：利用本征函数展开的l a n c z o s 矢量光束传播方法， 以及使用p a d 6 近似算子的光束传播方法等等。在以上各种束传播方法之中，本 征函数展开法是较为复杂的一种，随着器件结构的变换，本征函数都需重新计算， 因而，编程代码复杂，运算量较大。l a n c z o $ 矢量光束传播方法已经被发现在收 敛性上存在着很大的问题。因而，现在最常用的是p a d 6 算子近似光束传播方法， 该方法简单、方便，同以往的傍轴光束传播方法相比，在几乎不增加计算量的情 况下，使计算的精度、稳定性大大的提高，打破了原先光束传播法共同遵守的应 用原则，即要求折射率沿传播的方向上缓变。 若波源是单一频率时间调和的，则二维波场的决定方程是以下h e l m h o l t z 方 程 + + | i 2 ( 而z ) ”= oi ( 1 1 ) 其中工表示传播方向( 设为水平方向) 距离变量，z 表示垂直于传播方向( 设为深度 方向) 的变量，讯矽表示波数。若以d 表示z 的变化宽度，彳表示所要求的 h e l m h o l t z 方程的解的传播方向的距离。典型地，我们有 a d 坛。2 z 即波长z 远小于z 的变化宽度，且传播距离也远大于z 的变化宽度。 对上述h e l m h o l t z 方程，由于k 在水平x 方向上是弱衰减的，故可以把x 方 向分成两部分，即0 石a 和石a ，其中这两部分由水平方向o u t g o i n g 条件连 接【9 】。当石a 时，设k 与工无关，可由分离变量的方法得到方程的近似解。本 文主要考虑的是0 x a 这部分的h e l m h o l t z 方程的求解问题。 当k 是常数，并且定解区域是矩形区域时，以上h e l m h o l t z 方程有精确平面 波解和球面波解【1 0 1 。但实际的大部分问题中，k 并不是常数，而是与水平距离x 和深度z 等有关的f l l 】，故一般无法直接给出如下h e l m h o l t z 方程 + + j i 2 “z 如= o 在0 z a 的解析解。 对于z 方向上的无限深度问题，以往的方法是简单地设立假边界条件将无界 的求解区域有界化，如设这些边界条件为第一、第二或第三类边界条件，这种做 法只是将原方程粗糙地近似化，所以无论是对此近似问题进行精确求解还是离散 方程后进行数值计算，解得逼近效果都较差【1 1 ，产生误差的主要原因在于这样的 4 第一章绪论 边界条件实际上并不符合无界区域中声波传播的特性。 在这种情况下，必须引入适当的吸收边界层条件把计算空间截断，并且应保 证在截断边界处只有向外传播的波而没有向内传播的反射波，吸收边界的效果直 接关系到数值计算的正确性和精确性，是影响计算品质的决定因素。1 9 9 4 年， j p b e r e n g c r 首先提出了高效的二维理想匹配层( p c r f e m z t l ym a t c h e dl a y e r , 简写为 p m l ) 吸收边界条件的概念【2 羽。后来，又在理论上证明了该方法可以完全吸收 来自各个方向、各个频率的声波，而不发生任何反射。现在完美匹配层作为边界 吸收条件的观点它被认为是目前最好的吸收边界条件。这样，在计算边界的周围， 都有完美匹配层作为吸收介质。波由区域内通过边界传播达到完美匹配层时，基 本不会发生反射。波在完美匹配层中传播时，也不会发生反射，并且按传播距离 的指数呈规律衰减。当波传播到完美匹配层的外边界时，波场近似为零，也不会 发生反射。实际上，b e r e n g e r 所假设的完美匹配层内的电磁场并不存在，它是一 种假想的能产生电损耗和磁损耗的各向异性介质。所以完美匹配层自身可以被看 作是一个狭窄的介质层，在这一介质层中，任何角度的入射光都被完全吸收。1 9 9 7 年f r a n c i sc o l l i n o 针对上述添加的p m l ，在数学上相当于对坐标做了一个复的伸 展变换圣= x + i 【c r ( r ) d r ，把h e l m h o l t z 方程转换成一个新的复方程，我们称之 椰 为改进的复h d m h o l t z 方程。 在有界区域上求解h e l m h o l t z 方程有许多直接的数值方法，如有限元和有限 差分法等【2 1 1 。但方程定解区域的水平距离与波长相比非常的大，而且波场的内部 振荡也很大，每个波长范围内部需要用很密的格点( 或基函数) 来表示波场， 对于这个水平距离比波场还要大得多的区域来说，若用有限差分或者有限元方法 来处理这样一个水平区域很大的h e l m h o l t z 方程时，产生的线性系统的阶数将非 常的大，导致相当大的存贮空间，计算的代价也很高昂，同时这些系统常常也是 不定的，或非对称的，这就使得方程的求解更加困难。 以有限差分法来求解以上的h e l m h o l t z 方程为例，离散时分别在x 方向和z 方向取个点和m 个点，若步长相同，即有m ( a d ) ，总的存贮空间的需求 就是o ( a d ) n 2 ，离散所得的系统是非对称的和不定的，并且它的带宽是，即 使用带状矩阵的处理方法来求解，运算所需空间的数量级将是o ( a d ) n 4 。随着 5 北京化工大学硕：l 学位论文 划分的加密，所需的存贮空间将是非常惊人的【2 0 】。 就海洋中声波导问题来说，一方面，在波传播方向上求解区域的长度，要比 纵向( 垂直于传播方向) 的要大得多( 都远比波长大) 【9 】；另一方面，波导介质 沿着传播方向有变化，即和区域是相关的，但是很细微。显然一个好的数值方法 应该充分利用声波导的这两方面的特性，利用波导对区域水平距离依赖很弱的特 性，有很多更合适的方法【2 1 。 c o u p l e dm o d em e t h o d ( c p m ) 是求解h d m h o l t z 方程的很好的方法【6 3 1 。c p m 方法用一系列的水平分片的与水平距离无关的波导来近似整个依水平距离变化 的波导，由原边界条件以及相邻的小片所需满足的一系列衔接条件可得一组耦合 方程。与直接有限差分和有限元方法相比，c p m 方法水平方向步长的选取不受 波长的限制，特别是波导在依区域的变化很小时，c p m 方法能使计算时保证精 度要求的前提下，可以取到比波长还大很多的步长。但当给定的边界是弯曲的时， c p m 方法使用阶梯折线来近似，这使得数值结果的误差非常的大，同时也要求 步长要取得很小，这就失去了大步长的优点【6 】。 一些基于e x a c to n e - w a yr e f o m u l a t i o n s ( e o r ) 变换的近似方法有着比c p m 方法更好的优点【1 让1 2 5 1 。广泛使用的抛物方程近似方法( p a r a b o l i ce q u a t i o n ) 2 4 , 2 6 是多种对h e l m h o l t z 方程单向( o n e - w a y ) 近似中的一种。它们可以看作是对以 下理想单向方程的近似， “善= f 后2 ( x ，z ) + a ：“ 式( 1 2 ) 若将变量x 视作是时间变量，则抛物方程可以用步进方法( m a r c h i n gm e t h o d ) 2 1 当作初值问题来解。当区域的介质是水平方向无关的，即k 与x 无关的话，方程 的解正是满足无穷远处辐射条件的h e l m h o l t z 方程的解【1 9 】。 当k 弱依赖于z 时，为改进( 1 2 ) 的逼近效果，m d c o l l i n s ( 1 9 9 3 ) 提出了 能量守恒改进方法。y y l ua n dp l h o ( 2 0 0 2 ，2 0 0 3 ) 1 5 , 1 6 给出了另一形式的 改进。 1 9 9 3 年f i s h m a n 提出了一个基于散射算子的h e l m h o l t z 方程的步进方法。1 9 9 6 年yyl u 和t r m c l a u g h l i n 提出并实现了基于d i r i c h l e t t o - n e u m a r m 映射( 简称 d t n 映射) q ( x ) 的r i c a t i 方法【4 1 。d t n 映射q ( x ) 将h e l m h o l t z 方程的解“映射成 它的导数略，即蚝= q ( x ) u 。算子q ( 功可以看作是以x 为参数的，作用于定义在 6 第一章绪论 区间( 0 ，d ) 上的z 的函数的算子，它须满足以下r i c c a t i 方程 j _ a i d = - 箧+ 七2 ( z ，z ) 卜q 2 式( 1 3 ) r i c c a t i 方法利用d t n 映射将原h e l m h o l t z 方程的边值问题，转化成一个算子的初 值问题( 将水平方向的坐标变量看作是时间变量) ，然后再进一步用水平方向上 的步进方法来求解。数值实现时，若深度z 仍用刀个点来离散，q 可以用刀刀的 矩阵来近似，总的内存的需求就是o ( n 2 ) 。更进一步可以用局部的特征函数系展 开的方法，只需保存开个特征值中的前个最大的特征值，通常可以取到1 0 n 以。这样原本用t l y 阶的矩阵近似算子q ，可用阶数要小得多的n x n 阶矩阵来 代替，所需的内存空间将更小【1 3 】。而且运算所需的存贮空间只和深度方向的分划 大小有关，和水平距离的大小及分划无判1 4 1 。因此，在海洋这种水平距离远比深 度及波长要大得多的声场环境中，r i c c a t i 方法与传统的直接差分法和有限元方法 相比，是一种非常实用的方法。 1 9 9 9 年，yyl 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 ) ，该方法进一步解决了大步长计算的问题【3 , 1 2 。o l r s m 方法也可 以利用d t n 映射将原h e l m h o l t z 方程的边值问题，转化为一个算子的初值问题( 将 水平方向的坐标变量看作是时间变量) ，然后再进一步用水平方向上的m a r c h i n g m e t h o d ( m m ) 方法来求解。 实际应用中的声场区域的底边，即水下声道的底面边界，通常是弯曲的【1 8 1 。 针对这种弯曲的底边，有两大类的处理方法，一是用“阶梯状 的折线段( 实际 上是不连续的平行线段) 来近似代替原弯曲的底边。如图f i g 1 所示，在每一段 小区间内( 水平方向的) ，视作是平坦的底边，这样原来整个大区域( 弯曲底边 的) 求解问题，就转化成一系列连续的平坦的小区间上的求解问题。这种近似方 法很容易实现，但这也同时带来了一对矛盾。首先相邻阶梯底边的不连续性会带 来误差，特别是水平距离很大的场合，这误差会积累得很大；另一方面，若要很 好地近似弯曲底边，减小误差，另一类方法是通过坐标变换，将带有弯曲的底面 边界的区域变成易于求解的平坦的矩形带状区域( 如图x y z 所示) 。有三种坐标 变换，1 ) 全局变换；2 ) 局部变换( 非正交) ；3 ) 局部正交变换【1 刀。 全局变换是用一种共形映射( c o n f o r m a lm a p ) 拉平原弯曲底边【2 6 1 ，同时也 保持了原方程的简单形式，但对于这种全局变换来说，当波导范围很大，且边界 7 北京化工大学硕士学位论文 很复杂的时候，计算起来非常困难。广泛使用的局部坐标变换相比之下计算起来 要容易得多【1 9 2 1 1 ，但对于含有法向导数的底面边界来说，变换后得到的边界条件 是水平方向和深度方向的偏导数的组合形式【5 1 。水平方向的偏导数将给数值计算 上带来很大难度。 图1 - 1 阶梯近似弯曲底边 f i g 1 - 1a p p r o x i m a t i n gc u r v e db o u n d a r yb yt h es t a i r w a y 1 2 研究的主要内容 本论文的主要工作就是要在前面所述的方法基础上，做适当的改进，来提 高解的精度，减少计算所需的存储空间。在激光的传输和光晶仪器的设计等问题 研究中会遇到波数对传播方向j 的依赖较强，甚至依赖于，的情况，目前的算法 尚需大的改进。我们特别研究了对算子p 的更好的逼近，并且得到了更有效的算 法。 第二章h e im h olt z 方程 我们讨论二维h e l m h o l t z 方程。 2 1 基本方程 8 第二章h e l m h o l t z 方程 考虑矩形区域上二维h e l m h o l t z 方程 l k + + 七2 ( 五力“= o 多e ( 2 1 ) 的定界问题。设其求解区域为0 x 彳，0 z x 式( 3 1 ) c 】，( 石；力= ，式( 3 2 ) 其中j 是单位算子。由关系 r ( x ；功= y ( ，；x + 血) h 石+ 如功式( 3 3 ) 和y ( x + a x ；x ) 的t a y l o r 展开 y ( x + a x ；x ) y ( x ；x ) + a x q ( x ) y ( x ；x ) = i + a x q ( x ) 式( 3 4 ) 我们得到 第三章算子的计算与逼近 !：!：!：i-!：!：-掣=：-】，(】c；】f-ax)q(x)式(35) 因此 掣= 一l r ( x ，；x ) q ( 功 当， k 式( 3 6 ) c ! x 当石 z 时，由固定的这一方程可以被用来计算y o ；力初始条件是 y 0 ；j ，) - i 。 设鹤 为0 到之间的一系列点，满足o = 岛 s 2 j s ，= 三，其中解 是需要得到的。 我们从x = l 到x = o 求解方程 a r ：一r q 式( 3 7 ) d x 和r i c e a t i 方程 睾： + 七z o ，z ) 卜q ：式( 3 8 ) d x i 灭_ i b - j ( s j ，+ 1 ) ，y 表示r ( s p l ；力。在点x = s j + l 从m 开始( 例如y ( j ，+ l ；s ，+ 1 ) ) ， 方程( 3 7 ) 得到的结果为x = 勺。结果被l ，( + 。；勺) 输出到二级存储设备，例如硬 盘。然后，我们在点工= 勺重置h ，在下一个区间( 0 - p s j ) 继续这一求解过程。 当对q 的方程( 3 8 ) 和对y 的方程( 3 7 ) 同时求解到x = o 时，我们得到 q ( o ) 和y ( s j + 。；勺) ，其中j = o ，1 ，p 一1 。则“的结果可以在网格鹤) 中由 y ( s j + i ；z ) = 】，( i ；0 弦鸱z ) 其中= o ，1 ，p 一1 式( 3 9 ) 很容易的计算出。 注意到我们的核心存储是存储q 和y 的空间。算子】，q + 。；_ ) 一旦计算出就 被写入二级存储设备，当计算“在点+ 的值时，才被重新写入核心存储中。如 1 3 北京化工大学硕士学位论文 果需要求解的位置总数不是很大，则我们的方法仍是实用的。 3 2 有理p a d 6 逼近 而的有理p a d 6 逼近是求解o n e - w a y 波传播模型的非常有用的方法1 6 1 ，使 用它我们可以很方便地计算平方根算子岛= 止乒忑丽。g ( x ) = 瓜的标 准 p p p a d 6 逼近为 其中 以渊嘻蔫= 珥p 前1 4 - r ( p ) r 引3 加， 垆2 ( 舞) ，c ；p ) = s i p _ 2 而j t t 坩= 羔加 上述函数列捃p ( x ) ) 对复平面上任意的石都收敛到g 似，除了k 一1 的实数。事实 上，当x - i ，g 俐是一个具有正虚部的纯虚数，而g p ( 功是实的。 下面我们说明在计算o n e - w a yh e l m h o l t z 方程( 2 2 ) 时，用上面标准的有理 p a d e 逼近去逼近平方根算子会导致计算误差。在此后，我们将方程( 2 2 ) 重 新写为 “，= l u = 玩汀酉轧式( 3 1 2 ) 即令：疋i 干i ，其中x ：o ：+ k 2 i ( x 广, z 一) - k 2 ，氏称为参考波数。由( 3 1 0 ) ，f ： 方程( 3 1 2 ) 可以被近似为 1 4 第四章改进算法及实现 蚝= 玩( g ，( x ) m ，式( 3 1 3 ) 现假设波导的波数在x 方向是不变的。若乃 o 是算子a ：+ 七2 ( z ) 的一个衰减 模的特征值，c j ( x ) 是对应于此特征值的特征函数。此时，我们可知方程( 3 1 2 ) 有衰减模解 u j ( x , z ) = 咖( 乃。 我们知道当它沿x 的正方向传播时，它是呈指数级衰减的。若此时我们应用 标准p a d 饿，可知“j ( x ，z ) 办( x 弦乃弦，其中，万，= ( 乃一屠) g - 1 。 由上面分析我们知道此时g ，( 旯，) 是实的，故上面的衰减模的逼近解不会随着它 沿工正方向的传播而衰减，从而造成计算误差。 类似地，我们可以对p ，( x 形式的算子函数关于x 给出p a d 6 逼近。 3 2 改进的算法 第四章改进的算法及实现 我们考虑一个二维h e l m h o l t z 方程模型： + + 七2 ( 易z ) “= 0 式( 4 1 ) 其中七是波数。求解区域为0 x a ，0 z 图5 - 3 u ( l ，矽在h = 2 ( 小圆圈) 和7 = 1 2 5 6 ( 实线) 的图像 f i g 5 - 3f i g u r eo f “佤矽f o rh = 2 ( 1 i t t l ec i r c l e s ) a n dh = - 1 2 5 6 ( s o l i dl i n e s ) 2 0 m 第五章结论 可以看出我们得到的结果还是比较接近的，说明这种改进是可行的。 第五章结论 对大范围波导计算，最重要的就是利用问题本身的特殊性来设计数值方法。 对于波导问题，传导距离比横向长度大得多。波导特征沿主要传播方向( 如水平 方向) 变化非常缓慢。双节点方法和大步长方法( 基于o n e - w a y 的h e l m h o l t z 方 程) 对这些问题都很有价值，因为离散传播方向的步长大小可不受波长限制。 本文主要探讨了描述波的传导的二维h e l m h o l t z 方程模型数值求解问题，探 讨了求解大范围波导问题的主要的现代方法。利用d t n 映射q 满足的r i c c a t i 方 程，给出一种新的更精确的b p m 方法的改进。我们新的算法对波数不用作限制， 能自动满足散射边界条件。新算法直接采用了指数型算子函数的有理帕德逼近来 得到近似解，可以大大减少计算量和存储量，得到较为精确的结果，是一种求解 此模型有效的方法。在数值求解时通过引入理想匹配层( p m l ) 新算法也可求解 无限域的波导问题。 2 1 北京化- t 大学硕士学位论文 1 l 【2 】 参考文献 b a y l m 气g o l d s t e i nci t u r k e le a ni t