第3章 线性方程组的数值解法_shopo

1、1第三章 线性方程组的数值解法医学技术学院医学信息工程教研室赖小波 博士，副教授，硕士研究生导师3.1 3.1 引言引言 在工程技术、自然科学和社会科学中，经常在工程技术、自然科学和社会科学中，经常遇到的许多问题（管道网络，大型输电网络，弹遇到的许多问题（管道网络，大型输电网络，弹性力学，种群繁殖等）都归结为解线性方程组，性力学，种群繁殖等）都归结为解线性方程组，再如用最小二乘法求实验数据的曲线拟合问题，再如用最小二乘法求实验数据的曲线拟合问题，工程中的三次样条函数的插值问题，经济运行中工程中的三次样条函数的插值问题，经济运行中的投入产出问题以及大地测量、机械与建筑结构的投入产出问题以及大地测

2、量、机械与建筑结构的设计计算问题等等，都需要求解某些线性方程的设计计算问题等等，都需要求解某些线性方程组。因此线性方程组的求解对于实际问题是极其组。因此线性方程组的求解对于实际问题是极其重要的。重要的。 第第3 3章章 线性方程组的数值解法线性方程组的数值解法 v为了简化我们的讨论为了简化我们的讨论, ,对于一般形式的线性方程组而言对于一般形式的线性方程组而言, ,我我们总可以假定它的系数矩阵行列式值不等于零们总可以假定它的系数矩阵行列式值不等于零. .v在上面的假定下在上面的假定下, ,我们有我们有: :如果如果mn,mn;mn;那么方程组无解那么方程组无解, ,此时作为实际问题还是有意此时

3、作为实际问题还是有意义的义的, ,我们可以转向求它的最小二乘解我们可以转向求它的最小二乘解, ,这是第这是第5 5章中讨章中讨论的问题论的问题. .如果如果A A是满秩的是满秩的n n阶方阵阶方阵, ,那么方程有唯一解那么方程有唯一解. .v注释：在这一章中，我们主要研究上面第三种情形的问题注释：在这一章中，我们主要研究上面第三种情形的问题基本结论基本结论6 . 5118.189.8 . 925. 0.688 . 38 . 0.176250212502125021xxxxxxxxx当方程组是大型的线性方程组时，怎么办？当方程组是大型的线性方程组时，怎么办？你会想到使用计算机吗？你会想到使用计算

4、机吗？ 问题引入：求解下列有问题引入：求解下列有250250个未知数的方程组个未知数的方程组这就是本章介绍的主要内容这就是本章介绍的主要内容 如果把方程组变成以下形式是不是好点？如果把方程组变成以下形式是不是好点？6 .32534.2868 .87541188 .29453.2 .668 . 38 . 02.176250250249250249225024921xxxxxxxxxx nnnnnnnnnnbxaxaxabxaxaxabxaxaxa.22112222212111212111 nnnnnnnnbbbBxxxXaaaaaaaaaA.,.,.2121212222111211 常见的线性方

5、程组是方程个数和未知量个数相同常见的线性方程组是方程个数和未知量个数相同的的n n阶线性方程组，一般形式为阶线性方程组，一般形式为 简记为简记为 Ax=bAx=b，其中，其中 ( 3.1 ) 线性方程组的数值解法一般有两类：线性方程组的数值解法一般有两类：1.1.直接法：就是经过有限步算术运算，可求得方程组直接法：就是经过有限步算术运算，可求得方程组 精确解的方法（若计算过程中没有舍入误差），如精确解的方法（若计算过程中没有舍入误差），如 克莱姆法则就是一种直接法，直接法中具有代表性克莱姆法则就是一种直接法，直接法中具有代表性 的算法是高斯的算法是高斯(Gauss)(Gauss)消元法。重要的

6、直接法全都消元法。重要的直接法全都 受到受到GaussGauss消去法的启发。计算代价高消去法的启发。计算代价高. .2. 2. 迭代法迭代法: : ( ( 第四章介绍）第四章介绍）：基于一定的递推格：基于一定的递推格 式，产生逼近方程组精确解的近似序列式，产生逼近方程组精确解的近似序列. .收敛性收敛性 是其为迭代法的前提，此外，存在收敛速度误差是其为迭代法的前提，此外，存在收敛速度误差 估计问题。估计问题。简单实用简单实用, ,诱人诱人高斯高斯 3.2 3.2 解线性方程组的直接法（高斯消去法）解线性方程组的直接法（高斯消去法） 3.2.1 3.2.1 高斯消去法的基本思想高斯消去法的基本

7、思想先用一个简单实例来说明先用一个简单实例来说明GaussGauss法的基本思想法的基本思想例例3.1 3.1 解线性方程组解线性方程组 72452413221321321xxxxxxxx解解: : 该方程组的求解过程实际上是将中学学过的该方程组的求解过程实际上是将中学学过的消元法标准化消元法标准化, ,将一个方程乘或除以某个常数将一个方程乘或除以某个常数, ,然后然后将两个方程相加减将两个方程相加减, ,逐步减少方程中的未知数逐步减少方程中的未知数, ,最终最终使每个方程只含有一个未知数使每个方程只含有一个未知数, ,从而得出所求的解。从而得出所求的解。整个过程分为消元和回代两个部分。整个过

8、程分为消元和回代两个部分。 niaamii, 3 , 2) 1 (11) 1 (11令令用用-m-mi1i1乘以乘以a a1111加到第加到第i i个方程上去就消去个方程上去就消去了第了第1 1列系数列系数4 4和和1 1，这就是第，这就是第1 1部消元，即部消元，即（1 1）消元过程）消元过程第第1 1步步: :将方程乘上将方程乘上(-2)(-2)加到方程加到方程 上去上去, ,将方程将方程 乘上乘上 加到方程加到方程 上去上去, ,这样就消去了第这样就消去了第2 2、3 3个方程的个方程的 项项, ,于是就得到等价方程组于是就得到等价方程组 211x2132325241323232321x

9、xxxxxx85425)2(22)2(3232aam令令用用-m-m3232乘以乘以a a2222加到第加到第3 3个个方程上去，即方程上去，即72452413221321321xxxxxxxx第第2 2步：将方程步：将方程 乘上乘上 加到方程加到方程 上去，这样上去，这样就消去了第就消去了第3 3个方程的个方程的 项，于是就得到等价方程项，于是就得到等价方程组组 852x4218724132332321xxxxxx这样，消元过程就是把原方程组化为上三角形方这样，消元过程就是把原方程组化为上三角形方程组，其系数矩阵是上三角矩阵。程组，其系数矩阵是上三角矩阵。 （2 2）回代过程）回代过程回代过

10、程是将上述三角形方程组自下而上求解，回代过程是将上述三角形方程组自下而上求解，从而求得原方程组的解：从而求得原方程组的解： 6, 1,9321xxx前述的消元过程相当于对原方程组前述的消元过程相当于对原方程组 741021524312321xxx702145241312bAA21323250214013121213)2()21(rrrr42187002140131223)85(rr同样可得到与原方程同样可得到与原方程组等价的方程组组等价的方程组 的增广矩阵进行下列变换的增广矩阵进行下列变换( ( 表示增广矩阵的第表示增广矩阵的第 行）行） iri 由此看出由此看出, ,高斯消去法解方程组基本思

11、想是设高斯消去法解方程组基本思想是设法消去方程组的系数矩阵法消去方程组的系数矩阵A A的主对角线下的元素的主对角线下的元素, ,而而将将Ax=bAx=b化为等价的上三角形方程组化为等价的上三角形方程组, ,然后再通过回代然后再通过回代过程便可获得方程组的解。过程便可获得方程组的解。原方程组原方程组 Ax bAx b等价上三角形方程组等价上三角形方程组 Ux bUx b从上到下依次从上到下依次消去主元素以消去主元素以下的下的 列向量列向量回代过程便可回代过程便可思思路路首先将首先将A化为上三角阵，再回代求解化为上三角阵，再回代求解 。=3.2.2 Gauss3.2.2 Gauss消去法的计算过程

12、消去法的计算过程), 2 , 1nj),(),() 1 () 1 (bAbA;, 2 , 1(niijijaa) 1 (线性方程组线性方程组(3.1)(3.1)用矩阵形式表示为用矩阵形式表示为nnnnnnnnbbbxxxaaaaaaaaa2121212222111211( 3.3 ) 这样这样, ,方程组方程组(3.1)(3.1)又可写成又可写成 。消元过。消元过程就是要按确定的计算过程对方程组进行初等行程就是要按确定的计算过程对方程组进行初等行变换变换, ,将方程组化为上三角方程组将方程组化为上三角方程组. . )1()1(bxA nnnnnnnbaaabaaabaaa21222221111

13、211(1)(1)(1)(1)(1)11121311(2)(2)(2)(2)222322(3)(3)(3)3333( )( )000000nnnnnnnnaaaabaaabaabab第第1 1步消元步消元: :假设假设 , ,作初等行变换运算作初等行变换运算0)1(11 a步骤如下：步骤如下：niiaai, 2,1111行第行第nnnnnnnbaaabaaabaaa21222221111211)2()2()2(2)2(2)2(2)2(2211121100nnnnnnbaabaabaaa), 3 , 2(,)1(11)1(11niaamiinjibmbbamaaiiijiijij,3,2,)1(

14、11)1()2()1(11)1()2(其中其中 令令)3()3()3(3)3(3)3(3)3(33)2(2)2(2)2(23)2(221113121100000nnnnnnnbaabaabaaabaaaa第第2 2步：步：niiaai, 3,2)2(22)2(2行第行第) 2() 2() 2(2) 2(2) 2(2) 2(2211121100nnnnnnbaabaabaaa第第k k步步knkknkknnknkkknkkknnbbbbxxxxaaaaaaaaa)()2(2) 1 (121)()()()()2(2)2(22) 1 (1) 1 (12) 1 (11nkiiaakkkkik, 1,k

15、)()(行第行第nkjibmbbamaakkikkikikkjikkijkij,1,)()()1()()()1(), 1()()(nkiaamkkkkikik记为记为 其中其中)()(kkbxA只要只要 , ,消元过程就可以进行下去消元过程就可以进行下去, ,直到经过直到经过n-1n-1次消元之后，消元过程结束，得到与原方程组等次消元之后，消元过程结束，得到与原方程组等价的上三角形方程组价的上三角形方程组, ,记为记为 0)(kkka)()(nnbxA或者写成或者写成 )()2(2)1(121)()2(2)2(22)1(1)1(12)1(11nnnnnnnnbbbxxxaaaaaa)()()2

16、(2)2(22)2(22)1(1)1(12)1(121)1(11nnnnnnnnnnbxabxaxabxaxaxa即即 ( (3.7) （2 2）回代过程）回代过程就是对上三角方程组（就是对上三角方程组（3.73.7）自下而上逐步回代解方）自下而上逐步回代解方程组计算，即程组计算，即 )1 ,2, 1(,)(1)()()()(niaxabxabxiiijnijiijiiinnnnnn 消元过程消元过程; ;设设 计算计算 1,2, 1,0)(nkakkk对nkkjibmbbamaaaamkkikkikikkjikkijkijkkkkikik,2,1,)()()1()()()1()()( 回代过

17、程回代过程 1,2,1)(1)()()()(nniaxabxabxiiijnijiijiiinnnnnn（3 3）高斯消去法的计算步骤：）高斯消去法的计算步骤：高高斯斯消消去去法法的的算算法法流流程程图图)()(kkkkikikaam)()()1(kkikkikibmbb)()()1(kkjikkijkijamaa)()(nnnnnnabx)(1)()(iiijnijiijiiiaxabx#include#include#include#include#define N 3 /#define N 3 /* * N N为方程组系数矩阵的阶数为方程组系数矩阵的阶数 * */ /int Gauss(

18、float aNN,float bN)int Gauss(float aNN,float bN) int i,j,k,flag=1; int i,j,k,flag=1; float t; float t; for(i=0;iN-1;i+) for(i=0;iN-1;i+) if(aii=0) if(aii=0) flag=0; break; flag=0; break; else else for(j=i+1;jN;j+) / for(j=i+1;jN;j+) /* *消元过程开始消元过程开始* */ / t=-aji/aii; t=-aji/aii; bj=bj+t bj=bj+t* *bi

19、;bi; for(k=i;kN;k+) for(k=i;kN;k+) ajk=ajk+t ajk=ajk+t* *aik;aik; return(flag); return(flag); void zg_matric(float aNN,float bN) /void zg_matric(float aNN,float bN) /* * 输出增广矩阵输出增广矩阵 * */ / int i,j; int i,j; for(i=0;iN;i+) for(i=0;iN;i+) for(j=0;jN;j+) for(j=0;j=0;i-) for(i=N-2;i=0;i-) xi=bi; xi=bi;

20、 for(j=i+1;jN;j+) for(j=i+1;jN;j+) xi=xi-aij xi=xi-aij* *xj;xj; xi=xi/aii; xi=xi/aii; for(i=0;iN;i+) / for(i=0;iN;i+) /* * 输出方程组的解输出方程组的解 * */ / printf( x%d=%11.7fn,i+1,xi); printf( x%d=%11.7fn,i+1,xi); 453113503010.A351b程序运行结果：程序运行结果： 0.100000 0.300000 0.500000 1.0000000.100000 0.300000 0.500000 1.

21、000000 3.000000 1.000000 -1.000000 5.000000 3.000000 1.000000 -1.000000 5.000000 3.000000 5.000000 4.000000 3.000000 3.000000 5.000000 4.000000 3.000000 x1= 5.4583359 x1= 5.4583359 x2= -6.5416689 x2= -6.5416689 x3= 4.8333344 x3= 4.8333344 3.2.3 3.2.3 高斯消去法的适用条件高斯消去法的适用条件 定理定理3.1 3.1 方程组系数矩阵的顺序主子式全不方

22、程组系数矩阵的顺序主子式全不 为零则高斯消去法能实现方程组的为零则高斯消去法能实现方程组的 求解。求解。证明：由于证明：由于设方程组系数矩阵设方程组系数矩阵 ，其顺序主子式，其顺序主子式 nijaA)(01111mmmmmaaaaA（m =1,2,m =1,2,，n n） 经变换得到的上三角形方程组的顺序主子式经变换得到的上三角形方程组的顺序主子式 0)()2(22)1(11)()2(2)2(22)1(1)1(12)1(11mmmmmmmmmaaaaaaaaaA所以能实现高斯消去法求解所以能实现高斯消去法求解 （m =1,2,m =1,2,，n n） 定义定义3.1 3.1 设矩阵设矩阵 每一

23、行对角元素的绝每一行对角元素的绝对值都大于同行其他元素绝对值之和对值都大于同行其他元素绝对值之和 nijaA)(niaanijjijii,2, 1,1则称则称A A为严格对角占优矩阵。为严格对角占优矩阵。 定理定理3.2 3.2 若方程组若方程组 的系数矩阵的系数矩阵A A为严格对为严格对角占优，则用高斯消去法求解时，角占优，则用高斯消去法求解时， 全不为零。全不为零。 bAx )(kkka 一般线性方程组使用高斯消去法求解时，在消元一般线性方程组使用高斯消去法求解时，在消元过程中可能会出现过程中可能会出现 的情况，这时消去法将无的情况，这时消去法将无法进行；即使法进行；即使 ，但它的绝对值很

24、小时，用其，但它的绝对值很小时，用其作除数，会导致其他元素数量级的严重增长和舍入误作除数，会导致其他元素数量级的严重增长和舍入误差的扩散，将严重影响计算结果的精度。实际计算时差的扩散，将严重影响计算结果的精度。实际计算时必须避免这类情况的发生。主元素消去法就可弥补这必须避免这类情况的发生。主元素消去法就可弥补这一缺陷。一缺陷。 0)(kkka0)(kkka32主元为零的可能性主元为零的可能性定理定理3.13.1的意义与更多讨论的意义与更多讨论p 根据顺序主子式判断高斯消去过程是否中断根据顺序主子式判断高斯消去过程是否中断20101320如何解决主元为零的问题如何解决主元为零的问题)()()()

25、()2(2)2(22)1(1)1(12)1(11knnknkkknkkknnaaaaaaaaav高斯消去法解线性方程组高斯消去法解线性方程组, ,要进行除法运算要进行除法运算. .当某个当某个a aKKKK的绝对值非常小时的绝对值非常小时, ,这种方法的数值稳定性可能这种方法的数值稳定性可能不好不好, ,为此为此. .我们可用选主元消去法我们可用选主元消去法. .v选列主元的程序相对说来简单得多，所以我们只推选列主元的程序相对说来简单得多，所以我们只推荐选列主元。荐选列主元。v选列主元的思想就是把选列主元的思想就是把a aKKKK,a,aK+1KK+1K, ,a,aMKMK中绝对值最中绝对值最

26、大的元素移到主对角线上来大的元素移到主对角线上来. .3.2.4 3.2.4 高斯主元素消去法高斯主元素消去法v交换原则：通过方程或变量次序的交换，使在对角交换原则：通过方程或变量次序的交换，使在对角线位置上获得绝对值尽可能大的系数作为线位置上获得绝对值尽可能大的系数作为a akkkk(k)(k)，称，称这样的这样的a akkkk(k)(k) 为为主元素主元素，并称使用主元素的消元法，并称使用主元素的消元法为主元素法为主元素法记笔记记笔记3.2.4 3.2.4 高斯主元素消去法高斯主元素消去法选主元的方法有选列主元和全主元之分，基本原选主元的方法有选列主元和全主元之分，基本原理是互换任意两个方

27、程的位置不影响方程组的解；理是互换任意两个方程的位置不影响方程组的解；互换任意两个变量的位置也不影响方程组的解互换任意两个变量的位置也不影响方程组的解 例例3.2 3.2 用高斯消去法求下列方程组的解用高斯消去法求下列方程组的解 211021219xxxx解：确定乘数解：确定乘数 , ,再计算系数再计算系数 9)2(29122122)2(22102101bamaa这时方程组的实际形式是这时方程组的实际形式是 9292191010110 xxx9911212110101aam单精度解方程组单精度解方程组9999101010000000001. 01019999101010000000002. 0

28、102 由此回代解出由此回代解出 , ,但这个解不满但这个解不满足原方程组足原方程组, ,解是错误的。这是因为所用的除数太解是错误的。这是因为所用的除数太小使得上式在消元过程中小使得上式在消元过程中“吃掉吃掉”了下式，解决这了下式，解决这个问题的方法之一就是采用列选主元高斯消元法。个问题的方法之一就是采用列选主元高斯消元法。即按列选绝对值大的系数作为主元素，则将方程组即按列选绝对值大的系数作为主元素，则将方程组中的两个方程相交换，原方程组变为中的两个方程相交换，原方程组变为1, 021xx110221921xxxx929211021)101(2xxx得到消元后的方程组得到消元后的方程组 用小主

29、元用小主元1010-9-9 作除数，致作除数，致使其它元素使其它元素的数量级大的数量级大大增加，舍大增加，舍入误差的扩入误差的扩散将准确解散将准确解淹没了。淹没了。因而方程组的实际形式是因而方程组的实际形式是 12221xxx由此回代解出由此回代解出 , ,这个结果是正确的这个结果是正确的 1, 121xx 可见用高斯消去法解方程组时可见用高斯消去法解方程组时, ,小主元可能导致小主元可能导致计算失败计算失败, ,因为用绝对值很小的数作除数因为用绝对值很小的数作除数, ,乘数很大乘数很大, ,引起约化中间结果数量级严重增长引起约化中间结果数量级严重增长, ,再舍入就使得计再舍入就使得计算结果不

30、可靠了算结果不可靠了, ,故避免采用绝对值很小的主元素。故避免采用绝对值很小的主元素。以便减少计算过程中舍入误差对计算解的影响。以便减少计算过程中舍入误差对计算解的影响。全主元素消去法全主元素消去法 是通过方程或变量次序的交换，使在对角线位是通过方程或变量次序的交换，使在对角线位置上获得绝对值尽可能大的系数作为置上获得绝对值尽可能大的系数作为 称这样的称这样的 为主元素。尽管它的算法更稳定，但计算为主元素。尽管它的算法更稳定，但计算量较大，实际应用中大多数使用列主元素消去法即量较大，实际应用中大多数使用列主元素消去法即可满足需要。可满足需要。 )(kkka)(kkka不是按列选主元素，而是在不

31、是按列选主元素，而是在全体待选系数中选取，则得全体待选系数中选取，则得全主元素法。全主元素法。 10 x1 - 19x2 - 2x3=3 (1)-20 x1 +40 x2 + x3 =4 (2) x1 + 4x2 + 5x3=5 (3) 解：选择所有系数中绝对值最大的解：选择所有系数中绝对值最大的4040作为作为主元素主元素，交换第一、二行和交换第一、二列使该主元素位，交换第一、二行和交换第一、二列使该主元素位于对角线的第一个位置上，得于对角线的第一个位置上，得40 x2 - 20 x1 + x3 =4 (4)-19x2+10 x1 - 2x3=3 (5) 4x2+ x1 +5x3=5 (6)

32、例例3.3 3.3 用用全主元素法解下列线组全主元素法解下列线组3121212111111940.4750.14040aammaa 计算计算1 . 0404113131aam(5)- m(5)- m2121(4), (4), (6)- m(6)- m3131(4) (4) 消去消去x x2 2 得得(5)- m(5)- m2121(4), (6)- m(4), (6)- m3131(4) (4) 消去消去x x2 2 得得 0.5x1 1.525 x3 =4.9 (7) 3x1 + 4.9 x3 =4.6 (8)选选4.94.9为主元素为主元素 4.9 x3 + 3x1=4.6 (9)-1.5

33、25 x3 +0.5x1=4.9 (10)记笔记记笔记计算计算(10)- m(10)- m3232(9)(9)消去消去x x2 2得得1.43366x1.43366x1 1=6. 33161 (11)=6. 33161 (11)3232221.5250.311224.9ama3121212111111940.4750.14040aammaa 计算计算1 . 0404113131aam保留有主元素的方程保留有主元素的方程40 x2 - 20 x1 + x3 =4 (4) 4.9x3 + 3x1=4.6 (9) 1.43366x1=6. 33161 (11)进行回代进行回代x1=4.41637 x

34、3 =-1.76511x2=2.352303.2.4.1 3.2.4.1 列主元素法列主元素法 列主元素法就是在待消元的所在列中选取主元列主元素法就是在待消元的所在列中选取主元, ,经方程的行交换经方程的行交换, ,置主元素于对角线位置后进行消元置主元素于对角线位置后进行消元的方法的方法 10 x1 - 19x2 - 2x3=3 (1)-20 x1 +40 x2 + x3 =4 (2) x1 + 4x2 + 5x3=5 (3)解：选择解：选择-20-20作为该列的作为该列的主元素，主元素，-20 x1 +40 x2 + x3 =4 (4) 10 x1 - 19x2 - 2x3=3 (5) x1

35、 + 4x2 + 5x3=5 (6)计算计算m m2121=10/-20=10/-20=-0.5=-0.5 m m3131=1/-20=1/-20=-0.05=-0.05(5)- m(5)- m2121(4), (6)- m(4), (6)- m3131(4)(4)得得 x2 1.5x3=5 (7)6x2 + 5.05x3=5.2 (8)选选6 6为主元素为主元素6x2 + 5.05x3=5.2 (9) x2 1.5x3=5 (10)计算计算m m3232=1/6=1/6=0.16667=0.16667， (10)- m(10)- m3232(9) (9) 得得-2.34168x-2.3416

36、8x3 3=4.13332 (11)=4.13332 (11)记笔记记笔记-20 x1 +40 x2 + x3 = 4 保留有主元素的方程保留有主元素的方程 -20 x1 +40 x2 + x3 =4 (4) 6x2 + 5.05x3 =5.2 (9) -2.34168x3=4.13332 (11)进行回代进行回代x3 =-1.76511x2=2.35230 x1=4.41637记笔记记笔记 列选主元素的计算方法与高斯消去法完全一样列选主元素的计算方法与高斯消去法完全一样, ,不同的是在每步消元之前要按列选出主元不同的是在每步消元之前要按列选出主元例例3.5 3.5 用矩阵的初等行变换求解方程

37、组用矩阵的初等行变换求解方程组 754217743322321321321xxxxxxxxx 解解: : 用矩阵的初等行变换求解用矩阵的初等行变换求解, ,对增广矩阵对增广矩阵 ( (下面带下划线元素为主元素下面带下划线元素为主元素) ) 2 . 12 . 1005 . 65 . 85 . 7017745 . 25 . 05 . 105 . 65 . 85 . 7017745 . 65 . 85 . 705 . 25 . 05 . 101774754233221774754217743322232313121251_2121_) 1 (rrrrrrrrrrbAA列选主元消去法的算法实现列选主元

38、消去法的算法实现（1 1）输入方程组的阶数）输入方程组的阶数n n，系数矩阵，系数矩阵A A和右端常数矩阵和右端常数矩阵b b 。（2 2）列主元素）列主元素 对对 ，从，从 中选出中选出绝对值最大的元素绝对值最大的元素 ，将，将 行和行和 行交换后，再作行交换后，再作第第 次消元操作。次消元操作。 （3 3）消元过程：）消元过程： 对对 , ,计算：计算： 1, 2 , 1nk)k(kn,)k(k,1kk)(kk,a,aa)(,kkmakm1, 2 , 1nk)()(kkkkkikiaam)()() 1(kjkkikjikjiamaa),2,1,(nkkji)()()1(kkkikikibm

39、bbk（4 4） 回代过程：回代过程：（5 5） 输出方程组的解。输出方程组的解。（6 6） 结束结束)()(nnnnnnabx )(1)()(iiinijjijiiiiaxabx)1, 2,1( ni例例 3.6 3.6 用列主元高斯消去法求解线性方程组用列主元高斯消去法求解线性方程组 ，式中，式中 1121332300441135087200609521329495625060455101123330504133300031.A7723136040750053.bbAx void ColGauss(float aNN,float bN)void ColGauss(float aNN,flo

40、at bN) float t,max_fab; float t,max_fab; int i,j,k,l; int i,j,k,l; for(i=0;iN-1;i+) for(i=0;iN-1;i+) l=i; l=i; max_fab=fabs(aii); max_fab=fabs(aii); for(j=i+1;jN;j+) / for(j=i+1;jmax_fab) if(fabs(ajimax_fab) ll=j; max_fab=fabs(aji); ll=j; max_fab=fabs(aji); if(il) / if(il) /* * 交换交换i i、l l两行两行 * */

41、/ t=bi; t=bi; bi=bl; bi=bl; bl=t; bl=t; for(j=i;jN;j+) for(j=i;jN;j+) t=aij; aij=alj; alj=t; t=aij; aij=alj; alj=t; for(j=i+1;jN;j+) / for(j=i+1;jN;j+) /* * 消元过程开始消元过程开始 * */ / t=-aji/aii; t=-aji/aii; bj=bj+t bj=bj+t* *bi;bi; for(k=i;kN;k+) for(k=i;kN;k+) ajk=ajk+t ajk=ajk+t* *aik; aik; void zg_matr

42、ic(float aNN,float bN) /void zg_matric(float aNN,float bN) /* * 输出增广矩阵输出增广矩阵 * */ / int i,j; int i,j; for(i=0;iN;i+) for(i=0;iN;i+) for(j=0;jN;j+) for(j=0;jN;j+) printf(%10f,aij); printf(%10f,aij); printf(%10f,bi); printf(%10f,bi); printf(n); printf(n); printf(n); printf(n); #include #include #defin

43、e N 4 /#define N 4 /* * N N为方程组系数矩阵的阶数为方程组系数矩阵的阶数 * */ /main()main() float aNN=1.003,0.333,1.504,-0.333,-2.011,1.455,0.506,2.956, float aNN=1.003,0.333,1.504,-0.333,-2.011,1.455,0.506,2.956, 4.329,-1.952,0.006,2.087, 4.329,-1.952,0.006,2.087,5.113,-4.004,3.332,-1.1125.113,-4.004,3.332,-1.112; float

44、bN=3.005,5.407,0.136,3.772; float bN=3.005,5.407,0.136,3.772; float xN=0,0,0,0; float xN=0,0,0,0; int i,j; int i,j; zg_matric(a,b); zg_matric(a,b); ColGauss(a,b); ColGauss(a,b); xN-1=bN-1/aN-1N-1; / xN-1=bN-1/aN-1N-1; /* * 回代过程开始回代过程开始 * */ / for(i=N-2;i=0;i-) for(i=N-2;i=0;i-) xi=bi; xi=bi; for(j=i

45、+1;jN;j+) for(j=i+1;jN;j+) xi=xi-aij xi=xi-aij* *xj;xj; xi=xi/aii; xi=xi/aii; for(i=0;iN;i+) / for(i=0;i 0 A 的顺序主子阵的顺序主子阵 Ak 亦对称正定亦对称正定 A 的特征值的特征值 i 0 A 的全部顺序主子式的全部顺序主子式 det ( Ak ) 0（充要条件）（充要条件）3.4 3.4 平方根法平方根法定理定理3.6 3.6 对对 阶线性方程组阶线性方程组 ，当系数矩阵，当系数矩阵 为对称正定时可以惟一地分解为为对称正定时可以惟一地分解为 ，其中，其中 为上三角矩阵。即写成矩阵元

46、素的形式如下：为上三角矩阵。即写成矩阵元素的形式如下： nbAx ATLLA Lnnnnnnnnnnnnnnllllllllllllaaaaaaaaa00000022211211212221112122221112111112111112112122221222221212nnnnnnnnnnnnnnaaallllaaallllaaallll 按矩阵乘法展开，可逐行求出分解矩阵按矩阵乘法展开，可逐行求出分解矩阵L L的元素，的元素，计算公式是对于计算公式是对于i=i=1,21,2, ,n,n21111al121121all1111nnall22222122all211121all2211222

47、nnnallll 移项便得移项便得Cholesky分解公式分解公式21112)(ikikiiiilaliiikikjkjijilllal11j=i+j=i+1 1,i+,i+2 2, ,n,n 因对称性因对称性无需存储无需存储11l21l31l1nlStep122l32l2nlStep233l3nlStep312a13a1na23a2na3nannlStepn的计算过程的计算过程: :L逐逐 列列 计计 算算1 2(;, , )iklik kn元素元素A仍然存放在矩阵仍然存放在矩阵 的相应位置上的相应位置上乔累斯基（乔累斯基（CholeskyCholesky）分解所需要的乘除次数约为）分解所需

48、要的乘除次数约为 数量级，比数量级，比LULU分解节省近一般的工作量分解节省近一般的工作量 361n其分解法的计算过程分为两个步骤：其分解法的计算过程分为两个步骤：1. 1. 将对称正定矩阵将对称正定矩阵A A分解，有分解，有A=A=LLLLT T, , 算出算出L L。 2. 2. 将将A A= = 代入线性方程组代入线性方程组Ax=bAx=b，有，有 x=bx=b。因。因 此求解两个三角形方程组此求解两个三角形方程组L Ly=by=b及及TLLTLLyxLT相应的求解公式相应的求解公式iikikikiilylby11i=1,2,n iiknikkiiilxlyx1i=n,n-1,2,1 例

49、例3.123.12用用Cholesky分分解法求解下列方程组解法求解下列方程组123123123464252750527535125.xxxxxxxxx 解：解：系数矩阵为系数矩阵为414 2512 753 5.A 20 54 250 52 753 5. Step120 520 51 53 5. Step220 520 51 51. Step320 520 51 51.L Lyb 求解方程组求解方程组20 50 521 51.TL TL xy 求解方程组求解方程组60 51 25(.)Tb 30 51(.)Ty 211()Tx 3.5 3.5 追赶法追赶法 常常会遇到三对角线性方程组常常会遇到

50、三对角线性方程组n nAxfAR 其中其中1b1c2a2b2c3a3b3cnanb1nc 1nb 1na A f 1f2f3f1nf nf11;nnbcba 02 31;, ,iiiiibaca cin 则则 非奇异，且非奇异，且A01 2, ,iuin A如果三对角矩阵如果三对角矩阵A 为对角占优矩阵，则可为对角占优矩阵，则可以采用追赶法求解。以采用追赶法求解。三对角矩阵三对角矩阵11222111nnnnnbcabcAabcab满足条件满足条件系数系数A矩阵可作三角分解矩阵可作三角分解：ALU 112221111nnnlualuLUual 11112122111122211nnnnnnnnu

51、uulalalbacbacbacb1,2, 111111niluabulclbiiiiiii按乘法展开按乘法展开 1,2,1/11111niuabllcubliiiiiii则可计算则可计算 nnluulul12211可依次计算可依次计算 112221111nnnlualuAual nnlalal221=12nfff12nyyy111121nuuu=12nxxx12nyyy分解为分解为“追追”的过程的过程方程组求解的计算公式：方程组求解的计算公式： 解方程组解方程组Lyf 解方程组解方程组Uxy nnxy 111,iiiixyu xin “赶赶”的过程的过程111/yfl 12 3() /, ,

52、iiiiiyfl ylin annlalal221=12nfff12nyyy111121nuuu=12nxxx12nyyy1,2, 1/11111niuabllcubliiiiiii则可计算则可计算 111152242231123214433221uuulalalal11112122111122211nnnnnnnnuuulalalbacbacbacb1651541211252512225125224223112例例3.11 3.11 用追赶法求解三对角方程组用追赶法求解三对角方程组 501352242231124321xxxx解解 21211111lcubl54252221222lcuuab

53、l565123332333lcuuabl253444uabl16515412112525122251252242231121, 5 . 12321222111lyafylfy1,654344432333lyafylyafy0, 1433344xuyxyx2, 121113222xuyxxuyx由由Ly=fLy=f 解出解出y y又由又由Ux=yUx=y解出解出x x 开始开始输入数据输入数据c1/ b1 r1 d1/ b1 y12 ibi ri-1aiq ci /q ri(di yi-1ai)/qyii0,0,使得对任意使得对任意 恒有恒有qpxx,nRnRxpqpxCxxC21（证（证: :

54、略）略） 定理定理3.93.9 其中其中 为向量中的任一种范数。为向量中的任一种范数。 *)(limxxkk 0lim*)(xxkk证证: :由于由于 *)(limxxkk而对于而对于 上的任一种上的任一种 范数范数, , 由定理由定理3.73.7知存在常数知存在常数C1C1，C2C2，使，使 nR*)(2*)(*)(1xxCxxxxCkkk于是可得于是可得 ( )*( )*lim0lim0kkkkxxxx从而定理得证。从而定理得证。 ( )*lim0kkxx定义定义3.53.5（矩阵的范数）如果矩阵（矩阵的范数）如果矩阵 的的某个非负的实值函数某个非负的实值函数 ，满足，满足正定性：正定性：

55、齐次性：齐次性：三角不等性：三角不等性：A0,n nAR 且且A00A ,n nAAARR ,n nABABA BR 则称则称 为为 中矩阵中矩阵 的范数。的范数。An nR A相容性相容性:,n nABA BA BR nnRAAAN)(复习复习向量范数向量范数：设设12(,)Tnxx xx 1- -范数：范数： 2- -范数：范数： - -范数：范数： 11niixx 12221()( , )niixxx x 1maxii nxx 上述上述3种向量范数统称为种向量范数统称为P- -范数范数111()nppipixxp 列范数：列范数：111maxnijj niAa 行范数：行范数：11max

56、niji njAa 其中其中 是是 的最大特征值的最大特征值1 TA A 谱范数：谱范数：12A 12()TA A ()ijn nAa 矩阵范数的计算公式矩阵范数的计算公式定理定理3.10 3.10 对对n n 阶方阵阶方阵则则定义定义3.73.7（矩阵的谱半径）设（矩阵的谱半径）设 的特征的特征 值为值为 , , 称称 为为A A的谱半径。的谱半径。nnRA), 2 , 1(niiiniA1max)(), 2 , 1(pAp420420001A的三种常用范数的三种常用范数 例例 3.13 3.13 计算方阵计算方阵谱半径谱半径1( )maxii nA 解：解： 420420001A 88 ,

57、4, 1maxmax31311iijjaA66,6,1maxmax3131jijiaA)(max2AAA先计算先计算 3200080001420420001440220001AA32)(maxAA24322A定理定理3.11 3.11 设设A A为为n n阶方阵阶方阵, , 则对任意矩阵范数则对任意矩阵范数 都有都有AA )(证证: : 设设 为为A A的特征值，的特征值，x是对应于是对应于 的特征向量的特征向量, ,则则 x=A=Ax。两端取范数并。两端取范数并 依据其性质得依据其性质得xAAxxx由于由于x0 x0，故，故 ，所以，所以AAA )(3.7 3.7 误差分析误差分析3.7.1

58、 3.7.1 方程组的性态方程组的性态 在建立方程组时，其系数往往含有误差在建立方程组时，其系数往往含有误差（如观测误差或计算误差），就是说，所要（如观测误差或计算误差），就是说，所要求解的运算是有扰动的方程组，因此需要研求解的运算是有扰动的方程组，因此需要研究扰动对解的影响。究扰动对解的影响。 例例3.14 3.14 考察方程组考察方程组 0001. 20001. 122121xxxx20001.122121xxxx和和 上述两个方程组尽管只是右端项有微小扰动上述两个方程组尽管只是右端项有微小扰动, , 但解大不相同但解大不相同, , 第第1 1个方程组的解是个方程组的解是第第2 2个方程组

59、的解是个方程组的解是 。这类方程组称为病态的这类方程组称为病态的。121xx0, 221xx例例3.153.15考察方程组：考察方程组：123143 0001 14 0001.xx 精确解为精确解为1 1()Tx 设方程组存在扰动设方程组存在扰动123142 999914 0002.xx 精确解为精确解为210()Tx 上例说明该方程组的解对初始元素的扰动非常敏感。上例说明该方程组的解对初始元素的扰动非常敏感。定义定义3.8 A3.8 A或或b b的微小变化的微小变化( (又称扰动或摄动又称扰动或摄动) )引起引起方程组方程组A Ax=b=b解的巨大变化，则称方程组为病态方程解的巨大变化，则称

60、方程组为病态方程组，矩阵组，矩阵A A称为病态矩阵。否则方程组是良态方程称为病态矩阵。否则方程组是良态方程组，矩阵组，矩阵A A也是良态矩阵也是良态矩阵 为了定量地刻画方程组为了定量地刻画方程组“病态病态”的程度，要对的程度，要对方程组方程组A Ax=b=b进行讨论，考察进行讨论，考察A A（或（或b b）微小误差对解）微小误差对解的影响。为此先引入矩阵条件数的概念。的影响。为此先引入矩阵条件数的概念。 定义定义3.9 3.9 （矩阵条件数）设（矩阵条件数）设A A为非奇异矩阵，为非奇异矩阵，称称 为矩阵为矩阵A A条件数。条件数。 1)(AAAcond设方程组为设方程组为Axb 系数矩阵系数