油藏数值模拟中大型稀疏矩阵算法综述教程文件.doc

油藏数值模拟中大型稀疏矩阵算法综述

精品文档
大型稀疏矩阵数值解法综述
前言
线性方程组的求解方法基本上可以分为两种：直接法和迭代法。直接法是以Gauss消去法为基础的方法，实际上这种方法是初等代数中解多元一次方程组的直接推广。其基础思路是：对原方程组经过一定的运算处理后，逐个消去部分变量，最后得到一个与原方程等价的、便于逐步求解的方程组，以解除各个变量的值。如果不考虑计算式可能产生的舍入误差，则可以认为直接法是一种精确的方法，它可以一次性求得原线性方程组的解。与直接法相比，迭代法是一种近似的解法，它的基本思想是：先估计一组变量的数值，作为原方程组的第一次近似，也称作迭代初值，然后利用某种迭代处理，逐次修改这组数值，得到解的第二次、第三次以至第k次近似值。在方程组满足一定条件的前提下，这些结果可以逐次逼近原方程组的真实解。进行有限次重复计算后，如果所得到的近似解满足所规定的误差范围，则认为最后一次迭代的近似值就是原方程的解。
随着石油行业的发展，对油藏的认识也不断地深入，通过应用大型计算机，油藏数值模拟的规模也越来越大，另外随着一些新的方法如：、全隐式、局部加密网格方法的使用，以及混合驱、热力驱、化学驱等油藏模型的出现，使方程组矩阵的性质更加复杂多变，甚至带有严重的病态性质。开发能够快速求解各种复杂的矩阵线性方程组的新方法就成为油藏数值模拟进一步发展的一个重要方向。
通过对IMPES、全隐式的油藏数值模拟常用方法的研究发现，油藏数值模拟中所研究的矩阵只是在少数元素上为非零值，大部分元素值为零值，数学上称之为稀疏矩阵。当研究对象为二维油藏数值模拟时，稀疏方程组的系数矩阵为对称格式，当研究对象为三维油藏数值模拟时，稀疏方程组的系数矩阵则变为非对称格式，直接法在求解这种稀疏矩阵时，会使原矩阵新添加大量的非零元素，增加计算量和存储量，对于求解大规模稀疏矩阵，速度会变得非常缓慢，甚至无法求解，因此直接法已经逐渐被迭代法所取代，不再作为求解大型油藏数值模拟线性方程组的方法。
迭代法包括古典迭代法和Krylov子空间迭代法。古典迭代法包括Jacobi.、Gauss Seidel、SOR、ssor等方法，因其收敛速度很慢，目前也已很少用于求解大型稀疏线性方程组，而是用于与一些效率更高的方法相结合使用。目前，Krylov子空间方法是求解大型稀疏矩阵线性方程组最流行和最有效的方法之一，也是当前研究的热点，其主要思想是为各迭代步递归地构造残差向量，即第n步的残差向量通过系数矩阵A的某个多项式与第一个残差向量相乘得到：

通常，迭代多项式的选取应使用所构造的残差向量在某种内积意义下相互正交，从而保证某种极小性（极小残差性），达到快速收敛的目的。Krylov子空间方法具有存储量少，计算量少且易于并行等优点，非常适合于求解大型稀疏线性方程组，结合预条件的Krylov子空间迭代法是目前求解大型稀疏线性方程组的最主要方法，
1 Krylov子空间概念
给定初值，求解稀疏线性方程组：

其中
设为m维子空间，一般投影方法是从m维仿射子空间中寻找近
似解使残差向量满足Petrov—Galerkin条件

其中为另一个m维子空间，为迭代初值。如果是Krylov子空间，
则上述投影方法就称作Krylov子空间方法。Krylov子空间定义为

理想的Krylov子空间方法求解稀疏线性方程组的过程具有以下特征：
1. 极小剩余性或极小误差性（保证收敛速度快）；
2. 每一步迭代的计算量少，存储量小，以保证计算高效性。
2 Krylov子空间方法分类
选择不同的和就得到了不同的Krylov子空间算法。
2.1 基于正交投影方法
取，这一类Krylov子空间方法称为正交投影方法。共轭梯度（CG）方法是其中最重要的方法，此时要求矩阵A是对称正定矩阵。CG该方法要求在子空间中的能量范数到达极小，且具有短的迭代计算公式。此后发展得到的与各种预条件相结合的方法成为求解大型对称正定稀疏线性方程组最主要的方法。
2.1.1 共轭梯度法（CG）
具体流程为：
1.任选初值，计算初值残差和内积（），并置,置j=0。
2.计算参数更新向量与残向量，若满足精度要求，则停机。
3.计算，置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的处理，是系数矩阵的条件数下降，从而减少迭代次数和计算量。
预处理方法是目前解决系数矩阵条件数问题最有效的方法之一，其原理是选取预处理矩阵，此矩阵具有以下特征： 为对称正定的，而且具有与A相似的稀疏性，不但易于求逆，还能改善矩阵的病态特性，,使得矩阵的特征值比矩阵集中。的选择尤为重要，因为稀疏矩阵可进行分裂式近似分解，,其中L矩阵是下三角矩阵，R矩阵称为剩余矩阵，则可将看做,用矩阵代替矩阵进行迭代。
预处理矩阵得到流程：
1 先处理预处理矩阵第一列元素

i=2 3….n
2 循环嵌套得到预处理矩阵其它元素的值
k-2.3……n
如果
\
则

否则

i=k+1…..n

得到了预处理矩阵M 就可以进行预处理共轭梯度法了
预处理共扼梯度法(ICCG) 具体迭代算法如下：
1 输入系数矩阵A，预处理矩阵M,右端向量b和初值，于是得到一系列迭代初值
置k=1;
2 开始迭代 K=2.3….n

3 如果达到精度要求，则输出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 正交化法
取，则这类Krylov子空间方法称为正交化方法。GMRES（Generalized Minimal RESidual algorithm）即广义极小余量算法，是正交化方法的典型代表，也是目前求解大型稀疏非对称线性方程组的主要方法之一。
GMRES算法是以Galerkin原理为基础的算法，具体地，设所求方程为：，其中为一大型不规则稀疏矩阵，为一给定向量，记和是中的两个维Krylov子空间，分别由和所张成。取为任一向量，令，则方程等价于：

其中。
求解的Galerkin原理可叙述为：在子空间中找的近似解，使得残余向量和中的所有向量正交。即对，有：

令和，其中和分别为和的基底。因此，可以把表示成，。这样上式就可以写成：

假设为非奇异矩阵，则不难得到近似解：

要求出的值必须先选择合适的和以及它们的基底，在GMRES算法中，我们取，而，简记。对于这组特定的和的取法，实际上求等价于在中极小化残余向量的2-范数。
利用Arnoldi过程，及，其中是一个上Hessenberg矩阵，，可以找到中的一组标准正交基，不难得：

=
=
=
其中，。由于是正交的，即，所以：

于是，在中极小化就等价于在中极小化，具体地，GMRES算法步骤如下：
1）选择初始近似值，计算，记，则；
2）定义阶矩阵，并设；
3）计算和的值：
4）求解最小二乘问题，得到的解为；
5）计算。
从理论上讲，如果线性无关，当时GMRES算法应当给出方程的准确解。但是，不难看出，当很大时，计算中需要保存所有的，对于大型问题（）时，这将引起存储空间过多的需求，因此是不现实的。克服这个困难的方法是采用如下形式的循环GMRES（m）算法，具体如下：
1）选择初始近似值，计算，记，则；
2）定义阶矩阵，并设；
3）完成Arnoldi过程（不再详述），得到和；
4）极小化得到；
5）计算；
6）计算，如果满意，则停止；
7），，转向2）。
在GMRES的具体计算过程中，要根据具体问题选择合适的m值，因为m值取得过大会增加存储的需求，然而太小又会减缓收敛的速度。
精度测试
维数为9*9矩阵：
通过直接法所求的的矩阵真实解和通过CG法迭代求的的近似解误差值的数量级在10-5—10-6。
维数为36*36矩阵：
通过直接法所求的的矩阵真实解和通过CG法迭代求的的近似解误差值的数量级在10-5。
维数为64*64矩阵：
通过直接法所求的的矩阵真实解和通过CG法迭代求的的近似解误差值的数量级在10-4—10-5。。
维数为169*169矩阵：
通过直接法所求的的矩阵真实解和通过CG法迭代求的的近似解误差值的数量级在10-4—10-5。
2.3 双正交化方法
在双正交化过程中取
,,

显然，A是非对称矩阵时，这一类方法就是双正交化方法,Lanczos提出的双共轭梯度法（BICG法）是其中最基础的方法。
双共轭梯度法BICG具体流程为：
1.计算选取使得置方向向量，并置j=0.
2.计算参数与向量，
如果满足精度，停机。
3.计算向量，置j=j+1.转2
BICG法在计算过程总中要使用到矩阵的转置，并且由于振荡剧烈收敛不好，为避免这些缺陷，许多研究者在BICG法的基础上发展了一系列收敛好的方法，如共轭梯度平方法（CGS法），广义共轭梯度平方法（GCGS法）、共轭梯度稳定性方法（BICGSTAB法）等。
BICGSTAB算法的主要思想是基于双边Lanczos,是基于残差正交子空间的迭代方法。对于n阶线性方程组，假设初始近似解为，第k次近似解为，相应的第k次的残差。是与不正交的非零向量。是分别与有关的k阶Krylov子空间，并且

BICGSTAB算法在Krylov子空间，选择序列，通过适当选择参数，使得按下式产生序列的满足且，同时计算的残差总是取最小的范数。

这种基于Krylov残差子空间的迭代方法收敛速度快，精度高，而且稳定性好。BICGSTAB算法常常和预处理技术相结合处理方程求解问题,如果预处理矩阵M是单位阵，则流程为一般的BICGSTAB方法，其流程为
1 给定系数矩阵A，向量b，初始值，最大迭代次数，最大容许误差以及预处理矩阵M，计算，并令k=1,
2 如果且，转3，否则终止，输出
3 ，如果或者，算法失败，终止，否则转4
4 如果K=1 ，则，否则计算

5 求解方程，计算

6 ，如果，，否则停止循环，输出
7 求解

转2

精度测试
维数为9*9矩阵：
通过直接法所求的的矩阵真实解和通过CG法迭代求的的近似解误差值的数量级在10-5—10-6。
维数为36*36矩阵：
通过直接法所求的的矩阵真实解和通过CG法迭代求的的近似解误差值的数量级在10-5。
维数为64*64矩阵：
通过直接法所求的的矩阵真实解和通过CG法迭代求的的近似解误差值的数量级在10-4—10-5。。
维数为169*169矩阵：
通过直接法所求的的矩阵真实解和通过CG法迭代求的的近似解误差值的数量级在10-4—10-5。

收集于网络，如有侵权请联系管理员删除