（计算数学专业论文）奇异线性方程组的gmres型方法.pdf

摘要 在这篇文章中，首先讲了g m r e s 方法和r r g m r e s 方法的算法，以及这些 方法能够得到奇异( 可能不相容) 线性方程组的最小二乘解所需的条件其次，用 r r g m r e s 方法解s t o k e s 方程在解s t o k e s 方程时，采用块l u 预条件，这种预条 件算法和其他类似算法的主要不同点是，较精确的s c h u r 补矩阵有效的计算出来， 随后使用它再计算s c h u r 补矩阵的预条件结果表明这种预条件是较好的 关键词：k r y l o v 子空间，r r g m r e s 方法，奇异线性系统，预条件方法，不完 全l u 分解，块l u 分解 a b s t r a c t i n t h i s p a p er ，f i r s t l y , t h e b e h a v i o ro fi t e r a t i v em e t h o d so fg m r e sa n d r r g m r e sw h e na p p l i e dt os i n g u l a r ，p o s s i b l yi n c o n s i s t e n ，l i n e a rs y s t e m si sd i s c u s s e d a n dc o n d i t i o n su n d e rw h i c ht h e s em e t h o d sc o n v e r g et ot h el e a s t - s q u a r e ss o l u t i o no f m i n i n a ln o r ma r ep r e s e n t e ds e c o n d l y ，w ed i s c u s s e dh o w t os o l v es t o k e se q u a t i o n w i t hr r g m r e sm e t h o da n dab l o c kl up r e c o n d i t i o n e ri sp r e s e n t e d t h em a i n d i f f e r e r i c eb e t w e e nt h ea p p r o a c hp r e s e n t e dh e r ea n dt h a to fo t h e rs t u d i e si st h a t a n e x p l i c i t ，a c c u r a t ea p p r o x i m a t i o n o ft h es c h u rc o m p l e m e n tm a t r i xi se f f i c i e n t l y c o m p u t e d t h i si su s e dt oc o m p u t eap r e c o n d i t i o n e rt ot h es c h u rc o m p l e m e n tm a t r i x t h a ti si nt u r nd e f i n e sap r e c o n d i t i o n e rf o rag l o b a li t e r a t i o n t h er e s u l t si n d i c a t et h a t t h i sp r e c o n d i t i o n e ri se f f e c t i v e k e yw o r d s ：k r y l o vs u b s p a c e ，r r g m r e s m e t h o d ，s i n g u l a r l i n e a r s y s t e m s p r e c o n d i t i o n i n gm e t h o d s ，i n c o m p l e t el u ，b l o c kl u 2 引言 当线性方程组a x = b 的系数矩阵a 是稀疏的且阶数1 3 很大时，用高斯消去法 显然不理想，因为它要至少进行o ( n 3 ) 次的操作步数，时间花费太大。那么怎样 解类似的线性系统昵? 1 9 8 6 年，s a a d 和s c h u l t z 提出了广义极小残量法( t h eg e n e r a l i z e dm i n i m a l r e s i d u a lm e t h o d ，简称为g m r e s ) f 1 6 1 。当a 是非奇矩阵时，采用g m r e s 方法 至多迭代n 步就能得到解；而当a 是奇异矩阵时，情况就比较复杂，也许不能 得到解。p e t e r n b r o w n 和h o m e rf w a l k e r 在 1 2 1 中详细讨论了这种情况，并且 分析了在何时g m r e s 能够得到方程的最小二乘解。 1 9 9 9 年，d c a l v e t t i ，b l e w i s 和l r e i c h e l 三人又提出另一种类似g m r e s 方法的算法：限制值域的广义极小残量法( t h er a n g er e s t r i c t e dg m r e s m e t h o d ，简 称为r r g m r e s 方法) 1 3 。这种方法是用k m ( a ，a r 0 ) = s p a n a r o ，a 2 r o ，a r o ) 替换原先g m r e s 方法使用的k r y l o v 子空间k 。( a ，r o ) = s p a n ( r o ，a r o ，a ”1 0 0 ) 。 这样，当取初始向量x o = 0 时，迭代产生的解x 。k 。匕r ( a ) ，即x 。r ( a ) 。 在f 1 3 1 中讨论了当r ( a ) = r ( a 1 ) 时，r r g m r e s 方法能够产生范数最小的最小二 乘解。本文在此基础上证明了当r ( a ) r ( a 1 ) ，在适当的情况下，r r g m r e s 方 法也能得到方程的最小二乘解。 由 1 3 】中的定理32 容易知道：当r ( a ) = r ( a 1 ) ，且线性方程组是相容的， 则用r r g m r e s 方法可产生原方程的解。正因为如此，在本文的第三章采用 r r g m r e s 方法解s t o k e s 方程。用m a c 差分格式离散方程后，系数矩阵a 是 对称的( 即满足r ( a ) = r ( a 1 ) ) ，另外线性方程组是相容的，经过适当步数迭代 后就得到方程组的近似解。 第一章广义极小残量( g m r e s ) 方法 1 算法的介绍 g m r e s 算法是对一股非对称线性方程组 a x = b 的非定常迭代方法，其具体步骤如下： ( 1 ) 选初始迭代向量x o ，计算残量r 0 ：b - a x o 和初始正交向量v f 2 ) 对j = 1 , 2 ，直至满意 对j = l ，2 ，j h ，= ( a v ，v ，) h 川、，= l l 川f | ： v h = i 川 川， ( 3 ) 形成近似解x k x 。= 粕+ z k , 其中气k 。满足最小二乘条件 ( 1 o ) i i “r l ：；l i b - 爿( 屹) 峥快一4 气i i ：= m i n 骶一4 z l l ：zs 甄) ( 1 1 ) 也也可表示为 = 确十v , y 。，其中= v ，v ：，。 是由残量产生的伽如v 子空间， 即瓯；甄( 爿，) ：印口 ，0 ，4 ，爿“，0 ) 白勺正交基张成的矩阵，而_ y t 极小化，( y ) = 鼽虬中也牡勘盯副“，耻 + 孙卜圳“o 玎 月是上产k s s p n 6 p 愕矩阵： h = h l lh l2 h 2 fh 2 2 h 3 2 - 一h 】e 一h 2 1 h 3 k h i 一lh 女j 5 2 g m r e s 方法的性质 当( 1o ) 式是非奇线性方程组时，g m r e s 方法有如下性质：g m r e s 方法在第i 步产生的解x 。是精确解的充要条件为下述诸等价条件： ( 1 ) 算法在第j 步中断； ( 2 ) h j + 1 3 2 o ； ( 3 ) v j + 1 = o ； ( 4 ) r o 的最小多项式次数为j 确解 当( 1o ) 式是奇异线性方程组，用g m r e s 方法解时，情况较复杂，未必能得到精 令k k = k k ( a ，r o ) 是由初始残量r o 产生的k r y l o v 子空间，显然 d i m a ( k k 、三d i mk k 三k 定义1 1 若d i m a ( k k ) = k ，则称g m r e s 方法在第k 步无中断；若d i m a ( k k ) d i mk k ，称g m r e s 方法在第k 步是秩亏损最小二乘问题中断；若d i m k k 2 n , j 、二乘解 定理1 4 隐含的另一层意思是当a 不满足r ( a ) = r ( a 1 ) 时，且r a n k ( a ) d ，( 其 中d 是由如下定义：d i m k k ( a ，r o ) = k ，当k 三d ：d i m k k ( a ，r o ) 。d ，当k 兰d ) ，那么 g m r e s 方法可能不能产生( 1 0 ) 式的最小二乘解 第二章限制值域的广义极小残量( r r g m r e s ) 方法 1 算法的介绍 r r g m r e s 方法与g m r e s 方法相似，但它所用的k r y l o v 子空间和g m r e s 不同，r r g m r e s 使用的k r y l o v 子空间按如下定义： k m 三k m + ( a ，a r o ) = s p a n a t o ，a 2 。o ，a mr o ) ， m = 1 ，2 因此它与g m r e s 略有不同，具体算法如下 r r g m r e s 方法： ( 1 ) 选初始迭代向量x 0 ，计算残量r o = b - a x 。和初始正交向量v l = a r o 8a r oi 2 f 2 ) 对j = l ，2 ，直至满意 对i = 1 ，2 ，j h 。= ( a v ，v ，) i ，+ l = a v ，- z h v h 。，= h ，ij ： v 川= ；川 川 f 3 ) 形成近似解x k xk2 x n + y k yk 其中圪= v 。，v ：，。提由k ，y f o v 子空间k ：的正交基张成的矩阵，而儿极小化，( y ) = r i b - 爿( v 。y ) l t ：= 卜a v 。y i ：= ij r o 一圪+ 。j 孑乩 要极小化 ，( y ) = j i , o k 。7 洮 ( 2 o ) 等价于解法方程膏；v 。r 圪+ f i 。y = 百j 吆，r o ，而吃，圪+ = ，所以上面的法方程可 写为 蠢：两k y = a ：v 己、r o 另一方面，若要极小化 m ) = 懈r 0 一膏佻 ( 21 ) 也等价于解法方程 秘：商k y = 西：z 。r o 由此可知，要极小化( 2 0 ) 等价于极小化( 21 ) 下面考虑算法的具体实施情况 对最小二乘问题 i y l l n) yer j 的求解可用j 呼。的q r 分解来求得，但这里要求当疗。的每一列出现时更新厅。的 d r 分解以便于获得近似解的残向量范数，而决定是否停止计算。 用f l 表示在( e j ，e j + i ) 平面的g i v e n s 变换： = c | 一s i s ?c i 设f i ，i = l ，j ，己应用于厅，产生( j + 1 ) j 的矩阵豆， f i f 、hj 兰rj = x 0xx 0 o 下一步处理西。的最后一列，为获得巨。必须先前乘最后一列前面的j 个旋转。 由止e 得一f 歹0 形状的( j + 2 ) f j + 1 ) 矩阵 f ? f t h f + t = 0 o 0 其中h = h j + 2 j 而f j + 1 由 c = r ( r2 + 州7 2 j = 一h ( r 2 + 硎7 2 确定。自然f j 也需同时应用于吆 总之，k 步以后达到疗。的下述分解 由此 or 0矗 彤) = 瞻，0 一鼠y j | ： = 阱嵫。，0 一鼠巩 = 慨一砭y ： 困砭的最后一列为0 7 ，故极小值t j ( y 。) 即为或的最后一个分量的绝对值，这就是 的钱向量丘的范数恢 r r g m r e s 方法当k 上升时其计算量和存储量随之增加，为克服这个缺点，可 采用下述重开始算法 麓1 篱羰蓑蒋残量r o ：b - a x o 棚舡交向量v i ：a ，0 ”i ia r o ”i i 2 f ) 选初始迭代向量x o ，计算残量 2 和初始正父i 司重 。 ( 1 ) 对j = l ，2 ，m 对i = 1 ，2 ，j “ r 吆g = 甑 中蜘， i h = a v ，一h ，v n 川，= h n v 川= ；川 川， ( 3 ) 形成近似解x 。 x 。，= x o 十k ，y 。 其中y 。极小 2 - j ( y ) = l l r o + 。疗。y 0 ： f 4 ) 重开始 计算_ = b - a x 。，若 h | 1 适当小则停止，否则：= x 。，v ：= a r 。l l a r j ： 再转向( 2 ) 2 r r g m r e s 方法的性质 对于r r g m r e s 方法也有类似定理1 3 的性质 定理2 1若r ( a ) = r ( a 7 ) ，则对( 1 o ) 式应用r r g m r e s 方法可产生范数最 小的最小二乘解 那么当r ( a ) r ( a x ) ，又会得到怎样的结果呢? 下面就来说明定理2 2 定理2 2 若r a n k ( a ) = m l ，d i m ( a k m - i * ) = m 一1 ，则对( 1 o ) 式应用r r g m r e s 方法，在m 步中断时产生最d , - - 乘解 证明：r r g m r e s 方法在m 步中断，即d i m ( km 1 ) = d i m ( k 。) = m - 1 若不存在最4 , - - 乘解x k m - i ，则a t b 诺a t a k m - i + 因为d i m ( a 1 a k m - i * ) = d i m ( a k m 1 ) = m l ， 所以d i m a ls p a n ( b ，a 2b ，a b ) = m 然而r a n k ( a 1 ) = r a n k ( a ) = m - l ， 所以d i m a l s p a n b ，a 2b ，a b ) 三r n 1 矛盾! 因此，r r g m r e s 方法在m 步中断时产生( 1 o ) 式最小二乘解x e k m 1 ( 0 ) 注：因为d i m ( k m - l * ) = m 1 ，条件d i m ( a k m - l * ) = m 1 等价于n ( a ) n k m - i + = 一二，a = 7 ，。= ( 1 ) 容易知道满足定理2 2 的所有要求，r r g m r e s 方法在第2 步中断，且产 生最小二乘解x ： i 3 ，1 3 1 由于k 。( a ，a b ) c _ k 。+ l ( a ，b ) ，所以能够应用r r g m r e s 方法的情况比应 用g m r e s 方法少 例2 3 ，肛| j j ，6 通过计算得到：k 。+ = s ”a n 仃 订铆州例 k 3 + = s p a n ， ； 删 因此k l + k2 + = k3 + ，r r g m r e s 方法在m 23 时中断而 肌= ( | f f 2 a 7 a k2 = s p a n l3 心1 因为a 7 b a t ak 2 + ，所以不存在x ek 2 * 满足法方程a t a x = a t b ，即在此例子 r r g m r e s 方法不能产生最小二乘解 但应用o m r e s 方法却能产生解由计算可得 即鲫俐， 2 k 2 = s p a n k 3 2s p a n j 势 封蛰 f t k 4 = s p a n l0l “3 | ， i l 2 j l j 耕 ； ， 耕 k 2 k3 = k4 ，g m r e s 方法在第4 步中断因为a k 2 = a k3 ( 即k2 + = k3 + ) ，我 们有a ak 2 = a ak 3 所以g m r e s 方法可产生最小二乘解x e k2 ( a ，b ) 计算 得到( 1 0 ) 式最小二乘解是x = 3 8 一5 ，9 ，一1 3 1 1 k2 ( a ，b ) ，而范数最小的最小二乘 解是x + = 3 2 ，0 ，3 2 7 k3 ( a ，b ) _ 3 “h 第三章s t o k e s 方程的一种块l u 预条件算法 1 问题的介绍 ( 30 ) 的鞍点问题，其中a 和b 分别是n a n a 和n a x n b 矩阵，且m b 三n a 。如果a 是对称正定矩阵，这种问题称为是古典问题；如果a 是不定或非对称的，则称是 一般问题。这样的系统方程来自于s t o k e s 和n a y i e f s t o k e s 方程，离散后的方 程通常是稀疏的大矩阵，那么直接解是不可行的。要得到方程( 3 0 ) 精确数值 解的主要困难存于这种系统是不定的，这可从以下l d l 分解得出，其中一s = 1 5 a 。b 是s c h u r 补矩陌：。 肾( 硝鼻嬲瑞爿- ，i 刁 b - ， 如果时这类系统不使用预条件，那么k r y l o v 子空间方法效果很差，如f 7 ，因为 这类是不定系统。对于方程( 3 0 ) 一般有如下三种主要方法： ( 1 ) c o n d e n s a t i o n 方法。从两个等式中消去变量u ，得到s c h u r 补系统 s p = b t a 。1 f g 从这儿精确解出p ，则u 可从以下等式得到 a u = f - b p 从f 5 ，9 1 两篇文章中可看出这种方法化费很大，因为得到s 先要精确解出a 的解， 再需乘以矩阵b 和b 7 当a 是l a p l a c i a no p e r a t o r 时，a 的解能够有效的得到 但! ja 足一般的形式，则解不容易得到。 ( 2 ) u z a w a 算法 2 。如果内迭代几乎是精确的，则称为精确的u z a w a 算法； 在文章【1 1 】中说明精确的内迭代通常是不必要的，那样的算法就称为不精确的 u z a w a 算法。对于对称的情况及不对称情况分别在文章【3 】， 4 】中有介绍。 f 3 ) 对这系统作为整体来做预条件。利用距阵的块性质来设计有效的预条件， 相对柬浇较简单，因为不必要得到矩阵a 的精确解。一般，预条件有形如以下 的块对角和块三角形式 姬 厂g眦h 需 “， 壬、l， 碰 口o 常经 。矿 中题司 j 一 多在 勒弓= 鼻 若a 是对称正定的，r u s t e n 和w i n t h e r 6 介绍了块对角预条件方法，对a 作 预条件采用不完全c h o l e s k y 分解，对s 则可取j = i 或s = b 1 b 对于稳定 s t o k e s 方程，文章 8 】和 1 0 采用了块对角预条件方法，先是使用对角矩阵来做预 条件，后又采用更一般的预条件如m u l t i g r i dv - c y c l e s 方法。 本文对于s t o k e s 问题，利用m a c 差分格式离散后得到形如( 30 ) 的方程， 再采用类似第( 3 ) 种的预条件方法分别用g m r e s 方法和r r g m r e s 方法解此方 程。 2 m a c 差分格式 这篇文章主要是计算s t o k e s 方程，也即计算如下方程： v n + v p = f i nq ( 3 ) v u = 0i nq 其中q 是实空间r 2 的有界区域，向量值函数u 表示在q 上的速度，它有两个分量 u ；( u ，v ) 7 ，标量值函数p 表示压力，常数v 是粘性系数，与雷诺数( r e y n o l d s n u m b e r ) 成倒数关系，l a p l a c e 算子表示方程扩散部分。且定义速度的第一个分量u 在边界y = 1 上为l ，其余边界上都等于o ；第二个分量v 在所有边界上值都为0 ， 如图l ： f i g u r e l ：b o u n d a r yc o n d i t i o n so f s p e e dv e c t o r 应用m a c 有限差分法的优点是：离散过程相对有限元离散简单，且无需稳 。 0 】 f0 l。 fh。t 定化能够广泛的用更多的迭代法去解离散后的线性方程组。下面就简单介绍 m a c 有限差分法。 假设q 是( o ，1 ) ( o1 ) 正方形区域，被平均分成n x n 个正方形单元网格，且其边 长h ：1 n 用以下方法选取速度解和压力解在网格上的离散点：速度的第一个分 量u 的离散点定义在每个单元格子的左右两条边的中点，第二个分量v 的离散点 定义在每个单元格子的上下两条边的中点，压力的离散点则在每个单元格子的中 心。如h = 1 4 的网格，速度和压力的结点如图2 所示( 其中符号表示速度分量u ， 符号。表示速度分量v ，符号表示压力p ) ： uuuu ： ： 一一一 一 ： 一 一一 一 ：，： 一 一一 、j ： ： ： ： f i g2 s t a g g e r e dg r i df o rt h em a cf i n i t e d i f f e r e n c ed i s c r e t i z a t i o n 速度u = ( u ，v ) ，可以把方程组( 3 ) 写成以下形式 一v ( 窘芬心= 一叫【萨+ 矿j + 瓦叫1 一v 侈豺考= 1 l 萨+ 矿j + 方叫2 塑一8 v 0 苏 砂 其中( f l ，f 2 ) = f ， 那么按速度和压力在网格上的排列位置，对( 3 3 ) 进行如下中心差分 一au 】j k 。 a l u j k l h 2 ( 4 u j k u j 1 ，k u j + i ，k u j , k 1 一u j ，k 1 ) p x j k 2 b l p j k 51 h ( p j + i 2 k p j 1 2 ，k ) u 。+ v y b k 【b t u j k 三一 1 h ( u j + l 2 ，k u j i 2 k ) + l h ( v j ，k + l 2 一v j , k i 2 ) 】 ( 3 3 ) 就得到下面的矩阵形式 旧蒯= 其中爿= f 鲁爿0 ， ，且是对称正定的，b = f 篁 3 块l u 预条件 考虑一般情况，整个矩阵是非对称的 呤( 0 ： ( 3 4 ) 对这矩阵采用l u 分解后得到： 肾( 瑞爿- ，1 勺 s ， 其中s = 一c v a b 有这分解暗示可采用如下形式的块l u 预条件： 拈b 瑞夤一 s ， 其中五和s 分别是a 和s 的近似矩阵。如果形如j x = b 和雪y = c 的系统能够 不费时的解出，就可采用块预条件( 3 , 6 ) 。 对于系统 砌= ( 习 可从以下算法解得。 算法i 块l u 解过程 ( 1 1 解jx = w ( 2 ) 解雪y = c t x z ( 3 ) 解- t = b y ( 4 ) 计算x = x t 以e 算法涉及到矩阵j 的两个解和矩阵雪的一个解。 媾于( 36 ) 式的块l u 预条件需要两个成分：f 1 ) 矩阵a 的近似一，且爿的 解容易得剑：( ! ) 矩阵s 的近似s ，对于s 同样也可方便得剑解。因为a 通常可得 到，所以它的近似矩阵可方便计算出，如用不完全c h o l e s k y 分解得到爿。这方 2 ：关键是怎样近似矩阵s ，旦这近似矩阵选好，再适当分解s 。对矩阵s 的合 适分解也可作为其他分解积的形式直接算出，而不必先形成蜃。 s c h u r 补矩阵的近似 因为a 一般是稠密矩阵，所以通过a 。来计算毒是不实际的。两个较易得 j , i j 的方法足： s l ；c 1 b ， 2 。c t d i a g ( a ) b 但是这两个近似矩阵不足够精确，下面介绍另一种方法来近似s c h u r 补矩阵。 对于对称鞍点问题即c = b 而且a 是对称e 定矩阵，则采用不完全的 c h o l e s k y 分解得到a ， 五= l l t 通过这个分解来计算近似s c h u r 补矩阵： 鸢3 = b t ( l l l ) 。1 b = b l 玎l - 】b = x t x 其中x = l - i b 若分别令l = i ，l = d i a g ( a ) ，则就变为叠l 和j2 形式。 i f 如算法1 所示，颁条件的应用需要解形如jx = w 和s y = c 的方程，而这 些义是用来定义全局的顸条件，所以不需要得到精确的解。这两个方程可以用迭 代法汁算，用j 作为j 的预条件，取j 作为j 的预条件。然而，用迭代法解这些 方程比较花时间。在一种极端情况即采用精确的不完全分解，不用迭代法，而是 直接在方程中使用顸条件算子。这个方法也节省了存贮空间，因为一旦季计算出， j 就可丢弃了。 ：计算x 因了 或许会误以为矩阵x 是稠密的订算它会很花时问。事实上，x 常常是一个 相当稀疏的矩阵，因此计算矩阵x 的近似是较快的。f 面介绍计算矩阵x 的 方法，其中x 。表示矩阵x 的第i 行，b 表示矩阵b 的第i 行，l i 表示矩阵l 的 第i 行。 算法2 计算矩阵x ( 1 ) f o r 滓1 ，n a d o ： ( 2 ) x i ：= b i f 3 1f o rk = l ，2 ，i - 1a n di f l i k 0 d o ( 4 )x ：= x l l i k + x k f 5 1 e n d d o ( 6 ) x i ：= x l i i f 7 1e n d d o 矩阵x 的第i 行也即通过原来第i 行减去前面i - 1 行线性组合，再除以矩阵l 对角元l 。但对于标准l u 分解，l 。= l ，所以就不必除以l i i 了。以上算法可用如 下图3 表示： l x f i g u r e3 ：r o w o r i e n t e dc o m p u t a t i o n o fl 1 b 在此需要指出的是：x 的前几行是稀疏的，随后变得越来越稠密，所以对 一些元素要丢弃。一般有两种方法： ( 1 ) 完成算法2 第3 5 行后，那些小于t o l 的数就置为0 ( 2 ) 这是一种比较好的方法，先介绍怎样定义矩阵元素的l e v e l 值。对矩 阵m ，如果m 。= 0 ，则l e v ( m i ) = o o ；如果m u 0 ，则l e v ( m , j ) 。0 。 在高斯消去法中，每个元素都以如下公式更新： m u2 m u m i k4 m q 那么每个元素的l e v e l 值也随即更新，公式如下： l e v ( m i j ) = m i n l e v ( m 0 ) ，l e v ( m - k ) + l e v ( m k j ) + 1 其中当i j ，m i j 代表矩阵l 的元素；当i 三j ，m 。代表矩阵u 的元素。 不完全l u 分解是由通常的高斯消去过程构成，当元素的l e v e l 值大于p 时就丢弃。 假设对a 进行不完全c h o l e s k y 分解后得到矩阵l ，且l 的每个元素已 定义了各自的l e v e l 值。下面是对算法2 修改，丢弃了l e v e l 值大于p 的元 素。 算法3 修改后的计算矩阵x 的算法 ( 1 ) f o ri - l ，n ad o ： ( 2 ) x i = b l ；l e v ( x o ) = ( 0i f b i j - 0 ；c 。i f b i j2 0 f o r j 2 1 ，2 ，m b ( 3 ) f o rk = 1 , 2 ，i - 1a n di f l i k od o ： ( 4 ) x i ：2 x i l i k + x k ( 5 ) l e v ( x 0 ) 2m i n l e v ( x i j ) ，i e v ( l 。k ) + i e v ( x 日) + 1 j 2 1 ，m b ( 6 ) e n d d o ( 7 ) d r o pe a c hx 0s u c ht h a tl e v ( x , j ) pf o r j 。1 ，m b ( 8 ) x i = x i l i i f 9 ) e n d d o 三对矩阵雪的预条件 在计算矩阵x 中，当p 越来越大时，则以上算法能产生一个较精确的s c h u r 补 矩阵。然而，在应用预条件时，还需计算矩阵s 的解。这时可对做预条件，方 法如下： ( 1 ) 不完全c h o l e s k y 分解 ( 2 ) 不完全l u 分解 ( 3 ) 不完全q r 分解，即x 一- q r ，而q 1 q i ，因此x 1 x z r 7 q 1 q r - - 一r t r 。 这个方法也节约了内存，因为当计算出矩阵r 后，q 就不再需要，可 以丢弃了。 2 0 综上所述，当a 是对称正定，m 。是对称的，算法可描述为： ( 1 ) 通过不完全c h o l e s k y 分解计算a = l l ( 2 ) 用算法3 计算l x = b 的解 ( 3 ) 通过x x 的不完全c h o l e s k y 分解或不完全l u 分解计算s ，或通过矩 阵x 的不完全q r 分解计算sm r r ( 4 ) 应用算法l ，分别用a ，s 代替a ，s 对于般鞍点问题，对a 不能用不完全c h o l e s k y 分解，这时就可采用不完全 l u 分解，则定义量3 = c 1 ( l u ) 。1 b = c 7 u “l b = y t x ，其中x = l j b ，y = u 。t c 第四章数值试验 f 面取不同的n 和v ( 粘性系数) 做试验，分别用g m r e s 方法和r r g m r e s 方 去、目都采用本文第三章的第3 节介绍的块l u 预条件。得出的结果如下两表 所吖：，其中表i 是g m r e s 乃法，袁：是r r g m r e s 力法( 逗号莳是得到近似 解所需的h i 问，逗号后是迭代的次数) ： n 2 8n = 1 6 n = 3 2 n = 6 4 f 1 7 6 阶1( 7 3 6 阶、 ( 3 0 0 8 阶1 ( 1 2 1 6 0 阶) v = l o 0 1 4 ls ，1 5 次o 3 6 0 s ，2 9 次2 8 2 8 s ，5 8 次4 6 4 3 8 s 、1 3 4 次 v ：1 0 ： 0 1 2 5 s ，1 7 次0 3 7 5 s ，3 1 次2 8 9 0 s ，5 9 次4 8 3 2 8 s ，13 3 次 v = 1 0 10 1 2 6 s i 9 次 0 4 0 5 s ，3 5 次 34 0 7 s ，7 0 次 5 3 2 3 5 s 】4 9 次 v ：1 0 1 0 1 4 1 s 2 2 砍o 4 3 8 s ，3 7 次3 ，7 8 1 s 7 8 次6 5 9 8 4 s 1 7 5 次 、1 = 1 0 1 0 1 4 2 s 2 4 次o 4 5 3 s 4 0 次3 9 5 3 s ，8 4 次7 2 9 0 6 s 18 8 次 表l n 。8n = 1 6n = 3 2n = 6 4 ( 17 6 阶)( 7 3 6 阶)( 3 0 0 8 阶1( 1 2 1 6 0 阶) v = 1 0 10 1 0 9 s 。18 次 0 3 5 3 4 次3 1 8 7 s ，7 j 次6 4 6 8 0 0 s 2 7 2 次 v = 1 0 20 1 0 9 s ，2 0 次o ，3 7 5 s 3 6 次 3 2 5 0 s ，7 2 次6 4 6 2 5 s ，2 6 2 次 v2 1 0 一 0l l o s ，2 2 次o 4 2 2 s 3 9 次3 6 7 2 s 、8 0 次 6 9 6 2 5 s ，2 8 8 次 v = l ( ) 10 1 0 9 s 。2 4 次o 4 5 3 s ，4 2 次4 0 6 2 s 、8 7 次 8 4 1 2 5 s 3 6 7 次 v l ( ) 10 1 1 0 s 2 6 次o 4 6 9 s ，4 5 次4 3 4 3 s 、9 3 次 8 8 4 6 8 s ，3 8 4 次 表2 从上面两张表得出：随着粘性系数的减少，r r o m r e s 方法和g m r e s 方法 两爵所需的迭代次数都会增加，相应的，花费的时间也就多。一般的，当n 和 v 幽定不变，r r g m r e s 力法比g m r e s 方法所需的迭代次数多，时问也多一些， 特别当n 较大时，两者之间的差别就越明显。 下面再分别对n = 3 2 ，6 4 ，以及v = 1 0 ，1 0 ，1 0 ，1 0 4 作图，其中x 轴表 示迭代步数，y 轴表示蠛向量范数，图中实线代表g m r e s 方法所得的结果，而 虚线代表r r g m r e s 方法所得的结果。n j _ 从以下的八张图上很清楚看剑， r r g m r e s 力法相对于g m r e s 方法收敛速度慢一些，但差别不是非常大。 藕 谰 删 坦 鬣 饕 型 鬣 蕤 理 躐 鬟 鐾 激 2 4 蓁 型 鬣 豢挺硎定暇 参考文献 【1 l e i g hl i t t l e ，y o u s e fs a a d ， b l o c kl up r e c o n d i t i o n e r sf o rs y m m e t r i ca n d n o n s y m m e t r i cs a d d l ep o i n tp r o b l e m s ， 【2 】k a r r o w ，l h u r x i c z ，h u z a w a ，s t u d i e si n n o n l i n e a rp r o g r a m m i n g ，s t a n f o r d u n i v e r s i t yp r e s ss t a n f o r d ，c a ，1 9 5 8 3 】j h b r a m b l e ，j e ? a s c i a k ，a t v a s s i l e v ，a n a l y s i so ft h ei n e x a c tu z a w aa l g o r i t h m f o rs a d d l ep o i n tp r o b l e m s ，s i a mjn u m e r , a n a l ，5 ( 19 9 7 ) ，p p10 7 2 10 9 2 4 jh b r a m b l e ，a tv a s s i l e v ，j e p a s c i a k ，u z a w at y p ea l g o r i t h m sf o rn o n s y m m e t r i c s a d d l ep o i n tp r o b l e m s ，t oa p p e a r , m a t h c o m p 5 n d y n ，w f e r g u s o n t h en u m e r i c a ls o l u t i o no fe q u a l i t y c o n s t r a i n e dq u a d r a t i c p r o g r a m m i n gp r o b l e m s ，m a t h c o m p u t ，4 1 ( 1 9 8 3 ) ，p p 1 6 5 1 7 0 【6 】t r u s t e n ，r w i n t h e r ， ap r e c o n d i t i o n e di t e r a t i v em e t h o df o rs a d d l e p i o n t p r o b l e m s ，s i a mj m a t r i xa n a l a p p l ，1 3 ( 1 9 9 2 ) ，p p 8 8 7 9 0 4 7 y s a a d ，i t e r a t i v em e t h o d sf o rs p a r s el i n e a rs y s t e m s ，p w s ，n e wj e r s e y , 19 9 6 【8 】dj s i l v e s t e r ，aj w a t h e n ，f a s t i t e r a t i v es o l u t i o no fs t a b i l i s e ds t o k e ss y s t e m s p a r t 1 i ：u s i n gg e n e r a lb l o c kp r e c o n d i t i o n e r s ，s i a mj n u m e r a n a l ，3 i ( 1 9 9 4 ) ，p p 1 3 5 2 一 1 3 6 7 【9 rv e r f u r t h ， ac o m b i n e dc o n j u g a t e g r a d i e n t m u l t i g r i da l g o r i t h m f o rt h e n u m e r i c a ls o l u t i o n o ft h es t o k e sp r o b l e m ，i m aj o fn u m a n a l ，4 ( 1 9 8 4 ) ， p p4 4 1 4 5 5 i o a j w a t h e n ，d j s i l v e s t e r , f a s ti t e r a t i v es o l u t i o no fs t a b i l i s e ds t o k e ss y s t e m s p a r ti ：u s i