8.1代數方程組求解方法概要AX=b(8-1)其中b是已知的N維求解矢量，X為待求的N維未知矢量，A為N×N矩陣。當A中各元素與X無關時，方程組是線性的；反之，當A中元素與X有關時，方程組是非線性的。通常，線性代數方程組解法分為直接解法和迭代解法兩種。所謂直接求解，是指在不法、求逆矩陣法等，其中Gauss消去法是最基本且最有效的直接方法，而矩陣乘積分解直接解法和迭代解法各有其優缺點。直接解法主要優點是精確，對于階數不太高且系數矩陣為滿秩的線性方程組，Gauss消去法是最佳選擇。但由于離散方程(8-1)的系數矩陣對于非線性方程組，如用迭代法求解其線化后的方程，則非線性方程中的迭代求解包由于迭代解法的重要性，近年來，在前述的基本迭代方法基礎上，發展了許多能夠提代過程法（StrongImplicitProcedure—SIP）,多重網格法（Multi-Grid8.2拓展的三對角陣算法方陣；每個節點上的求解函數和方程組的右端元素也不再是標量而是維數為變數個數的矢采用由標量三對角陣算法拓展得到的塊三角陣算法（BlocAjUj?1+BjUj+CjUj+1=Dj(j=1,2,...,J?1,J)(8-2)其中Aj,Βj,Cj均為n×n的系數塊矩陣，Uj?1,Uj,Uj+1為n維矢量函數，顯見Α當j=1和j=2時，按照(8-2)分別有BD1A2U1+B2U2+C2U3=D2E2U2+F2U3=G2如此對j=3,4,...,j逐個向后，可得第j行方程消元Uj?1后的一般關系式為EjUj+FjUj+1=Gj(8-9)Ej=E1Fj?1?A1Bj,Fj=?A1Cj,Gj=E1Gj?1?A1Dj(8-10)(j=2,3,...,J)而,FUj=E1Gj?E1FjUj+1(j=J?1,J?2,...,2,1)(8-13)從前至后計算消元遞推系數Ej,Fj,Gj(j=1,2,...,J)；再按(8-12)和(8-13)從后至前進行回代計算所求矢量函數Uj(j=J,J?1,...2,1)。利用三對角陣算法的求解思路，而發展出所謂的環形三對角陣算法（CTDMA設周期性計算區域離散為N個節點，第1個同，實際需求解的節點數為N?1。對任意節點j，三點式離散代數方程為+Dj(j=1,2,...,N?1)(8-14)Ejφj+1+FjφN?1+Gj這是一個把任意φj既與它的鄰點j+1上的φj+1相聯系，又與最后一個求解節點N?1上的對j=1，由方程(8-14)，(8-15a)，(8-16)，顯然有)φN?1+D1A1將方程(8-16)在j?1點寫出，再代入(8-14)，合并整理有上式兩邊同除以(Aj?CjEj?1)，再與式((j=2,3,...,N?2)進而看如何求出φN?1。對j=N?1，利用(8-15b),方程(8-14)可寫成DN1(8-19)令p1φN?1+F2=r將方程(8-16)在j?1點寫出，并代入方程(8-24)，可以得到其中，系數pj,qjrj的遞推關系是pj=pj?1?qj?1Fj?1,qj=qj?1Ej?1,rj=rj?1+qj?1Gj?1(j=2,3,...N?2)將式(8-25)對j=N?2寫出，有pN?2φN?1=(qN?2+CN?1)φN?2+rN?2(8-27a)將式(8-16)對j=N?2點寫出，有EN?2+FN?2)φN?1+GN?2(8-27b)（1）按式(8-17)和(8-18)計算系數Ej,Fj,Gj(j=1,2,...,N?2)（2）按式(8-20)和(8-26)計算系數pj,qj,rj(j=1,2,...,N?2)（4）按式(8-16)依次回代求出φj(j=N?2,N?3的方程，系數矩陣為五對角陣。二維問題則圍繞離散點(i,j)分別在兩個方向上各有4個鄰點(i±1,j),(i±2,j)及(i,j±1),i1,2,...N-1,N（8-29）過消元后只剩下1個可以直接解出，則消元步需要進行兩次：第一次消去所有的ej，第二次消去變化后的aj'。類似三對角陣算法推導消元系數遞推公式的做法，可以得到五對角陣1）消去系數矩陣中全部的ej元素，矩陣形式方程變為8.3迭代解法的收斂性A=M?Ν(8-36)MX=NX+b(8-37)MX(n+1)=NX(n)+b(8-38)式中X(n+1)和X(n)分別為第n+1次和第n迭代的近似解。上式又可寫成X(n+1)=M?1NX(n)+M?1b=HX(n)+d(8-39)其中H=M?1N稱為迭代矩陣。定義收斂解X與迭代解之間的偏差為迭代誤差ε,則第n次迭代解的迭代誤差為ε(n)=X?X(n)(8-40)定義兩次迭代之間的迭代函數差為δ,則第n+1次與第n迭代之間函數差為δ(n+1)=X(n+1)?X(n)(8-41)定義第n次迭代后原方程的余量為r(n)，則r(n)=b?AX(n)(8-42)將式(8-37)減去(8-38)，并利用迭代誤差函數Mε(n+1)=Nε(n)(8-43)即ε(n+1)=M?1Nε(n)=Hε(n)(8-44)迭代收斂要求。顯然迭代收斂速度與迭代矩陣H的結構特性相關，而這些特性體現在它的特征值和特征矢量上。令H的特征值和對應的特征矢量分別為λk和gk，則有Hgk=λkgk(k=1,2,...,N)(8-45)這里的N為代數方程組未知數的個數。假設gk在N維空間是完備的，則迭代初始誤差矢量ε(0)可通過特征矢量表成為ε(n)=Hε(n?1)=H2ε(n?2)=...=Hnε(0)(8-47)λλ特征值中絕對值最大者稱為矩陣的譜半徑，故迭代收斂的充要條件是迭代矩陣H的譜半徑小于1。記譜半徑為ρ(H)=λ1含λε(n)≈α1(λ1)ng1λλn度會很慢。通常稱(?lnρ(H))為收斂速率。δ(n+1)=X(n+1)?X(n)=(X(n+1)?X)?(X(n)?X)=?ε(n+1)+ε(n)(8-51))(λ?λ)(λ利用(8-53)和(8-54)預估出的ρ(H)和ε(n)ε(n)ε=X?X(n)≤ε0時，則稱迭代達到收斂，并取X(n)作為方程組的迭代近似按照(8-54),迭代誤差既與相鄰兩次迭代解之差δ(n+1)=X(n+1)?X(n)有關，還與迭代矩陣的譜半徑相關，因此，判斷收斂與否的嚴格做法，應先利用(8-53)預估出ρ(H)，再利用(8-54)求出ε(n)，進而看ε(n)是否小于或等于ε0，從而做出相應判斷。2）規定離散方程的迭代余量r(n)=b?AX(n)小于一定的值。由迭代余量定義Aε(n)=r(n)(8-55)代收斂，余量下降應達到迭代誤差的精度要求；而內迭代則可降低1~2個數量級，以便既不同迭代方法有不同的收斂速度，而影響迭代收斂廣、越快，迭代收斂越快。如Jacobi點迭代，每完成一輪迭代，邊界與邊界相鄰的那批節點，傳入范圍僅限一個網格，且收斂速度與掃描方向無關。若采用Gauss-Seidel點迭代，設掃描方向為從區域的傳入以及增加迭代計算中的直接求解成分，可以提高收斂速度。近年來，發展了比ADI迭代方法隱式關聯性更強，直接解法比例更高的強隱過程（StrongImplSIP）迭代法[6]來求解非對稱稀疏矩陣代數方程組；還發展了不需選擇松弛因子參數，而迭代收斂速度不低于松弛迭代，每次迭代都是一個直接求解過程的共軛梯度算法（ConjugateGradient—CG）[7]。原始CG算法可解常物性均分網格上系數矩陣為正定對稱型題，而發展的CG系列算法可以解非均分網格、非對稱系數矩陣的擴散和對3）網格尺度對迭代誤差矢量ε(n)中不同8.4加速迭代收斂的塊修正技術前面說過，促使計算區域內物理量的守恒，可以加速迭代解法的收斂速度。基于這一它的絕對值大小可以反映溫度場的收斂程度。為要求塊上各節點有一個共同的修正值Ti'（或Tj'即在做列掃描時，每列的修正值滿足TT可化為三對角陣方程組，可以用三對角陣解法求解。具體做法是，在做完一輪通常的ADI 迭代后，先在x方向按照方程(8-59)逐列掃描，求出修正值Ti'，再以校正后的值作為新的T值，在y方向逐行掃描，求出修正值Tj'，加到列掃描后的T上，塊修正過程，就是要讓修正的塊在平均意義上滿足方程，符合能量守恒。這樣雖然在一輪ADI迭代中增加了一些計算量，但由于它促8.5時間相關法所謂時間相關法，是一種把穩態問題化為非穩態問題求解的方法。這種方法由于數學我們知道，物理中的穩態問題是相應的非穩態問題在一定的初始條件和邊界條件下經應著穩態問題迭代求解進行了一輪迭代。如二維問題的渦流函數解法中，其流函數Ψ與渦采用交替方向隱式格式離散并迭代求解，則迭如用時間相關法求解方程(8-61)，則在方程中添加一非穩態項A?Ψ?t，使其為采用交替方向隱式格式對該方程進行離散，可以得到與(8-62)形值相當迭代的超松弛，松弛作用的大小取決于2A/Δt。因此可通過調整2A/Δt而使非穩態合參數2A/Δt的取值，使其起到亞松弛的效應，以保證非線性迭代計算的收斂性。8.6強隱過程迭代法非對稱的大型稀疏矩陣，采用通常適用于求解滿秩或近滿秩系數矩陣的完全LU分解法是不經濟的。為此，Stone于1968年提出了一種基于矩陣不完全LU分解法的求解技術，文設矩陣形式的離散方程組為式(8-1),將矩陣A按式(8-36)作和分解。要求分解的兩個矩陣，其M是矩陣A的一個很好近似，其N是一個小矩陣。所謂“小”矩陣，是指N應滿足條件NX≈0。所謂不完全分解，是指將矩陣M而不是矩陣A作上、下三角陣的LU乘積分解，且要求確定L和U的方法簡單。此即M=A+N=LU(8-67)(A+N)X=(A+N)X+(b?AX)(8-68)或MX=MX+(b?AX)(8-69)如果A的近似矩陣M能夠得到，則可以構成如下迭代格式MX(n+1)=MX(n)+(b?AX(n))(8-70)即M(X(n+1)?X(n))=(b?AX(n))(8-71)LUδ(n+1)=r(n)(8-72)Z=Uδ(8-73)迭代解出δ(n+1)后，由δ(n+1)=X(n+1)?X(n)即可從前一次的迭代值X(n)算出本次迭代值X(n+1)。由上分析顯見，只要能找到M作乘積分解的下三角陣L和上三角陣U，迭代計算確定L和U的簡單方法是[5]：讓L和U均為稀疏三角陣，且非零元素所在位置與矩陣A一致。進而根據矩陣乘積運算法則，推知矩陣M的稀疏結構形態以及各元素與L和U中元素的對應關系。再假定矩陣N具有與M一樣的稀疏結構形態，并按照N為小矩陣的要求，另加上一些合理的簡化近似假設[6]，就可以先得到矩陣N各元素的表達式。將N的元素和A的對應元素相加，即得M的元素表達式。有了M的元素表達式及其這些元素與L和U中元素的對應關系，通過對照比較，即可最終得到L和U各元素的值。的導熱問題和對流換熱問題。但這種方法的一輪迭代計算中，確定中間計算量Z各分量和8.7多重網格法多重網格方法是基于迭代誤差矢量的衰減過程與迭代求解方程的網格尺度密切相關的在均勻網格上對節點j作中心差分離散，得如下代數方程x2fj(8-76)令迭代誤差函數為ε,則有+εj代以上三式入式(8-76)，并減去式(8-77εj態問題作格式穩定性分析時，我們采用vonNeumann方法分析誤化特性，同樣，亦可用該方法分析迭代誤差函數在迭代過=eikxj(8-79)εj(n+1)=geikxj,εj式中k稱波數，g表示誤差函數相應波數為k的分量的振幅在一次迭代中衰減程度，稱為衰減因子。波數k和它對應的波長λ之積為2π。代式(8-80)入式(8-78)，得geik?2geikxj+eik(8-81)這說明，隨著θ增大，g逐漸減小，即迭代衰減加快。連續5別很大。因為θ=kΔx=2πΔx/λ,網格步長Δx一定，較大的θ值對應波長較短的高頻分分量由于具有較大θ一開始就迅速衰減，但長波分量卻因只有較小θ而衰減緩慢，致使迭使迭代在高θ條件下進行。多重網格方法就是為了解決上述在固定網格上實施迭代所存在的問題而發展起來的[8,9]，目的在于使迭代計算中，不同波長的誤差分量能得到比較均勻一致的衰減，促進迭代收斂。這種方法將迭代在步長尺度不同的多套網格上進行：先在步長小的細網格上迭也可采用直接求解。以下針對求解方程(8-1),即AX=b，討論實施多種網格計算的主要相格上求解未知量X自身，其它網格上求解未知量的迭代誤差量或曰修正量εk，稱之為修正定算子；而從疏網格到密網格的傳遞稱為延拓（或插），該法僅在最密網格上以未知數X作為求解對象。令在某個初值下迭代數次后得到近似 k值為X，此時方程的迭代余量rk k令Xk為該重上的準確解，則該重上迭代誤差εk是0=bk-AkXk(8-86)對于線性方程，A不隨Xk改變，因此有Akεk=rk(8-88)速度逐漸變慢時，轉入第k-1重計算，相應于(8-88)的計算方程為Ak-1εk-1=I-1rk(8-89) kk其中Xnew和Xold分別為第k重網格上改進的新值和上一迭代循環在第k重網格上迭代所如果網格重數多于兩重，則按式(8-89)迭代數次得到εk-1后，隨之計算第k-1重上的余量rk-1=I-1rk-Ak-1εk-1，再按式(8-89)往下計算εk-2，如此下去，直到最疏網格層，得到εk-1,εk-2,...ε1；進而按(8-90)一重重返回，依次得到更新的ε2,ε3,...εk-1，當達到 k最密網格時，就求出未知量在多重網格一個循環中的更新值Xnew。這一循環過程常需進行較粗的第k-1重網格繼續迭代計算。此時在第k-1重網格上所求解的任然是第k重網格上 k-1正確解的一個近似值，記為X，所求解的代數方程的矩陣形式為對線性代數方程，Ak-1和bk-1各元素是常數，但對非線性方程，Ak-1各元素要根據第k重網格已得到的解Xk之值確定是第k重網格上的計算余量。由于第k-1重上所求解的是第k重網格上的解，因此要把該余量從第k重網格傳遞到第k-1重網格上。傳遞方法通過限定算子I-1來確定。如果只有兩重網格，則在第k-1重網格上迭代進行到一定程度時，再把所得解 k-1X傳遞回第k重網格上，以改進密網格上的解。令第k重網格限定過程之前迭代計算出的值為即中的第k重網上的改進解為X，則其中，表示通過限定算子將第k重網格算得的解轉移到第k-1重網格上，而表示第k-1重網格上得到的解的修正值，這個修正值再通過延拓算子I-1傳遞到第k重網格上，以修正第k重網格上前一輪迭代所算的老值，從而獲得第k重網但此時k-1重網格上的迭代計算次數不必過多，即可將解向k-2重網格傳遞，在k-2重網格上形成k重網格解的一組代數方程，其形式同式（8-91）,只要把相應的上標指數各減1就是。如此下去，一直到最疏那重網格時（相應的k-1=1可以將代數方程迭代求解3．確定限定算子I-1和延拓算子I-1的常用方法V循環是最常用的形式。圖8－3(a)所示的這8.8共軛梯度法共軛梯度（ConjugateGradient－CG）法，是一類迭代解法，但在每一輪迭代中是直接求解的。它先是用來求解系數矩陣A是對稱正定型的代數方程組，進而推廣到A為非對稱型的情況，并發展成一種系列算法。這類算法，對于一個N階的線性代數方程組，最多不一個系數矩陣A為對稱正定的N階線性代數方程組AX=b的求解問題，等價于一個按A所作的關于X的二次函數的極小值求解問題。式中上標T于是在極值點X*的鄰域t處，有即f(X*+t)-f(X*)>0(8-98)數f(X)的極小值點，也就求得了代數方程組AX=b的解。欲對f(X)直接求在定義域RN上的極小值點是困難的，一般可通過在RN的某些給定搜索方向p(k)上，逐次求得向f(X)的極小值過渡的極小值序列{X(k)}，要求該序列的極限就是f(X)的極小值。在RN上取X(0)為f(X)最小值的初始近似，并給出從X(0)點出發的一個方向p(0)，進而在方向（直線）X(0)+αp(0)上求出f(X)最小值的下一個近似點X(1)。再取方向p(1)，在方向（直線）X(1)+αp(1)上求出f(X)最小值的又下一個近似點X(2)。如此下去，就得到f(X)極小化的近似序列{X(k)}，該序列的極限就是f(X)取極小值的點。實施這個過程要解決兩個問題，即如何定搜索方向p(k)，又如何由X(k)和p(k)求得f(X)在X(k)+αp(k)方向上的極小化序列點X(k+1)。先討論第二個問題。對確定的X(k)和p(k)，X(k+1)將只與α相關。令?(α)=f(X(k)+αp(k))則f(X)在方向X(k)+αp(k)上的極小值點應滿足令g(k)=AX(k)-b=▽f(X(k))(8-99)即g(k)為f(X)在X(k)上的梯度（或稱斜量），于是有而X(k+1)=X(k)+αkp(k)(8-101)且由于X(k+1)是X(k)+αkp(k)方向上f(X)的最小值近似點，故總有f(X(k+1))<f(X(k))由此可知，對于給定的X(0)和方向序列{p(k)}，由(8-100)和(8-101)就可得到使f(X)逐漸變小的序列{X(k)}。只要p(k)取得恰當，{X(k)}能夠很快收斂于f(X)在RN上的極小值點X*，即方程AX=b的解。進而討論如何選取搜索方向序列{p(k)}的問題。最初人們是取f(X)的負梯度方向-▽f(X(k))作為p(k)，構成了所謂的梯度法，或曰最速下降法。但梯度法只有線性收斂速率，且其收斂因子依賴f(X)的二階導數矩陣H(X)的特征值λ,當H(X)的條件數ρ=λmaxλmin很大成為病態矩陣時，會嚴重影響此法的收斂速度[11]。為此將f(X)的負梯度方向改為矩陣A的一組共軛向量系，發展成共軛梯度方法。設{p(k)}是RN上的非零向量系，A為正定矩陣，如果{p(k)}滿足(p(i){p(k)}為A的共軛向量系。共軛向量系具有以下性質1）在RN上線性無關2）其m≤N?1，即共軛向量系中向量個數不會超過N3）以為迭代方向作出的迭代系列{X(k)}滿足(g(k)即f(X)在X(k)上的梯度（或斜量）g(k)與共軛向量p(1),p(2),...p(k?1)正交。按照第三條性質，如能找到矩陣A的包含有N個向量組成的共軛向量系{p(k)}?1，則以該組向量作為迭代方向所得到的最后一個點X(N)的梯度g(N)正交于共軛向量系中所有向量，而這些向量是線性無關的，因此只能g(N)=AX(N)?b=0。這就是說至多迭代N次，理論上講方程p(0)和g(1)導出p(1)。如此下去，設由p(k?1)和g(k)可推得的p(k)具有以下遞推關系p(k)=?g(k)p(k?1)(k)p(k?1)(8-104)的N階代數方程組AX=b的求解問題等價于用A構成的一個二次函數的極小值的求解問題；而為了求解該二次函數的極小值，采用正定矩陣A的共軛向量系作為搜索極小值的迭代方向，在最多不超過方程組階數N的迭代