（计算数学专业论文）大型稀疏线性方程组的krylov子空间方法.pdf

中文摘要 很多科学工程计算问题都转化为求解大型稀疏的线性方程组血= b ，如偏微 分方程组的差分格式，有限元方法离散得到的刚度矩阵等等，系数矩阵a 都具 有大而稀疏的特点。由于问题的规模往往非常大，因此迭代法成为求解大型稀疏 的线性方程组a x = b 的最常用的方法之一迭代法的基本思想是：从解的某个近 似值出发，通过构造一个无穷列去逼近精确解( 一般有限步内得不到精确解) 。 与直接法相比，迭代法能保持矩阵的稀疏性，具有计算简单，编制程序容易的优 点，并在许多情况下收敛较快，因而能有效地解大型稀疏的方程组。 本论文将回顾一些迭代法，包括i c c g 、c g s 、b i c g s t a b 等k r y l o v 子空间方法， 并探讨是否会发生中断的现象。再从中挑选几种算法来求解一些利用有限差分法 离散化偏微分方程式所化为的大型稀疏线性方程组，并比较其迭代的时间与收敛 的速度。 第一章，我们简单介绍研究k r y l o v 子空问法的背景，动机及目的。 第二章，介绍k r y l o v 子空间的定义、定理以及一些相关名词，并探讨子空 间与矩阵的关系来作为各种k r y l o v 子空间法发展的基础。 第三章，介绍最近几年比较有名的k r y l o v 子空间法。 第四章，数值实验。 第五章，总结本文工作。 关键词：大型稀疏线性方程组，k r y l o v 子空间法，迭代法。 a b s t r a c t m a n ys c i e n t i f i ca n de n g i n e e r i n gc o m p u t a t i o n a lp r o b l e m sh a v eb e e n c h a n g e di n t oal a r g es p a r s el i n e a re q u a t i o na x = b ，s u c ha sd i f f e r e n c e f o r m a to fp a r t i a ld i f f e r e n t i a le q u a t i o n s ，s t i f f n e s sm a t r i xo b t a i n e df r o m d i s c r e t i o nb yu s i n gf i n i t ee l e m e n tm e t h o de t c d u et ot h es c a l eo ft h e p r o b l e ma r eo f t e nv e r yl a r g e ，t h ei t e r a t i v em e t h o db e c o m e so n eo ft h em o s t c o m m o n l yu s e dm e t h o d sf o rs o lr i n gl a r g es p a r s eli n e a re q u a t i o n sa x = b t h e b a s i ci d e a o fi t e r a t i v em e t h o di sa sf o l l o w s ：s t a r t i n gf r o m s o m e a p p r o x i m a t i o n o ft h es o l u t i o n ，w ec o n s t r u c ta ni n f i n i t e s e r i e st o a p p r o x i m a t et h ee x a c ts o l u t i o n s ( g e n e r a ll ys o l u t i o nw o u l dn o tb et h ee x a c t s o l u t i o nw i t h i ns e v e r a ls t e p ) c o m p a r e dw i t ht h ed i r e c tm e t h o d ，i t e r a t i v e m e t h o dc a nm a i n t a i nt h es p a r s i t yo ft h em a t r i x ，s i m p l ec a l c u l a t i o n ，t h e a d v a n t a g e so fe a s yp r o g r a m m i n g ，a n di nm a n yc a s e s ，r a p idc o n v e r g e n c e ，t h u s c a ne f f e c t i v e l ys o l v el a r g es p a r s ee q u a t i o n s i nt h i sp a p e r ，w ew i l lr e v i e ws o m ef a m o u si t e r a t i v em e t h o d so fk r y l o v s u b s p a c es u c ha si g c g ，c g s ，a n db i c g s t a ba n dd is c u s sa tw h a t k i n d o f s i t u a t i o nt h eb r e a k d o w nw i l lo c c u r ，t h e ns e l e c ts e v e r a la l g o r i t h m st o s o l v el a r g e - s c a l es p a r s e1i n e a re q u a t i o n s ，o b t a i n e df r o md i s c r e t i z a t i o n o fp a r t i a ld i f f e r e n t i a le q u a ti o n sb yu s i n gt h em e t h o do ff i n it ed if f e r e n c e ， a n dc o m p a r et h e i ri t e r a t i o nt i m e sa n ds p e e do fc o n v e r g e n c e f i v ec h a p t e r sa r ef o l d e di nt h isp a p e r i nc h a p t e ro n e ，w es i m p l yi n t r o d u c et h eb a c k g r o u n d ，m o t i v a t i o na n d p u r p o s eo fk r y l o v s u b s p a c em e t h o d s i nc h a p t e rt w o ，w eg i v et h ed e f i n i t i o na n ds o m et h e o r e m s o fk r y l o v s u b s p a c e a n dr e v i e ws o m er e l a t i o n sb e t w e e ns u b s p a c ea n d m a t r i xt od i s c u s s t h ed e v elo p m e n to fk r ylo vs u b s p a c e i nc h a p t e rt h r e e ，w ei n t r o d u c es o m ew e l 卜k n o w nm e t h o d so fk r y l o v s u b s p a c e i nc h a p t e rf o u r ，w eg i v es o m en u m e r i c a le x p e r i m e n tt od e m o n s t r a t et h e a d v a n t a g e so ft h e s em e t h o d s i nc h a p t e rf i v e ，w eg i v eac o n c l u s i o no ft h i sp a p e r k e yw o r d s ：l a r g es p a r s es y s t e m so fe q u a t i o n s ；k r y l o vs u b s p a c em e t h o d s ；i t e r a t i v e m e t h o d 2 厦门大学学位论文原创性声明 本人呈交的学位论文是本人在导师指导下，独立完成的研究成 果。本人在论文写作中参考其他个人或集体已经发表的研究成果，均 在文中以适当方式明确标明，并符合法律规范和厦门大学研究生学 术活动规范( 试行) 。 另外，该学位论文为() 课题( 组) 的研究成果，获得() 课题( 组) 经费或实验室的 资助，在() 实验室完成。( 请在以上括号内填写课 题或课题组负责人或实验室名称，未有此项声明内容的，可以不作特 别声明。) 声明人( 签名) ：邗锄孝 矿哗 二只q e t 厦门大学学位论文著作权使用声明 本人同意厦f - j ) k 学根据中华人民共和国学位条例暂行实施办 法等规定保留和使用此学位论文，并向主管部门或其指定机构送交 学位论文( 包括纸质版和电子版) ，允许学位论文进入厦门大学图书 馆及其数据库被查阅、借阅。本人同意厦门大学将学位论文加入全国 博士、硕士学位论文共建单位数据库进行检索，将学位论文的标题和 摘要汇编出版，采用影印、缩印或者其它方式合理复制学位论文。 本学位论文属于： () 1 经厦门大学保密委员会审查核定的保密学位论文， 于年月日解密，解密后适用上述授权。 () 2 不保密，适用上述授权。 ( 请在以上相应括号内打“”或填上相应内容。保密学位论文 应是已经厦门大学保密委员会审定过的学位论文，未经厦门大学保密 委员会审定的学位论文均为公开学位论文。此声明栏不填写的，默认 为公开学位论文，均适用上述授权。) 声明人( 签名) ： 年月日 第一章绪论 1 1 引言 很多科学工程计算问题都转化为求解大型稀疏的线性方程组a x = b ，如偏微 分方程组的差分格式，有限元方法离散得到的刚度矩阵等等，系数矩阵a 都具 有大而稀疏的特点。由于问题的规模往往非常大，因此迭代法成为求解大型稀疏 的线性方程组a x = b 的最常用的方法之一迭代法的基本思想是：从解的某个近 似值出发，通过构造一个无穷列去逼近精确解( 一般有限步内得不到精确解) 。 与直接法相比，迭代法能保持矩阵的稀疏性，具有计算简单，编制程序容易的优 点，并在许多情况下收敛较快，因而能有效地解大型稀疏的方程组。 在许多科学工程上的应用问题，如飞行器模拟等，在像模拟这样的连续事件， 将之离散化之后会得到一个线性方程组a x = 6 ，其中系数矩阵通常很大且稀疏。 因此如何快速且有效率的求解这样的线性方程组变得十分重要。 而这样的线性方程组一般都利用迭代法( 或称非直接解法) 来求解。早在十 九世纪初，就有许多学者对用迭代法来求解线性方程组很感兴趣。这类型的迭代 法最基本的想法就是将给定的线性方程组a x = b 中的系数矩阵a 分解成两个矩 阵的和，其中一个矩阵将会使得给定的线性方程组容易求得其真解。例如，最简 单的分法是a = i 一( j a ) ，这个分法得到的迭代法就是著名的r i c h a r d s o n 法： 给定初始值z ( o x = ( i a ) x h + 乃i = 1 ，2 ，j ，3 一 此类型的迭代法都充分利用稀疏矩阵的特性来避免执行零元素的加乘运算 及浪费电脑存储空间去存储零元素。但当方程组很大时，这些方法就显得较没效 率。因为并无法保证在有限的迭代内求得线性方程组的真解，且常花很长的时间 来逼近真解。以上述的r i c h a r d s o n 法来说，能保证对于任给初始值要收敛的必 要条件是矩阵，一彳中的最大特征值的模必须小于l ，而这个模的值越小，收敛的 速度越快。e w a l db o d e w i g 在他的书 1 中提到：除非系数矩阵a 是接近对角矩 阵，否则这类型的迭代法并不实用 在1 9 5 2 年左右，h e s t e n e s 与s t i e f e l 2 2 提出了另一种形式的迭代法( 共轭 梯度法) 。虽然这种方法仅限于求解对称j 下定的矩阵但其收敛的速度却快于之前 的基本迭代法。此法在迭代过程中构建了一序列的子空间，称之为k r y l o v 子空 间。不久，l u e n b e r g e r 3 最先提出求解对称但非j 下定方程组，之后又有p a i g e 与 s a u n d e r s 【4 】，b u n c h 与k a u f m a n 5 】，f f i d m a n 【6 】及s t o e r 与f r e u n d 【7 】陆续提 出有关求解对称但非正定方程组的方法。 之后许多学者发明了一些方法去扩展c g 法，最早如c o n c u s 与g o l u b 8 】 跟w i d l u n d 9 】给了一种广义的共轭梯度法，其中系数矩阵是实的萨定矩阵。而 v i n s o m e 1 0 提出一种n o r t h o m i n 的方法与y o u n g 和j e a 1 l 】 1 2 】提出广义的共 轭梯度法( i g c g ) 中的o r t h o m i n 中的方法是相似的。而其他和o r t h o m i n 方法类似 的有a x e l s o n 【1 3 】以及由e l m a n 1 4 提出的迭代法。 近几年来又有很多新的迭代法陆续地被提出，这些方法的确比早期的基本迭 代法收敛要来得快( 见表一) 。本篇论文将回顾k r y l o v 子空间法的构建方式，并探 讨各类型的k r y l o v 子空间法会发生中断的可能性。 方法迭代次数 g a u s s s e i d e l超过1 0 0 0 0 0 k r y l o v 子空间法 大约3 0 0 10 0 0 表一：求解第四章的偏微分方程式p 1 所需的迭代次数 4 1 2 研究目的 在当今科技上总是希望获得偏微分方程之高精确度的数值解而且求解所需 的时间短，如此才能获得最大的经济效益。本文将利用几种类型的k r y l o v 子空 间法求解3 个边界值问题的偏微分方程并比较其差异、解题速度以及收敛快慢。 第二章构建k r y lo v 子空间法 2 1k r y l o v 子空间与矩阵的基本性质 本章将介绍k r y l o v 子空间的定义、定理以及一些相关名词，并探讨子空间与 矩阵的关系来作为各种k r y l o v 子空间法发展的基础。 定义2 1 令v r ，由v ，a v ，a ”1 v 所生成的子空间称之为由1 ，与 彳所生成的第m 个k r y l o v 子空间，并记k ( v ，么) 。 借由k r y l o v 子空间的定义可以得知对于任意的，l 大于1 ，可以得 到e ( 1 ，a ) 互e + 。似彳) 。因为尺空间中最多仅有个线性无关量，因此最后由 ，跟彳所生成的向量集最后会成线性相关，将1 ，彳v ，a 4 v 成线性无关的最大j 下 整数记为d ，称d 为v 对于彳的指数。由指数的定义可以知道d 必定小于或等 于一l 。借由基本迭代法可以推得对于任意的n 大于零， k ( r ，彳) = s p a n x n x m ，x 孙一石，工一工o ) 。因此得知x 一工o 会属于 k n ( ，彳) 中，因此也会属于k a + i ( ，们，彳) ，但是 x x o 是否也会属于 k d + ( ，叭，彳) ，下面这个定理将给予这个保证。 定理2 2 令，对于a 的指数为d ，若i - - x o k 。( r 们，a ) ，则m d + l 。 证明见 1 1 。 既然x - - x o 巧+ l ( r 们，彳) ， 必然存在口r “1 满足 ；= x ( o + ( r ( 们，a r 们，a d ，o ) 云 ( 2 1 ) 也就是；可以借由求取磊得到，可惜，( o 对于a 的指数d 难以得知，所以必须借 由一个矩阵z ，并经由迭代的方法，从超平面x o + k 。( ，么) 中求取满足 g a l e r k i n 条件的向量；也就是从x o + k 。( ，们，a ) 中求取近似解x ”使得其残向量 6 ，”与子空间z k 。( ，们，彳) 垂直。所求的z 须满足对于任何的，l = l ，2 ，d + l ，线性 方程组 ( ，o ) ，彳，们，a ”1 ，o ) z r o = ( ，m ，爿厂们，a ”一1 ，) z 4 ( ，们，彳，们，a ”一1 ，o ) 口“ ( 2 2 ) 有唯一解。因为，= b - a x = 0 且x = x o + 髟+ l ( 厂0 ，彳) ，显然z 是 工o + 巧+ i ( ，a ) 中满g a l e r k i n 条件的一个向量，而在什么条件下才能使 x ( n 1 ) = 工成立，将在后面再予讨论。 为了从k r y l o v 子空间求得较好的解，因此需要找出一个较合适的基来作为 k r y l o v 子空间的基。根据数值的观点用，们，彳，们，a 川r o 作为k 。( ，彳) 的基是 相当不适合的，因为当n 够大时a ”，0 会趋近于系数矩阵中最大的固有值所对应 的固有向量。因此可以在迭代的过程中建立一组利于求取近似解x 伽的基 仞们，p n ，p ”) 来取代原先的基 厂们，彳，们，a ”1 ，l o ) ，在【1 5 中提到两种方 向向量的取法常常被拿来作为k r y l o v 子空间的基。 由于定理2 2 保证x 一i f o 髟+ l ( 厂们，彳) ，但是指数d 不知道，为了确保n = d + l 时x ”= x 会成立，因此要求对于每个刀= 0 ，l ，d ，近似解 工斛1 = x o + ( p ，p n ，p ”) 口” ( 2 3 ) 唯一确定。由于p 们，p n ，p ”是一组线性无关的向量，所以满足( 2 3 ) 的是1 2 “ 唯一的。 前面提到利用一个矩阵z 来求取口枷，到底z 在什么条件下能保证在满足 g a l e r k i n 条件下，所导出的线性方程组 ( p 们，p n ，p ”) z r o = ( p 们，p ，p ”) z 4 ( p 们，p ”，p ”) 口” ( 2 4 ) 有唯一解，且在n = d 时，( 2 4 ) 式的唯一解口。能够满足 工= 工o + ( p 们，p ，p ”皿4 也是唯一的，因此下面介绍矩阵与子空间的关系， 首先先介绍两个定义： 定义2 3 - 令日为n x n 的矩阵且嵋，屹，是r 空间的向量。假如 ( k ，_ ) = o l f 甩 ( 2 5 ) 则v l , ，k 为r 中的日- 半正交向量。 定义2 4 给定一个n xn 的矩阵日，若能使得对于任意非零的n x l 的实向量 z 都能满足( x ，x ) 0 ，则称此矩阵h 为非零实矩阵。 下面这个定理说明了对于任意的线性无关组都可以利用广义的 g r a m s c h m i d t 过程将它转换成日半正交的线性无关组。 定理2 5 令日为n xn 的非零实矩阵，而且 崩，见) 是尺中的线性无关 组。此外定义 m = a ，唯+ 。= 以+ ，- z ，_ ，k = 1 ，2 ，n - 1 ，其中 ( 踟，眺) 一i - i 。，( _ ，如) ，= 弋考茜习_ 则“，) 为日- 半正交并依旧是线性无关 组。 由于，( o 对于a 的指数d 常小于，因此并不需要限定矩阵z 满足z a 为非零 实的，即对于矩阵z 的条件可以放宽。 定义2 6 令s 为r v 中的子空间，而且令h 为n xn 的实矩阵。若对所有子 空间s 中的向量石0 都能满足( 工，x ) 0 ，则称日对于子空间s 是非零实的。 定理2 7 假设s 为尺空间中的n 维子空间，h 为n xn 矩阵且对于子空间s 是非零实的，则s 中仅有零向量垂直子空间h s ，并脚中唯有零向量垂直s 。 由定理2 7 得知矩阵z 若满足烈对于髟( 厂们，彳) 是非零实的，则对于任意 的刀n ，子空间z a k 。( ，们，彳) 中唯有零向量垂直子空间k 。( 厂们，彳) 。是否线性方 程组( 2 4 ) 能因为z a 对于k d ( ，j 们，彳) 为非零实的而保有唯一解，我们将介绍下面 两个定理： 定理2 8 令s 为尺空间中的n 维子空间，且日对于s 是非零实的，若 a ，见) 是子空间s 的一组基，则矩阵曰= ( 蝴，h p 2 ，地) ( p 1 p 2 ，p n ) 以及 c = ( p 。p 2 ，岛) ( 地，h p 2 ，地) 均是非奇异的。 定理2 9 令s 为尺空间中的n 维子空间，a 和z 均为n xn 的矩阵，其中 a 为非奇异矩阵，v 是尺中的一个向量。若z 满足别对于子空间s 是非零实 的，则超平面v + a s 中存在唯一向量垂直子空间z s 。 定理2 9 证明了若z ，4 对于畅( ，0 ，爿) 为非零实的，则线性方程组( 2 4 ) 有唯一 解，也就是在超平面x o + k 。( ，彳) 中仅存在一个向量满j 已g a l e r k i n 条件。而在 定理2 2 中说明x 是x + k ( ，彳) 中满足g a l e r k i n 条件的向量，再由上面这个定 理证明能满足g a l e r k i n 条件是唯一的，因此可以确信当刀= d - t - 1 次迭代时，所得 到的近似迭代解x ( ”必定会和解析解x 相同。 k r y l o v 子空间法是借l i t k r y l o v 子空间的性质来求取真解x ，k r y l o v 子空间法 和选用的矩阵有相当密切的关系，也就是在n d - i - 1 次迭代时，从超平面 x + k 。( r 们，a ) 求取满足g a l e r k i n 条件的迭代近似解工，并且用 k s m ( z ，【风，1 ，】) 来表示矩阵为z 的k r y l o v 子空间法。而以【饥，1 ，】表示方向 向量的准则，其中方向向量p 满足 其中日为n xn 的实矩阵，k 是非负整数或无限大，而v 是 o ，l ，) 一尺的函 数。较常见到的函数，有下列两种： ( 1 ) v ( o ) = ，叭，1 ，( 刀) = a p 川，n = 1 ，2 ，d ( 2 ) v ( n ) = 厂，n = 0 ，l ，2 ，d 将具有以上两种函数的向量准则分别表示为 饥，p 】， 饥，尺 。虽然 v ( n ) = a ”厂0 在理论上是可以做为方向向量准则，但在之前提到根据数值的观点 9 c j，厶 一 d 聆 ，b ，b q + = 尼 n 一聆，尼棚肛 = r ， k 仉 + = 1叫坳 町 p 0 ) 2 ，、【 并不合适，因此以后将不讨论方向向量为彳”，( o 的类型，并将在第四章中以实际 例子说明此类型的k r y l o v 子空间法并不合适。 2 2 停止条件 要决定何时停1 1 2 k s m ( z ， 1 - i , ，v ) 法，必须经由一些情况来判定停止的条件。 k s m ( z ， 饥， ) 法将会因为下列的理由而停止： 1 外部的错误； 2 迭代次数超过； 3 不够充分的财力或资源去继续程序运行； 4 收敛结果满意； 5 程序发生中断。 外部的错误是迭代法本身之外的问题，有可能是由于求解的问题无解。使 用者必须去发现这样的问题并且停止继续迭代下去。如果使用者不停止，迭代 法继续做下去后可能会停止，或是不收敛，甚至收敛到一个错误的解。迭代次 数的超过是由于k r y l o v 子空问法本身所生成的子空间的维度将不会超过系数 矩阵的维度，因此所设的迭代次数不需要超过系数矩阵的维度。而不够充分的 财力或资源去继续程序运行是没有足够的存储空间去执行；如g m r e s ( m ) 法， 若将m 取作无限大，此时所需要的存储空间将会越来越大，最后导致存储空间 不足而无法继续下去。而收敛结果的满意也就是收敛测试。测试的条件通常有 下列两种： ( 1 ) 残向量的范数够小 ( 2 ) 估计误差的范数够小 一般所使用的范数有1 范数，2 范数，以及无穷大一范数。而中断发生的情 况有下列三种： 情况一：若k s m ( z ， 以，川) 在刀d + 1 次迭代时，超平面x o + 民( 厂们，么) 中 不存在或是不只一个向量满足g a l e r k i n 条件。 情况二：若k s m ( z ，【以，1 ，】) 在n d + 1 次迭代时，方向向量p “k 。( ，m ，么) 。 情况三：若k s m ( z ， 巩，川) 在，z d + 1 次迭代时，无法从，( 胛) + 疋( r 们，爿) 中 1 0 找出方向向量p ”使得能与p ( n - o , p ”劲，p ”都是衅正交的向量。 下面给出三个实例来阐述上面三种情况： 例i 给定一4 x 4 的线性方程组a x = b ，其中彳= 用k s m ( i ， 风，p 】) 法求解，且选取x o = p ( o ) = ，( o ) = ol 1o 0o o0 0o 0 0 0l 10 ，b = ，则，( o ) = b 一彳石( o ) = 。若 ，且令 。因为( p ，z r o ) = ( p ，0 ) = 2 且( p 们，z 4 r 。) = o ，所以线性 方程组( p 0 ，z a r ) 口5 ( p ，z r ) 无解。即无法从超平面一+ k 。( ，彳) 中求 出满足g a l e r k i n 条件的向量。 例2 给定一2 2 的线性方程组彳石= 6 ，其中爿= ( ? 三) ，6 = ( 芋) 。假设初 始值x ( o ) = ( ：) 。：g i 罚k s m ( 彳【4 ，尺 ) 法求解，也就是选取矩阵z = 彳，而且要求 方向向量p 佃取自厂( o ) + k ( 厂( o ) ，彳) 。令初始方向向量p ( o ) = ，( o ) = ( 习，并假设 x ”= x + a p o 是x o + k i ( ，0 1 ，彳) 中满足g a l e r k i n条件的向量。 则 ( p 0 ，z 厂i l ) = ( p ，z ，一a z a p ) = 一a p 0 ，z a p ) = 0 ，显然( p ，z a p ) o ， 所以得到口= 0 。因此x 1 = x 叭，且厂1 = ，o k l ( ，彳) 。所以不论p 1 无论怎么 从r o + k l ( r 们，彳) 中选取，一定会落在k l ( ，们，么) 。 例3 所要解的线性方程组和例2 相同，但却利用k s m ( z ， 4 ，p 】) 法求解并取 起始值x ( o ) = ( ： 并且令p ( o ) = ，( 0 ) = ( 三 。贝u4 p ( o ) = ( ? ) 、c 4 p ( o ) ，4 p ( o ，= 且 ( p ( 们，却o ) = 0 。由于k s m ( z ， 4 ，尸】) 法要从彳p ( o ) + k ，( ，( 0 1 ，彳) 中求出一个能和 p 。彳一半正交的向量。但是线性方程组( p 们，却。) 口= ( 却们，却。) 显然无解。 所以在却+ k l ( ，们，彳) 中不存在能与p a 一半正交的向量。 在前面提到当矩阵z 满足黝是非零实的，则对任何的nsd + 1 ，则超平面 x o + k ( ，们，彳) 中存在唯一向量满足g a l c r k i n 条件，也就是当z能满足z 4 对 + 。( ，彳) 是非零实的，则情况一将不会发生。 若对于n = 0 ，1 ，d ，将方向向量p “取自k n + 。( 厂们，a ) - k ( 厂们，彳) ，则 仞们，p “1 ) 是秘+ 。( ，们，彳) 的一组基，因此可以避免情况二发生。可是在 k + ，( ，们，彳) 的基确定前，很难从k 川( r 们，a ) - 疋( ，们，彳) 中求出方向向量p 。因 此f 面将介绍构建叫的方式。 定理2 1 3 令p o = ，们。若将p “取自却川+ k ( 厂们，么) 1 n d ，则 p 们，p n ，p ”是e + i ( 厂们，彳) 的一组基。 证明见 “】。 定理2 1 3j , 正明若强制p ”取自却”1 + k 。( ，们，彳) ，则情况二将不会发生。 c h u a n g 1 5 】提到若将p ”取自厂”+ k ( ，彳) ， 且x 肿1 落在 x o + e + l ( ，们，彳) 一蟛( ，们，彳) ，因此就可以确保 ，川 落在 k 棚( ，们，爿) 一k 川( 厂们，么) 。因此p ”取自，肘”+ k + l ( ，| 们，彳) 将不会发生情况二的 中断。若x 川工o + k 。( 厂，彳) ，由于x + k 。( 厂叭，彳) 中满f f z g a l e r k i n 条件的向量 唯 一 存在，即x ”) = 工( ” 。 因此 册胛 p 们，p ，p ”) = s p a n p t o ) , p n ，p ，p 川) ，所以情况二就发生了。 情况 二 是由 于乃o 垂直 k 。( 厂，a ) 而导致 ( p ，p ，p ”) 7 z r o = ( p 们，p ，p ”) 7 别( p 0 1 ，p ，p ”) 口”呈齐次线性 方程组，因此可以利用矩阵z 来避免。 定理 2 1 4 令z 对于k d + ，( ，0 ，么) 是非零实的，对于任意的 m = 0 ，l ，刀d ， z r ”垂直k m ( r 们，彳) 且p ”，埘+ k ( 厂们，彳) ， 则 ( p ，p ，p ”) z r 。o ，0 k n d 证明见 1 1 】 1 2 定理2 1 5 若日为非零实矩阵，则对于任意的矩阵z 以及任意的k ， k s m ( z ，【巩，尺 ) 、k s m ( z ， 风，尸】) 两种准则都不会发生情况三的现象。 证明见 1 5 】 利用k s m ( z ， 吼，1 ，】) 来求解线性方程组血= b 可能发生中断的三种情况都 可借由矩阵z 和的选取来避免。要求矩阵z 满足别是非零实的来避免情况一 的发生，而要求日为非零实的来避免情况三的发生。至于情况二的现象若方向 向量条件为尸型，则不会发生中断，但是r 型的方向向量条件则必须要求矩阵z 对于子空间是非零实的才不会发生情况二的中断。 第三章k r y l o v 子空间法 由于k r y l o v 子空间法在近几年有相当多种的类型，在本章将介绍最近几年 比较有名的k r y l o v 子空间法，其余的在表三中将介绍其为k r y l o v 子空间法的哪 种类型。 c gm e t h o d ： c g 法是由h e s t e n e s 与s t i e f e l 【2 】在1 9 5 2 所提出的，这是k s m ( i ， ，r ) 的 k r y l o v 子空间法。c g 法用来解对称且正定方程组十分有效，但是若拿来解非对 称或是非正定的线性方程组则会发生中断。它是借由迭代的方式产生一序列的方 向向量用来更新迭代解以及残向量，虽然产生的序列会越来越大，但是却只需要 存储少数的向量。当系数矩阵a 相当大而且稀疏，此时c g 法几乎就是高斯消去法。 c g 法理论上虽然保证最多n 步能解出线性方程组a x = b 的解，但是若系数矩阵是 病态的，此时累进误差会让c g 法在n 步后无法求得充分精确的近似解。 给一个初始值z ( 叭，计算，o ) = 6 一a x o 当i = l ，2 ， 如果f = l p ( 1 ) = ，( o ) 否则 屈= r ( i - i ) , r ( i - 1 ) ) 伦卜孙，) p ”= 厂卜1 + 层一l p 卜” 结束，如果 a ，= r ( i - i ) , r ( i - i ) ) ( a p “) x ( ) = x 卜1 ) + a g p ( ) t h ec o n j u g a t eg r a d i e n tm e t h o d cg ne andcgnr ： 利用c g 法仅能求解对称且正定的方程组，而非对称的方程组就无法求解。 1 4 h e s t e n e s 与s t i e f e l 2 提到求解线性方程组a x = 6 和解正规方程有相同解，且彳r 么 是一个对称且半正定的矩阵，因此新的线性方程组就可以利用c g 法求解。c g n e 求解的正规方程是( 州r ) y = b ，其中z = a r y 。c g n r 求解的正规方程是 ( a7 a ) x ：占，其中占= a r b 。这两种方法对于非正定的阵有可能达到极快速度的 收敛；例如系数矩阵a 若为么j 下矩阵，此时一般的k r y l o v 子空间法收敛缓慢，但 是c g n e 跟c g n r 就能在一步求解出线性方程组之解。但是一般由偏微分方程离 散化所得到的线性方程组的系数矩阵彳不可能如此完美，所以a x e l s s o n 1 3 在 1 9 8 0 年提到如此的解法常常会因为a7 a 或州7 的条件数大于a 的条件数而导致 收敛速度比求原来的方程组要来的缓慢。因此这两种方法的实用性并不是很强。 g e n e r a l i z e dm i n i m a lr e s i d u a l ( m ) ( gmres ( m ) ) g m r e s ( m ) 是由s a a d 与s c h u l t z 1 6 所提出的，这是k s m ( a t , l ，尸】) 的 k r y l o v 子空间法。g m r e s ( m ) 是m i n r e s ( 仅能求解对称方程组) 的推广，且可 以求解非对称的线性方程组。这种方法是在m 步迭代中将残向量变成最小的。 g m r e s ( m ) 的主要缺点在于每多一次迭代时，所需要的计算次数以及存储空间 都成线性的增加，如果不是很幸运的得到一个快速的收敛，此时去计算线性方程 组将会变得成本很高。因此最常用来克服上述缺点的就是利用一种叫做重置的方 法去控制存储空间；也就是当选取好迭代的次数m 后，将之前所存储的资料清 除并且将最后的迭代解用来当作接下来新的m 步迭代的初始值。将这重置的过 程一直持续下去直到收敛到满意的值。如果不用重置的方法，虽然收敛的结果最 后不超过n 步。但是n 很大时，这种方法就没什么实用价值，因为存储空间以及 运算的要求将会非常的大。因此g m r e s ( m ) 要成功的运用的重点在于何时去 重置，也就是m 的选取。因此如何去选取m 是最困难的，如果m 太小，则 g m r e s ( m ) 将会收敛的很慢或是根本不收敛；如果m 选取的太大，则需要额外 的运算以及更多的存储空间。不幸的是并没有什么法则告诉我们如何来选取这样 的m ，一切都靠经验来决定这样的m 。由于矩阵z 取的是a ，显然z a 是非零实 的，h = i 对于k 叶。( r ，a ) 是非零实的，且方向向量的条件是p 型，因此不会 发生中断。 给初始值x ( m ，计算，= b 一血o 当_ ，= l ，2 ， v n = r l l r l l ：；s = l l r i ：e l 当i = l ，2 ，m w = 彳1 ，7 ) 当k = 1 ，2 ，i 以广( w y 耻) ；w = w - h k f v 似 v i + 。= w 1 1w i i ： 当k = 2 ，i u = _ l - f 纸- 1 f = c k l u + s k l 魂f 喽f = 一s k - l u + q i f 万= ；c i = h t , i 6 ；s f = h i 一1 2 6 ，；f = q ，+ j f h i “f b i + i = 一s f 岛；b i = c i b t p = 忙l l l 2 ：i i b - 血州叫i ： 如果p 是足够小，则( n r = f ；转到s o l ( x 歹) ) n r = m ；= 占雌h l l r ，刀， s o l ： 当k = n ，一1 ，1 x 一= y + + ， 1 6 七嚷 、， 乃纯 雌小 一 饥 ，- = 七y 如果；是一个够精确的近似值，则存在 否贝u 工( o ) ：； thegene r al iz edminimumresiduaimeth0d i g c g 法是f 1 3 y o u n g 与j e a 【1l 】【1 2 所发明出来的，其中包括t - - 种k r y l o v 子 空间法：o r t h o d i r 、o r t h o m l n 以及o r t h o r e 。这三种方法分别属于 k s m ( z ， 巩，尸】) 、k s m ( z ， 巩，尺】) 以及k s m ( z ，【z 如，尺】) 子空间法。矩阵z 通 常选为a 。倘若矩阵不选为a ，则我们很难选取其他的矩阵使得z 4 为对称且正 定。因此o r t h o m i n以及o r t h o r e 则会发生中断。这三种方法由于所选取 的方向向量需要之前的方向向量做运算，因此当迭代次数增加时，就会发生和 g m r e s 法相同的问题。因此i g c g 法利用了一种叫截断的方法，也就是将某些 之前的方向向量舍弃，因此i g c g 法的三种k r y l o v 子空间法就变成 k s m ( z ， 弛，p ) 、k s m ( z ， 戤，尺 ) 以及k s m ( z ， 戤，尺】) 的子空间法。 p ( o ) = ，( o ) u ”+ 1 = 材”+ 以p 4 以= 黼 p ”= 印川+ 尾 l p 川+ + 尾o p o ( 别z p ，p ) + i - i 成z a pj ) , p ) 玑，=一_z云虿歹占专j万了一，f=。，1，l一1 p ( o ) = ，( o ) 1 4 ”+ 1 = “”+ l p ” 丸= p ”= ，”+ 尾n i p ”一1 + + 尾。o p o ( 别2 ) ，p ) + i - i 尾，z a p ) p ) p n , i2 - - 弋荔易丁一。0 1 川。 o r t h o r e so ”( 肿1 ) = r ( ”) + 口n + 1 h ( “) + + 口玎+ 1 o ”( 0 ) 五”= ( 尾+ 1 o + 成+ i i + + 尾+ i ，。) - 1 吒= 乃尾“f ( 别2 ，p ) + i - ! 川，( z ，) 成+ - j 2 一- j ；毒i _ 了万歹一，f = 。1 甩一1 i g c gm e t h o d s b i c o n ju g a t eg r a d i e n t ( b i c g ) ： b i c g 法是i 士i f l e t c h e r 【1 7 所提出的，这是k s m ( i ， 1 。，r ) 的k r y l o v 子空间法。 倘若系数矩阵彳是一个对称正定阵，贝j j b i c g 法就和c g 法一致。但是我们所花的 计算成本将是c g 法的两倍。但是若系数矩阵彳是个非对称的矩阵，由于要借由短 的过程来构建互相垂直的方向向量是不可能的( 证明见f a b e r 与m a n t e u f f e l 2 ) 。因 此b i c g 采用另一种方法的逼近，也就是利用两组互相垂直的方向向量来取代一 般用一组相互垂直的方向向量，因此不再如同之前的方法一样将残向量最小化。 由于b i c g 不是将残向量最小化，所以他的收敛情况相当不规则，而且这种方法 有可能会发生中断。发生中断的情况有可能是仁( i - 1 ) 1 o 或是( g ( f ) ，；m ) o 。 、 ，、， 而要避免发生中断可以在快发生中段时借由重置的方法来避免，或是转换成一种 较健全的方法( 如g m r e s ( m ) )