（应用数学专业论文）求解bvms离散线性系统的块二对角预处理算子.pdf

摘要 摘要 现今，微分方程在各学科和工程领域有着广泛的应用。b o u n d a r yv a l u e m e t h o d s ( b v m s ) 作为一种无条件稳定的常微分方程离散方法也得到了广泛应用， 而且引起了许多学者的关注。因为由b v m s 离散的常微分方程线性系统通常是非对 称的，大型而且稀疏的，如果使用直接法来求解将耗费大量资源，所以迫切需要 寻找快速求解的方法。前些年，c i r c u l a n t 类预处理技术被提出来用于求解这类 线性系统，而且也取得了良好的效果。这类预处理技术的一大特点是使用了快速 傅立叶变换( f f t ) 技术，该技术能够有效提高计算效率，缩减计算时间。本文 在国内外关于这类线性系统预处理方法研究状况的基础上，从另一角度考虑新的 构造方法，提出了一类块二对角预处理算子。该算子在不使用f f t 技术的前提下， 也能有效提高计算的效率，本文将围绕其构造以及性质进行相关讨论。 首先，分别从基于c r a n k n i c o l s o n 方程、最小条件数和最小二乘法这三个方 面提出这类块二对角预处理算子的三种构造方法。针对每类型，探讨了它们关 于收敛性的问题以及给出了一些相应的定理证明。 其次，对它们进行了g m r e s 迭代计算量的分析，得出这类块二对角预处理技术 在不使用f f t 技术条件下，也能降低迭代的计算量的结论。在与传统的c i r c u l a n t 预处理技术需要o ( m n l o g n ) 的计算量相比较时，其所需计算量仅为o ( m n ) 。 最后，对这三种构造方法得到的预处理算子进行数值实验，并将它们与 t c h a n 、s t r a n g t y p e 和p c i r c u l a n t 这几类较典型的c i r c u l a n t 预处理算子作比 较，验证了新预处理算子的可行性和有效性。 理论分析和数值结果表明，本文提出的块二对角预处理算子对b v m s 离散的微 分方程线性系统有较好的预处理作用。结合g m r e s 方法，能有效加速方程求解的 过程，并在一定条件下使得加速效果更加明显，是计算上行之有效的预处理技术。 关键词：边界值方法；预处理算子；广义极小残量法；c r a n k n i c o l s o n ；条件数； 最小二乘法 广东工业大学理学硕士学位论文 a b s t r a c t t h en u m e r i c a la p p r o x i m a t i o no fs o l u t i o n so fd i f f e r e n t i a le q u a t i o n sh a s b e e na n dc o n t i n u e st ob eo n eo ft h ep r i n c i p a lc o n c e r n so fn u m e r i c a l a n a l y s i s b o u n d a r yv a l u em e t h o d s ( b v m s ) ，t h ed i s c r e t i z a t i o nm e t h o d s ，h a v e b e e nw i l d l yu s e di no r d i n a r yd i f f e r e n t i a le q u a t i o n s ( o d e s ) t h ea d v a n t a g e i nu s i n gb v m so v e rc l a s s i c a li n i t i a lv a l u em e t h o d s ( i v m s ) c o m e sf r o mt h e i r u n c o n d i t i o n a ls t a b i l i t yp r o p e r t i e s i nt h i sp a p e r w ec o n s i d e rt h el i n e a r s y s t e m so fo d e sb yu sin gb v m s t h e s es y s t e m sa r eu s u a ll yu n s y m m e t r ic ，l a r g e a n ds p a r s e r e c e n t l y ， ac l a s so fc i r c u l a n tp r e c o n d i t i o n e r sh a sb e e n p r o p o s e dt os o l v et h e s es y s t e m s t h ec i r c u l a n tp r e c o n d i t i o n e r su s et h e f f tt e c h n o l o g yt oa c c e l e r a t et h ep e r f o r m a n c ea n dt h e yh a v eb e e n s u c c e s s f u l l yu t i l i z e di ns o l v i n gm a n ys y s t e m s i nt h i sa r t i c l e ，w ep r o p o s e an e wc l a s so fb l o c kl o w e rd i a g o n a lp r e c o n d i t i o n e r sf o rs o l v i n gt h e s e li n e a rs y s t e m s f i r s t l y ，b a s e do nt h ec r a n k n i c o l s o n ，l e a s tc o n d i t i o nn u m b e ra n dl e a s t s q u a r e sm e t h o d ，t h eb l o c k l o w e rd i a g o n a lp r e c o n d i t i o n e r sh a v e b e e n c o n s t r u c t e dt oa c c e l e r a t et h eg m r e sm e t h o d ，a sw e lla ss o m eo ft h e i r c o n v e r g e n tp r o p e r tie sh a v eb e e nd is c u s s e d s e c o n d l y ，t h ec o m p u t a t i o n a lc o s t so ft h eb l o c kl o w e rd i a g o n a l p r e c o n d i t i o n e r sa n dt h ec i r c u l a n tp r e c o n d i t i o n e r sh a v eb e e na n a l y z e d c o m p a r ew i t ht h ec i r c u l a n tp r e c o n d i t i o n e r s ，t h ec o s t so fs o l v i n gt h e s e p r e c o n d i t i o n e ds y s t e m sw i t ho u rm e t h o d sa r er e d u c e df r o mo ( m n l o g 刀) t o o ( m n ) f i n a l l y ， n u m e r i c a l e x p e r i m e n t sp e r f o r m e d w i t ht h e s en e w p r e c o n d i t i o n e r sa n dc o m p a r a t i v es t u d yb e t w e e nt h e ma n dt h ec i r c u l a n t p r e c o n d i t i o n e r sd e m o n s t r a t et h ee f f i c i e n c yo ft h e s en e wm e t h o d s t h i sc l a s so fp r e c o n d i t i o n e r sh a sb e e np r o v e dt ob eh i g h l ye f f i c i e n c y f o rt h e1 i n e a rs y s t e m so fo d e sb yu s i n gb v mo nt h eg r o u n d so ft h e o r e t i c a l a n a l y s i sa n dn u m e r i c a lr e s u l t s i ti sr e a s o n a b l eo nt h e o r ya n de f f e c t i r e i i o nc o m p u t a tio n a lp r a c tic e a b s t r a c t k e y w o r d s ：b o u n d a r yv a l u em e t h o d s ( b v m s ) ；p r e c o n d i t i o n e r ；g m r e s ； c r a n k n i c o l s o n ：c o n d i t i o nn u m b e r ；l e a s ts q u a r e sm e t h o d i i i 独创性声明 独创性声明 秉承学校严谨的学风与优良的科学道德，本人声明所呈交的论文是我个人在 导师的指导下进行的研究工作及取得的研究成果。尽我所知，除了文中特别糯以 标注和致谢的地方外，本论文中不包含其他人已经发表或撰写过的研究成果，不 包含本人或其他用途使雳过的成果。与我一同工作的同志对本研究所徽的贡献均 已在论文中作了明确的说明，并表示了谢意。 本学位论文成果是本人在广东工业大学读书期间在导师的指导下取得的，论 文成果归广东工业大学所有。 申请学位论文与资料若有不实之处，本人承担一切相关责任，特此声明。 论文作者签字； 指导教师签字： 年月日 4 1 第一章绪论 第一章绪论 1 1 研究背景 微分方程在应用学科、工程领域，如流体力学、空气动力学、机构力学、航 天航空工业、电子工程、建筑领域等等，有着极其广泛的应用，并且求解大型离 散的微分方程组成为了这些领域解决问题的核心。但往往求解这些问题都需要耗 费大量的计算，从而这些领域迫切寻找使用更简便，消耗更少的数值方法来求解 问题的近似解。 本文主要研究在使用b v m s n l 方法离散如下形式常微分方程初值问题( i v p ) 后 所得线性系统的求解方法， 1y ( f ) = j 脚y ( 于) + g ) ，t o o ，t 】 1y ( t o ) = y 。， ，( 1 1 ) 其中y ( t ) ，g ( t ) ：贸一贸“，y o 孵”，并且j 。吼“”。在通过使用阶的b v m s 离散该问题后，将得到如下线性系统： m y = ( 彳 i m h b o 厶) y = b ， 。 ( 1 2 ) 其中么、召是依赖线性多步法( l m f ) 的0 + 1 ) ( ”+ 1 ) 的矩阵，具体形式将在后 面给出，。是m m 的单位阵，h = ( t - t 。) in 是步长，圆是张量乘积瞳3 。 b v m s 是一类基于l m f ，并且无条件稳定的常微分方程离散方法n 5 6 7 l 。其在微 分方程的初值、边界值等问题中有着广泛的应用聃9 1 m 川。为说明b v m s 的离散形 式，首先考虑在等步长的格点f ，= + 乃，j = 0 ，疗，h = ( t - t o ) i n 上使用隐式 步l m f 作用于问题( 1 1 ) ，则有如下n 一+ 1 个方程组 p - v 一v 口m y ( t 川) = 办氏y ( f 州) s = 1 ，n - ，t + v ， ( 1 3 ) i = - v i = - v 其中v 根据不同b v m ，取不同值，具体参考随后各种b s 的类型介绍或文献 1 ，口，、 尼：是l m f 系数，并满足 兰口，：o ，和兰屈：1 。 i = 0i = 0 ( 1 4 ) 因为初值问题( 1 1 ) 仅提供了y 。初始值，而b 还需要一1 个条件方程组才能完 广东工业大学理学硕1 ：学位论文 整，这就需要通过使用称为“a d d i t i o n a lm e t h o d s ”的方法n 1 进行添加。于是有 卢 a ；川y ( t ，) = 乃罗j 7 y ，) 歹= 1 ，v 一1 ， ( 1 5 ) j = o1 = 0 , u - i一i 口；门y ( t 。一，) = 是夕j 7 y 。一，) 歹= n - r + v + l ，嚣。 ( 1 6 ) i = 0j = o 上两式中的参数口j 川、7 需要被选取为能使( 1 5 ) 、( 1 6 ) 中方程与( 1 3 ) 中方程 具有相同的阶的系数。于是综合( 1 3 ) 、( 1 5 ) 和( 1 6 ) ，将得到( 1 2 ) 式，更详细 的式子如下： m y = ( aoi 。一h b 厶) y = e 1 y o + h ( boi ，) g = b ( 1 7 ) 其中e l = ( 1 ，o ，0 ) 丁贸川，y = ( y o ，y 。) r 贸州咖，g = ( 9 0 ，g 。) r 贸川m ， 彳、召具有如下形式： a = b = 1 口5 ” ： 口 a o o 毹” 尉”1 9 。 口“ - a 0 口“ 0 口5 ”声+ 蚪n 口? 一声“ 口5 口? 0 雕 u o p n pu 麒”肿什? 叫”“ ； p 窖。 p 紫 2 一f：“：乍 ；雌口 第一章绪论 通常彳、b 可以分解为一个t o e p l i t z 矩阵徽1 未找粥旧谭。加上一个秩很小的扰动矩阵， 可以将它们看作是准t o e p l i t z ( q u a s i t o e p l i t z ) 矩阵。根据不同的b v m s ，将有 不同的参数，这里简要归纳一下b v m s 的几种类型，具体参考 4 中描述： ( 1 ) b d f ( b a c k w a r dd i f f e r e n t i a t i o nf o r m u l a e ) 和g b d f ( g e n e r a l i z e d b a c k w a r dd i f f e r e n t i a t i o nf o r m u l a e ) ，它们的通用公式为： “一v 口卅y ( t 阳) = b y ( f ，) ，其中1 ， j = - - v 这一类型很适合刚性问题。 ( 2 ) g a m ( g e n e r a li z e da d a m sm e t h o d s ) 和e t r ( e x t e n d e dt r a p e z o i d a l 肿肿。母j-vy(t)-y(t ，y v 川鸺惟= 烽篡三。 肿 肿一，) = 办以，饥，) ，其中1 ，1 ，= 。，；1 当为偶数时称g a m ，它们适用- y n u 性问题；而当为奇数时称为e t r ，适用于 h a m i l t o n i a n 问题和连续的边界值问题( b v p ) 。 ( 3 ) e t r ，s ( t h ee x t e n d e dt r a p e z o i d a lr u l e so fs e c o n dk i n d ) 它们的通 用公式为： “一v 口m y ( t 州) = h ( f l y ( f ，) + ( 1 - f 1 ) y ( f 州) ，1 ，v = f = v 等，为偶数， _ p 厂+ l ，为奇数。 这一类型与g a m s 有相似的性质，同样适用于刚性问题。 当使用了b v m s 后，如果通过直接法解( 1 7 ) 这一线性系统，将耗费大量的计算。 而 1 中作者虽然给出了b v m s 的许多性质，但并未就( 1 7 ) 这类问题给出高效的解 决方法。这类问题实际是类t o e p l i t z 系统问题，在 1 7 ，2 0 ，2 2 中提出了许多有 关求解t o e p l i t z 系统的技术。在2 0 0 0 年，r c h a n 、d b e r t a c c i n i 分别在文章 3 ，4 中提出使用“循环矩阵n 2 ( c i r c u l a n tm a t r i c e s ，它是指一种每一行最后一 个元素和该行下一行第一个元素相同的t o e p l i t z 矩阵) 作为预处理算子并结合 g m r e s n 3 _ 方法来解决这类型问题。自从1 9 8 6 年开始将循环矩阵作为预处理算子 以来，这类算子被经过深入的研究，得出了许多很好的性质，参考 2 3 ，2 4 ，2 5 3 。 数 数 偶 奇 为 伪 2 0 i 譬字 广东t 业大学理学硕士学位论文 而且它们被证明对许多t o e p l i t z 系统的求解有很好的加速效果，参考 1 9 ，2 6 ， 3 2 。 在文献 3 中，r c h a n 提出使用s t r a n g t y p e 型的循环预处理算子，他证明了 预处理算子是可逆的并且预处理系统的所有特征值除了最多2 r a p 外全为l ，并且 g m r e s 方法最多迭代2 r a p + 1 次。在 4 中，b e r t a c c i n i 介绍了两种循环预处理算 子，第一种是由t c h a n 首先提出的最佳循环预处理算子丁( 其性质参考文献 2 7 ) ， 第二种是他自己研究的p c i r c u l a n t 型的。他证明了在某些特殊的b v m s 系统中 t c h a n 的预处理算子可能是奇异的或者是病态的，而p c i r c u l a n t 型却是可逆的。 以上三种是c i r c u l a n t 类预处理算子中的经典代表，下面给出它们各自的形式。 设循环预处理算子为q ，通常它们都有如下类似的形式： q = a i 。一h bo ，。 ( 1 8 ) 其中才、百是循环矩阵，它们的元素由( 1 - 3 ) 所生成： a = c i r c ( - 瓦j ) ，巧2 c j ，1 ( ，? ) 口一+ c ，2 口川州 f 19 、 b = c i r c ( j ) ，岛= c ( 玎) 色+ ，+ c 肼岛一( 州) ，j = o ，2 ， 这里c 讹( 万，) 表示由系数 万，) 生成的循环矩阵， 而系数 c ( ，z ) ，o = 1 ，4 ；j = 0 ，刀) 是根据不同的预处理算子而不同。另需要注意的是， 由于采用的是阶的b s ，所以当j 时，口，、，都为0 。前文提到的 三种主要循环预处理算子的系数c 研) 如下： ( 1 ) t c h a n 的循环预处理算子t ： “雄) = 气s ( 诈) = 1 一矗， c 川( 功= c ( 刀) 2 鬲j ，= o ，刀 ( 2 ) b e r t a c c i n i 的p c i r c u l a n t 循环预处理算子p ： “疗) 司小( h ) = 1 + 矗， c 川( ，z ) 2 勺，( 聆) 2 焘，= o ，聆，1 r1 ( 3 ) r c h a n 的s t r a n g t y p e 循环预处理算子s ： 4 第一章绪论 c ，l ( ，z ) = c j , 3 ( 力) = 1 ， c j , 2 ( ，2 ) = c j , 4 ( 玎) = 1 ， c f ，j ( 甩) = 0 ， 舢，降j ， 斗爿扎膨 其他情况。 使用这些循环预处理算子的其中一大优势就是可以通过快速傅里叶变换 ( f f t ) 来加速计算过程。这将能大大减少g m r e s 迭代过程中的运算时间，从而让 整个计算过程简化。 1 2 研究内容 对于求解由b v m s 离散微分方程所得形如( 1 7 ) 的线性系统，虽然c i r c u l a n t 类预处理技术拥有很好的效率和稳定性，而且f f t 技术也能很好的加速循环预处 理系统，但计算量将被限定在o ( n l o g n ) 的水平上。于是本文尝试不使用循环矩阵 作为算子而从另一角度构造新的预处理算子，即通过块二对角矩阵算子来逼近这 些线性系统的主矩阵，从而获得高效的加速效果。 本文将从基于c r a n k - n i c o l s o n 方程、最小条件数和最小二乘法这三个角度 去讨论新预处理算子的构造，并结合g m r e s 方法分别给出相关的收敛性分析。这 一类算子在计算过程中可不使用f f t 技术也能使每一次g m r e s 迭代大大减少所需 的计算量。 本文结构如下：第二章将从总体上给出块二对角矩阵类的算子的构造针对各 个不同的算子，并给出它们一些共有的和特有的性质。第三章对这类块二对角矩 阵算子和c i r c u l a n t 类算子的进行了计算量分析，并比较了它们之间的差异。第 四章给出了针对三种方法构造所得算子的数值实验结果，并将它们和c i r c u l a n t 类算子的结果进行比较，得出与理论相符合的结论。最后，第五章是总结和展望。 广东工业大学理学硕上学位论文 第二章预处理算子的构造及性质 本文提出的块二对角矩阵预处理算子如下形式： c = a i m h b o 厶， 其中叠和雪具有如下形式： a =，b2 o q b ( 2 1 ) ( 2 2 ) 这里口、移为待定实数，本文随后从三个不同角度讨论它们的取值。 2 1 理论背景 当使用g m r e s 方法来解决( 1 7 ) 这类主矩阵m 为非对称的线性方程时，通常有 三种形式的预处理方式来加速求解过程：左、右和s p l i t 预处理。但只有前面两 种预处理方式是适用于本文所提出的方法。又因为相较于左预处理算子，使用右 预处理方式所得残差的2 一范数是线性方程的真实范数，不像使用左预处理方式时 所得残差的2 一范数与预处理算子有关，所以使用右预处理算子的g m r e g 方法每次 迭代所得的残量范数不会受到该预处理算子的影响n 3 m9 | 。因此本文将使用结合右 预处理算子的g m r e s 进行随后的分析和计算，该算法如下： 算法2 1g m r e sw i t hr i g h tp r e c o n d i t i o n i n g ( 1 )计算：= 6 一坻，= r o2 和v 1 = r o 8 ( 2 ) f o r j = 1 ，m d o ： ( 3 )计算w ：= m p v ， ( 4 ) f o r f = 1 ，j d o ： ( 5 ) h 打：= ( w ，1 ，f ) ( 6 ) w ：= w 一办，v f ( 7 ) e n d d o ( 8 ) 计算哆+ l ，：= 0 w | i ：和v j + ，：= w h i “j ( 9 ) e n d d o 6 第二章预处理算子的构造及性质 ( 1 0 ) 设圪：= 【v l ，一，1 ，。】，以；魄) l 螂，+ 1 ；l s j 。 ( 1 1 ) 乇- i g y 聊：= a r g m i n y 0 尾一以y | 1 2 ，和x 。：= + p - 1 圪虼 ( 1 2 )如果满足条件就停止，否则设x 。：= x 。并返回第( 1 ) 步继续运算 通过观察上面的算法2 1 ，可以知道如果假设线性方程中：= l 并且迭代算法 的初始值为0 ，也不会影响对于收敛率的分析，所以本文将把这一假设运用到之 后的分析之中。 为保证预处理算子c 的有效性，我们需要对它的可逆性进行研究。假设( 1 1 ) 中的矩阵，。可被对角化，则对于预处理算子c 有如下可逆性定理。 定理2 1如果柳女( ，。) l ，其中，九( ，。) ，( 后= 1 ，m ) 分别是，。的m 个特 征值，则预处理算子c ( 2 1 ) 是非奇异的。 证明：由( 2 2 ) 中的a 和雪都是二对角矩阵这一性质，容易得知它们的特征值 分别为： 九，= 1 ，= o ，”， 九o = 0 ，旯亩= ，j = 1 ，刀， 于是得到c 的特征值旯腑( c ) 为 fk ( c ) = 1 ，= 0 ，k = 1 ，m 【( c ) = l h f l a k ( j 。) ，歹= 1 ，n k = 1 ，m 所以当蟛魄( ，。) 1 时，见政( c ) 0 ，= o ，刀，且七= 1 ，m 。定理成立。 因为，。的特征值一般都分布在左半轴，所以定理的假设是合理的。由以上定 理可以知道预处理算子c 的可逆性，然后我们需要对具体的构造方法分别进行研 究。 2 2 基于c r a n k - nic ois o n 的构造方法 从使用低精度b v m 矩阵作为高精度b v m 线性系统的预处理算子这一思想出发， 于是我们想到使用c r a n k - n i c o l s o n 方法将非常适合于构造块二对角矩阵预处理 算子。首先考虑2 阶l m f 方法用于前文提到的等步长的格点f ，= + 办， = 0 ，n ，h = ( t - t o ) i n 上，则有 y j + i 一：昙( y + l - i - y y ，)(23)j ，一 2 i ( y + ly j ) 3 ) 7 广东工业大学理学硕十学位论文 其中y 是y ( t ) 的近似估计，y j = j m y + 毋，g j = g ( t ，) 。上式( 2 3 ) 也称为 c r a n k n i c 。l s 。n 方法 3 1 ，如将这一方法运用到( 2 ，2 ) 中则有a = 声= 圭，即得该 预处理算子形式为 c = l - 0 11 l 圆，。一h r l2 2 11 i l 22 为分析的简便起见，这里仅考虑问题( 1 1 ) 的齐次情况， 到较为简单的形式： f 少( f ) = f ( t ，y ( f ) ) = j m y ( f ) ，f ( ，丁 【y ( ) = y o ， 假设其中的函数足够光滑。当使用b v m 离散( 2 5 ) 后得 即g ( t ) = 0 ，于是得 ( 2 5 ) m y = ( 么 i 。一柚 。) y = e l y o = b 。 ( 2 6 ) 如果使用书 1 中c h a p t e r l 0 7 的假设，则我们将可以导出随后的结果。首先如 书中所设夕= ( y ( t 。) ，y ( t 。) ) 丁是在前文所述等距格点 t ，) 上的严格连续解，并且 使用阶的b v m 方法离散问题( 2 5 ) ，于是可得 m 夕一b = 以，忆k = o ( 办”1 ) 。 ( 2 7 ) 从书 1 中还得到如下十分重要的定理，它将对随后的定理其关键作用。 定理2 2 n 3 成立等式m m ；1 f p = f p + o ( h p + 2 ) ，其中a p 。 因为这一节( 2 4 ) 中预处理算子c 是基于c r a n k - n i c o l s o n 方法，它同时又是 2 阶l m f ，于是有 c 夕一6 = f c ，忆忆= o ( h p + 1 ) 。 ( 2 8 ) 需要注意的是( 2 7 ) 中的右端项b 并未改变，与( 2 6 ) 中的相同。接着我们给出如 下的定理来说明该方法的有效性。 定理2 3 通过使用带右预处理的g m r e s 算法求解由阶的b v m 方法离散问题 ( 2 5 ) 所得线性方程组后，将有如下关系式： 0 ( m c i l - i ) b | i 。= d ( 乃3 ) 。 证明：根据定理2 2 并结合( 2 7 ) ，( 2 8 ) ，可得 8 第二章预处理算子的构造及性质 ( m 卢c 了一，) 6 =m 卢c 二1 b b = m 。多+ m 。c ：t c b = b + r 口+ r c + d ( j 1 2 3 ) 一b = f “+ f c + d ( 乃3 ) 于是，定理得证。 对使用这一构造方法所得预处理算子，我们对于m 固定而”一( 或h 0 0 ) 的情况更感兴趣，而且也得出了下面的一些结果。通过使用定理2 3 ，我们将可 能估计使用带预处理g m r e s 算法一次迭代后的残量范数和收敛率。因为如果这一 收敛率足够小，整个计算过程将可能因此在随后较少的迭代次数中达到要求的停 止标准，并且数值结果也展示了这一效果。为分析这一收敛率，首先需要再次参 考算法2 1 的过程。在这个使用g r a m - s c h m id t 正交化的g m r e s 中，一般通过使 用平面旋转矩阵将正交化后所得h e s s e n b e r g 矩阵转变为上三角矩阵，具体参考 - 1 3 ，c h 6 5 3 。该平面旋转矩阵形式如下： c j s j s j c j 按照书中所说，我们假设s j ，c ，是第个平面旋转矩阵的系数： ( 2 9 ) 且它们满足c ；+ s ；= 1 ，其中y ，是歹次迭代后h e s s e n b e r g 矩阵对应的右端项的最 后一个非零元。根据 1 3 中p r o p o s i t i o n6 9 ，i y ，l 等于歹次迭代后残量模，而且 这一性质对于使用右预处理算子的g m r e s 方法也适用。另从 1 3 ，c h 6 5 3 中可以 得到如下关于y ，的关系式： ? j2 一s3 y j - 1 于是b l 可看成是次迭代后的收敛率。从s j 的公式( 2 9 ) 和前面所列算法2 1 可 知，第一次迭代后收敛率为 9 燕蔬 广东工业人学理学硕士学位论文 卅焘2 弦h 2 丽1 。 ( 2 l o ) 定理2 4通过使用带右预处理的g m r e s 方法结合基于c r a n k n ic o ls o n 的 预处理算子求解( 2 6 ) ，则第一次迭代后得收敛率蚓为o ( 4 - 磊h 3 ) 。 证明：由算法2 1 可知， h ：，= lw 。0 ：= im cn 。v 。一( m cn 叫v 。，v 。) v 。 从几何角度上看，因为v ：垂直于1 ，。( v 2 上v 。) ，所以向量脚m - ( m c n l m ，v 1 ) 的 模比任意其它形如m c v 1 一册1 的向量的模都小，其中巧是常数。于是根据定理2 3 ，可 得 | i 加 一( 埘；l v l ，v 1 ) v 。小i l ( 脚i 1 - i ) v 。0 ：= o ( 4 7 a h 3 ) 。 又因为| i v 。8 ：= 1 ，则有 0 埘刘l ：i ：一l l ( 加二1 一咖，1 1 2 = l o ( 4 - 荔h 3 ) 。 最后，由( 2 1 0 ) 可得 s l 2 南怒 蚓一o ( 4 - 磊h 3 ) ，( 办一o ) 。 所以，定理得证。 2 3 基于最小条件数的构造方法 和2 1 中基于c r a n k n ie o ls o n 的构造方法不同，这里的方法是通过确定( 2 2 ) 中口、的数值使得条件数t c ( m c 。1 ) 尽可能小。由 3 7 中定理1 和 3 3 中各收敛 定理可知道，g m r e s 方法的收敛率和主矩阵的条件数以及其谱分布有密切关系， 即其谱分布越是集中在“1 的附近，g m r e s 方法将收敛越快。所以这节的方法就 是基于这一结论，选取适合的口、作为参数。然而，直接求取这样的口、是 比较困难得。我们这里针对使用基于3 阶l m f 的g a m 方法离散的问题( 1 1 ) 进行说 明，其得到的形如( 1 7 ) 的线性系统，其中a 、b 的形式为： 1 0 第二章预处理算了的构造及性质 a = ，b = oo 58 l 1 21 21 2 581 1 21 21 2 581 1 21 21 2 185 1 21 21 2 这里先给出本节所得的较优口、值结果，并在随后讨论它们的选取过程祁 预处理系统的收敛性分析。当聊”时，取口= 主，= 西8 ；当聊 0 ，其中m i n 是厂的最小值，则有歹( ) t ( 1 f ) 。 引理2 2 n 2 3 如果丁( ) 是h e r m i t i a n 矩阵，则有p ( 刮i ： _ m a x ，其中朋缸，是厂 的最大值。 由上面的两个引理，可以得到以下定理，其揭示了瞥l 与k ( 脚- 1 ) 关系： 定理2 5 1 i mk ( m c 。1 ) m a x i ，q i m i n i m l i 。 _ o1，1 证明：因c 日c 是h e r m i t i a n 矩阵，且易知h 2 是大于零的实函数，所以通过引 理2 1 和引理2 2 可得，当n 寸o o 时， ( m c _ 1 ) ( 埘- 1 ) = m ( c 叼。1 m 川， 于是有忙( i 厂g j 2 ) l l ：朋q ，吲：，即得! 觋j | 脚。14 ：朋巩i ，吲。 同理有! 觋l l ( 脚。1 ) 。1 ：m a x i 。卅= 1 1 m i # 1 ，。l 。所以有 ! i m x ( m c 一) m a x i s q i m i n i 厂q i ，于是定理成立。 有了定理2 5 ，我们将可以通过计算旨i 的界来估计条件数k ( 脚一) ，并对这一 方法构造的预处理算子进行收敛性证明。由( 2 1 2 ) 、( 2 1 3 ) 可得 阱 ( 可1 + 1 ) + 告_ + 8 - z ) ( 一z 一1 + 1 ) + h j ( a z 一1 + 7 ) ( 一1 + 百5 m ) z - 1 + ( 1 + 百8 h j ) 一型1 2 z ( 一1 + a m ) z 。+ ( 1 + 村) 直观地，当令口= 主时，将一定程度的简化上述方程： 阱 1 + ( ( 丧卅一五1 z ) h j ( 一1 + l z z 。1 + ( 1 + 舢) + ：w ) 如前面提到的，我们希望选取适当的，使矩阵m c 一的谱分布集中在“1 附近。 1 2 第二章预处理算予的构造及性质 这也就是需要选取适当的卢，使得! ! 毒二竺! 盘三! 竺一的值分布在零附近。由 这也就是需要选取适当的，使得争j 互一的值分布在零附近。由 ( - 1 + 云彬) z - 1 + ( 1 + 卢：柑) 于过大或者过小都可能使得分布远离零，所以这里我们选取= 西7 和= 熹进 行讨论。 当= 三1 2 时， 阱 1 + ( 1 - z ) 1 篙2 - ( 一l + 瓦5 z - i + ( 1 + 老彬) ，因分子( 1 一z ) 等实部大 于等于零、分母( 一1 + 言z _ 1 + ( 1 + 云彬) 实部大于零，所以z 分别取“l ”、“一l 时旨i 有最小、最大值，即，i 吾l 1 + 1 2 五+ ，h j 。 当= 扣阱 1 + h j z 1 2 ( - l + 5 l zh j ) z 一1 + ( 1 + 8 1 2 h j ) ，同理可知，z 分别取“1 ”、 “一，”时瞥i 有最小、最大值，即一i 1 ：1 号1 + 一h j 。 从上面的分析可以看出，当m _ 。时， 取= 西7 ，l 号l 有较窄的界，相应得 条件数磐妥k ( 脚。1 ) = 1 。而当彬_ 时，取= 西8 ，l 号l 有较窄的界，相应的条 嗍驰删，：莲：詈。所大致可愀定，抛砒取 口= 主，= 西8 ；当聊 甩时，取口= 西5 ，= 云。 要注意的是由于上面的分析都是针对简化的问题( 2 1 1 ) 进行的，与问题( 1 1 ) 有一定差距，所以得出的结果也只是大致的。对于本节中所给出的参数还不是最 优的参数，而是较好的参数，但其已经能大大提高g m r e s 方法的计算速度，这在 广东t 业大学理学硕上学位论文 后面的数值结果中可以看到。而如何取得最佳的口、参数值还有待以后的进一 步研究。 2 4 基于最小二乘法的构造方法 这一节里，我们延续上一节中构造方法的思路，也是将通过确定( 2 2 ) 中a 、 罗的数值使得条件数t c ( m c 1 ) 尽可能小作为目标，而所采用的方法将较为简便， 即使用基于最小二乘法来构造预处理算子。它的适用性较广，可以适用于多种不 同的b v m 方法上，但所取的口、罗值也只是较优的而不是最优的。 为说明该构造方法的过程，我们采用4 阶e t r ，s 对问题( 1 1 ) 进行离散后的所 得形如( 1 7 ) 的线性系统为例进行说明，其中么、b 矩阵的形式为： a = 1 1 79 2 42 4 19 1 21 2 991 1 21 21 2 991 7 1 2 42 42 4 ，b = 同样地，因当1 7o 时，m 、c 都可看作t o e p l i t z 矩阵，所以它们的生成函 数 3 2 ，3 4 ( g e n e r a t i n gf u n c t i o n ) 可分别定义为 厂= z + 六= _ z - 2 - - 9 z - 1 + 9 + z ) + _ h z j - 1 + 1 ) ， q = 吼+ 9 2 = ( 一z 一1 + 1 ) + h j ( a z 一1 + ) 。 由定理2 5 可以看出当瞥l 的界越窄，预处理系统的收敛性越好，又因 其中 。 g9 1 ( f l q 2 一一厶) 9 1 ( q l + q 2 ) 1 4 ( 2 1 2 ) ( 2 1 3 ) ( 2 1 4 ) 一2，一4 ，一2 3 4 o 3412 0一4，一2 ，一抖。一眩。，一压。一他。 9 一m 9一屹 ，一住一弭，一住。一抖 第二章预处理算子的构造及性质 妻鲨茅=z-1+lo+zq 1 21 ) 一= 一一= - 1 ( 一z 叫) 。 则容易得卧如。为取得蚓较窄的界，直观得需要鲁一 越小越好，因此， 我们想通过调整( 2 2 ) 中参数口、来达到这一目的。由于 即有 同理有 鲁= 彬芒1 2z - 2 + ( 坐1 2 拦1 2 矿1 + ( 旦1 2 + 坐1 2 ) + 旦1 2 m c 7 1 、。 。 丝：m ( z - 2 , z - l , 1 ，z ) 9 1 = 办，( z - 2 ，z - 1 , 1 ，z ) 1 1 2 l o 1 2 1 1 2 0 o 1 1 2 1 0 1 2 1 1 2 协 于是可以想到使用最小二乘法求解如下方程： 降1 衢。搋。对1 衢1 2 | l 0 _ 1 黝酬喾 ( 2 1 5 ) 并得出口= = 百3 3 。 以上是该方法对于使用4 阶e 豫：s 时的构造过程，实际上对于任意b v m s ，只要 1 5 o 1，j 2 2 0 0 ，l 口 。l 1j 2 眩2 o 儿儿o吖叭 广东工业人学理学硕士学位论文 按照上面的步骤，得出类似( 2 1 5 ) 的方程，就可求出相应的口和值。由于2 3 节中基于最小条件数的构造方法较难求出所需的口和值，相比较而言，本方法 能更为简便的求出所需的口和值，达到使用较小条件数的目的。 1 6 第三章计算量分析 第三章计算量分析 3 1 块二对角类预处理算子计算量分析 由于本文给出的三种构造方法所得预处理算子的结构都是块二对角矩阵，所 以它们的计算量是一样的。对于使用右预处理算子的g m r e s 方法中的每一次迭代， 主要计算量集中在这一过程中： m c y = ( 4 i m h b 圆厶) ( 彳 i 。一h b 圆d 。) - 1y m c 爿y 的计算可分解为两步。首先，计算c y ，因其形式如下： c _ y = i m c l c 2 c l c 1 c 2 其中c 1 = 一l 一砌y 。，c ：= 一，。一届u 。，所以这一过程可通过使用块的直接算法 求解，实际需要进行胛次阶为m 的矩阵向量乘法和解玎个阶为m 的线性方程组。 由于问题( 1 1 ) 中。的不确定性，在解这露个线性方程组时参照 3 ，4 采取直接 法，则所需的计算量为o ( m n l ) 口0 | ，其中，为，。中非零元的个数。 然后第二步是进行关于m 的矩阵向量乘法，而m 又是一个块的对角矩阵， 于是一次迭代所需总计算量为 o ( m n l ) + d ( 刀( + 1 ) z ) 。 ( 3 1 ) 3 2cir c uia n t 类预处理算子计算量分析 类似的，对于c i r c u l a n t 类预处理算子，主要计算量也是集中在如下过程 m q y = ( 4 i m h b o 厶) (