油藏?cái)?shù)值模擬中大型稀疏矩陣算法綜述

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

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

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

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

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

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

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

8、陣具有以下特征： 為對(duì)稱(chēng)正定的，而且具有與A相似的稀疏性，不但易于求逆，還能改善矩陣的病態(tài)特性，,使得矩陣的特征值比矩陣集中。的選擇尤為重要，因?yàn)橄∈杈仃嚳蛇M(jìn)行分裂式近似分解，,其中L矩陣是下三角矩陣，R矩陣稱(chēng)為剩余矩陣，則可將看做,用矩陣代替矩陣進(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 開(kāi)始迭代

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

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

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

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

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

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