第4章 解線性方程組的迭代法用迭代法求解線性方程組與第4章非線性方程求根的方法相似，對方程組進行等價變換，構造同解方程組(對可構造各種等價方程組，如分解，可逆，則由得到)，以此構造迭代關系式 （4.1）任取初始向量，代入迭代式中，經計算得到迭代序列。若迭代序列收斂，設的極限為，對迭代式兩邊取極限即是方程組的解，此時稱迭代法收斂，否則稱迭代法發散。我們將看到，不同于非線性方程的迭代方法，解線性方程組的迭代收斂與否完全決定于迭代矩陣的性質，與迭代初始值的選取無關。迭代法的優點是占有存儲空間少，程序實現簡單，尤其適用于大型稀疏矩陣；不盡人意之處是要面對判斷迭代是否收斂和收斂速度的問題。可以證明迭代矩陣的與譜半徑是迭代收斂的充分必要條件，其中是矩陣的特征根。事實上，若為方程組的解，則有再由迭代式可得到由線性代數定理，的充分必要條件。因此對迭代法(4.1)的收斂性有以下兩個定理成立。定理4.1 迭代法 收斂的充要條件是。定理4.2 迭代法收斂的充要條件是迭代矩陣的譜半徑因此，稱譜半徑小于1的矩陣為收斂矩陣。計算矩陣的譜半徑，需要求解矩陣的特征值才能得到，通常這是較為繁重的工作。但是可以通過計算矩陣的范數等方法簡化判斷收斂的工作。前面已經提到過，若|A|p矩陣的范數，則總有。因此，若，則必為收斂矩陣。計算矩陣的1范數和范數的方法比較簡單，其中于是，只要迭代矩陣滿足或，就可以判斷迭代序列是收斂的。要注意的是，當或時，可以有，因此不能判斷迭代序列發散。在計算中當相鄰兩次的向量誤差的某種范數小于給定精度時，則停止迭代計算，視為方程組的近似解（有關范數的詳細定義請看3.3節。）4.1雅可比（Jacobi）迭代法4.1.1 雅可比迭代格式雅可比迭代計算元線性方程組（4.2）寫成矩陣形式為。若將式（ 4.2）中每個方程的留在方程左邊，其余各項移到方程右邊；方程兩邊除以則得到下列同解方程組：記，構造迭代形式或 （4.3）迭代計算式（4.3）稱為簡單迭代或雅可比迭代。任取初始向量 ，由式（4.3）可得到迭代向量序列雅可比迭代矩陣設由，得到等價方程：記不難看出，正是迭代式（4.3）的迭代矩陣，是常數項向量。于是式（4.3）可寫成矩陣形式：（4.4）其中：雅可比迭代算法下面描述解線性方程組的雅可比迭代算法，為了簡單起見，在算法中假定矩陣滿足雅可比迭代要求，即，并設由系數矩陣構造迭代矩陣是收斂的。1定義和輸入系數矩陣與常數項向量的元素。2FOR i：=1，2，n/假定，形成常數項向量 FOR j:=1,2,n /形成迭代矩陣元素3 / 賦初始值，x1和x2分別表示和4WHILE x1:=x2 x2:=B*x1+g/ FOR u:=1,2,n/ s:= gu;/FOR v:=1,2,n s:=s+buv*x1v;/ x2u:=s； ENDWHILE5輸出方程組的解例4.1 用雅可比方法解下列方程組：解：方程的迭代格式：或雅可比迭代收斂。取初始值，計算結果由表4.1所示。表4.1 計算結果0111 1-1.51.60.90.252-1.252.081.090.483-0.9152.0681.0170.3354-0.95751.98640.98470.08165-1.014451.988440.997110.056956-1.007222.002311.00260.013877-0.9975432.001971.000490.009687方程組的準確解是4.1.2 雅可比迭代收斂條件對于方程組，構造雅可比迭代格式其中，。當迭代矩陣的譜半徑時，迭代收斂，這是收斂的充分必要條件。迭代矩陣的某范數時，迭代收斂。要注意的是范數小于1只是判斷迭代矩陣收斂的充分條件，當迭代矩陣的一種范數|B|1，并不能確定迭代矩陣是收斂還是發散。例如，則，但它的特征值是0.9和0.8。是收斂矩陣。當方程組的系數矩陣具有某些性質時，可直接判定由它生成的雅可比迭代矩陣是收斂的。定理4.3 若方程組的系數矩陣，滿足下列條件之一，則其雅可比迭代法是收斂的。（1）為行對角占優陣，即（2）為列對角占優陣，即證明：（1）雅可比迭代矩陣其中（2）為列對角優陣，故為行對角占優陣，由系數矩陣構造的迭代矩陣為行對角占優陣，則有又得到而，得由系數矩陣構造的雅可比迭代矩陣收斂。（如矩陣既是行對角占優陣，也是列對角占優陣）定理4.4若方程組系數矩陣 為對稱正定陣，并且也為對稱正定，則雅可比迭代收斂。4.2 高斯-塞德爾（Gauss-Seidel）迭代法高斯-賽德爾迭代計算在雅可比迭代中，用的值代入方程（4.2）中計算出的值，的計算公式是事實上，在計算前，已經得到的值，不妨將已算出的分量直接代入迭代式中，及時使用最新計算出的分量值。因此的計算公式可改為：即用向量計算出的值，用向量計算出的值，用向量計算出的值，這種迭代格式稱為高斯塞德爾迭代。對于方程組AX=y ，如果由它構造高斯-塞德爾迭代和雅可比迭代都收斂，那么，多數情況下高斯塞德爾迭代比雅可比迭代的收斂效果要好，但是情況并非總是如此。構造方程組的高斯-塞德爾迭代格式的步驟與雅可比類似，設將式（4.1）中每個方程的留在方程的左邊，其余各項都移到方程的右邊；方程兩邊除以，得到下列同解方程組：記，對方程組對角線以上的取第步迭代的數值，對角線以下的取第步迭代的數值，構造高斯塞德爾迭代形式：（4.5）例4.2 用高斯-塞德爾方法解方程組：解：方程的迭代格式：取初始值有時，時，計算結果如表4.2所示。表4.2 計算結果012340002.52.11.142.50.88 2.0040.98761.621.0042 1.9984 1.00060.12421.0005 2.0002 1.00000.0037高斯塞德爾迭代矩陣設寫成等價矩陣表達式：構造迭代形式：有 （4.6）則高斯-塞德爾迭代式（4.4）為 （4.7） 稱為高斯-塞德爾迭代矩陣。高斯-塞德爾迭代算法高斯塞德爾迭代的程序實現與雅可比迭代步驟大致相同，對于方程組，在前面的雅可比算法中，假定雅可比迭代矩陣為表示表示,其迭代核心部分是。下面的算法給出由和計算的過程，省略了形成迭代矩陣和對初始化的部分。雅可比迭代的核心部分：WHILEfor(u:=1;u=n,u+) x1u:=x2ufor(i:=1;j=n;j+) s:=gi; for(j:=1;i=n;i+) s:=s+bix2j /注意x2jENDWHILE在高斯-賽德爾迭代計算中并不需要形成迭代矩陣，由式（4.5）可知在計算中只要形成矩陣在程序中令向量。它的核心部分是計算迭代式，計算中只需及時將放到的位置上。高斯-塞德爾迭代的核心部分：WHILE for(u:=1;u=n;u+) x1u:=x2u for(i:=1;j=n;j+) s:=gi; for(j:=1;j=n;j+) s:=s+bix2 j /注意x2jx2i:=s ENDWHILE上列算法是在假定迭代收斂的前提下，使用當型（WHILE）結構控制循環。更一般地，可將上列算法中WHILE循環改為FOR循環，通過控制循環次數和觀測計算誤差終止循環。屆將上列算法中WHILE語句改為WHILE 循環次數這時在程序中要增加循環變量的設定和運算。判斷高斯塞德爾迭代收斂的方法與判斷雅可比迭代收斂類似，一方面從高斯-塞德爾迭代矩陣獲取信息，當或的某種范數時，迭代收斂；另一方面，直接根據方程組系數矩陣的特點作出判斷。定理4.5 若方程組系數矩陣A為列或行對角優時，則高斯塞德爾迭代收斂。定理4.6 若方程組系數矩陣A為對稱正定陣，則高斯塞德爾迭代收斂。例4.3 方程組中，,證明當 時Gauss-Seidel法收斂，而Jacobi迭代法只在時才收斂。解：對法，因為是對稱矩陣，因此只要證時正定即可，由順序主子得而 得于是得到時故對稱正定，法收斂。對 法，可根據定理4.2，由于迭代矩陣 即 是 法收斂的充要條件，故法只在時才收斂。當時，法收斂，而，法不收斂，此時不是正定的。4.3 逐次超松弛（SOR）迭代法逐次超松弛迭代計算方程組的雅可比迭代形式記其中是下三角矩陣，是上三角矩陣。得高斯-塞德爾的迭代形式：記，有這樣可視為的修正量，而恰是由加修正量而得，如果將改為)加上修正量乘一個因子，迭代格改為：整理得 (4.8)這里為修正因子，稱為松弛因子，而式（4.8）稱為松弛迭代。迭代的分量形式為(4.9)式（4.9）稱為松弛迭代的計算格式。例4.4 給定方程組用SOR法求解,取解：用SOR迭代公式可得取初始值：如果用高斯-賽德爾迭代法迭代72次得：用SOR迭代法，只須迭代25次即得： 逐次超松弛迭代算法下列算法假定迭代矩陣收斂，否則要在WHILE循環中增加判斷條件。1定義和輸入系數矩陣與常數項向量的元素，輸入松弛因子的值。2FOR i:=1,2,n / 假定，形成常數項向量 FOR 34. WHILE ; ENDWHILE5輸出.松弛迭代矩陣將式（4.9）中含有與的項分別放在方程兩邊： 用 代入上式，有 則松弛因子為的迭代矩陣為 其中 定理4.7 逐次超松弛迭代收斂的必要條件。定理4.8 若為正定矩陣，則當時，逐次超松弛迭代恒收斂。以上定理給出了逐次超松弛迭代因子的范圍。對于每個給定的方程組，確定究竟取多少為最佳，這是比較困難的問題，對某些特定的方程組，我們可以得到一些理論結果。通常，把的迭代稱為亞松弛迭代，把的迭代稱為高斯-塞德爾迭代，而把的迭代稱為松弛迭代。4.4逆矩陣計算在線性代數中逆矩陣是按其伴隨矩陣定義的，若則方陣可逆，且，其中為的伴隨矩陣。要計算個階的列式才能得到一個伴隨矩陣，在數值計算中因其計算工作量大而不被采用。通常對做行的初等的效換，在將化成的過程中得到。在數值計算中，這仍然是一種行之有效的方法。由逆矩陣的定義令，有化為個方程組j是第個分量為1，其余分量為0的維向量?；蛴洖椋?。用直接法或迭代法算出也就完成了逆矩陣計算。如果依次對用高斯若爾當消元法，組合起來看有（當然也能組合起來做）：這正是在線性代數中用初等變換計算逆矩陣的方法。由此可見，計算一個階逆矩陣的工作量相當于解個線性方程組。在數值計算中常常將計算矩陣逆的問題轉化為解線性方程組的問題。例如，已知方陣和向量有迭代關系式，在計算中不是先算出，再作與的乘積得到；而將作為線性方程組系數矩陣，求解方程組作為常駐數項解出。4.5程序示例程序4.1用雅可比迭代求解方程組：算法描述輸入系數矩陣及常數項向量c。按雅可比迭代公式：求解。程序源碼/ Purpose:雅可比迭代求解線性方程組 /#include#include#define MAX-N 20 /方程的最大維數#define MAXREPT 100#define epsilon 0.00001 /求解精度int main( ) int n; int I, j, k; double err; static double a MAX-N MAX-N,bMAX-N MAX-N,cMAX-N,g MAX-N; static double xMAX-N,nxMAX-N; printf(nInput n value(dim of Ax=c):); /輸入方程的維數 scanf(%d,&n); if(nMAX-N) printf(The input n is larger then MAX-N,please redefine the MAX-N.n); if(n=0) printf(please input a number between 1 and %d.n), MAX-N;return 1;/輸入Ax=c的系數矩陣A.PRINTF(Now input the matrix a(I,j,)=0,，%d:n,n-1);for(i=0; in; i+) /形成迭代矩陣b for(j=0; jn; j+)bij=aij/aii;gi=ci/aij; /為了簡化程序，假設aii!0,/否則要附加對aii的判斷及其相應的處理for(i = 0; i MAXREPT; i+) for(j=0; jn; j+ ) nxj=gj; for(j=0, jn; j+ )for(k=0, kn; k+ )if (j=k)continue;nxj+ =bjk*xk; /迭代 err = 0 ;for(j=0, jn; j+ )if (errfabs(nxj-xj)err=fabs(nxj-xj); /誤差for(j=0, jn; j+ )xj=nxj;if(errepsilon)printf(“Solvex_I=n”); /輸出for(i=0, in; i+ )printf(“%fn”,xi);return 0; printf(“After% d repeat, no result n”, MAXREPT ); /輸出return 1; 計算實例用雅可比方法解下列方程組： 程序輸入輸出Inprt n value(dim of Ax=c):3Now input the matrix a(j, j), i, j =0,2:64 -13 -1 2 90 1 1 1 40Now input the matrix c(i), i= 0, 2:14 -5 20Solvex_i=0.2295470.0661300.492608程序4.2用逐次超松弛迭代（ 作參數）求解方程組： 算法描述輸入矩陣及列向量c。按因子為的超松弛迭代公式：求解。程序源碼/Purpose:超松弛迭代求解線性方程組 /# include# include# define MAX_N 20 / /方程的最大維數# define MAXREPT 100# define epsilon 0.00001 / /求解精度int main( ) int n;int i, j, k ;double err, w;static double aMAX_NMAX_N, bMAX_NMAX_N,cMAX_N,cMAX_N,gMAX_N ;static double aMAX_N , nxMAX_N;PRINTF(“nlnput n value(dim of Ax=c ):”); /輸入方程的維數Scanf(“% d”, &n);If (nMAX_N) printf(“The input n-is larger then MAX_N, please redefine the MAX_N.n”); return 1;if (n=0) printf(“Please input a number between 1 and % d. n”, MAX_N); return 1； / /輸入的A矩陣printf(“Now input the matrix a(i, j), i, j =0,% d: n”, n-1);for(i=0; in; i+)for(j=0; jn; j+) scanf(“%1f” , &sij);printf(“Now input the matrix c(i), i =0,% d: n”, n-1);for(i=0; in; i+) scanf (“%1f” , &ci);printf(