（计算数学专业论文）一类奇异线性系统的解法.pdf

一类奇异线性系统的解法 摘要 本文主要目的在于研究任意奇异线性系统的一类迭代解法，这类解法主要是用 来求解系统的d r a z i n 逆解。在第一章中，我们简要地叙述了有关这方面研究的现 状和最新进展及本文所解决的一些问题。在第二章中，我们给出了一些基本定理， 这些定理保证了本文中的一些算法的收敛性。在第三章，我们给出了本文的主要算 法。在第四章，我们对d g m r e s 算法，d b c g 算法，d f o m 算法做了误差分析， 并对d f o m 算法和d g m r e s 算法做了比较。在第五章，我们对这类问题做了扰动 分析。在第六章，我们从一个比较新的角度去讨论了d g m r e s 算法的停滞问题。在 第七章，我们给出了数值实验的结果。并且利用这些结果对算法进行了比较研究， 尤其对算法4 ，我们通过数值实验对它进行了比较详细的研究。 关键词：奇异线性系统k r y l o v 子空间指标n a v i e r - s t o c k s 方程迭代法 g m r e s 方法d f o m 算法p o l s s o n 方程 ac l a s so fi t e r a t i v em e t h o df o rl a r g es i n g u l a rl i n e a rs y s t e m a b s t r a c t t h ep r e s e n tp h d d i s s e r t a t i o ni sc o n c e r n e dw i t hs o m ea l g o r i t h m sw h i c ha r eu s e d t o g e t t h ed r a z i ni n v e r s es o l u t i o nf o rs i n g u l a rl i n e a rs y s t e m i nc h a p t e ri ，t h es t a t e q u t o ，n e w l ye n m r g i n ga c h i e v e m e n ta n ds o m ep r o b l e m sw h i c h a r es o l v e db yt h i sp a p e ra r e c h i e f l yp r e s e n t e d i nc h a p t e ri i ，s o r t i et h e o r e m sw h i c h e n s u r et h ec o n v e r g e n c eo fs o m eo f t h e s ea l g o r i t h m sa r ep r e s e n t e d i nc h a p t e ri i i ，w ep r e s e n tt h em a i na l g o r i t h m so ft h i s p a p e r i nc h a p t e ri v ，w ed e l i v e rt h e e r r o ra n a l y s i so ft h ed g m r e sa l g o r i t h m ，t h ed f o m a l g o r i t h ma n d t h ed b i c g a l g o r i t h m a n dt h e n ，w ec o m p a r e t h ed g m r e s a l g o r i t h mw i t h t h ed f o m a l g o r i t h m h ic h a l ) e rv ，w em a k es o m ep e r t u r b a t i o na n l y s i sf o rt h e s ep r o b l e m i nc h a p e rv i ，w ed i s c u s st h es t a g n a t i o no ft h ed g m r e s a l g o r i t h m a san e wm e t h o d i n c h a p t e rv i i w ep r e s e n tt h ec o m p u t a t i o n a lr e s u l t s a n dt h e nw ee x a m i n et h em e r i ta n d t h es h o r t a g eo ft h e s ea l g o r i t h m sa tl a s t ，w ee x a m i n et h ea l g o r i t t m l4c a r e f u l l y k e y w o r d s ：s i n g u l a rs y s t e mk r y l o vs u b s p a c e i n d e xn a v i e r s t o c k se q u a t i o n i t e r a t i v em e t h o dg m r e sm e t h o dd f o m a l g o r i t h m p o i s s o ne q u a t i o n 第一章绪论 求解奇异线性系统的方法很多，诸如矩阵分裂方法【1 , 2 ，5 2 】，半迭代方法【3 ，4 ，5 ，6 ，5 4 这方面的综叙可以看文 4 ，最近人们更多的开始注意k r y l o v 子空间方法，这种方法 的主要特点是只涉及矩阵向量乘积。特别地，在一些研究中，涉及到的矩阵向量乘 积是a 和a + 与向量的乘积。g m r e s 7 ，8 ，q m r 和t f q m r 9 能用来求解在m a r k o v 链模型问题中出现的相容的奇异线陛系统。当系数矩阵是厄米特正半定的，并且系 统是相容的时候，在文 1 0 中，使用了共轭梯度法。在文 1 1 中指出当系数矩阵不 是厄米特但相容，且指标是1 的时候，a r n o l d i 方法 1 2 1 ，e i s e n s t a t 在文 1 3 中的 广义共厄残量法( g c r ) 和l a n e z o s 方法 1 4 】都是有效的方法。并且在文f 1 1 , 2 5 】中， s i d i ，魏益民和吴和兵还给出了所有这些方法的收敛理论。对于不相容的系统，求 解是一件很困难的事情，在c a l v e t t i 等人的文 1 5 中，对于指标为1 不相容的情 况，使用了共厄梯度类型方法。在 1 6 中，f i s c h e r 对于同样的情况，在文 1 5 的 基础上做了些推广。必须指出，当指标大于l 的时候，所有的这些方法在求解时都 变得非常困难。因此这些方法有着一定的局限性，从而寻找新的方法就很有必要了。 因为d r a z i n 逆在指标为1 的时候就是群逆，也就是方程的解逆，这样求解系统 的d r a z i n 逆的方法就比以上的方法更具有广泛性了，并且d r a z i n 逆有着广泛的应 用，例如：在奇异微分，差分方程，m a r k o v 链，迭代法，多体动力学和数值分析等 方面，具体可以见文献【3 , 4 ，1 7 - 2 8 ，因此我们可以从求解系统的d r a z i n 逆解入手来 寻求新的方法，而且求解奇异线性系统的d r a z i n 逆解也有着重要的意义。在文 2 9 中，s i d i 给出了用k r y l o v 子空间方法求解任意指标的奇异线性系统d r a z i n 逆解的 统一框架，在文 3 0 ，3 1 1 中，s i d i 又给出了在这个框架下的d g m i ：t e s 算法和d b i c g 算法，本文给出了在同样框架下的d f o m 算法和收敛性分析，并证明在效果上与 d g m r e s 方法等效， 所有的这些用k r y l o v 子空间来计算d r a z i n 逆解的方法的最大优势是不受指 标限制，但这些方法要求系数矩阵有比较好的性质，这点在数值实验中计算摄动 后的n a v i e r s t o c k s 方程时表现得很明显，尤其是这些方法不能加预条件更是他们 的缺陷。所以，本文针对这点还给出了两次迭代的方法，这种方法的最大优势是可 以加预条件，因此可以计算n a v i e r - s t o c k s 方程这样的病态问题，数值实验也很好的 说明了这点。因为算法4 实际上就是两次计算，所以在选择预条件矩阵的时候， 我们可以参考经典的鞍点问题的预条件 3 2 4 0 ，在本文中，我们选择了上三角预条 件，因为上三角预条件在实际计算中效果是比较好的，对于这点我们可以参看【4 5 。 最后本文给出了有关这些算法的数值实验，分析了这些结果，并且从多方面比 较了这些算法的优缺点和算法在不同条件下的表现，希望能为以后的研究工作提供 一些思路。 2 第二章一些基本定理 2 1引言 在这一节，我们将给出用k r y i o v 子空间方法进行半迭代的一些基本定理。 如果令z 。是迭代的初始向量，用我们这里所讨论的迭代方法计算出来的迭代 序列为：z t ，x z ，它们具有如下的一般形式： # ，n = 。o + q m - l ( a ) r o ；r o = b a x o ，( 21 ) 这里口。一1 ( a ) 是a 的次数最多是”。一1 的多项式。下面我们定义 p 。( a ) = 1 一a 。一1 ( ) 这称为是第m 次迭代残向量多项式，因为： 注意到 r r 。= b a x ，r = p ”l ( j 4 ) 7 、0 p 。( 0 ) = 1 正如文【4 ，所有象上面那样产生的迭代序列收敛的充要条件是 1 i mp f ：) ( o ) = 0 ，i = l ，2 ，o ；= i n d ( a ) m - - 4 且 ( 2 2 ) ( 23 ) ( 2 4 ) ( 2 5 ) l i m p ( , ) ( a j ) = 0 ，i = 1 ，2 ，一，b 一1 ( 26 ) 这里是a 的非零特征值，k j = i n d ( a a i i ) ，i n d ( a ) 表示矩阵a 的指标。 矩阵a 的指标是指满足，a n k ( a 。) = m n k ( a 。+ 1 ) 的最小非负数n ，它也等于矩 阵a 的零特征值的最大j o r d a n 块的维数。如果i n d ( a ) = 1 ，a d 也是a 的群逆， 通常用a l 来表示。 3 如果 p 聊( o ) = 0 ，i = 1 ，2 ，对于所有的m = 0 ，1 ，( 27 ) 成立，那么条件( 2 5 ) 是自然成立的，所以我们这里讨论的迭代序列的残向量多项 式满足( 2 7 ) 和( 2 4 ) 为了方便起见，我们在下面介绍一些记号 ，。表示所有i n 次多项式的集合。我们定义： n o = ( p e i i ，。：p ( o ) = l ，p ( 2 ( o ) = o ，i = 1 ，2 ，o ) ( 2 , 8 ) 因此我们所讨论的残向量多项式全属于n o 。注意到当m = 0 ，l ，a 时，p 。( a ) = 1 是曼中唯一的成员，当m a 时，袅中p m ( a ) 有这样的形式： m a p m ( a ) = 1 一c 。a 叶4 ， t = 1 让我们用g 来表示对应于a 的非零特征值的不变子空间的直和，用s 来表示对 应于零特征值的不变子空间的直和。因此，占就是r ( a “) ，即 n 的值域，而s 是 n ( a o ) ，即1 4 0 的零空间。c “中的每个向量都能被写成一个在g 中的向量和一个 在s 中向量的直和。 显然，a a d 是到r ( a “) 上的斜投影，所以我们能分解b ：b = 6 + 5 ，这里 6 s = r ( a “) ，5 g = n ( a 。) ，且a x = b 的d r a z i n 逆解a d b 是在s 中唯一满足 相容系统a x = b 的向量。 让我们也把x o 同样分解，那么z o = 岳o + i o ，这里童oes ，童o s ，并且正如 t h e o r e m4 1 o fc l i m e n te t a 3 j ： z 。一a 。b = 只。( a ) ( o a 口b ) 十i o 也就是说，t 。= i ，。十i 。，i 。= i o = ( i a a d ) 。o 占对于所有的m 成立。显然， 我们能取x o = 0 或者z o = a 8 f 对于任意的f c “。在这两种情况下i o = 0 。 我们定义。= 。一a d b ，”一0 ，l ，简单的讲，我们想让趋向于零。而 在文 2 9 】中，s i d i 指出只要a “ 0 就可以了。 4 在本文中，我们使用投影法来构造算法，使俨r 。+ 0 。下面，我们给出投影 法的一般形式： 因为p 。( a ) = 1 一：兰。c ；a 。“且。l ( a ) = ；m - - 1 。c f ”“，我们有： z 。= ：c 0 + q + 如果a 是n 阶方阵，我们定义n ( m a ) 矩阵v 和m a 维向量c 如下 那么 v = a 。”o ，a “+ 1 r o ，- ，a “一1 r 0 1 ，c = c 1 ，一，一。 t ( 21 0 ) z m = o o + v c ，r m = ”0 一a v c ( 2 1 1 ) 我们再用n ( r n a ) 矩阵w 表示子空间w 。然后，再让向量a 。r 。与子空间w 正交，即+ a o r 。= 0 ，因此有 w + a “+ l y c = w + a “r o 所以，我们只要选取不同的w ，就可以得到不同的算法。 2 2收敛定理 ( 2 1 2 ) 在本节，我们将要给出一些收敛定理。这些定理都是来自于文 2 9 ，因此这里 只引用，而不给出证明了： 定义1 1 2 9 a 奇异，i n d ( a ) = n 。我们将p ( a ) 称为a 对于向量o s = r ( a n ) 的a 一不完全最小多项式，如果p 并且m 是可能满足p ( a ) n = 0 的最小值。 引理1 _ 2 2 9 p ( a ) 存在并且唯一。 里q 是a 对于i 的最小多项式的次数， 在r i 。o + 。中满足p ( a ) f i = 0 唯一的成员。 并且，它的次数m 满足q m q + a ，这 因而q d i m ssn 一“。实际上，p ( a ) 是 这个定理给出了a 不完全最小多项式的存在唯一性，在这个基础上，s i d i 得到 了有关投影方法的基本定理。 5 9 + a “ 一 一 = m 定理1 3 2 9 】令p ( ) 是a 对于向量南= 2 0 a d b 的。不完全最小多项式，m 是 次数。那么p ( a ) = 1 一“i = 一i 。a 。舻押，蠡唯一，并且a d b + 洳= x o + m ；- - i 。e t 钟+ “r o 。 定理1 4 2 9 】令m 是a 对于向量o 的次数。而且z ，。是由投影法产生的向量。 z m = a d b + j o 。 。一a d b 的a 不完全最小多项式p ( a ) 那么，只要d e t ( w + a a + l v ) 0 ，我们有 从这个定理，我们可以发现，如果用投影法来求解奇异线性系统a x = b ，只要 满足d e t ( w + a 。+ 1 矿) 0 ，我们总可以使由投影法构造的迭代法不中断地收敛到系 统的d r a z i n 逆解再加上这个迭代初值在n ( a 8 ) 中的斜投影，因此我们可以适当的 选择迭代初始值，从而使迭代解收敛到d r a z i n 逆解。 6 成了 第三章算法 3 1d g m r e s 算法 在本节里，我们选取= a ”1 v ，这里的v 是( 2l o ) 中的v ，那么( 2 1 2 ) 变 ( a “+ 1 y ) + ( a “十1 y ) c = ( a 。+ 1 v ) + a “r 0 ， ( 3 1 ) 通过这样方法计算出来的z 。满足 | | a 。r m l l = 。+ r a ，i 。n ( ； 。) 1 1 ( a “( 6 一a 。) 1 1 ， ( 3 2 ) 这里的k r y l o v 子空间( a ；6 ) = s p a n b ，a b ，a m - 1 b ) 。 在 2 9 】中的定理3l 保证了d g m r e s 算法产生的迭代序列z m 的存在唯一陛 同时给出了d g m r e s 算法： 算法1 1 ，选取z u 并且计算r o = b a ：c o 和4 “r o 2 ，令卢= l i a “r o | | 且u l = 卢_ 1 ( a 8 r o ) 3 ，用g r a m - s c h m i d t 过程来正交化k r y l o v 向量a o r o ，a a + l r o 对于i = 1 ，2 ，一循环做如下运算 计算h j i = ( v j ，a v i ) ，j = l ，2 ，z 计算魂= a v i 一：1v j h j 令h 川。= j 1 啦| j 并且v i = 。 + 向量虬，就是这样得到的正交组。 4 ，对于 ：= 1 2 ，形成矩阵亿c “。和凰e ( 1 ) 柚 7 亿= v l l w z 卜h ，凰= ，0 1 1h 1 2 h 2 lh 2 2 0 h 3 2 0 h l k h 2 k h k k 0 h k + l ， 5 ，形成矩阵豆。= 鼠。豆。一1 豆。一。鼠。e ( m + 1 ) “( m - a ) 有如下形式 = 7 u 2 h 2 z n + 2 1 ： 0 h a + 3 ，2 0 1 m n 7 0 2 m n 九m n 一1 仃p 。n n m 一8 m n ( 3 - 3 ) ( 3 4 ) 6 ，计算膏，。的q r 分解：直。= q 。r 。；q 。g ( “+ 1 ) 。( 一n ) 并且r 。g ( m n ) ( m - ) 的上三角矩阵。 7 ，解上三角方程组r 。f 。= 1 3 ( q 未e 1 ) ，这里e l = f 1 ，0 ，o 】r 。 8 ，计算x 。= z o + e f m 。( 那么，显然i i a 。r m l l = 卢r 邗百磊了酽) 从上面我们可以发现，计算向量z 。需要存储m a + 1 个向量v l ，”。一。， 和z 。因此随着迭代的增加，所需要的内存也越来越多，因此我们在实际中大多 采用了重开始的方法，这也就是重开始的d g m r e s 方法： 1 ，选择适当大的整数m 0 。 选择迭代初始值u ( o ) ，并且令z o = ( o ) 且i = 1 。 8 由迭代初始值z o 开始用d g m r e s 算法计算迭代序列。- ，2 ：2 ，并且 m 3 ，如果u ( ；) 达到所要的精度就停止算法，否则，令z o = “( ，用i + l 取代i 重复第2 步。 由“( 一1 ) 计算出“( 2 ) 的过程构成了算法的第i 个循环。象这样用d g m r e s 算法 产生迭代序列( u ( ) ) 罂。就被称之为m 步重开始的d g m r e s 算法，与i n 步重开始 的g m r e s 算法类似，被记为d g m r e s ( i n ) 。这样我们就可以克服内存无限制地增 长的困难了。 3 2d b i c g 算法 在这个算法中，我们使用l a n c z o s 类型的方法来求解d r a z i n 逆解。这里的迭代序 列z m z o + k m ，这里= s p a n a 。t o ，a “+ 1 r o ，a m 一1 ”o ) ，并且r m = b a z m 上 l m ， l 。= s p a n ( a + ) 。+ 1 f o ，( a + ) “十2 而，- ，( a + ) ”f o ) ， 这里的= b a + o ，而j o 是任意的。在这里我们选取的w 是 s p a n a + f o ，( a + ) 2 f o ，( ) “一。怕) ， 并使小r 。正交于它。文 3 0 中的定理3 1 保证了这一正交性质，从而也保证它的 收敛性。 这个算法是短递推的，因而它要存储的向量是固定的，独立于系数矩阵的指标 a ，所以它在计算速度上可能要比d g m r e s 算法快，但它的残量性质没有d g m r e s 好。下面我们给出d b i c g 算法的具体形式。 算法2 1 ，任意选取。o 和圣o ，计算r o ；b a z o ，9 0 = b - a 圣o ，并让。q = 9 2 0 ，= ”o ，v d l = a o r o ，仉一1 = ( a 4 ) 。怕，u 口一1 = 1 ； n = a ： 9 2 ，如果| | 4 。一l i i t o l 如果n 三“+ 1 ，那么5 。：= 如果n2a + 2 ，那么伽：= ，重复下面过程： 一垮暑竽旱，否则如：= 0 l d 一1 ， n lj 7 一 ” 一镳童划，否则：= 0 ( “一2 ， 一2 ) d n ：= u n 一1 ( v n 一1 + 西l d n l + 7 n d n 一2 ) u n ：= u w l ( a n 一1 + “ n 一1 + r n v n 一2 ) 西n ：= o n l ( a + 面nl + 眈西n l + 何t 西n 一2 ) 一! ! 竺! 1 2 l 。“：2 丽； r n + l ：= n l u n n ； z n + l ：= z n + u n d n ； n ：= n + l ： 很显然，我们可以不用o ，但是我们要选取o 。 从d g m r e s 算法中发现，这种算法在( o 。，u 。) = 0 或者( o 。，r n ) = 0 时会导致u 。没 定义或者= 0 ，从而会导致算法失效。但是在( o 。，v 。) z0 ( 并且) 或者( ，_ 。) z0 时会引起算法停滞。对于这一点，我们可以用l o o k a h e a d 技巧来克服。 3 3d f o m 算法 在这个算法里，我们选取的w 是v ，这里的v 是( 2 1 0 ) 中的v ，那么( 2 1 2 ) 变成了 ( y ) ( a 8 + 1 v ) c = ( y ) + a 。”o ，( 3 5 ) 通过这样方法计算出来的与d g m r e s 算法有类似性质，最多迭代n a 步就可 以得到精确解，但它有可能会发生中断。 下面我们来推导算法： 因为。，。= z o + i = 1 c 。a 叶”1 r 。o ，首先用a r n o l d i g r a m s c h m i d t 这样的过程来 正交化k r y l o v 向量a 8 t 0 ，a ”1 t o ，。我们具体执行的时候用和d g m r e s 算法第3 步一样的修改g r a m s c h m i d t 方法来执行这一过程。 同样我们把正交化后的向量记为v l ，v 2 ，当i q 一1 的时候，这里q 是矩 1 0 阵a 关于小珊的最小多项式的次数，因而也是关于v l 的最小多项式的次数，满足 而且，对于每个k i + 1 a = v j h j i ，i = l ，2 8 p a r o l v l ，”2 ) +j 啦 2 = 如果我们定义n k 矩阵亿如下： s p a n a 。r o ，a 。+ 1 r o ，a “+ “一1 o ) k k ( a ：a o r o ) 垓= v l ，v 2 ，冁 ，k = 1 ，2 ( 3 6 ) 那么对于m m o ，有。= z o + 亿，对于某个ec m 。因此我们必须求 出岛。首先，由r 。= r o a 一。f 。，我们有 a 。r m = a 8 r o 一4 “+ 1 一n f m = 卢 1 一a 。+ 1 一n 矗n 其次，只要k 茎q 一1 ，从( 3 6 ) ，我们有 a k = + l 王k ；h k = uh i 2 h 2 1 2 2 0 h 3 2 h l k h 2 k ： + h k k 0一0 h k + l ( 37 ) 注意到鼠eg ( + 1 ) “。因为当ksg 一1 时，向量a 1 ，a 2 ，4 是线性独立的， 所以r a n k ( a 咴) = k 。又因为r a n k ( 仇+ t ) = k + l 且 r a n k ( a v e ) 曼m i n r a n k ( 记k + 1 ) ，r a n k ( 凰) ) ， 我们有r a n k ( i k ) = k 。也就是说凰是满秩的。下面把( 3 7 ) 用于a 8 + 1 一。，只要 m g i ，我们有 a 卅1 吒一。= a “玩。+ 1 忘， ；a ”1 吨讲2 凰叫l 曰m 一。一= + l 颤； - i m ；瓯_ 1 砜一 1 1 因此，这里的宣。就是d g m r e s 算法中的岛。 只要7 t i , q l ，我们有 a 。r 。= 卢w 1 一+ l 矗。 。，( 3 8 ) 并且因为”- = + l e l 和昧一。+ l = ) 。( m - a ) ，o ( ma ) 。( 1 ) 由上式，最终我们 有 0 = ( 优。一。) 4 a 。r 。= ( e 。一。) + ( 3 e 1 一直k f 。) = 卢e 。一直南一。f 。( 39 ) 这里 庙k 一。= i c m - a ) 。( m - a ) ，o ( m - a ) 。( 。+ 1 ) 西- m h l lh 1 2 h 2 lh 2 2 n + 2 ，1 ： 0 h n + 3 2 00 所以我们为了求只需要求解下面的线性系统 f k 。= 卢e 1 l m 2 、m n f 3 1 0 1 f 3 1 1 1 显然，我们发现当或。一奇异时，算法会中断。当直。一。非奇异时，我们能 将鼠。一。进行q r 分解，即或。一。= q 这里q 。一。g ( m - a ) 。( m “) 是酉矩 阵，r 。一。g ( m - a ) ( m - a ) 是非奇异的上三角矩阵。所以系统( 3 1 1 ) 就变成了 r 岛= 卢( q 麓一。8 1 ) f 3 1 2 1 在具体计算时，我们发现d f o m 算法比d g m r e s 算法要少做n + 1 次正交化。 在形成矩阵或。一。之前，要首先形成矩阵宣。，关于矩阵盘。的形成我们可以参见 文 3 0 】。下面我们把算法具体写出来： 算法3 i ，选取z o 并且计算”o = b a x o 和a “r o 1 2 2 ，令卢= i i a “r o | | 且 l = 卢“( a o r o ) 3 用g r a m s c h m i d t 过程来正交化k r y l o v 向量a 。r 。o ，a 外1 7 0 对于i = l ，2 ，循环做如下运算 计算h j i = ( v j 、a v i ) ，j = l ，2 ，i 计算砚= a v i 一：lv j h j i 令h i 扎。= 慨| | 并且v g = 讥 h 向量钆v 2 ，就是这样得到的正交组。 4 ，对于k = 1 ，2 ，形成矩阵亿c 。2 和凰g ( 。+ 1 ) 玩= p t 池卜h ，盈 h i th t 2 2 lh 2 2 0 h 3 2 0 h l k h 2 k h k k 0 h k + l ，k 5 ，形成矩阵岛。= 露。蟊。一1 鼠。立。g ( ”+ 1 ) ( m 一“，然后由这个矩阵 形成豆。鱼。一。c ( m 一“】( 一。) 有如下形式： 一。= h l lh i 2 h 2 th 2 2 。2 l 0 7 乜+ 32 00 l m n 2 m 一 6 ，计算茸。一。的q r 分解：藏。一。= q m - ，, r m - , 。；q 。一。e ( m - a ) 一n ) 并且 兄c ( m - a ) 。( m o ) 的上三角矩阵。 7 ，解上三角方程组j h 一。= 卢( 0 未一。e 1 ) ，这里e l = 【1 ，0 ，o t 1 3 8 计算z ，。= x o + 以 从上面我们可以发现，与d g m r e s 算法一样，计算向量z 。需要存储m a + 1 个向量叽，一，”。和z o ，因此随着迭代的增长，和d g m r e s 算法一样，所需 要的内存也增长。所以我们也可以用重开始方法来解决这个问题，这就是是重开始 的d f o m 算法，它的具体实行过程与重开始d g m r e s 算法一样，在这里就不再重 复了。 3 4对于指标1 的两次算法 前面的三个算法可以计算系数矩阵有任意指标的奇异线性系统，但是它们不能 加预条件，因而在解病态问题时，效果将会很差。并且因为很大一类病态问题的系 数矩阵的指标是1 ，所以我们这里使用了两次算法来克服这个困难。 假设线性系统a x = b 的系数矩阵指标为1 ，由d r a z i n 逆的性质( 参看 1 7 ，4 1 ) ，求 d r a z i n 逆解问题的本质就是求解下面的线性系统： a a x = a b ( 3 1 3 ) 因为系统( 3 1 3 ) 总是在r ( a ) 里，而且是相容的。所以我们可以加预条件。而在具 体计算时，我们可以分两步来计算：( z ) ，a y = a b ；( 吲，a z = y 。在具体计算时，我 们可以选择各种k r y l o v 子空间算法，如：g m r e s 算法，b i c g 算法，c g s 算法 等，并对各种不同的问题加具体的预条件。 1 4 第四章误差分析 在本节，我们将考虑对d g m r e s 算法，d b i c g 算法，d f o m 算法进行误差 分析。 我们首先假设d f o m ，d b i c g 可以不中断的迭代下去。在文 3 1 中，s i d i 等人对 d 1 3 i c g 算法做了误差分析，实际上这种误差分析也可用于其他的两种算法， 由本章第一节，我们知道这三个算法本质上都是投影法，只不过是选取了不同正交 的矩阵w 。在开始讨论前，我们先介绍一些记号，和前面一样，我们用r o 表示初 始残向量，但这里矩阵一。定义为 k n 一。= a 。r o ，a 。+ 1 ”o ，- 一，a m 一1 r 0 一。的列张成的线性子空间用圮。一。表示。而在投影法中，用来正交的子空间w 的列张成的线性子空间用工表示。 当我们用算法l ，算法2 ，算法3 进行迭代的时候，我们有z 。x 0 + k m ，且 4 “r 。上，所以可以得到w + 小r 。= 0 ，即有 w a “( 7 0 a u 。) = 0 ，“。j 矗，( 4 1 ) 这里“。是奇异线性系统 a u = r o( 4 2 ) 的近似d i a z i i - 逆解。因为我们可以将u 。= yy 满足 w + a 。( r 0 一a c m 一。) 分= 0 由假设，我们知道a n + - 一。非奇异，所以 y = ( w + a 。+ 1 吒一。) _ 1 w + r o( 4 3 ) 注意到加= + 拍，这里f o r ( a “) ，而n ( a 8 ) ，并且这样的分解是唯一的。又因 为a 。r o = a n f o ，所以( 4 3 ) 就变成了 口= ( w + a 。+ 1 一n ) 一1 w + 。f o 令z 】= 孰+ 钆，这里o r ( a “) ，o n ( a “) ，且令n 是系统( 4 2 ) 的d r a z i n 逆解， 令j = a d b ，因此n = 一2 0 。 1 5 令只。是到子空间凰。上的正交投影，我们这里将按照= f l ( s 一玮。) 训的尺 度来分析算法的误差。这里_ | | | 表示欧氏范数。 下面的引理来自于文 3 1 中的引理4 3 。 引理4 1n 到k r y l o v 子空间一。的距离i l ( 一f k ) 训满足 i i ( s p r o ) a l l2 。m 。i n 。l l p ( a ) ( i o i ) n( 4 - 4 ) p 1 1 i 我们这里要用算子方程来解释这种误差分析方法，我们首先要定义到o + 墨。并且 正交于子空间l 的算子q 。： 且 引理4 2 假设 q m o f o + 墨。 a 。( 。一q ，。z ) 上j 0 。 d e t ( w + a o v ，m 一。) 0 那么算子q 。在a m 的标准基下的矩阵表示是 q 。z = f o + 窿。一。( w a 一。a “蟓一。) 一1 w + a 。( z f o ) 证明：由q 。的定义可以知道 且 q 。z = f o + k n n y a 8 ( z q 。o ) = o ( 45 ) 因此 w + a “卫= h aq 。z = w + a “( f o + v a - a y ) 所以我们有 y = ( w + a 。一。) 。1 + a 。( z f o )( 4 6 ) 然后将( 4 6 ) 代入q 。z 的表达式就可以得到结论了。 口 1 r 引理4 ，3 假设( 45 ) 成立，则q 。是到o + 玛n 上的投影。 证明：对于任意的z ，我们有 q 纛z = 拍+ 记。一。( + a “e 。一。) 一1 w + a 。( q m z f o ) = f o + e 。一。( w + a “e 。一。) 一1 w + a 。( 拍+ 识。一。( 仉7 4 4 。谚n 一。) 一1 w + a “( z 一庐o ) 一声o ) ：托+ 亿( + a “谚。) 。w + a “( 。一o ) = o 邢z 口 注意到( 4 5 ) ，所以q 。z 是唯一的。 下面定义a 。： a 。= q 。a p m 那么我们有 引理4 4 当w + a n + 1 姝一。与w + a “吒一。非奇异时，线性系统 f o a 。u = 0 ，u 一。 ( 47 ) 有唯一解“。 证明：因为u 识，。，所以u = 一。y ，而且p m u = “。因此由怕一a m “= 0 有 f o a 。u = f o q m a 只n “ = 矿。一q m a u ：f o f o e 。一。( + a o e 。一。) 一1 w + a 。( a u 一怕) 所以 一。( + a “吒一。) 。w + a 。( a u f o ) = 0 又因为瓯一。是列满秩的，所以 w + a 8 f o = w + a 。+ 1 一n 9 ， 从而有 y = ( + j 4 叶1 一。) 。1 w + a “f o ， 即意味着( 4 7 ) 有唯一解。 口 17 这里( 4 7 ) 是原问题的近似，下面给出( 4 7 ) 的误差分析： 引理4 5 令 0 。= l i q ，。a ( i p ”圳， 那么a u = f o 的d r a z i n 逆解n 满足： | | f o 一4 n | | o m z 。 证明：因为而f o + 蜀。，和q 。是到f o + k 。上的投影，在加上= q 。f o ，a 也= ， 我们有 i o a m b , = q m e o q m a p m d = q 。( f o a j ) m n ) = 0 。a ( i 一晶。) n = q 。a ( i p m ) ( ，一j ) m ) 血 所以结论成立。 口 引理4 6 设( 4 5 ) 成立，a 。k 。表示a 。在子空间啊。上的限制，那么a 。i k 。可 逆。 证明：我们只需要证明 ( a ml k ) “= 而+ ， w k 二( 48 ) 有唯一解u k m 一。，叉因为“，”墨。，所以我们有 “= 一n z ， = 一y ，y ，z c ”1 由a 。与q 。的定义，那么( 4 8 ) 成为 f o + 一。( + a 。一。) 一1 w + a “( a j “一内) = i o + 谚n - a y 因为一。的列是线性无关的，所以上面的线性系统就变成了 ( + a “e ，。) 一1 w 4 a “( 4 只。u f o ) = y 两边同时乘上w + 俨一。并注意到p m u = “，我们得到 + a “十1 一z = w + 产o + w p a 8 v m 一。y 1 8 因为缈+ a ”t 一。是非奇异的，所以唯一解存在。 口 定理4 7 令引理3 5 那样定义，并且令一。= i i ( a 。i k 。) 。| | 。那么误差| | z m 一 ( + 2 0 ) | | 满足 i z ，。一( j + i o ) i is ( 1 + 皖k 毛) 1 2 。 证明：由引理4 4 和引理4 6 ，我们有u 。= ( 4 。j n ，。) 一1 f o ，因为t z m 墨n ，我们有 p k t 。，。= “。，从而有 p m ( “。一血) = “m 一f ) m i = ( a ，。1 。) 一1 嗡一( a 。f 。) 只。训 其次，因为焉矗墨。并且焉= 只。，我们有( a 。k 。) 只。= a 。p m n = a m n ，因 此 f ( “m o ) = ( a mi k m ) 一1 ( 庐。一a ，。血) 由引理4 5 ，我们有 1 1 f ) m ( u 。一血) | | 一m p m e m 因为 也一“。= n 一u 。= ( 。一只。) 也十j 陋一u ，。) 并且因为f 和k 一马。是互相正交的投影算子，所以 i | “。一二i j 2 = j j ( ，一p ，。) i i l 2 + l i p , 。( 也一u 。) 1 1 2 f 2 + k ，2 。2t 。2 又由w 。= 。o + t m ，因此 j i z m 一0 + 量o ) = | l z o + ”m s 一童d i = i l u ，。一( 一o ) l l = n 。一训 s ( 1 十一未嚷) 1 2 e 。 ( 49 ) 口 还有些有关结论，可以参看文 2 9 第6 节和文 3 0 的定理3 2 。 d f o m 计算出来的第z n 步迭代解我们用z 景”表示，它的残量用螺9 表示；d g m r e s 1 9 计算出来的第。步迭代解用。铲表示，它的残量用r 。d g 表示。同时用m 未表示 宜三直。，那么我们有 峨：穰砒= m 矾一 ：直：一。曾。一。+ b t b ：由：一。( 一。+ 碟一。只。一n ) 爿k n 由假设，d f o m 不中断，所以或。一。是非奇异的，这里有f ，。= b 直三! 。下面的 定理是有关d f o m 和d g m r e s 的比较。 定理4 8 不失一般性，我们令迭代的初始值z o = 0 ，直m n 非奇异，那么： z 。d = 1 i a 1 j v r m n h m - n e l ； i i a “r 驴| | = 刚f m 一。e l l l ； 。d g ：i i a n 圳吒一。烁2 藤e 1 _ b 1 1 一。蛎2 h m - a e l ； r 船i i ：卢1 ( ( 一。e 1 ) 7 1 ( k + 一。碟一。) 一1 ( 一。e 1 ) ) 1 胆； i i a 8 r 。d g i isi i a 。螺f m 证明；首先计算a 。r 驴： i i a n r 。d | | = t l a 。b i i a 。b l l a 。+ 1 识。一。曾i ! 。e 1 | ：k l a n b i i a 。b 1 1 + ，鼠厩! 。e 1 i i 钏椭w 酬( 鼍。) 昧一1 | ：”a a “一l l a ”+ ，( i , n - 。) + ( z 一。) ) e - “ = u 肌圳肌m ( 。) e - i i a i i ( 0 。) e - = n 删m ( 0 。) 蚓i 钏删u ( 只训i 嘲| f ( e 曼。) 圳i = 卢1 1 1 ，一n e l m 下面计算a 。r 寮g ： i i a n r 宝g | i = l i d 。b i i a 。b l l a 。+ 1 1 一。直高e lf = i i a 。b i i a 。圳+ 1 直。2 1 。+ e 1 这里2 1 + 表示成。的m o o r e p e - - r o s e 逆，所以直。矗志是投影阵，令z 为丘。的规范 正交的左零空间组成的矩阵：矿鱼。= 0 ，所以岛。2 1 + = 一z z t 。我们有： i 【a 。r 是g l i = i i a “b ir a “b 1 1 + l ( 。一z z r ) e i j j = i i a “b l l l l z z 丁e l |