第3章 線性方程組直接求解

1、第3章 線性方程組直接求解v順序高斯消元法順序高斯消元法v列主元高斯消元法列主元高斯消元法v全主元高斯消元法全主元高斯消元法v高斯約當消元法高斯約當消元法v消元形式的追趕法消元形式的追趕法vLULU分解法分解法v矩陣形式的追趕法矩陣形式的追趕法v平方根法平方根法3.1 引言引言n階線性方程組的一般形式為：階線性方程組的一般形式為： a11x1+a12x2+a1nxn=b1 a21x1+a22x2+a2nxn=b2 an1x1+an2x2+annxn=bnn階線性方程組一般形式的矩陣形式為：階線性方程組一般形式的矩陣形式為：bAx 階數(shù)在階數(shù)在100150以上的線性方程組為高階線性方程組以上的線

2、性方程組為高階線性方程組階數(shù)在階數(shù)在100150以下的線性方程組為低階線性方程組以下的線性方程組為低階線性方程組實際應(yīng)用中，經(jīng)常見到高階稀疏線性方程組，低階稠密線性方程組，實際應(yīng)用中，經(jīng)常見到高階稀疏線性方程組，低階稠密線性方程組，對稱正定線性方程組，帶狀線性方程組等等。對稱正定線性方程組，帶狀線性方程組等等。系數(shù)矩陣的大部分元素為零元素的線性方程組為稀疏線性方程組系數(shù)矩陣的大部分元素為零元素的線性方程組為稀疏線性方程組系數(shù)矩陣的大部分元素為非零元素的線性方程組為稠密線性方程組系數(shù)矩陣的大部分元素為非零元素的線性方程組為稠密線性方程組克萊姆法則并不實用。常用的數(shù)值解法主要分為兩類：克萊姆法則并

3、不實用。常用的數(shù)值解法主要分為兩類： 直接求解方法是指經(jīng)過有限次四則運算，求出線性方程組精確解的方法。直接求解方法是指經(jīng)過有限次四則運算，求出線性方程組精確解的方法。 迭代求解方法是指構(gòu)造一種迭代方法，由某個（套）迭代初值（粗略解），迭代求解方法是指構(gòu)造一種迭代方法，由某個（套）迭代初值（粗略解），得到近似解序列，用序列極限逐步逼近線性方程組精確解的方法。得到近似解序列，用序列極限逐步逼近線性方程組精確解的方法。用順序高斯消元法求解線性方程組的過程分為消元過程和回代過程。消元用順序高斯消元法求解線性方程組的過程分為消元過程和回代過程。消元過程是自上而下，把原線性方程組化為上三角方程組的過程；回

4、代過程是過程是自上而下，把原線性方程組化為上三角方程組的過程；回代過程是自下而上，對這個上三角方程組進行求解的過程。自下而上，對這個上三角方程組進行求解的過程。 3.2 順序高斯消元法順序高斯消元法消元的過程是對線性方程組做同解變換，把原線性方程組變?yōu)樯先欠匠滔倪^程是對線性方程組做同解變換，把原線性方程組變?yōu)樯先欠匠探M的過程。組的過程。 消元過程消元過程n階線性方程組在消元過程一開始，可以表示為：階線性方程組在消元過程一開始，可以表示為： )1(1,)1(,2)1(2,1)1(1 ,)1(1, 2)1(, 22)1(221)1(21)1(1, 1)1(, 12)1(121)1(11nn

5、nnnnnnnnnnnaxaxaxaaxaxaxaaxaxaxa外層循環(huán)每循環(huán)一輪，就用涉及到的矩形區(qū)域中上面的第外層循環(huán)每循環(huán)一輪，就用涉及到的矩形區(qū)域中上面的第1行消去它下面行消去它下面各行左面的第各行左面的第1列，此矩形區(qū)域最上列，此矩形區(qū)域最上1行和最左行和最左1列在之后的消元過程中不列在之后的消元過程中不再涉及（上標不再改變），除此之外的右下矩形區(qū)域中所有的系數(shù)被更新再涉及（上標不再改變），除此之外的右下矩形區(qū)域中所有的系數(shù)被更新一遍，且上標一遍，且上標k自增自增1。第。第k次消元的初始狀態(tài)如左圖所示。次消元的初始狀態(tài)如左圖所示。 3.2 順序高斯消元法順序高斯消元法在第在第k次更新

6、的初始，涉及到的待更新矩形區(qū)域中，系數(shù)的上標都為次更新的初始，涉及到的待更新矩形區(qū)域中，系數(shù)的上標都為k，此矩形區(qū)域中最左上角的系數(shù)的行標、列標也為此矩形區(qū)域中最左上角的系數(shù)的行標、列標也為k。更新的過程，就是用。更新的過程，就是用第第k行消去第行消去第k1行至第行至第n行的第行的第k列，并且讓第列，并且讓第k1行至第行至第n行，第行，第k1列至第列至第n1列的矩形區(qū)域的系數(shù)更新，上標變?yōu)榱械木匦螀^(qū)域的系數(shù)更新，上標變?yōu)閗1。在此后的消元過。在此后的消元過程中，第程中，第k行及之上和第行及之上和第k列及其左邊的系數(shù)不再變化。列及其左邊的系數(shù)不再變化。)(1,)(,)(,)(1,)(,)(,)2

7、(1, 2)2(, 2)2(, 22)2(22)1(1, 1)1(, 1)1(, 12)1(121)1(11knnnknnkkknknknknkkkkknnnkknnnkkaxaxaaxaxaaxaxaxaaxaxaxaxa)(1,)(,)2(1, 2)2(, 22)2(22)1(1, 1)1(, 12)1(121)1(11nnnnnnnnnnnnnaxaaxaxaaxaxaxa消元結(jié)束之后得到的上三角方程組中，每向下一行，上標消元結(jié)束之后得到的上三角方程組中，每向下一行，上標k增增1，即上標，即上標k=行標行標i，如右圖所示。，如右圖所示。3.2 順序高斯消元法順序高斯消元法順序高斯消元法的

8、算法順序高斯消元法的算法輸入原方程組的階數(shù)輸入原方程組的階數(shù)n。輸入原方程組的增廣矩陣輸入原方程組的增廣矩陣ann+1。for(k=0;k=n-2;k+)，循環(huán)，循環(huán)1次消去次消去1列。列。 （消元過程）（消元過程） for(i=k+1;i=n-1;i+)，循環(huán)，循環(huán)1次更新次更新1行。行。 行乘子行乘子aik/=-akk。 for(j=k+1;j=0;k-) （回代過程）（回代過程） s=0; for(j=k+1;j=n-1;j+) s+=akj*xj; xk=(akn-s)/akk;輸出原方程組的解輸出原方程組的解xn。3.2 順序高斯消元法順序高斯消元法順序高斯消元法對應(yīng)的程序（順序高斯

9、消元法對應(yīng)的程序（1/2）#include #include #define MAXSIZE 50void input(double aMAXSIZEMAXSIZE+1,long n);void output(double xMAXSIZE,long n);void main(void)double aMAXSIZEMAXSIZE+1,xMAXSIZE,s;long n,i,j,k;printf(n請輸入原方程組的階數(shù)：請輸入原方程組的階數(shù)：);scanf(%ld,&n);input(a,n);for(k=0;k=n-2;k+)for(i=k+1;i=n-1;i+)aik/=-akk;

10、for(j=k+1;j=0;k-)s=0;for(j=k+1;j=n-1;j+)s+=akj*xj;xk=(akn-s)/akk; output(x,n);3.2 順序高斯消元法順序高斯消元法順序高斯消元法對應(yīng)的程序（順序高斯消元法對應(yīng)的程序（2/2）void input(double aMAXSIZEMAXSIZE+1,long n)long i,j;printf(n請輸入原方程組的增廣矩陣：請輸入原方程組的增廣矩陣：n);for(i=1;i=n;i+)for(j=1;j=n+1;j+)scanf(%lf,&ai-1j-1);void output(double xMAXSIZE,l

11、ong n)long k;printf(n原方程組的解為：原方程組的解為：n);for(k=1;k=n;k+)printf( %lf,xk-1);3.2 順序高斯消元法順序高斯消元法當方程組的階數(shù)增長時，順序高斯消元法的運算量增長速度，比克萊姆法當方程組的階數(shù)增長時，順序高斯消元法的運算量增長速度，比克萊姆法則求解線性方程組的運算量增長速度要慢得多。則求解線性方程組的運算量增長速度要慢得多。 順序高斯消元法的乘除次數(shù)順序高斯消元法的乘除次數(shù)MD 。 33n由克萊姆法則，可知當線性方程組的系數(shù)行列式由克萊姆法則，可知當線性方程組的系數(shù)行列式不等于不等于0時時，則線性方程，則線性方程組有唯一解，否

12、則線性方程組可能無解（如組有唯一解，否則線性方程組可能無解（如0 x0y=1），或有無窮多解），或有無窮多解（如（如0 x0y=0）。）。 性質(zhì)：當線性方程組有唯一解時，要想讓順序高斯消元法能求出解，需要性質(zhì)：當線性方程組有唯一解時，要想讓順序高斯消元法能求出解，需要 （k=1,2,n）都不為）都不為0。 )(,kkka定義：定義： 對求解起突出作用，稱為主元素。對求解起突出作用，稱為主元素。 )(,kkka定理：對定理：對n階線性方程組階線性方程組 ，順序高斯消元法能求出解的充要條件是：，順序高斯消元法能求出解的充要條件是：系數(shù)矩陣系數(shù)矩陣 的所有順序主子式的所有順序主子式i(i=1,2,n

13、)均不為均不為0。 bAx A3.3 列主元高斯消元法列主元高斯消元法順序高斯消元法的缺陷有：順序高斯消元法的缺陷有： 當線性方程組的系數(shù)行列式不等于當線性方程組的系數(shù)行列式不等于0，且存在，且存在1個主元素等于個主元素等于0時，時，線性方程組有唯一解，但順序高斯消元法不能求出解。線性方程組有唯一解，但順序高斯消元法不能求出解。 交換增廣矩陣的交換增廣矩陣的2行（行（2個方程的位置互換）或交換個方程的位置互換）或交換2列（列（2個變元的位置互個變元的位置互換），得到的方程組與原方程組同解。通過交換行、列，來避免主元為換），得到的方程組與原方程組同解。通過交換行、列，來避免主元為0和小主元，然后

14、再消元求解的方法，稱為選主元高斯消元法。和小主元，然后再消元求解的方法，稱為選主元高斯消元法。 一、列主元高斯消元法的主要思想一、列主元高斯消元法的主要思想列主元高斯消元法是一種常見的選主元高斯消元法，與順序高斯消元法相列主元高斯消元法是一種常見的選主元高斯消元法，與順序高斯消元法相比，它的改進之處是：在用主元素消去第比，它的改進之處是：在用主元素消去第k1行至第行至第n行的變元行的變元xk之前，之前，從第從第k行至第行至第n行所有的行所有的xk系數(shù)中，尋找最大系數(shù)中，尋找最大的系數(shù)的系數(shù)，若，若maxik，則交換，則交換第第k行和第行和第maxi行，使最大的系數(shù)成為主元。行，使最大的系數(shù)成為

15、主元。 當線性方程組有唯一解時，列主元高斯消元法必能求出解。當線性方程組有唯一解時，列主元高斯消元法必能求出解。 列主元高斯消元法能避免小主元時舍入誤差被放大的情況，提高解的精度。列主元高斯消元法能避免小主元時舍入誤差被放大的情況，提高解的精度。 小主元（主元素遠小于其他元素）時能求出解，但求出的解誤差很大。小主元（主元素遠小于其他元素）時能求出解，但求出的解誤差很大。 列列主主元元高高斯斯消消元元法法的的算算法法輸入原方程組的階數(shù)輸入原方程組的階數(shù)n。輸入原方程組的增廣矩陣輸入原方程組的增廣矩陣ann+1。for(k=0;k=n-2;k+)，循環(huán)，循環(huán)1次消去次消去1列。列。 （消元過程）（

16、消元過程） 暫設(shè)最大系數(shù)暫設(shè)最大系數(shù)max=akk，最大系數(shù)所在行的行號，最大系數(shù)所在行的行號maxi=k。 for(i=k+1;i|max| Y N 令令max=aik;maxi=i; 系數(shù)最大值系數(shù)最大值max為為0 Y N break; maxik Y N for(j=k;j=n;j+)，交換，交換2行。行。 t=akj; akj=amaxij; amaxij=t; for(i=k+1;i=n-1;i+) 行乘子行乘子aik/=-akk。 for(j=k+1;j=0;k-) （回代過程）（回代過程） s=0; for(j=k+1;j=n-1;j+) s+=akj*xj; xk=(akn-

17、s)/akk; 輸出原方程組的解輸出原方程組的解xn。3.3 列主元高斯消元法列主元高斯消元法列主元高斯消元法對應(yīng)的程序（列主元高斯消元法對應(yīng)的程序（1/2）#include #include #define MAXSIZE 50void input(double aMAXSIZEMAXSIZE+1,long n);void output(double xMAXSIZE,long n);void main(void)double aMAXSIZEMAXSIZE+1,xMAXSIZE,s,max,t;long n,i,j,k,maxi;printf(n請輸入原方程組的階數(shù)：請輸入原方程組的階數(shù)：

18、);scanf(%ld,&n);input(a,n);for(k=0;k=n-2;k+)max=akk;maxi=k;for(i=k+1;ifabs(max)max=aik;maxi=i;if(max=0)break;if(maxi!=k)for(j=k;j=n;j+)t=akj;akj=amaxij;amaxij=t; for(i=k+1;i=n-1;i+)aik/=-akk;for(j=k+1;j=0;k-)s=0;for(j=k+1;j=n-1;j+)s+=akj*xj;xk=(akn-s)/akk; output(x,n);void input(double aMAXSIZEM

19、AXSIZE+1,long n)long i,j;printf(n請輸入原方程組的增廣矩陣：請輸入原方程組的增廣矩陣：n);for(i=1;i=n;i+)for(j=1;j=n+1;j+)scanf(%lf,&ai-1j-1);void output(double xMAXSIZE,long n)long k;printf(n原方程組的解為：原方程組的解為：n);for(k=1;k=n;k+)printf( %lf,xk-1);3.4 全主元高斯消元法全主元高斯消元法全主元高斯消元法是另一種選主元高斯消元法。全主元高斯消元法是另一種選主元高斯消元法。 全主元高斯消元法的主要思想全主元高

20、斯消元法的主要思想全主元高斯消元法與順序高斯消元法和列主元高斯消元法相比，它的全主元高斯消元法與順序高斯消元法和列主元高斯消元法相比，它的改進之處是：在用主元素消去第改進之處是：在用主元素消去第k1行至第行至第n行的變元行的變元xk之前，從主元之前，從主元素右下方的矩形區(qū)域所有系數(shù)中，尋找最大者，設(shè)素右下方的矩形區(qū)域所有系數(shù)中，尋找最大者，設(shè)最大的系數(shù)在最大的系數(shù)在maxi 行行maxj列列，若，若maxik，則交換第，則交換第k行和第行和第maxi行；若行；若maxjk，則交換，則交換第第k列和第列和第maxj列，使最大的系數(shù)成為主元。列，使最大的系數(shù)成為主元。 在程序中不存儲變元的名稱，確

21、定變元的依據(jù)是變元的位置，即變元在程序中不存儲變元的名稱，確定變元的依據(jù)是變元的位置，即變元位于哪一列。因此交換列時，需要記錄交換后變元的位置（位于哪一位于哪一列。因此交換列時，需要記錄交換后變元的位置（位于哪一列），求解完畢，按照原方程組各變元的次序輸出各變元的值。列），求解完畢，按照原方程組各變元的次序輸出各變元的值。 當線性方程組當線性方程組 有唯一解時，全主元高斯消元法必能求出解。有唯一解時，全主元高斯消元法必能求出解。與列主元高斯消元法相比，全主元高斯消元法也能夠避免小主元時舍與列主元高斯消元法相比，全主元高斯消元法也能夠避免小主元時舍入誤差被放大的情況，而且解的精度一般比列主元高斯

22、消元法高。入誤差被放大的情況，而且解的精度一般比列主元高斯消元法高。 總之，與全主元高斯消元法增加的運算量相比，列主元高斯消元法增總之，與全主元高斯消元法增加的運算量相比，列主元高斯消元法增加的運算量可以忽略。全主元高斯消元法的程序結(jié)構(gòu)比列主元高斯消加的運算量可以忽略。全主元高斯消元法的程序結(jié)構(gòu)比列主元高斯消元法要復(fù)雜得多。只要原方程組有解，這元法要復(fù)雜得多。只要原方程組有解，這2種方法都能求出解。一般而種方法都能求出解。一般而言，全主元高斯消元法比列主元高斯消元法精度高。言，全主元高斯消元法比列主元高斯消元法精度高。 3.5 高斯約當消元法高斯約當消元法與順序高斯消元法的求解過程相比，高斯約

23、當消元法求解過程與之不與順序高斯消元法的求解過程相比，高斯約當消元法求解過程與之不同的地方有：同的地方有： 在消去變元在消去變元xk那一列時，不僅把主元之下的那一列時，不僅把主元之下的xk消去，也把主元之上消去，也把主元之上的的xk消去。也就是把系數(shù)矩陣除主對角線之外的元素都化為消去。也就是把系數(shù)矩陣除主對角線之外的元素都化為0。 在消去變元在消去變元xk那一列時，先把主元化為那一列時，先把主元化為1，再用它去消其它方程的，再用它去消其它方程的變元變元xk。也就是把系數(shù)矩陣主對角線上的元素都化為。也就是把系數(shù)矩陣主對角線上的元素都化為1。 一、高斯約當消元法的主要思想一、高斯約當消元法的主要思

24、想消元完畢，原線性方程組化為下面的形式：消元完畢，原線性方程組化為下面的形式： )(1,)(,)2(1, 22)2(22)1(1, 11)1(11nnnnnnnnnaxaaxaaxa這時系數(shù)矩陣這時系數(shù)矩陣 被化為單位陣，右端向量就是解向量，不再需要回代過被化為單位陣，右端向量就是解向量，不再需要回代過程。因此高斯約當消元法又稱為無回代過程消元法。程。因此高斯約當消元法又稱為無回代過程消元法。3.5 高斯約當消元法高斯約當消元法二、高斯約當消元法的算法二、高斯約當消元法的算法輸入原方程組的階數(shù)輸入原方程組的階數(shù)n。輸入原方程組的增廣矩陣輸入原方程組的增廣矩陣ann+1。for(k=0;k=n-

25、1;k+)，循環(huán)，循環(huán)1次消去次消去1列。列。 for(j=k+1;j=n;j+)，把主元化為，把主元化為1。 akj/=akk。 for(i=0;i=n-1;i+)，循環(huán)，循環(huán)1次更新次更新1行。行。 ik Y N for(j=k+1;j=n;j+) aij+=(-aik)*akj;輸出原方程組的解（此時右端向量就是解向量）。輸出原方程組的解（此時右端向量就是解向量）。3.5 高斯約當消元法高斯約當消元法高斯約當消元法對應(yīng)的程序（高斯約當消元法對應(yīng)的程序（1/2）#include #include #define MAXSIZE 50void input(double aMAXSIZEMAX

26、SIZE+1,long n);void output(double aMAXSIZEMAXSIZE+1,long n);void main(void)double aMAXSIZEMAXSIZE+1;long n,i,j,k;printf(n請輸入原方程組的階數(shù)：請輸入原方程組的階數(shù)：);scanf(%ld,&n);input(a,n);for(k=0;k=n-1;k+)for(j=k+1;j=n;j+)akj/=akk;for(i=0;i=n-1;i+)if(i!=k)for(j=k+1;j=n;j+)aij+=(-aik)*akj;output(a,n);3.5 高斯約當消元法高斯

27、約當消元法高斯約當高斯約當消元法對應(yīng)的程序（消元法對應(yīng)的程序（2/2）void input(double aMAXSIZEMAXSIZE+1,long n)long i,j;printf(n請輸入原方程組的增廣矩陣：請輸入原方程組的增廣矩陣：n);for(i=1;i=n;i+)for(j=1;j=n+1;j+)scanf(%lf,&ai-1j-1);void output(double aMAXSIZEMAXSIZE+1,long n)long k;printf(n原方程組的解為：原方程組的解為：n);for(k=1;k=n;k+)printf( %lf,ak-1n);3.6 消元形式

28、的追趕法消元形式的追趕法帶狀對角矩陣是一種常見的高階稀疏矩陣。帶狀對角矩陣是一種常見的高階稀疏矩陣。3對角陣是帶狀對角矩陣的對角陣是帶狀對角矩陣的一種。在做三次樣條插值、求解微分方程等問題時，經(jīng)常需要解一種。在做三次樣條插值、求解微分方程等問題時，經(jīng)常需要解3對角對角方程組。方程組。 一、消元形式的追趕法的主要思想一、消元形式的追趕法的主要思想定義：如果在方陣定義：如果在方陣 中，除主對角線和中，除主對角線和2條次對角線之外的元素都為條次對角線之外的元素都為0，那么那么 是是3對角陣。即對角陣。即 可以表示為下面的形式：可以表示為下面的形式：AAA=nnnnncacbacbacb1112221

29、1A定義：如果線性方程組定義：如果線性方程組 的系數(shù)矩陣的系數(shù)矩陣 是是3對角陣，那么對角陣，那么 是是3對角方程組。對角方程組。 bAx bAx A3.6 消元形式的追趕法消元形式的追趕法追趕法用于求解追趕法用于求解3對角方程組。在求解對角方程組。在求解3對角方程組時，原來的求解過程對角方程組時，原來的求解過程被大大簡化。要求解的被大大簡化。要求解的3對角方程組的一般形式為：對角方程組的一般形式為：nnnnnnnnnnnndxbxadxcxbxadxcxbxadxcxb1111121232221212111這里按照順序高斯消元法的求解思路來求解這個這里按照順序高斯消元法的求解思路來求解這個3

30、對角方程組。對角方程組。在消去在消去xi一列時，行乘子一列時，行乘子 = 。更新第更新第i1行中變元系數(shù)的公式為：行中變元系數(shù)的公式為：bi+1=bi+1ci ，di+1=di+1di其中其中i=1,2,n1。 iliiba1ilil3.6 消元形式的追趕法消元形式的追趕法消元完畢，原方程組化為下面的形式：消元完畢，原方程組化為下面的形式： nnnnnnnndxbdxcxbdxcxbdxcxb11112322212111回代過程，仍然是自下而上求解。回代的公式為：回代過程，仍然是自下而上求解?；卮墓綖椋簒n=dn/bnxi=(dicixi+1)/bi 其中其中i=n1,n2,1。 采用這一

31、思路，可以求解帶狀對角線性方程組，運算量和存儲需求也會采用這一思路，可以求解帶狀對角線性方程組，運算量和存儲需求也會大大簡化。大大簡化。 3.6 消元形式的追趕法消元形式的追趕法二、二、消元形式的追趕法的消元形式的追趕法的算法算法輸入原方程組的階數(shù)輸入原方程組的階數(shù)n。輸入增廣矩陣對應(yīng)的數(shù)組輸入增廣矩陣對應(yīng)的數(shù)組a,b,c,d。for(i=1;i=0;i-) xi=(di-ci*xi+1)/bi;輸出原方程組的解輸出原方程組的解xn。3.6 消元形式的追趕法消元形式的追趕法消元形式的追趕法消元形式的追趕法的程序（的程序（1/2）#include #include #define MAXSIZE

32、 50void input(double a,double b,double c,double d,long n);void output(double x,long n);void main(void)double aMAXSIZE,bMAXSIZE,cMAXSIZE,dMAXSIZE,xMAXSIZE;long n,i;printf(n請輸入原方程組的階數(shù)：請輸入原方程組的階數(shù)：);scanf(%ld,&n);input(a,b,c,d,n);for(i=1;i=0;i-)xi=(di-ci*xi+1)/bi;output(x,n);3.6 消元形式的追趕法消元形式的追趕法消元形式

33、的追趕法的程序消元形式的追趕法的程序 （2/2）void input(double a,double b,double c,double d,long n)long i;printf(n請輸入原方程組的增廣矩陣：請輸入原方程組的增廣矩陣：n);printf(nb1,c1,d1:);scanf(%lf,%lf,%lf,&b0,&c0,&d0);for(i=2;i=n-1;i+)printf(na%ld,b%ld,c%ld,d%ld:,i,i,i,i);scanf(%lf,%lf,%lf,%lf,&ai-1,&bi-1,&ci-1,&di-1

34、);printf(na%ld,b%ld,d%ld:,n,n,n);scanf(%lf,%lf,%lf,&an-1,&bn-1,&dn-1);void output(double x,long n)long i;printf(n方程組的解為：方程組的解為：n);for(i=1;i=n;i+)printf( %lf,xi-1);3.7 LU分解法分解法三角分解法求解線性方程組的方法，是指把線性方程組的系數(shù)矩陣進行三角分解法求解線性方程組的方法，是指把線性方程組的系數(shù)矩陣進行三角分解，然后按照線性方程組的矩陣形式進行回代求解的方法。三角分解，然后按照線性方程組的矩陣形式進行回

35、代求解的方法。 三角分解法主要有三角分解法主要有LU分解法（又稱直接三角分解法、分解法（又稱直接三角分解法、Doolittle分解分解法）、法）、LDR分解法、分解法、Crout分解法等。分解法等。 LU分解法把系數(shù)矩陣分解為單位下三角陣分解法把系數(shù)矩陣分解為單位下三角陣L和上三角陣和上三角陣U的積；的積；LDR分解法把系數(shù)矩陣分解為單位下三角陣分解法把系數(shù)矩陣分解為單位下三角陣L、對角陣、對角陣D和單位上三和單位上三角陣角陣R的積，其中的積，其中DR=U；Crout分解法把系數(shù)矩陣分解為下三角陣分解法把系數(shù)矩陣分解為下三角陣 和單位上三角陣和單位上三角陣R的積，的積，其中其中 =LD。 LL

36、除此之外，還有求解對稱正定線性方程組的除此之外，還有求解對稱正定線性方程組的LLT分解法（又稱平方根分解法（又稱平方根法），求解帶狀對角線性方程組的矩陣形式的追趕法等方法。法），求解帶狀對角線性方程組的矩陣形式的追趕法等方法。 3.7 LU分解法分解法對任意矩陣對任意矩陣Amn做初等變換的結(jié)果，等于做初等變換的結(jié)果，等于Amn與對應(yīng)初等方陣的積。與對應(yīng)初等方陣的積。 Em(i（k）,j)= Em(i（k）)= Em(i,j)= 由Em得到的初等方陣10110111k1111k數(shù)k乘以Em的j列，加到i列上去，得到Em(i（k）,j)。Em(i（k）,j)右乘Bnm的積，等于數(shù)k乘以Bnm的j列

37、，加到i列上去的結(jié)果。數(shù)k乘以Em的i行，加到j(luò)行上去，得到Em(i（k）,j)。Em(i（k）,j)左乘Amn的積，等于數(shù)k乘以Amn的i行，加到j(luò)行上去的結(jié)果。數(shù)k乘以Em的i列得到Em(i（k）)。Em(i（k）)右乘Bnm的積，等于數(shù)k乘以Bnm的i列的結(jié)果。數(shù)k乘以Em的i行得到Em(i（k）)。Em(i（k）)左乘Amn的積，等于數(shù)k乘以Amn的i行的結(jié)果。Em交換i列和j列得到Em(i,j)。Em(i,j)右乘Bnm的積，等于交換Bnm的i列和j列的結(jié)果。Em交換i行和j行得到Em(i,j)。Em(i,j)左乘Amn的積，等于交換Amn的i行和j行的結(jié)果。對應(yīng)初等方陣右乘Bnm對

38、應(yīng)初等方陣左乘Amn一、相關(guān)的初等方陣性質(zhì)一、相關(guān)的初等方陣性質(zhì)3.7 LU分解法分解法方陣的方陣的LU分解是指把方陣分解為單位下三角陣分解是指把方陣分解為單位下三角陣L和上三角陣和上三角陣U的積過程。的積過程。 二、二、LU分解與順序高斯消元的聯(lián)系分解與順序高斯消元的聯(lián)系定理：若方陣存在定理：若方陣存在LU分解，那么分解是唯一的。分解，那么分解是唯一的。 用順序高斯消元法求解線性方程組，消元過程對應(yīng)于對增廣矩陣進行一用順序高斯消元法求解線性方程組，消元過程對應(yīng)于對增廣矩陣進行一系列倍加變換，直到系數(shù)矩陣化為上三角陣為止。系列倍加變換，直到系數(shù)矩陣化為上三角陣為止。 定理：如果系數(shù)矩陣定理：如

39、果系數(shù)矩陣A存在存在LU分解分解A=LU，那么，那么U就是消元后的系數(shù)矩陣就是消元后的系數(shù)矩陣A，即：即：U=，L= )(,)2(, 2)2(22)1(, 1)1(12)1(11nnnnnaaaaaa111113121nlll1111232nll11111,nnl上式中上式中l(wèi)i,k= ，(k1,n1,ik1,n)為消元過程中全部的行乘子，為消元過程中全部的行乘子，)(,)(,kkkkkiaa3.7 LU分解法分解法對方陣對方陣A進行進行LU分解分解A=LU的通項公式為：的通項公式為： 三、對方陣進行三、對方陣進行LU分解的過程分解的過程= ，（j=i,i+1,n） ijuija11ikkji

40、kul=( )/ ，（i=j+1,j+2,n） ijlija11jkkjikuljju在用通項公式對方陣在用通項公式對方陣A進行進行LU分解分解時，一個較好的計算次序是：時，一個較好的計算次序是： U第第1行，行， L第第1列，列， U第第2行，行， L第第2列，列， U第第3行，行， L第第3列，列，3.7 LU分解法分解法 對系數(shù)矩陣對系數(shù)矩陣A進行進行LU分解，得到分解，得到L和和U。四、四、LU分解法求解線性方程組的過程分解法求解線性方程組的過程方程組方程組Ax=b(LU)x=bL(Ux)=b。 令向量令向量y=Ux，代入上式得：，代入上式得：Ly=b?；卮?，求解回代，求解Ly=b，得

41、到，得到y(tǒng)。通項公式為：通項公式為：yi=bi yk，i=1,2,n。 11ikikl 回代，求解回代，求解Ux=y，得到解向量，得到解向量x。 通項公式為：通項公式為：xi=(yi xk)/uii，i=n,n-1,1。 nikiku13.7 LU分解法分解法五、五、LU分解法的算法。分解法的算法。輸入原方程組的階數(shù)n。輸入原方程組的系數(shù)矩陣ann，右端向量xn。for(k=0;k=n-2;k+) for(i=k+1;i=n-1;i+)，循環(huán)1次求出L的1個元素。 s=0; for(j=0;j=k-1;j+) s+=aij*ajk; aik=(aik-s)/akk; for(j=k+1;j=n

42、-1;j+) s=0; for(i=0;i=k;i+) s+=ak+1i*aij; ak+1j-=s;for(i=1;i=n-1;i+) s=0; for(j=0;j=0;i-) s=0; for(j=i+1;j=n-1;j+) s+=aij*xj; xi=(xi-s)/aii;輸出原方程組的解xn。3.7 LU分解法分解法LU分解法分解法的程序（的程序（1/2）#include #include #define MAXSIZE 50void input(double aMAXSIZEMAXSIZE,double xMAXSIZE,long n);void output(double xMAX

43、SIZE,long n);void main(void)double aMAXSIZEMAXSIZE,xMAXSIZE,s;long n,i,j,k;printf(n請輸入原方程組的階數(shù)：請輸入原方程組的階數(shù)：);scanf(%ld,&n);input(a,x,n);for(k=0;k=n-2;k+)for(i=k+1;i=n-1;i+)s=0;for(j=0;j=k-1;j+)s+=aij*ajk;aik=(aik-s)/akk;for(j=k+1;j=n-1;j+)s=0;for(i=0;i=k;i+)s+=ak+1i*aij;ak+1j-=s;3.7 LU分解法分解法LU分解法分

44、解法的程序的程序 （2/2）for(i=1;i=n-1;i+)s=0;for(j=0;j=0;i-)s=0;for(j=i+1;j=n-1;j+)s+=aij*xj;xi=(xi-s)/aii;output(x,n);void input(double aMAXSIZEMAXSIZE,double xMAXSIZE,long n)long i,j;printf(n請輸入原方程組的增廣矩陣：請輸入原方程組的增廣矩陣：n);for(i=0;i=n-1;i+)for(j=0;j=n-1;j+)scanf(%lf,&aij);scanf(%lf,&xi);void output(dou

45、ble xMAXSIZE,long n)long i;printf(n原方程組的解為：原方程組的解為：n);for(i=0;i=n-1;i+)printf( %lf,xi);3.8 矩陣形式的追趕法矩陣形式的追趕法用三角分解法求解用三角分解法求解3對角方程組，求解過程同樣要大大簡化，這就是追對角方程組，求解過程同樣要大大簡化，這就是追趕法的矩陣形式。趕法的矩陣形式。 3對角方程組對角方程組Ax=d的矩陣形式為：的矩陣形式為：nnnnncacbacbacb11122211 = nxxx21nddd21用矩陣形式的追趕法求解用矩陣形式的追趕法求解3對角方程組，需要對系數(shù)矩陣進行三角分解。對角方程組

46、，需要對系數(shù)矩陣進行三角分解。三角分解的方法同樣分為三角分解的方法同樣分為LU分解法、分解法、LDR分解法、分解法、Crout分解法。用這分解法。用這3種方法構(gòu)造的追趕法的原理、求解過程和能求出解的條件相似。這里種方法構(gòu)造的追趕法的原理、求解過程和能求出解的條件相似。這里以以Crout分解法為例，介紹矩陣形式的追趕法。分解法為例，介紹矩陣形式的追趕法。Crout分解的存在、唯一分解的存在、唯一性同性同LU分解。分解。 3.8 矩陣形式的追趕法矩陣形式的追趕法一、一、3對角陣對角陣Crout分解的過程分解的過程存儲存儲3對角方程組增廣矩陣可以用一維數(shù)組對角方程組增廣矩陣可以用一維數(shù)組an、bn、

47、cn和和dn，不存，不存儲系數(shù)矩陣中除儲系數(shù)矩陣中除3條對角線之外的零元素，以節(jié)省存儲存儲空間。條對角線之外的零元素，以節(jié)省存儲存儲空間。 Crout分解法把系數(shù)矩陣分解法把系數(shù)矩陣A分解為下三角陣分解為下三角陣 和單位上三角陣和單位上三角陣R的積，的積，即即A= R： LLnnnnncacbacbacb11122211=nnnnlalalal1122111111121nrrr對對3對角系數(shù)矩陣進行對角系數(shù)矩陣進行Crout分解的通項公式為：分解的通項公式為： 的左下次對角線元素與的左下次對角線元素對應(yīng)相等。的左下次對角線元素與的左下次對角線元素對應(yīng)相等。 l1=b1， li=bi-airi-

48、1，（，（i=2,3,n） ri=ci/li，（，（i=1,2,n-1）在用上述通項公式對在用上述通項公式對3對角系數(shù)矩陣進行對角系數(shù)矩陣進行Crout分解時，合邏輯的計算次分解時，合邏輯的計算次序是：序是：l1,r1,l2,r2,ln-1,rn-1,ln。按照這一次序計算，能夠保證在用到某一。按照這一次序計算，能夠保證在用到某一元素之前，這一元素已被求出。元素之前，這一元素已被求出。 3.8 矩陣形式的追趕法矩陣形式的追趕法用矩陣形式的追趕法求解用矩陣形式的追趕法求解3對角線性方程組對角線性方程組Ax=d的步驟如下：的步驟如下： 二、二、矩陣形式的追趕法的求解步驟矩陣形式的追趕法的求解步驟

49、對對3對角系數(shù)矩陣對角系數(shù)矩陣A進行進行Crout分解，得到下三角陣分解，得到下三角陣 和單位上三角陣和單位上三角陣R。 L方程組方程組Ax=d( R)x=d (Rx)=d。 LL 令向量令向量y=Rx，代入上式得：，代入上式得： y=d。 L回代，求解回代，求解 y=d，得到，得到y(tǒng)。 L通項公式為：通項公式為：y1=d1/l1，yi=(di-aiyi-1)/li，i=2,3,n。 回代，求解回代，求解Rx=y，得到解向量，得到解向量x。 通項公式為：通項公式為：xn=yn，xi=yi-riyi+1，i=n-1,n-2,1。步驟步驟的回代過程自上而下進行，下標的回代過程自上而下進行，下標 ，

50、稱為，稱為“追追”，步驟，步驟的回代過的回代過程自下而上進行，下標程自下而上進行，下標 ，稱為，稱為“趕趕”，因此這一求解方法稱為，因此這一求解方法稱為“追趕追趕法法”。 3.8 矩陣形式的追趕法矩陣形式的追趕法三、矩陣形式追趕法的算法。三、矩陣形式追趕法的算法。輸入原方程組的階數(shù)輸入原方程組的階數(shù)n。輸入增廣矩陣對應(yīng)的數(shù)組輸入增廣矩陣對應(yīng)的數(shù)組a,b,c,d。for(i=0;i=n-2;i+) （Crout分解）分解） ci/=bi; bi+1-=ai+1*ci;x0=d0/b0; （回代，求解）（回代，求解）for(i=1;i=0;i-) （回代，求解）（回代，求解） xi-=ci*xi+

51、1;輸出原方程組的解輸出原方程組的解xn。3.8 矩陣形式的追趕法矩陣形式的追趕法矩陣形式追趕法的程序（矩陣形式追趕法的程序（1/2）#include #include #define MAXSIZE 50void input(double a,double b,double c,double d,long n);void output(double x,long n);void main(void)double aMAXSIZE,bMAXSIZE,cMAXSIZE,dMAXSIZE,xMAXSIZE;long n,i;printf(n請輸入原方程組的階數(shù)：請輸入原方程組的階數(shù)：);scanf

52、(%ld,&n);input(a,b,c,d,n);for(i=0;i=n-2;i+)ci/=bi;bi+1-=ai+1*ci;x0=d0/b0;for(i=1;i=0;i-)xi-=ci*xi+1;output(x,n);3.8 矩陣形式的追趕法矩陣形式的追趕法矩陣形式追趕法矩陣形式追趕法的程序的程序 （2/2）void input(double a,double b,double c,double d,long n)long i;printf(n請輸入原方程組的增廣矩陣：請輸入原方程組的增廣矩陣：n);printf(nb1,c1,d1:);scanf(%lf,%lf,%lf,&am

53、p;b0,&c0,&d0);for(i=2;i=n-1;i+)printf(na%ld,b%ld,c%ld,d%ld:,i,i,i,i);scanf(%lf,%lf,%lf,%lf,&ai-1,&bi-1,&ci-1,&di-1);printf(na%ld,b%ld,d%ld:,n,n,n);scanf(%lf,%lf,%lf,&an-1,&bn-1,&dn-1);void output(double x,long n)long i;printf(n方程組的解為：方程組的解為：n);for(i=1;i=n;i+)print

54、f( %lf,xi-1);3.9 平方根法平方根法一、基礎(chǔ)知識一、基礎(chǔ)知識平方根法（又稱平方根法（又稱LLT分解法）用來求解對稱正定線性方程組。分解法）用來求解對稱正定線性方程組。 定義：如果線性方程組定義：如果線性方程組Ax=b的系數(shù)矩陣的系數(shù)矩陣A是對稱正定矩陣，那么是對稱正定矩陣，那么Ax=b是是對稱正定線性方程組。對稱正定線性方程組。 定義：如果方陣定義：如果方陣A滿足滿足A=AT，那么，那么A是對稱陣。是對稱陣。 定理：定理：(AB)T=BTAT定義：如果定義：如果n行行n列的實對稱陣列的實對稱陣A，對任意非零，對任意非零n維實向量維實向量x，恒有，恒有xTAx0，則稱則稱A為對稱正

55、定矩陣。為對稱正定矩陣。 定理：定理：判定實對稱陣判定實對稱陣A A是對稱正定矩陣的充要條件為：是對稱正定矩陣的充要條件為：A A的各階順序主子式的各階順序主子式都大于都大于0。 定理：定理：判定實對稱陣判定實對稱陣A A是對稱正定矩陣的充要條件為：存在可逆實矩陣是對稱正定矩陣的充要條件為：存在可逆實矩陣P，使使A A=PTP。 定理：如果定理：如果A是對稱正定矩陣，那么是對稱正定矩陣，那么A滿足：滿足： A的主對角線元素都大于的主對角線元素都大于0。 A的順序主子式正定。的順序主子式正定。 A-1正定。正定。 3.9 平方根法平方根法一、基礎(chǔ)知識（續(xù)）一、基礎(chǔ)知識（續(xù)）定理：如果定理：如果A是對稱正定矩陣，那么是對稱正定矩陣，那么A滿足：滿足： A必定存在必定存在LLT分解，分解，A= ，其中，其中 為下三角陣。為下三角陣。 的主對角線元素都大于的主對角線元素都大于0。 A的的LLT分解唯一。分