（计算数学专业论文）求解大型方程组的对称反对称分裂预条件迭代法.pdf

复旦大学硕士毕业论文 摘要 文章首先介绍了当今在求解非对称、正定的大型稀疏线性系统a x = b 中常用的g m r e s 算法，以及实际计算中必不可少的预条件技术，之后分析了基于系数矩阵对称反对称分裂 的预条件( h s s ) ，最后通过一个具体问题的求解，验证了预条件的效果，并且实验证明， 与其他相关预条件相比，无论对于对称部分占优，还是反对称部分占优的情况，对称反对 称分裂预条件都有较好的效果 关键词；对称，反对称，分裂，g m r e s ，迭代，预条件 量里盔堂堕坐些堡窒 2 a b s t r a c t w ei n t r o d u c eg m r e s ，t h em o s te f f i c i e n ti t e r a t i v em e t h o df o rt h e l a r g es p a r s eu n s y m m e t r i c p o s i t v ed e f i n i t es y s t e mo fl i n e a re q u a t i o n s ，a n dt h e n e c e s s a r yt e c h n o l o g y , p r e c o n d i t i o n t h e n w es t u d yt h ep r e c o n d i t i o nb a s e do nt h eh e m i t i a na n ds k e w h e r n f i t i a ns p l i t t i n go ft h ec o e f f i c i e n t m a t r i x ( h s sp r e c o n d i t i o n e r ) ，i e ，t h es y m m e t r i ca n du n s y m m e t r i c s p l i t t i n gi nr e a lf e t i d s ，a tl a s t 。 n u m e r i c a le x a m p l e sa r ep r e s e n t e dt oi l l u s t r a t et h e e f f e c t i v e n e s so fh s s p r e c o n d i t i o n e ri na d d i t i o n ， c o m p a r e dw i t ho t h e r s ，h s sp r e c o n d i t i o n e ri se f f i c m n tw h e t h e rt h es y m m e t r i c p a r to f t h ec o e 毋c i e n t m a t r i xi sd o m i n a n to rt h eu n s y m m e t r i c p a r td o k e y w o r d s ：s y m m t r i e ，u n s y m m e t r i c ，s p l i t t i n g ，g m r e 8 ，i t e r a t i v e ，p r e c o d i t i o n 复旦大学硕士毕业论文 引言 随着科学和技术的发展，越来越多的计算问题需要求解形如 a o = b 4 的大型线性方程组的问题。这时的矩阵不一定是对称的，但在大多数情况下，a 是稀疏矩 阵，即a 的大多数元素为零，且是正定的。求解这类问题已成为数值代数界研究热点之一。 由于y s a a d 等人一系列的工作( 参看文献 1 6 1 9 1 ) ，在求解大型非对称线性方程组的理论 和算法上都得到了重要的结果。从g a r l e k i n 原理出发，得到的a r n o l d i 算法，g m r e s 算 法，成为目前发展的诸多算法的基础。 由于实际情况中，方程组一般是病态的，迭代法收敛很慢，因此在具体求解这样的大 型方程组时，预条件技术是必不可少的，近年来，针对不同的问题与模型，已提出了多种多 样的预条件。本文把自中治、g h g o l u b 和m i c h a e l k n o 提出的对称反对称分裂迭代 1 这一思想，作为预条件应用到g m r e s 算法中去，通过对下述具体问题的求解， 一p u + f w 9 r a d ) u = f i nq 验证了h s s 预条件的效果。 文章的章节是这样安排的，第一章简要介绍了g m r e s 算法和预条件技术；第二章首 先讲了白中治等提出的h s s 迭代，并作为预条件实现到g m r e s 算法中去；第三章用五点 中心差分格式离散化所要求解的微分方程；第四章数值实验，验证了h s s 预条件的效果及 与其他预条件的比较。 复旦大学硕士毕业论文 第一章g m r e s 方法及预条件技术 我们知道，对一个大型稀疏的对称的正定线性系统来讲 a z = b 5 最有效的方法就是共轭梯度法，再结合有效的预条件方法，而对于更一般的对称系统，p a i g e 和s a u n d e r 4 中提出了m i n r e s 方法，即利用对称l a n c z o s 基来计算一个方程的近似解， 并使得残向量最小( 在k r y l o v 子空间上) ，1 9 8 6 年，y s a a d 和m h s c h u l t z 进一步推广 了此方法 3 ，提出了g m r e s ，用以解决非对称线性系统，类似于m i n r e s 方法所基于的 对称la i l c z o s 方法，在g m r e s 中则是基于a n o r l d i 方法。本章分别介绍g m r e s 方法和 预条件技术。 9 1 i g m r e s 方法介绍 1 1 1a r n o l d i 方法 a r n o l d i 方法【5 】是通过g r a m s c h m i d t 方法计算k r y l o v 子空间兰s p a n ( v 1 ，a v 小。u 1 的正交基， u 1 ，7 3 2 ，u k ) ，其方法如下， 算法1 a r n o l d i 方法 1 、首先选取一个范数为j 的初始向量”。，即1 m 1 1 1 = 1 2 、迭代计算，f o rj = 1 ，2 ，d o j = ( a v j ，魄) ，i = 1 ，2 ，j j = a v j 一h i ，j ， t = 1 = jj o j + l 扎 = o j + l h j + 1 ，j 1 1 2f o m 方法 若记k 是一个n x 七的矩阵，其每列是正交基 u ，啦，魄) ，即u = u 1 ， 。，仇) ， 那么，h k = 昭a k 是一个上h e s s e n b e r g 矩阵，其每一个元素h i ，j 就是a r n o l d i 方法中每步 k 阱帆 复旦大学硕士毕业论文 所产生的。为了用正交基k ，通过g a l e r k i n 方法来求解系统( 1 ) ，我们寻找一个如下形式 的近似解，z b = z o + ，其中，z o 是z 的初始解，z k = 8 p a n l ，a v l ，a 1 u l ， 此时残向量r o = 6 a z o ，假设经过以_ u 1 = t o l l r o l l 为初始向量的k 步a r n o l d i 过程后， 由g a l e r k i n 条件，即“= b a x t 正交于k ，这样，= k 玑，y k = h f l i l r o 慨，其中 e 1 = ( 1 ，o ，0 ) 7 ，这样就有了f o m 算法，现描述如下， 算法2 f o m 算法 1 、首先选取z o ，并计算r o = b a x o ，v l = r 0 l l r 0 1 l 2 、迭代计算，f o rj = 1 ，2 ，kd o ： k ，= ( a v j ，) ，i = 1 ，2 ，j ， 0 升l h i + l 。 v j + l = a v ，一h i ，j u ， = lj 吗十1 1 t ， = _ d j + t h j + l # 3 、得到解c = x 0 + k 珊，其中，玑= h f l r o 忙l f o m 算法的缺点是当步数k 增大时，运算量以及存储量都显著增大，且当凰奇异时， z 女不能计算。 1 1 3g m r e s 方法 现在再由a r n o l d i 方法，在经过k 步后，产生一组正交基k + i 和矩阵凰r ( + 1 ) ”，玩 与上h e s s e n b e r g 矩阵日k 相比，除了最后一行( 0 ，0 ，十1 ，k ) 以外，都是一样的，向量u ； 以及玩满足以下一个重要的关系， a k = k + l 凰( 2 ) g m r e s 的迭代解是求以下最小二乘问题的解， 翼骢怕一a 陋。+ z 1 1 1 = m i n 。l i t o a z l l ( 3 ) 设z = u ，那么以上将要最小化的向量的范数就是如下函数， j ( y ) = l l 肋- 一a y k ( 4 ) 复旦大学硕士毕业论文 为了方便起见，总是令卢= ll r o l l ，这样 g ( y ) = l l k + 1 卢e l t l k y e 1 = ( 1 ，0 ，o ) t 由于u + 的列向量是正交向量基，所以， j ( y ) = 怕e 1 一反目 这样最小二乘问题的解就是 = x o + k y k 其中，弧r 使得j ( y ) 取得最小值 以上算法就是g m r e s ( t h eg e n e r a l i z e dm i n i m a lr e s i d u a lm e t h o d ) 算法，描述如下 算法3 g m r e s 算法 、首先选取z o ，并计算t o = b a x o ， 1 = r o l l r o l i 2 、迭代计算，f o rj = 1 ，2 ， ，u n t i la a t i s f i e dd o ： v j + l = o 件l b 札。 7 ( 5 ) 3 、得到解t 。k = 2 7 0 + k 讥，其中，弧使( 5 ) 式达到最小。 以上所描述的g m r e s 算法与h o m 算法具有相同的缺憾，即当k 增大时，运算量显 著增大，为了解决这个问题，可以采用重开始的技术，表示为g m r e s ( m ) ，每计算m 步 后，就重新开始计算，算法描述如下， 算法4 g m r e s ( m ) 算法 、首先选取z o ，并计算r o = b a x o ，= r o l l r o l l 2 、迭代计算，f o rj = 1 ，2 ，md o ： = 咄，试， 一 一 川吩懈 | | = = 吼 卅 复旦大学硕士毕业论文 = ( a v j ，地) ，i = 1 ，2 ，j 3 、得到解： z k = 。o + v k y ，其中，y ker ”使( 5 ) 式达到最小。 4 、重新开始： 计算= b a x 。i 若满足条件，迭代结束，反之， z o ：= x m t v l ：= t m j l r m | lg dt o2 5 1 2 预条件技术 8 共轭梯度法和g m r e s 方法等在理论上都是直接法。若方程阶数是n ，则x 。必为精确 解，但实际上由于舍入误差的影响，x 。往往不能满足原方程。另外，对某些迭代法，如常 用的共轭梯度法，当系统的系数矩阵a 的条件数较大时，它收敛的将很慢，这也是为什么 自1 9 5 2 年共轭梯度法诞生直到7 0 年代一直无法用于实际计算的原因。改进这一问题的方 法是预条件技术近年来，各种新的预条件不断提出和完善，预条件技术已成为了解大型 稀疏方程组的必要技术，但是到目前为止，预条件作用的机理仍未能全部为人所知，构成 了当今科学计算领域的一个难点。很多预条件技术都与方程的物理背景或结构有密切的关 系本文后面所提出的预条件也是针对某一特定模型的。 对于系统( 1 ) ，预条件技术就是找一个“近似于”a 的矩阵q ，要求q 是非异的， 用q _ 1 左乘原方程两边，得到其等价的方程， q 。a x = q - 1 b( 6 ) 如果q 。a 的条件数很小( 对于共轭梯度法) ，或特征值分布在复平面内尽可能小的区域内 ( 对于k r y l o v 子空间方法) ，那么将大大改进迭代的收敛性。 预条件技术结合迭代法就形成了新的迭代方法，如预条件共轭梯度法。 ， ，：l， 一 “ 0 什m 忪 | | = “ 奶 塞里盔堂塑主望些篓塞 9 算法5 预条件共轭梯度法 1 、首先选取3 ；0 ，并计算7 0 = b a x o ，2 0 = q 1 7 0 ，p d3z o 2 、迭代计算，f o rj = 1 ，2 ，直到收敛d o ： n j = ( q ，z j ) ( a p j ，p j ) ， x j + z = z 3 + a i p y ， r j + l = r j a j a p j ， z j + 12 口_ 。r j + 1 ， 岛+ 1 = ( z j 十1 ，z j + 1 ) ( t j ，) ， p j + i 2 z j + l + 83 一l p j 相对于标准共轭梯度法，在预条件共轭梯度法中，每步都要求解形如 q z = r ( 7 ) 的方程因此除要求q 尽可能靠近a 之外，还要求方程( 7 ) “容易”求解。另外，实现 一个预条件关键就是如何快速求解方程( 7 ) 。下面介绍在求解方程( 7 ) 时常用的两种技 术：不完全的c h o l e s k y 分解和不完全的l u 分解。当预条件矩阵q 是对称正定矩阵时，可 以用不完全的c h o l e s k y 分解， 0 = r r 7 它的算法是， 算法6 ( 15 ) 不完全c h o l e s k y 分解 m = q ：f 2 f o r i = 1 ：n 如r j = 1 ：n 一1 i ，q ”= 0t h e n i j = 0 e l s e n j = ( g 巧一高r i k q k ，) q j e n d n ，= ( 吼。一jr 丧) 1 2 复旦大学硕士毕业论文 1 0 e n d 利用上述不完全的c h o l e s k y 分解技术得到的r 和q 有完全相同的稀疏性结构，无论 从储存空间的需求和计算量的角度，这种分解对求解方程( 7 ) 都十分有益。 不完全的l u 分解i l u ( i n c o m p l e t el ud e c o m p o s i t i o n ) 若a 是一个n xn 的矩阵，而 且他的顺序主子矩阵都是非奇异的，则一定可以表示成a = l u ，其中l 是主对角线元素 均为1 的下三角矩阵，u 是一个上三角矩阵，如a 是个大型稀疏矩阵，那么l 和u 的一 般来说都是满的矩阵，所以这样分解理论上虽好，但实际上却行不通。而不完全的l u 分 解是 q l u ，三和d 与l 、u 有相同的形状，但是他们与a 的下上三角部分有完全相同的非零结构。 下面给出用m a t l a b 语言写成的i l u 算法。 算法7 ( 【1 5 ) 不完全l u 分解儿u 知ri = 2 ：n ，0 r k = 1 ：n 一1 i f a ( i ，) 0 a ( i ，) = a ( i ，k ) a ( k ，k ) ，0 r j = k + 1 ：礼 i f a ( i ，j ) 0 a ( i ，j ) = a ( i ，j ) 一a ( i ，) a ( k ，j ) j e n d e n d e n d e n d e n d 分解后的d 占用了a 的上三角部分，三占用了下三角部分，因为的主对角线上的 元素均为1 ，所以无需存储。 总之，预条件技术是求解大型线形系统的必不可少的手段，但又是十分地不成熟，往 往需要结合所研究问题的具体背景加以分析。 复旦大学硕士毕业论文 第二章对称反对称预条件 2 1 分裂预条件 在工程计算中经常出现这样一种线性系统， a z = b ( 8 ) a 卯“是一个大型稀疏不对称，但正定的矩阵。若对此系统进行迭代求解，往往需要对 a 作一个有效的分裂，比如，j a c o b i 和g a u s s s e i d e l 迭代 6 】把系数矩阵分裂成对角部分 与非对角部分( 分别是严格上三角，下三角部分) ，而广义共轭梯度法【8 和广义l a n c z o s 方法 9 】是把a 分裂成他的对称与反对称部分具体可见f l o 一1 3 】很自然地，a 具有对称 与反对称分裂 8 ， a = 抒+ s ( 9 ) ，其中， 日= ；( a + a t ) ，s = ；( a - a t ) ( 1 0 ) 在本章以及以后的章节里，我f f 探讨的都将是基于这种分裂的预条 牛方法 由以上分裂，我们有 a = 日( j + 日_ 1s ) ，那么， a = ( i + 日。s ) 。1 h 。 如果我们把( + 日。s ) - 1 用他的一阶近似j 一日- 1 s 来替换，那么 p l = ( j i - i _ s ) h 。( 1 1 ) 就可以作为原系统的预条件，当然预条件的性能将取决于日。s 的谱半径的分布，并且当 反对称h 部分占优时，效果是明显的( 1 3 】。 那么当反对称部分s 占优时，g o l u b 和v a n d e r s t r a e t e n 类似地提出了一个相应的预条 件。其基本想法就是，为了保证反对称部分s 可逆，对其增加了一个位移，s + 。，这样 复旦大学硕士毕业论文 分裂相应的变为： a = h a ij rs + 旺l 那么， a = ( s + q ) ( ，+ ( s + n ) 一1 ( 日一q ，) ) 这样， a 一1 = u + ( s + a ，) 一1 ( 日一0 f ，) ) 一1 ( s + a ，) 一1 同样，用一阶近似来替换( ，+ ( s + a ，) - 1 ( 日一q 川。 岛= ( i 一( s + o ) 一1 ( 日一q d ) ( s + ，) 一1 p 2 就可以作为当反对称部分占优时的预条件从p 2 的表达式可以看出 乜有直接的关系 5 2 2 h s s 迭代 1 2 ( 1 2 ) ( 1 3 ) 预条件的效果与 z h o n g - z h ib a i 、g e n eh g o l u b 和m i c h a e lk n g 在 1 】文中，基于以上分裂，对系统 ( 8 ) 给出了一种新的迭代方法，h s s 迭代， 算法8 ( 1 ) h s s 迭代( h e r m i t i a n s k e y h e r m i t i a ns p l i t t i n gi t e r a t o n ) 。 j 、首先给出一个初始解。o ， 2 、f o rj = 0 ，1 ，2 ，u n t i l 收敛，d o ： ( a 日) z k + j = ( a i s ) z k + b ( a - 5 ) z k + l = ( 。j 日) z + + 6 其中，是给定的正的常数。 从以上算法中可以看出，h s s 迭代是以对称部分与反对称部分为主而交替进行的两步 迭代。在迭代过程中，每步都要求解分别以h + c d 和s + a i 为系数矩阵的方程组，除特别情 形外，精确求解将是不可行的。为了提高计算的效率，可以用迭代法求近似解，如使用共轭梯 复旦大学硕士毕业论文 1 3 度法( c g ) 来解以日+ n 为系数的方程组，使用k r y l o v 子空间的方法求解以s + a ，为系 数的方程组， 1 中称这种具体的求解迭代过程为i h s s ( i n e x a c th e r m i t i a n s k e y h e r m i t i a n s p l i t t i n gi t e r a t o n ) 。 f 1 1 文中给出了h s s 迭代解的收敛性分析( 引理2 1 和定理2 2 ) ，下面列出其结果。 引理1 ( 1 ) 假设a c “，a = m i 一批( i = 1 ，2 ) 是a 的两个分裂，并且给定初始向量 x ( o ) ，若z ( ) 是由以下迭代产生的序列， 女= 0 ，l ，2 ，那么 蜴+ ；= 1 2 9 , k + 6 茁女+ l = n 2 x k + ；+ b z ( + 1 ) = 蚜1 n 2 m ( 1 n i x ( + 屿1 ( ，+ 2 m f l ) 6 ，= 0 ，1 ，2 ，( 1 4 ) 另外，若迭代矩阵m f l 2 m f l n l 的谱半径p 【m f l n 2 m i - 1 n 1 ) 0 ，p ( m ( d ) ) 1 ，即可得出一m ( o ) 非奇异。另外，由假设知 a 也是非奇异，因此 k e r ( a ) = k e r ( 1 一m ( ) ) = o ， 再由 7 】， q 一1 = ( s + ，) 一1 + ( a i 一日) ( 日+ 。，) 一1 直接计算就得到 q = 去( 日+ 。州s + a ) 因为常数对加速迭代算法是无关紧要的，所以 p s s = ( h + q ，) ( s + a i )( 1 9 ) 塑旦盔堂亟圭坐些监塞 1 5 就是由h s s 迭代所推出的预条件矩阵，称为h s s 预条件。同时我们也可以单独用对称部 分或反对称部分作预条件，分别记为， p h = h + i p s = j + q 在后面的数值实验里将讨论他们各自的效果。 5 2 4 预条件的实现 我们知道在g m r e s 算法等迭代法中，具体实现预条件时，总是要求解 q z = u( 2 0 ) 这样一个线性方程组，又知道预条件矩阵是可逆的，所以方程( 2 0 ) 存在唯一解。 根据预条件的意义，对于方程( 2 0 ) 的解，并不要求精确解，只要求一定精度就可以 了。求解方程( 2 0 ) 可以采用直接法和迭代法。 关于直接法，可以采用第一章介绍的不完全分解的方式来近似求解。对预条件p h s s = ( 日+ 州s + ，) ，用不完全的c h o l e s k y 分解和不完全的l u 分解分别来分解对称部分 日+ a j 和反对称部分j 5 _ + a ， 日+ a i r t r 日+ a i l u 实现方法描述如下： 算法9 对称反对称预条件实现 1 、h = h + a i ； 2 雪= s + a i ； 3 、r = c h o l i n c ( f t ，d r o p t 0 1 ) ；不完全c h o l e s k y 分解 4 、 l ，u = l u i n c ( s ，d r o p t 0 1 ) ；不完全l u 分解 5 、y = n n 7 u ； 6 、z = u l y 复旦大学硕士毕业论文 在具体的实现不完全的c h o l e s k y 分解与不完全的l u 分解时，关系到一个丢弃宽度的 问题( d r o p t 0 1 ) ，应当适当地选取参数d r o p t o l ，它控制了不完全分解的“近似性”，太小 了，会使得每步迭代的时间过长，反之，若太大了，又会失去预条件的意义。 另外，还可以通过迭代的方法求方程( 2 0 ) 的解。如对日+ a ，部分，可以采用共轭 梯度法( c g ) 来求解，而对于s + n ，部分，可以采用一些k r y l o v 子空间的方法来解，与 直接法一样，也不要求太精确的解，可只简单地迭代一定步数，或设定一个较低的指标。 最后在g m r e s 算法的基础上，加入预条件算法，形成后面做数值实验时使用的h s s g m r e s 算法，描述如下， 算法i 0 h s s - g m r e s 算法 、赋初植，取x o = 0 ，7 0 = b ，并计算 v h = 尸i ；s r o预条件处理 ? j l = v h v h l i 2 、预条件处理( 焉；s u ) 疗= h + a i ； 雪= s + a i ； i fk i n d = = 1 由直接法实现预条件 r = c h o l i n c ( h ，d r o p t 0 1 ) ； bu - l u i n c ( s ，d r o p t 0 1 ) ； y = r r 7 u ； z = u l 可； e l s e i fk i n d = = 2由迭代法实现预条件 y = i t e r a t e l ( h ，u ) j z = i t e r a t e l ( s ，) i e n d 3 、迭代计算，f o rj = 1 ，2 ，k ，u n t i ls a t i s f i e dd o ： “2 = a v j “= 尸晶s u 2 ，j = ( u ，u 。) ，i = 1 ，2 ，j 复旦大学硕士毕业论文1 7 o j 十1 = u 一；= 1 4 ，j t j h j + t ，= l i 吗+ l i i ， v j + l = _ o j + l h j + l ，j 3 、得到解： x k = x o + k y k ，其中，y 女使( 5 ) 式达到最小。 算法中i t e r a t e l ( ) 和i t e r a t e 2 ( ) 指求解系数矩阵分别为日+ q 和s + a i 的迭代 法，具体迭代法可根据具体情况而选择。同样可以得到以f 和p s 为预条件的g m r e s 算 法，此处不再赘述 复旦大学硕士毕业论文 第三章应用实例 3 1 问题描述 本文中所考虑的模型是定常n a v i e r - s t o k e r s 方程，也就是o s s e n 方程 1 8 一u + ( 伽。g r a d ) u + g r a d p = | ，i n n ( 2 1 ) 同时压力p = 0 的情形，即对流扩散方程 一p u + f 叫g r a d ) u = f 打。n v 是r e y n o l d s 数的倒数，在实验中取，= 0 ，n 为一个方的区域 0 z 1 ，0 y 茎1 ( 2 2 ) 取边界条件u i 。：o = uv ；o = o , u r 。：1 = “1 9 ：1 = 1 铷分别取常数情形叫。= 1 ，w g = 2 和函数 情形 两种情形 叫。= 8 x ( x 一1 ) ( 1 2 y ) = 8 ( 2 x 一1 ) ( 2 y 一1 ) 5 3 2 离散模型 对于上一节所介绍的模型，用五点中心差分格式进行离散，在建立差分格式之前，我 们先剖分网格，为方便起见，取 h ：j - _ 扎+ 1 这样容易得到n 及a ，网格剖分见图一： 复旦大学硕士毕业论文 1 9 图【一) 设u ，k 表示未知函数在点0 ，k h ) 而的值，指标( j ，k ) 代表网格剖分的交点，j 代表横 向，k 代表纵向，当w 。= a ，”。= b 为常数时，得到以下差分格式： 一u j k a u j k 三( 4 u j ，k 一一1 ，k u j + i ，k 一乱，女十l u j k 一1 ) 【a u 。 ，k 【( 。) 钆b k 三去( 吩+ ，k u j + l ，k o ，一 ，k 一l ，k ) b u v 北 n ( 9 ) u j k 三击( ，+ ，k + 1 一。j ，k 一 钆似一1 ) 这样，对于模型( 2 2 ) 的离散矩阵的形式为 f = u a + n 其中n = ( 。) + ( 驯，由离散格式， ( a ) ”。= 萨1 二j ；1 ( 2 3 ) ， 星里盔堂塑主望些堡塞 2 0 这里n = 礼2 ，i 是礼n 的单位矩阵，t 是三对角矩阵， 所以a 是对称矩阵 其中 t= 胪) - 去 ( g ) ：一b 2 九 d = d o， 一， 。1 10 i 仍为n n 的单位矩阵。 当w 。= a ( x ，y ) ，”，= b ( x ，y ) 为函数时，离散矩阵的整体形式不变，仍为( 2 3 ) 式 a 也不变，但( “，( ) 分别变为， 胪) _ 去 、，lllll i一 4 1 一 一 1 4 一 、， 巩 复旦大学硕士毕业论文 其中 d 2 厶= ( ) ：l 2 h 0 i i - i l 0 o ( 1 5 ，i ) - a ( 1 5 ，i ) - a ( n 一0 5 ，i ) b ( 1 ，i + o 5 ) 。 厶一1 一厶一1 0 5 ( n o5 ，i ) 0 6 ( n ，i + 0 5 ) i = 1 ，2 ，n i = 1 ，2 ，n 一1 2 l 这里用a ( i ，j ) ，5 ( i ，j ) 分别表示。( i + ，j + ) ，o ( i + 九，j + 危) ( ，( g 都是反对称矩阵，所以n = ( 。) + ( 是系统( 2 3 ) 的反对称部分， u a 是的对称部分相应的预条件岛s s ，p 日，b 就分别为： ( u a + 口，) ( + a i ) a + q ， + a i 銎助 b 复旦大学硕士毕业论文 第四章数值实验 在这一章，我们将对称反对称预条件方法用到g m r e s 中，求解系统( 1 0 ) ，通过具体 的数值实验来测试对称反对称预条件的效果本章所有实验都在c p u 为1 g ，内存2 5 6 m 的p c 机上进行，初始向量都取为0 向量，所用程序都由机器精度为1 0 _ 1 6 的m a t l a b 编 写，迭代的停止标准均定为： 燃10-6l1r ( 0 ) l l 。 其中，r ( ) 表示第k 步迭代的残向量在没有明确说明的情况下，我们所讨论的预条件都是 通过直接法实现的，即不完全的c h o l e s k y 分解与不完全的l u 分解，并且在缺省情况下， 不完全分解的丢弃宽度d r o p t o l = 0 0 0 5 ，并针对叫。= 1 ，w = 2 为常数的情形 4 1 数值例子 首先看下面一组实验数据表( 一) 这一实验是在固定剖分网格大小( = 志，即系统矩 阵大小为4 0 9 6 4 0 9 6 ) ，及= o 0 0 1 的情况下，而采用不同的参数q ，迭代收敛所用的步 数从表中可以看出，a 对预条件的效果有明显的影响选取适当的q ，如表中取o = o ，0 5 ， 可以显著地加快收敛速度，若选取不当的话，则可能会没有效果，有时甚至会减缓了原收 敛速度， 表( 一) 未使用预条件使用h s s 预条件 取值 o 。7 6 9 2o 0 10 0 5o 0 6o ，0 70 0 8o 1o 91 5 迭代步数 4 74 33 l2 12 22 42 52 74 34 4 注1 ：h = 矗，即系统矩阵大小为4 0 9 6 4 0 9 6 ) ，扩= 0 0 0 1 注2 ：表中第一个q 是a = 刍 比较迭代收敛的过程，图( 二) ，可以得出同样的结论，不同的有不同的效果。 复旦大学硕士毕业论文 c o n v e r g e n c eh i s t o r i e sf o ri t e r a t i o nw i t hd i f f e r e n ta l p h a 图( 二) 迭代收敛比较；h = 矗p = 0 0 0 1 从上面实验中看出，适当地选取。可以得到明显的效果，那么如何选取呢? 在 1 中提出，取。= 丘= 嘉是一个较好的选择，看表( 二) 与表( 三) 表( 二) 是取。= 丘= 丽h 时，不同网格大小及取不同值时的迭代步数表( 三) 所列 出的是在不使用预条件时g m r e s 方法的迭代步数 从以上两表中可以看出，当用作预条件时，取a = 丘= 岛并非都是较好的选择 表( 二) 网格 = 1= 0 1= 0 0 1= 0 0 0 1 85 01 25 86 4 1 67 71 86 02 2 8 3 28 42 94 12 6 4 6 49 84 44 32 0 9 ，p!=o e8c芑elj日口9一芒nlc m 兰 星星盔堂塑生些望塞 2 4 表( 三) 网格 = 1v = 0 1= 0 0 1= 0 0 0 1 82 41 76 06 4 1 64 23 06 22 2 8 3 28 15 84 52 6 4 6 41 5 21 1 04 72 0 9 表( 四) 给出了一组经过实验所得出的几乎最优的o ，方便起见，定义为q 。 表( 四) 网格 工，= 1= 0 1 = 0 0 1y = 0 0 0 1 q 迭代步数o 迭代步数a t迭代步数0 ：t迭代步数 811 30 41 20 0 51 1 0 0 0 69 1 60 41 7o 21 8o 0 51 2 0 0 0 61 0 3 2o 22 30 12 80 0 41 3 0 0 0 61 2 6 40 13 40 23 30 0 31 80 0 0 51 3 表( 四) 与表( 二) ，表( 三) 相比，h s s 预条件有非常明显的效果从表的纵向看，对某 个，当剖分网格变化时，。t 变化不大依此，我们可以根据粗网格的q 。推测细网格的数 据这一结论也从以下数据，表( 五) ，表( 六) ，中得到证实 表( 五) m = 1 2 8 ，p = 0 0 1 时不同。取值的迭代步数 迭代步数l 5 4 l 3 1 l 4 3 l 7 0 l 8 4 表( 六) m = 1 2 8 ，v = 0 0 0 1 时不同。取值的迭代步数 a l 0 0 0 0 6l0 0 0 6 1 0 0 0 1 f0 0 8l 1 迭代步数l 7 0 l1 3 1 9l 7 2 l 9 5 复旦大学硕士毕业论文 表( 七) 网格 = 1= 0 1= 0 0 1= 0 0 0 1 o t迭代步数口t迭代步数n 迭代步数a t迭代步数 81921 60 0 590 0 0 45 1 694 23 3 0o 0 41 30 0 0 57 3 258 115 60 0 41 9 0 0 0 41 0 6 451 5 221 0 60 0 43 50 0 0 41 5 注j ：表格中所列出的迭代步数都指实验最优时的 表( 八) 网格 = 1= o 1= 0 0 1= 0 0 0 1 o z t迭代步数o t迭代步数a t迭代步数o t迭代步数 80 0 0 1811 60 0 0 66 00 0 0 66 4 1 60 0 11 00 21 70 26 2o 22 2 8 3 20 0 11 80o l1 8o 0 54 30 12 6 4 6 4o 0 13 2o 0 0 12 5o 0 13 6o 0 12 0 1 注i ：表格中所列出的迭代步数都指实验最优时的 从表( 七) 中迭代收敛的步数与表( 四) 中数据相比，可以发现，多数情况下，使用预 条件砌s s = ( a i + 日) ( q j + s ) 比使用预条件b = n ，4 - s 效果要好，但当比较小时，使 用预条件岛有时会取得更好的效果。这是因为v 比较小时，相对于对称部分，反对称部分 占优，因此这时b 会有较好的效果，当”稍大时，与表( 三) 相比，会发现，用b 作预条 件几乎没有效果。反之，从表( 八) 中，可以看出，当”稍大时，使用预条件嘞= 口，4 - h 有时会取得更好的效果。从以上分析中，我们可以总结出，预条件晶s s 综合了p 鼻和b 的优点，无论当”取何值时，p u s s 都有比较好的效果，但是当v 很小时，有时不如b ， 反之，当u 较大时，有时又不如p 0 ，但从整体上讲，使用尸h s s 作预条件还是一个明智的 选择。 下面我们从时间及残向量收敛的角度来比较一下预条件的效果 看表( 九) ，表中列出了对p = 0 0 0 1 ，取0 1 = 画时，从时间、迭代步数、残向量的三 复旦大学硕士毕业论文 个角度对使用与不使用h s s 预条件所做的比较从表中可以明显看出，无论从收敛时间、 迭代步数、残向量，哪一个角度来比较，在使用h s s 预条件后的收敛速度都明显比不使用 预条件时快。 表( 九) 阿格不使用h s s 预条件使用h s s 预条件 收敛时间迭代步数残向量收敛时间迭代步数残向量 80 2 86 41 6 6 3 3 1 0 一l o0 1 797 4 6 1 5x1 0 一 1 63 1 42 2 88 6 8 8 7 x1 0 70 2 21 09 7 1 7 4 x1 0 7 3 21 2 7 92 6 49 5 4 8 7 1 0 一，0 5 51 26 3 4 7 7 1 0 7 6 42 5 72 0 99 8 5 4 7x1 0 71 4 21 38 4 6 5 7 x1 0 7 注1 ：= o 0 0 1 ，口= o l $ 同样对= o 0 1 ，o = o l t = 0 0 4 ，可见表( 十) 。 表( 十) 网格不使用h s s 预条件使用h s s 预条件 收敛时间迭代步数残向量收敛时间迭代步数残向量 80 2 86 06 8 0 4 9 1 0 7o 1 71 12 1 6 0 1 1 0 7 1 60 3 36 26 7 3 3 5x1 0 7o 1 61 23 。6 4 0 7 x1 0 7 3 2o 6 64 57 6 6 2 8x1 0 70 3 81 3 8 9 5 8 2x 1 0 一o 6 41 7 64 78 1 9 0 7 x1 0 713 7 7t 85 7 7 1 5 1 0 7 注1 ：2 = 0 0 0 1 ，o = n t = o 0 0 4 图( - - ) 给出了对同一个系统，使用h s s 预条件和不使用预条件的残向量收敛比较从 这个图中我们可以看到，在迭代过程中，使用了预条件后的残向量收敛的明显比不使用时 的快 复旦大学硕士毕业论文 c o n v e r g e n c eh i s t o r i e sf o ri t e r a t i o nw i t hh s sp r e c o n d i t i o n e ra n dw i t h o u tp r e c o n d i t i o n e r 图( 三) 丢弃宽度对预条件的效果，可以从表( 十一) 中看出宽度太小，迭代步数可能会减 少，但相应的时间就会变多；宽度太大，预条件的效果将变得不明显。 最后看一下，当。= 8 x ( x 一1 ) ( 1 2 y ) ，w ，= 8 ( 2 z 一1 ) ( 2 y 一1 ) 为函数时，预条件的 效果。比较表( 十二) 和表( 十三) ，可以看出预条件处理后，收敛速度明显加快，迭代步 数减少。 5 4 2 结论 对于非对称、正定的大型稀疏线性系统，基于系数矩阵的对称反对称分裂的对称反对 称预条件( h s s ) ，结合g m r e s 算法，无论对于对称部分占优，还是反对称部分