第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 阶数在阶数在100150以上的线性方程组为高阶线性方程组以上的线

2、性方程组为高阶线性方程组阶数在阶数在100150以下的线性方程组为低阶线性方程组以下的线性方程组为低阶线性方程组实际应用中，经常见到高阶稀疏线性方程组，低阶稠密线性方程组，实际应用中，经常见到高阶稀疏线性方程组，低阶稠密线性方程组，对称正定线性方程组，带状线性方程组等等。对称正定线性方程组，带状线性方程组等等。系数矩阵的大部分元素为零元素的线性方程组为稀疏线性方程组系数矩阵的大部分元素为零元素的线性方程组为稀疏线性方程组系数矩阵的大部分元素为非零元素的线性方程组为稠密线性方程组系数矩阵的大部分元素为非零元素的线性方程组为稠密线性方程组克莱姆法则并不实用。常用的数值解法主要分为两类：克莱姆法则并

3、不实用。常用的数值解法主要分为两类： 直接求解方法是指经过有限次四则运算，求出线性方程组精确解的方法。直接求解方法是指经过有限次四则运算，求出线性方程组精确解的方法。 迭代求解方法是指构造一种迭代方法，由某个（套）迭代初值（粗略解），迭代求解方法是指构造一种迭代方法，由某个（套）迭代初值（粗略解），得到近似解序列，用序列极限逐步逼近线性方程组精确解的方法。得到近似解序列，用序列极限逐步逼近线性方程组精确解的方法。用顺序高斯消元法求解线性方程组的过程分为消元过程和回代过程。消元用顺序高斯消元法求解线性方程组的过程分为消元过程和回代过程。消元过程是自上而下，把原线性方程组化为上三角方程组的过程；回

4、代过程是过程是自上而下，把原线性方程组化为上三角方程组的过程；回代过程是自下而上，对这个上三角方程组进行求解的过程。自下而上，对这个上三角方程组进行求解的过程。 3.2 顺序高斯消元法顺序高斯消元法消元的过程是对线性方程组做同解变换，把原线性方程组变为上三角方程消元的过程是对线性方程组做同解变换，把原线性方程组变为上三角方程组的过程。组的过程。 消元过程消元过程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外层循环每循环一轮，就用涉及到的矩形区域中上面的第外层循环每循环一轮，就用涉及到的矩形区域中上面的第1行消去它下面行消去它下面各行左面的第各行左面的第1列，此矩形区域最上列，此矩形区域最上1行和最左行和最左1列在之后的消元过程中不列在之后的消元过程中不再涉及（上标不再改变），除此之外的右下矩形区域中所有的系数被更新再涉及（上标不再改变），除此之外的右下矩形区域中所有的系数被更新一遍，且上标一遍，且上标k自增自增1。第。第k次消元的初始状态如左图所示。次消元的初始状态如左图所示。 3.2 顺序高斯消元法顺序高斯消元法在第在第k次更新

6、的初始，涉及到的待更新矩形区域中，系数的上标都为次更新的初始，涉及到的待更新矩形区域中，系数的上标都为k，此矩形区域中最左上角的系数的行标、列标也为此矩形区域中最左上角的系数的行标、列标也为k。更新的过程，就是用。更新的过程，就是用第第k行消去第行消去第k1行至第行至第n行的第行的第k列，并且让第列，并且让第k1行至第行至第n行，第行，第k1列至第列至第n1列的矩形区域的系数更新，上标变为列的矩形区域的系数更新，上标变为k1。在此后的消元过。在此后的消元过程中，第程中，第k行及之上和第行及之上和第k列及其左边的系数不再变化。列及其左边的系数不再变化。)(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消元结束之后得到的上三角方程组中，每向下一行，上标消元结束之后得到的上三角方程组中，每向下一行，上标k增增1，即上标，即上标k=行标行标i，如右图所示。，如右图所示。3.2 顺序高斯消元法顺序高斯消元法顺序高斯消元法的

8、算法顺序高斯消元法的算法输入原方程组的阶数输入原方程组的阶数n。输入原方程组的增广矩阵输入原方程组的增广矩阵ann+1。for(k=0;k=n-2;k+)，循环，循环1次消去次消去1列。列。 （消元过程）（消元过程） for(i=k+1;i=n-1;i+)，循环，循环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 顺序高斯消元法顺序高斯消元法顺序高斯消元法对应的程序（顺序高斯

9、消元法对应的程序（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请输入原方程组的阶数：请输入原方程组的阶数：);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 顺序高斯消元法顺序高斯消元法顺序高斯消元法对应的程序（顺序高斯消元法对应的程序（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 顺序高斯消元法顺序高斯消元法当方程组的阶数增长时，顺序高斯消元法的运算量增长速度，比克莱姆法当方程组的阶数增长时，顺序高斯消元法的运算量增长速度，比克莱姆法则求解线性方程组的运算量增长速度要慢得多。则求解线性方程组的运算量增长速度要慢得多。 顺序高斯消元法的乘除次数顺序高斯消元法的乘除次数MD 。 33n由克莱姆法则，可知当线性方程组的系数行列式由克莱姆法则，可知当线性方程组的系数行列式不等于不等于0时时，则线性方程，则线性方程组有唯一解，否

12、则线性方程组可能无解（如组有唯一解，否则线性方程组可能无解（如0 x0y=1），或有无穷多解），或有无穷多解（如（如0 x0y=0）。）。 性质：当线性方程组有唯一解时，要想让顺序高斯消元法能求出解，需要性质：当线性方程组有唯一解时，要想让顺序高斯消元法能求出解，需要 （k=1,2,n）都不为）都不为0。 )(,kkka定义：定义： 对求解起突出作用，称为主元素。对求解起突出作用，称为主元素。 )(,kkka定理：对定理：对n阶线性方程组阶线性方程组 ，顺序高斯消元法能求出解的充要条件是：，顺序高斯消元法能求出解的充要条件是：系数矩阵系数矩阵 的所有顺序主子式的所有顺序主子式i(i=1,2,n

13、)均不为均不为0。 bAx A3.3 列主元高斯消元法列主元高斯消元法顺序高斯消元法的缺陷有：顺序高斯消元法的缺陷有： 当线性方程组的系数行列式不等于当线性方程组的系数行列式不等于0，且存在，且存在1个主元素等于个主元素等于0时，时，线性方程组有唯一解，但顺序高斯消元法不能求出解。线性方程组有唯一解，但顺序高斯消元法不能求出解。 交换增广矩阵的交换增广矩阵的2行（行（2个方程的位置互换）或交换个方程的位置互换）或交换2列（列（2个变元的位置互个变元的位置互换），得到的方程组与原方程组同解。通过交换行、列，来避免主元为换），得到的方程组与原方程组同解。通过交换行、列，来避免主元为0和小主元，然后

14、再消元求解的方法，称为选主元高斯消元法。和小主元，然后再消元求解的方法，称为选主元高斯消元法。 一、列主元高斯消元法的主要思想一、列主元高斯消元法的主要思想列主元高斯消元法是一种常见的选主元高斯消元法，与顺序高斯消元法相列主元高斯消元法是一种常见的选主元高斯消元法，与顺序高斯消元法相比，它的改进之处是：在用主元素消去第比，它的改进之处是：在用主元素消去第k1行至第行至第n行的变元行的变元xk之前，之前，从第从第k行至第行至第n行所有的行所有的xk系数中，寻找最大系数中，寻找最大的系数的系数，若，若maxik，则交换，则交换第第k行和第行和第maxi行，使最大的系数成为主元。行，使最大的系数成为

15、主元。 当线性方程组有唯一解时，列主元高斯消元法必能求出解。当线性方程组有唯一解时，列主元高斯消元法必能求出解。 列主元高斯消元法能避免小主元时舍入误差被放大的情况，提高解的精度。列主元高斯消元法能避免小主元时舍入误差被放大的情况，提高解的精度。 小主元（主元素远小于其他元素）时能求出解，但求出的解误差很大。小主元（主元素远小于其他元素）时能求出解，但求出的解误差很大。 列列主主元元高高斯斯消消元元法法的的算算法法输入原方程组的阶数输入原方程组的阶数n。输入原方程组的增广矩阵输入原方程组的增广矩阵ann+1。for(k=0;k=n-2;k+)，循环，循环1次消去次消去1列。列。 （消元过程）（

16、消元过程） 暂设最大系数暂设最大系数max=akk，最大系数所在行的行号，最大系数所在行的行号maxi=k。 for(i=k+1;i|max| Y N 令令max=aik;maxi=i; 系数最大值系数最大值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 列主元高斯消元法列主元高斯消元法列主元高斯消元法对应的程序（列主元高斯消元法对应的程序（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请输入原方程组的阶数：请输入原方程组的阶数：

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之前，从主元之前，从主元素右下方的矩形区域所有系数中，寻找最大者，设素右下方的矩形区域所有系数中，寻找最大者，设最大的系数在最大的系数在maxi 行行maxj列列，若，若maxik，则交换第，则交换第k行和第行和第maxi行；若行；若maxjk，则交换，则交换第第k列和第列和第maxj列，使最大的系数成为主元。列，使最大的系数成为主元。 在程序中不存储变元的名称，确

21、定变元的依据是变元的位置，即变元在程序中不存储变元的名称，确定变元的依据是变元的位置，即变元位于哪一列。因此交换列时，需要记录交换后变元的位置（位于哪一位于哪一列。因此交换列时，需要记录交换后变元的位置（位于哪一列），求解完毕，按照原方程组各变元的次序输出各变元的值。列），求解完毕，按照原方程组各变元的次序输出各变元的值。 当线性方程组当线性方程组 有唯一解时，全主元高斯消元法必能求出解。有唯一解时，全主元高斯消元法必能求出解。与列主元高斯消元法相比，全主元高斯消元法也能够避免小主元时舍与列主元高斯消元法相比，全主元高斯消元法也能够避免小主元时舍入误差被放大的情况，而且解的精度一般比列主元高斯

22、消元法高。入误差被放大的情况，而且解的精度一般比列主元高斯消元法高。 总之，与全主元高斯消元法增加的运算量相比，列主元高斯消元法增总之，与全主元高斯消元法增加的运算量相比，列主元高斯消元法增加的运算量可以忽略。全主元高斯消元法的程序结构比列主元高斯消加的运算量可以忽略。全主元高斯消元法的程序结构比列主元高斯消元法要复杂得多。只要原方程组有解，这元法要复杂得多。只要原方程组有解，这2种方法都能求出解。一般而种方法都能求出解。一般而言，全主元高斯消元法比列主元高斯消元法精度高。言，全主元高斯消元法比列主元高斯消元法精度高。 3.5 高斯约当消元法高斯约当消元法与顺序高斯消元法的求解过程相比，高斯约

23、当消元法求解过程与之不与顺序高斯消元法的求解过程相比，高斯约当消元法求解过程与之不同的地方有：同的地方有： 在消去变元在消去变元xk那一列时，不仅把主元之下的那一列时，不仅把主元之下的xk消去，也把主元之上消去，也把主元之上的的xk消去。也就是把系数矩阵除主对角线之外的元素都化为消去。也就是把系数矩阵除主对角线之外的元素都化为0。 在消去变元在消去变元xk那一列时，先把主元化为那一列时，先把主元化为1，再用它去消其它方程的，再用它去消其它方程的变元变元xk。也就是把系数矩阵主对角线上的元素都化为。也就是把系数矩阵主对角线上的元素都化为1。 一、高斯约当消元法的主要思想一、高斯约当消元法的主要思

24、想消元完毕，原线性方程组化为下面的形式：消元完毕，原线性方程组化为下面的形式： )(1,)(,)2(1, 22)2(22)1(1, 11)1(11nnnnnnnnnaxaaxaaxa这时系数矩阵这时系数矩阵 被化为单位阵，右端向量就是解向量，不再需要回代过被化为单位阵，右端向量就是解向量，不再需要回代过程。因此高斯约当消元法又称为无回代过程消元法。程。因此高斯约当消元法又称为无回代过程消元法。3.5 高斯约当消元法高斯约当消元法二、高斯约当消元法的算法二、高斯约当消元法的算法输入原方程组的阶数输入原方程组的阶数n。输入原方程组的增广矩阵输入原方程组的增广矩阵ann+1。for(k=0;k=n-

25、1;k+)，循环，循环1次消去次消去1列。列。 for(j=k+1;j=n;j+)，把主元化为，把主元化为1。 akj/=akk。 for(i=0;i=n-1;i+)，循环，循环1次更新次更新1行。行。 ik Y N for(j=k+1;j=n;j+) aij+=(-aik)*akj;输出原方程组的解（此时右端向量就是解向量）。输出原方程组的解（此时右端向量就是解向量）。3.5 高斯约当消元法高斯约当消元法高斯约当消元法对应的程序（高斯约当消元法对应的程序（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请输入原方程组的阶数：请输入原方程组的阶数：);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、约当消元法高斯约当高斯约当消元法对应的程序（消元法对应的程序（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对角阵是带状对角矩阵的对角阵是带状对角矩阵的一种。在做三次样条插值、求解微分方程等问题时，经常需要解一种。在做三次样条插值、求解微分方程等问题时，经常需要解3对角对角方程组。方程组。 一、消元形式的追赶法的主要思想一、消元形式的追赶法的主要思想定义：如果在方阵定义：如果在方阵 中，除主对角线和中，除主对角线和2条次对角线之外的元素都为条次对角线之外的元素都为0，那么那么 是是3对角阵。即对角阵。即 可以表示为下面的形式：可以表示为下面的形式：AAA=nnnnncacbacbacb1112221

29、1A定义：如果线性方程组定义：如果线性方程组 的系数矩阵的系数矩阵 是是3对角阵，那么对角阵，那么 是是3对角方程组。对角方程组。 bAx bAx A3.6 消元形式的追赶法消元形式的追赶法追赶法用于求解追赶法用于求解3对角方程组。在求解对角方程组。在求解3对角方程组时，原来的求解过程对角方程组时，原来的求解过程被大大简化。要求解的被大大简化。要求解的3对角方程组的一般形式为：对角方程组的一般形式为：nnnnnnnnnnnndxbxadxcxbxadxcxbxadxcxb1111121232221212111这里按照顺序高斯消元法的求解思路来求解这个这里按照顺序高斯消元法的求解思路来求解这个3

30、对角方程组。对角方程组。在消去在消去xi一列时，行乘子一列时，行乘子 = 。更新第更新第i1行中变元系数的公式为：行中变元系数的公式为：bi+1=bi+1ci ，di+1=di+1di其中其中i=1,2,n1。 iliiba1ilil3.6 消元形式的追赶法消元形式的追赶法消元完毕，原方程组化为下面的形式：消元完毕，原方程组化为下面的形式： nnnnnnnndxbdxcxbdxcxbdxcxb11112322212111回代过程，仍然是自下而上求解。回代的公式为：回代过程，仍然是自下而上求解。回代的公式为：xn=dn/bnxi=(dicixi+1)/bi 其中其中i=n1,n2,1。 采用这一

31、思路，可以求解带状对角线性方程组，运算量和存储需求也会采用这一思路，可以求解带状对角线性方程组，运算量和存储需求也会大大简化。大大简化。 3.6 消元形式的追赶法消元形式的追赶法二、二、消元形式的追赶法的消元形式的追赶法的算法算法输入原方程组的阶数输入原方程组的阶数n。输入增广矩阵对应的数组输入增广矩阵对应的数组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请输入原方程组的阶数：请输入原方程组的阶数：);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分解法分解法三角分解法求解线性方程组的方法，是指把线性方程组的系数矩阵进行三角分解法求解线性方程组的方法，是指把线性方程组的系数矩阵进行三角分解，然后按照线性方程组的矩阵形式进行回代求解的方法。三角分解，然后按照线性方程组的矩阵形式进行回

35、代求解的方法。 三角分解法主要有三角分解法主要有LU分解法（又称直接三角分解法、分解法（又称直接三角分解法、Doolittle分解分解法）、法）、LDR分解法、分解法、Crout分解法等。分解法等。 LU分解法把系数矩阵分解为单位下三角阵分解法把系数矩阵分解为单位下三角阵L和上三角阵和上三角阵U的积；的积；LDR分解法把系数矩阵分解为单位下三角阵分解法把系数矩阵分解为单位下三角阵L、对角阵、对角阵D和单位上三和单位上三角阵角阵R的积，其中的积，其中DR=U；Crout分解法把系数矩阵分解为下三角阵分解法把系数矩阵分解为下三角阵 和单位上三角阵和单位上三角阵R的积，的积，其中其中 =LD。 LL

36、除此之外，还有求解对称正定线性方程组的除此之外，还有求解对称正定线性方程组的LLT分解法（又称平方根分解法（又称平方根法），求解带状对角线性方程组的矩阵形式的追赶法等方法。法），求解带状对角线性方程组的矩阵形式的追赶法等方法。 3.7 LU分解法分解法对任意矩阵对任意矩阵Amn做初等变换的结果，等于做初等变换的结果，等于Amn与对应初等方阵的积。与对应初等方阵的积。 Em(i（k）,j)= Em(i（k）)= Em(i,j)= 由Em得到的初等方阵10110111k1111k数k乘以Em的j列，加到i列上去，得到Em(i（k）,j)。Em(i（k）,j)右乘Bnm的积，等于数k乘以Bnm的j列

37、，加到i列上去的结果。数k乘以Em的i行，加到j行上去，得到Em(i（k）,j)。Em(i（k）,j)左乘Amn的积，等于数k乘以Amn的i行，加到j行上去的结果。数k乘以Em的i列得到Em(i（k）)。Em(i（k）)右乘Bnm的积，等于数k乘以Bnm的i列的结果。数k乘以Em的i行得到Em(i（k）)。Em(i（k）)左乘Amn的积，等于数k乘以Amn的i行的结果。Em交换i列和j列得到Em(i,j)。Em(i,j)右乘Bnm的积，等于交换Bnm的i列和j列的结果。Em交换i行和j行得到Em(i,j)。Em(i,j)左乘Amn的积，等于交换Amn的i行和j行的结果。对应初等方阵右乘Bnm对

38、应初等方阵左乘Amn一、相关的初等方阵性质一、相关的初等方阵性质3.7 LU分解法分解法方阵的方阵的LU分解是指把方阵分解为单位下三角阵分解是指把方阵分解为单位下三角阵L和上三角阵和上三角阵U的积过程。的积过程。 二、二、LU分解与顺序高斯消元的联系分解与顺序高斯消元的联系定理：若方阵存在定理：若方阵存在LU分解，那么分解是唯一的。分解，那么分解是唯一的。 用顺序高斯消元法求解线性方程组，消元过程对应于对增广矩阵进行一用顺序高斯消元法求解线性方程组，消元过程对应于对增广矩阵进行一系列倍加变换，直到系数矩阵化为上三角阵为止。系列倍加变换，直到系数矩阵化为上三角阵为止。 定理：如果系数矩阵定理：如

39、果系数矩阵A存在存在LU分解分解A=LU，那么，那么U就是消元后的系数矩阵就是消元后的系数矩阵A，即：即：U=，L= )(,)2(, 2)2(22)1(, 1)1(12)1(11nnnnnaaaaaa111113121nlll1111232nll11111,nnl上式中上式中li,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分解法分解法 对系数矩阵对系数矩阵A进行进行LU分解，得到分解，得到L和和U。四、四、LU分解法求解线性方程组的过程分解法求解线性方程组的过程方程组方程组Ax=b(LU)x=bL(Ux)=b。 令向量令向量y=Ux，代入上式得：，代入上式得：Ly=b。回代，求解回代，求解Ly=b，得

41、到，得到y。通项公式为：通项公式为：yi=bi yk，i=1,2,n。 11ikikl 回代，求解回代，求解Ux=y，得到解向量，得到解向量x。 通项公式为：通项公式为：xi=(yi xk)/uii，i=n,n-1,1。 nikiku13.7 LU分解法分解法五、五、LU分解法的算法。分解法的算法。输入原方程组的阶数n。输入原方程组的系数矩阵ann，右端向量xn。for(k=0;k=n-2;k+) for(i=k+1;i=n-1;i+)，循环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请输入原方程组的阶数：请输入原方程组的阶数：);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对角方程组，需要对系数矩阵进行三角分解。对角方程组

46、，需要对系数矩阵进行三角分解。三角分解的方法同样分为三角分解的方法同样分为LU分解法、分解法、LDR分解法、分解法、Crout分解法。用这分解法。用这3种方法构造的追赶法的原理、求解过程和能求出解的条件相似。这里种方法构造的追赶法的原理、求解过程和能求出解的条件相似。这里以以Crout分解法为例，介绍矩阵形式的追赶法。分解法为例，介绍矩阵形式的追赶法。Crout分解的存在、唯一分解的存在、唯一性同性同LU分解。分解。 3.8 矩阵形式的追赶法矩阵形式的追赶法一、一、3对角阵对角阵Crout分解的过程分解的过程存储存储3对角方程组增广矩阵可以用一维数组对角方程组增广矩阵可以用一维数组an、bn、

47、cn和和dn，不存，不存储系数矩阵中除储系数矩阵中除3条对角线之外的零元素，以节省存储存储空间。条对角线之外的零元素，以节省存储存储空间。 Crout分解法把系数矩阵分解法把系数矩阵A分解为下三角阵分解为下三角阵 和单位上三角阵和单位上三角阵R的积，的积，即即A= R： LLnnnnncacbacbacb11122211=nnnnlalalal1122111111121nrrr对对3对角系数矩阵进行对角系数矩阵进行Crout分解的通项公式为：分解的通项公式为： 的左下次对角线元素与的左下次对角线元素对应相等。的左下次对角线元素与的左下次对角线元素对应相等。 l1=b1， li=bi-airi-

48、1，（，（i=2,3,n） ri=ci/li，（，（i=1,2,n-1）在用上述通项公式对在用上述通项公式对3对角系数矩阵进行对角系数矩阵进行Crout分解时，合逻辑的计算次分解时，合逻辑的计算次序是：序是：l1,r1,l2,r2,ln-1,rn-1,ln。按照这一次序计算，能够保证在用到某一。按照这一次序计算，能够保证在用到某一元素之前，这一元素已被求出。元素之前，这一元素已被求出。 3.8 矩阵形式的追赶法矩阵形式的追赶法用矩阵形式的追赶法求解用矩阵形式的追赶法求解3对角线性方程组对角线性方程组Ax=d的步骤如下：的步骤如下： 二、二、矩阵形式的追赶法的求解步骤矩阵形式的追赶法的求解步骤

49、对对3对角系数矩阵对角系数矩阵A进行进行Crout分解，得到下三角阵分解，得到下三角阵 和单位上三角阵和单位上三角阵R。 L方程组方程组Ax=d( R)x=d (Rx)=d。 LL 令向量令向量y=Rx，代入上式得：，代入上式得： y=d。 L回代，求解回代，求解 y=d，得到，得到y。 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 矩阵形式的追赶法矩阵形式的追赶法三、矩阵形式追赶法的算法。三、矩阵形式追赶法的算法。输入原方程组的阶数输入原方程组的阶数n。输入增广矩阵对应的数组输入增广矩阵对应的数组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请输入原方程组的阶数：请输入原方程组的阶数：);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 平方根法平方根法一、基础知识一、基础知识平方根法（又称平方根法（又称LLT分解法）用来求解对称正定线性方程组。分解法）用来求解对称正定线性方程组。 定义：如果线性方程组定义：如果线性方程组Ax=b的系数矩阵的系数矩阵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 平方根法平方根法一、基础知识（续）一、基础知识（续）定理：如果定理：如果A是对称正定矩阵，那么是对称正定矩阵，那么A满足：满足： A必定存在必定存在LLT分解，分解，A= ，其中，其中 为下三角阵。为下三角阵。 的主对角线元素都大于的主对角线元素都大于0。 A的的LLT分解唯一。分