内容发布更新时间 : 2024/12/23 17:28:15星期一 下面是文章的全部内容请认真阅读。
大型稀疏矩阵数值解法综述
前言
线性方程组的求解方法基本上可以分为两种：直接法和迭代法。直接法是以Gauss消去法为基础的方法，实际上这种方法是初等代数中解多元一次方程组的直接推广。其基础思路是：对原方程组经过一定的运算处理后，逐个消去部分变量，最后得到一个与原方程等价的、便于逐步求解的方程组，以解除各个变量的值。如果不考虑计算式可能产生的舍入误差，则可以认为直接法是一种精确的方法，它可以一次性求得原线性方程组的解。与直接法相比，迭代法是一种近似的解法，它的基本思想是：先估计一组变量的数值，作为原方程组的第一次近似，也称作迭代初值，然后利用某种迭代处理，逐次修改这组数值，得到解的第二次、第三次以至第k次近似值。在方程组满足一定条件的前提下，这些结果可以逐次逼近原方程组的真实解。进行有限次重复计算后，如果所得到的近似解满足所规定的误差范围，则认为最后一次迭代的近似值就是原方程的解。
随着石油行业的发展，对油藏的认识也不断地深入，通过应用大型计算机，油藏数值模拟的规模也越来越大，另外随着一些新的方法如：、全隐式、局部加密网格方法的使用，以及混合驱、热力驱、化学驱等油藏模型的出现，使方程组矩阵的性质更加复杂多变，甚至带有严重的病态性质。开发能够快速求解各种复杂的矩阵线性方程组的新方法就成为油藏数值模拟进一步发展的一个重要方向。
通过对IMPES、全隐式的油藏数值模拟常用方法的研究发现，油藏数值模拟中所研究的矩阵只是在少数元素上为非零值，大部分元素值为零值，数学上称之为稀疏矩阵。当研究对象为二维油藏数值模拟时，稀疏方程组的系数矩阵为对称格式，当研究对象为三维油藏数值模拟时，稀疏方程组的系数矩阵则变为非对称格式，直接法在求解这种稀疏矩阵时，会使原矩阵新添加大量的非零元素，增加计算量和存储量，对于求解大规模稀疏矩阵，速度会变得非常缓慢，甚至无法求解，因此直接法已经逐渐被迭代法所取代，不再作为求解大型油藏数值模拟线性方程组的方法。
迭代法包括古典迭代法和Krylov子空间迭代法。古典迭代法包括Jacobi.、Gauss Seidel、SOR、ssor等方法，因其收敛速度很慢，目前也已很少用于求解大型稀疏线性方程组，而是用于与一些效率更高的方法相结合使用。目前，Krylov子空间方法是求解大型稀疏矩阵线性方程组最流行和最有效的方法之一，也是当
前研究的热点，其主要思想是为各迭代步递归地构造残差向量，即第n步的残差向量rn通过系数矩阵A的某个多项式与第一个残差向量r1相乘得到：
rn?Pn?1(A)r1
通常，迭代多项式的选取应使用所构造的残差向量在某种内积意义下相互正交，从而保证某种极小性（极小残差性），达到快速收敛的目的。Krylov子空间方法具有存储量少，计算量少且易于并行等优点，非常适合于求解大型稀疏线性方程组，结合预条件的Krylov子空间迭代法是目前求解大型稀疏线性方程组的最主要方法，
1 Krylov子空间概念
给定初值x0，求解稀疏线性方程组：
Ax?b
其中A?Rn?n,x,b?R
设Km为m维子空间，一般投影方法是从m维仿射子空间x0?Km中寻找近 似解xm使残差向量满足Petrov—Galerkin条件
rm?b?Axm?Lm
其中Lm为另一个m维子空间，x0为迭代初值。如果Km是Krylov子空间， 则上述投影方法就称作Krylov子空间方法。Krylov子空间Km(A,r0)定义为
Km(A,r0)?span{r0,Ar0,A2r0.?,Am?1r0}
理想的Krylov子空间方法求解稀疏线性方程组Ax?b的过程具有以下特征： 1. 极小剩余性或极小误差性（保证收敛速度快）；
2. 每一步迭代的计算量少，存储量小，以保证计算高效性。
2 Krylov子空间方法分类
选择不同的Km和Lm就得到了不同的Krylov子空间算法。
2.1 基于正交投影方法
取Lm?Km?Km(A,r)，这一类Krylov子空间方法称为正交投影方法。共轭梯度（CG）方法是其中最重要的方法，此时要求矩阵A是对称正定矩阵。CG该方法要求x?xm在子空间x0?Km中的能量范数到达极小，且具有短的迭代计算
公式。此后发展得到的与各种预条件相结合的方法成为求解大型对称正定稀疏线性方程组最主要的方法。 2.1.1 共轭梯度法（CG） 具体流程为：
1.任选初值x0，计算初值残差b?Ax0?r0和内积（r0,r0），并置s0?r0,置j=0。 2.计算参数aj?(rj,rj)/(Ajs,js更新向量xj?1?xj?ajsj与残向量
jrj?1?rj?jaA，若srj满足精度要求，则停机。
3.计算?j?1?(rj?1,rj?!)/(rj,rj),sj?1?rj?1??j?1sj，置j=j+1,转（2）
精度测试
维数为9*9矩阵：
通过直接法所求的的矩阵真实解和通过CG法迭代求的的近似解误差值的数量级在10-12。
维数为36*36矩阵：
通过直接法所求的的矩阵真实解和通过CG法迭代求的的近似解误差值的数量级在10-8—10-12。
维数为64*64矩阵：
通过直接法所求的的矩阵真实解和通过CG法迭代求的的近似解误差值的数量级在10-8—10-9。
维数为169*169矩阵：
通过直接法所求的的矩阵真实解和通过CG法迭代求的的近似解误差值的数量级在10-9—10-10。 2.1.2 预条件共轭梯度法
共轭梯度法在解决大型线性方程组本身是有不足之处的，当系数矩阵A接近单位阵，共轭梯度法收敛较快，而当系数矩阵A的条件数大于102时，共轭梯度法就非常缓慢了，因此需要通过对矩阵A的处理，是系数矩阵的条件数下降，从而减少迭代次数和计算量。
预处理方法是目前解决系数矩阵条件数问题最有效的方法之一，其原理是选取预处理矩阵M，此矩阵具有以下特征： M为对称正定的，而且具有与A相似的稀疏性，不但易于求逆，还能改善A矩阵的病态特性，,使得M矩阵的特征值比A矩阵集中。M的选择尤为重要，因为稀疏矩阵A可进行分裂式近似分解，
TA?LLT?R,其中L矩阵是下三角矩阵，R矩阵称为剩余矩阵，则可将LL看做M,
用M矩阵代替A矩阵进行迭代。 预处理矩阵得到流程：
1 先处理预处理矩阵第一列元素
m1,1?a1,1
mi,1?ai,1/m1,1 i=2 3….n
2 循环嵌套得到预处理矩阵其它元素的值
k-2.3……n
如果
k?1t?1
则
ak,k??mk,t2\\
mk,k?ak,k??mk,tt?1k?12
否则
mk,k?ak,k i=k+1…..n
k?11mi,k?(ai,k??mi,tmk,t)mk,kt?1
得到了预处理矩阵M 就可以进行预处理共轭梯度法了
预处理共扼梯度法(ICCG) 具体迭代算法如下：
1 输入系数矩阵A，预处理矩阵M,右端向量b和初值x0，于是得到一系列迭代初值
r0?b?Ax0z0?M?1r0p1?z0置k=1;
?0??r0,z0?,2 开始迭代 K=2.3….n
w?Apk,ak??k?1/?pk,w?,xk?xk?1?akpk,
rk?rk?1?akwzk?Mrk,?1
?k??rk,zk?,?k??k/?k?1,pk?1?zk??kpk3 如果达到精度要求，则输出x值，结束；否则k=k+1,转2； 精度测试
维数为9*9矩阵：
通过直接法所求的的矩阵真实解和通过CG法迭代求的的近似解误差值的数量级在10-12。
维数为36*36矩阵：
通过直接法所求的的矩阵真实解和通过CG法迭代求的的近似解误差值的数量级在10-12。
维数为64*64矩阵：
通过直接法所求的的矩阵真实解和通过CG法迭代求的的近似解误差值的数量级在10-12。
维数为169*169矩阵：
通过直接法所求的的矩阵真实解和通过CG法迭代求的的近似解误差值的数量级在10-12。
2.2 正交化法
取Lm?AKm(A,r0)，则这类Krylov子空间方法称为正交化方法。GMRES
（Generalized Minimal RESidual algorithm）即广义极小余量算法，是正交化方法的典型代表，也是目前求解大型稀疏非对称线性方程组的主要方法之一。
GMRES算法是以Galerkin原理为基础的算法，具体地，设所求方程为：
Ax?b，b?Rn为一给定向量，其中A为一大型不规则稀疏矩阵，记Km和Lm是Rnmn中的两个m维Krylov子空间，分别由{vi}im和所张成。取为任一向{w}x?R?1ii?1o量，令x?xo?z，则方程Ax?b等价于：