油藏數(shù)值模擬中大型稀疏矩陣算法綜述

1、 大型稀疏矩陣數(shù)值解法綜述前言線性方程組的求解方法基本上可以分為兩種：直接法和迭代法。直接法是以Gauss消去法為基礎(chǔ)的方法，實際上這種方法是初等代數(shù)中解多元一次方程組的直接推廣。其基礎(chǔ)思路是：對原方程組經(jīng)過一定的運(yùn)算處理后，逐個消去部分變量，最后得到一個與原方程等價的、便于逐步求解的方程組，以解除各個變量的值。如果不考慮計算式可能產(chǎn)生的舍入誤差，則可以認(rèn)為直接法是一種精確的方法，它可以一次性求得原線性方程組的解。與直接法相比，迭代法是一種近似的解法，它的基本思想是：先估計一組變量的數(shù)值，作為原方程組的第一次近似，也稱作迭代初值，然后利用某種迭代處理，逐次修改這組數(shù)值，得到解的第二次、第三次以

2、至第k次近似值。在方程組滿足一定條件的前提下，這些結(jié)果可以逐次逼近原方程組的真實解。進(jìn)行有限次重復(fù)計算后，如果所得到的近似解滿足所規(guī)定的誤差范圍，則認(rèn)為最后一次迭代的近似值就是原方程的解。隨著石油行業(yè)的發(fā)展，對油藏的認(rèn)識也不斷地深入，通過應(yīng)用大型計算機(jī)，油藏數(shù)值模擬的規(guī)模也越來越大，另外隨著一些新的方法如：、全隱式、局部加密網(wǎng)格方法的使用，以及混合驅(qū)、熱力驅(qū)、化學(xué)驅(qū)等油藏模型的出現(xiàn)，使方程組矩陣的性質(zhì)更加復(fù)雜多變，甚至帶有嚴(yán)重的病態(tài)性質(zhì)。開發(fā)能夠快速求解各種復(fù)雜的矩陣線性方程組的新方法就成為油藏數(shù)值模擬進(jìn)一步發(fā)展的一個重要方向。通過對IMPES、全隱式的油藏數(shù)值模擬常用方法的研究發(fā)現(xiàn)，油藏數(shù)值

3、模擬中所研究的矩陣只是在少數(shù)元素上為非零值，大部分元素值為零值，數(shù)學(xué)上稱之為稀疏矩陣。當(dāng)研究對象為二維油藏數(shù)值模擬時，稀疏方程組的系數(shù)矩陣為對稱格式，當(dāng)研究對象為三維油藏數(shù)值模擬時，稀疏方程組的系數(shù)矩陣則變?yōu)榉菍ΨQ格式，直接法在求解這種稀疏矩陣時，會使原矩陣新添加大量的非零元素，增加計算量和存儲量，對于求解大規(guī)模稀疏矩陣，速度會變得非常緩慢，甚至無法求解，因此直接法已經(jīng)逐漸被迭代法所取代，不再作為求解大型油藏數(shù)值模擬線性方程組的方法。迭代法包括古典迭代法和Krylov子空間迭代法。古典迭代法包括Jacobi.、Gauss Seidel、SOR、ssor等方法，因其收斂速度很慢，目前也已很少用于

4、求解大型稀疏線性方程組，而是用于與一些效率更高的方法相結(jié)合使用。目前，Krylov子空間方法是求解大型稀疏矩陣線性方程組最流行和最有效的方法之一，也是當(dāng)前研究的熱點，其主要思想是為各迭代步遞歸地構(gòu)造殘差向量，即第n步的殘差向量通過系數(shù)矩陣A的某個多項式與第一個殘差向量相乘得到：通常，迭代多項式的選取應(yīng)使用所構(gòu)造的殘差向量在某種內(nèi)積意義下相互正交，從而保證某種極小性（極小殘差性），達(dá)到快速收斂的目的。Krylov子空間方法具有存儲量少，計算量少且易于并行等優(yōu)點，非常適合于求解大型稀疏線性方程組，結(jié)合預(yù)條件的Krylov子空間迭代法是目前求解大型稀疏線性方程組的最主要方法，1 Krylov子空間概

5、念給定初值，求解稀疏線性方程組：其中設(shè)為m維子空間，一般投影方法是從m維仿射子空間中尋找近似解使殘差向量滿足PetrovGalerkin條件其中為另一個m維子空間，為迭代初值。如果是Krylov子空間，則上述投影方法就稱作Krylov子空間方法。Krylov子空間定義為理想的Krylov子空間方法求解稀疏線性方程組的過程具有以下特征：1. 極小剩余性或極小誤差性（保證收斂速度快）；2. 每一步迭代的計算量少，存儲量小，以保證計算高效性。2 Krylov子空間方法分類選擇不同的和就得到了不同的Krylov子空間算法。2.1 基于正交投影方法取，這一類Krylov子空間方法稱為正交投影方法。共軛梯

6、度（CG）方法是其中最重要的方法，此時要求矩陣A是對稱正定矩陣。CG該方法要求在子空間中的能量范數(shù)到達(dá)極小，且具有短的迭代計算公式。此后發(fā)展得到的與各種預(yù)條件相結(jié)合的方法成為求解大型對稱正定稀疏線性方程組最主要的方法。2.1.1 共軛梯度法（CG） 具體流程為：1.任選初值，計算初值殘差和內(nèi)積（），并置,置j=0。2.計算參數(shù)更新向量與殘向量，若滿足精度要求，則停機(jī)。3.計算，置j=j+1,轉(zhuǎn)（2）精度測試維數(shù)為9*9矩陣：通過直接法所求的的矩陣真實解和通過CG法迭代求的的近似解誤差值的數(shù)量級在10-12。維數(shù)為36*36矩陣：通過直接法所求的的矩陣真實解和通過CG法迭代求的的近似解誤差值的數(shù)

7、量級在10-810-12。維數(shù)為64*64矩陣：通過直接法所求的的矩陣真實解和通過CG法迭代求的的近似解誤差值的數(shù)量級在10-810-9。維數(shù)為169*169矩陣：通過直接法所求的的矩陣真實解和通過CG法迭代求的的近似解誤差值的數(shù)量級在10-910-10。2.1.2 預(yù)條件共軛梯度法共軛梯度法在解決大型線性方程組本身是有不足之處的，當(dāng)系數(shù)矩陣A接近單位陣，共軛梯度法收斂較快，而當(dāng)系數(shù)矩陣A的條件數(shù)大于102時，共軛梯度法就非常緩慢了，因此需要通過對矩陣A的處理，是系數(shù)矩陣的條件數(shù)下降，從而減少迭代次數(shù)和計算量。預(yù)處理方法是目前解決系數(shù)矩陣條件數(shù)問題最有效的方法之一，其原理是選取預(yù)處理矩陣，此矩

8、陣具有以下特征： 為對稱正定的，而且具有與A相似的稀疏性，不但易于求逆，還能改善矩陣的病態(tài)特性，,使得矩陣的特征值比矩陣集中。的選擇尤為重要，因為稀疏矩陣可進(jìn)行分裂式近似分解，,其中L矩陣是下三角矩陣，R矩陣稱為剩余矩陣，則可將看做,用矩陣代替矩陣進(jìn)行迭代。預(yù)處理矩陣得到流程：1 先處理預(yù)處理矩陣第一列元素 i=2 3.n2 循環(huán)嵌套得到預(yù)處理矩陣其它元素的值 k-2.3n如果 則否則 i=k+1.n 得到了預(yù)處理矩陣M 就可以進(jìn)行預(yù)處理共軛梯度法了預(yù)處理共扼梯度法(ICCG) 具體迭代算法如下：1 輸入系數(shù)矩陣A，預(yù)處理矩陣M,右端向量b和初值，于是得到一系列迭代初值置k=1;2 開始迭代

9、K=2.3.n 3 如果達(dá)到精度要求，則輸出x值，結(jié)束；否則k=k+1,轉(zhuǎn)2；精度測試維數(shù)為9*9矩陣：通過直接法所求的的矩陣真實解和通過CG法迭代求的的近似解誤差值的數(shù)量級在10-12。維數(shù)為36*36矩陣：通過直接法所求的的矩陣真實解和通過CG法迭代求的的近似解誤差值的數(shù)量級在10-12。維數(shù)為64*64矩陣：通過直接法所求的的矩陣真實解和通過CG法迭代求的的近似解誤差值的數(shù)量級在10-12。維數(shù)為169*169矩陣：通過直接法所求的的矩陣真實解和通過CG法迭代求的的近似解誤差值的數(shù)量級在10-12。2.2 正交化法 取，則這類Krylov子空間方法稱為正交化方法。GMRES（Genera

10、lized Minimal RESidual algorithm）即廣義極小余量算法，是正交化方法的典型代表，也是目前求解大型稀疏非對稱線性方程組的主要方法之一。GMRES算法是以Galerkin原理為基礎(chǔ)的算法，具體地，設(shè)所求方程為：，其中為一大型不規(guī)則稀疏矩陣，為一給定向量，記和是中的兩個維Krylov子空間，分別由和所張成。取為任一向量，令，則方程等價于：其中。求解的Galerkin原理可敘述為：在子空間中找的近似解，使得殘余向量和中的所有向量正交。即對，有： 令和，其中和分別為和的基底。因此，可以把表示成，。這樣上式就可以寫成：假設(shè)為非奇異矩陣，則不難得到近似解：要求出的值必須先選擇合

11、適的和以及它們的基底，在GMRES算法中，我們?nèi)?，而，簡記。對于這組特定的和的取法，實際上求等價于在中極小化殘余向量的2-范數(shù)。利用Arnoldi過程，及，其中是一個上Hessenberg矩陣，可以找到中的一組標(biāo)準(zhǔn)正交基，不難得： = = = 其中，。由于是正交的，即，所以：于是，在中極小化就等價于在中極小化，具體地，GMRES算法步驟如下：1）選擇初始近似值，計算，記，則；2）定義階矩陣，并設(shè)； 3）計算和的值：4）求解最小二乘問題，得到的解為；5）計算。從理論上講，如果線性無關(guān)，當(dāng)時GMRES算法應(yīng)當(dāng)給出方程的準(zhǔn)確解。但是，不難看出，當(dāng)很大時，計算中需要保存所有的，對于大型問題（）時，這將

12、引起存儲空間過多的需求，因此是不現(xiàn)實的?？朔@個困難的方法是采用如下形式的循環(huán)GMRES（m）算法，具體如下：1）選擇初始近似值，計算，記，則； 2）定義階矩陣，并設(shè)；3）完成Arnoldi過程（不再詳述），得到和；4）極小化得到；5）計算；6）計算，如果滿意，則停止；7），轉(zhuǎn)向2）。在GMRES的具體計算過程中，要根據(jù)具體問題選擇合適的m值，因為m值取得過大會增加存儲的需求，然而太小又會減緩收斂的速度。精度測試維數(shù)為9*9矩陣：通過直接法所求的的矩陣真實解和通過CG法迭代求的的近似解誤差值的數(shù)量級在10-510-6。維數(shù)為36*36矩陣：通過直接法所求的的矩陣真實解和通過CG法迭代求的的近似

13、解誤差值的數(shù)量級在10-5。維數(shù)為64*64矩陣：通過直接法所求的的矩陣真實解和通過CG法迭代求的的近似解誤差值的數(shù)量級在10-410-5。維數(shù)為169*169矩陣：通過直接法所求的的矩陣真實解和通過CG法迭代求的的近似解誤差值的數(shù)量級在10-410-5。2.3 雙正交化方法 在雙正交化過程中取,顯然，A是非對稱矩陣時，這一類方法就是雙正交化方法,Lanczos提出的雙共軛梯度法（BICG法）是其中最基礎(chǔ)的方法。雙共軛梯度法BICG具體流程為：1.計算選取使得置方向向量，并置j=0.2.計算參數(shù)與向量，如果滿足精度，停機(jī)。3.計算向量，置j=j+1.轉(zhuǎn)2BICG法在計算過程總中要使用到矩陣的轉(zhuǎn)

14、置，并且由于振蕩劇烈收斂不好，為避免這些缺陷，許多研究者在BICG法的基礎(chǔ)上發(fā)展了一系列收斂好的方法，如共軛梯度平方法（CGS法），廣義共軛梯度平方法（GCGS法）、共軛梯度穩(wěn)定性方法（BICGSTAB法）等。BICGSTAB算法的主要思想是基于雙邊Lanczos,是基于殘差正交子空間的迭代方法。對于n階線性方程組，假設(shè)初始近似解為，第k次近似解為，相應(yīng)的第k次的殘差。是與不正交的非零向量。是分別與有關(guān)的k階Krylov子空間，并且BICGSTAB算法在Krylov子空間，選擇序列，通過適當(dāng)選擇參數(shù)，使得按下式產(chǎn)生序列的滿足且，同時計算的殘差總是取最小的范數(shù)。這種基于Krylov殘差子空間的迭代方法收斂速度快，精度高，而且穩(wěn)定性好。BICGSTAB算法常常和預(yù)處理技術(shù)相結(jié)合處理方程求解問題,如果預(yù)處理矩陣M是單位陣，則流程為一般的BICGSTAB方法，其流程為1 給定系數(shù)矩陣A，向量b，初始值，最大迭代次數(shù)，最大容許誤差以及預(yù)處理矩陣M，計算，并令k=1,2 如果且，轉(zhuǎn)3，否則終止，輸出3 ，如果或者，算法失敗，終止，否則轉(zhuǎn)44 如果K=1 ，則，否則計算5 求解方程，計算 6 ，如果，否則停止循環(huán)，輸出