第4章线性方程组的数值解法精选PPT.ppt

《第4章线性方程组的数值解法精选PPT.ppt》由会员分享，可在线阅读，更多相关《第4章线性方程组的数值解法精选PPT.ppt（119页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第4章线性方程组的数值解法第1页，本讲稿共119页4.1引言引言在工程技术和科学研究中，很多科学计算的问题往往直接在工程技术和科学研究中，很多科学计算的问题往往直接或间接地归结为求解线性代数方程组。常见的线性代数方程或间接地归结为求解线性代数方程组。常见的线性代数方程组是方程个数和未知量个数相同的阶线性方程组，它的一般组是方程个数和未知量个数相同的阶线性方程组，它的一般形式是形式是其中，其中，、为常数，为常数，为待求的未知量为待求的未知量(4.1)第2页，本讲稿共119页用矩阵形式表示是用矩阵形式表示是其中：其中：一般一般。当系数矩阵。当系数矩阵非奇异（行列式不为零），即非奇异（行列式不为零）

2、，即时，方程组时，方程组(4.1)有有惟一解。惟一解。第3页，本讲稿共119页对于线性代数方程组的解法大致可分为两类，即直接解对于线性代数方程组的解法大致可分为两类，即直接解法和迭代解法。直接解法是指在假设没有舍入误差的条件下，法和迭代解法。直接解法是指在假设没有舍入误差的条件下，经过有限次算数运算就能求得方程组精确解的方法。经过有限次算数运算就能求得方程组精确解的方法。然而实际计算中由于舍入误差的影响，这类方法也只能求然而实际计算中由于舍入误差的影响，这类方法也只能求得近似解；迭代解法就是从一个已知的初始近似值开始，按一得近似解；迭代解法就是从一个已知的初始近似值开始，按一定的法则逐步求出解

3、的各个更准确的近似值的方法，它是用某定的法则逐步求出解的各个更准确的近似值的方法，它是用某种极限过程去逐步逼近精确解的方法。种极限过程去逐步逼近精确解的方法。在本章中主要介绍高斯消去法、列主元高斯消去法、约在本章中主要介绍高斯消去法、列主元高斯消去法、约当消去法、三角分解法等直接解法以及雅可比、高斯当消去法、三角分解法等直接解法以及雅可比、高斯赛德赛德尔和超松弛等迭代解法。尔和超松弛等迭代解法。第4页，本讲稿共119页4.2高斯（高斯（Gauss）消去法）消去法4.2.1 4.2.1 高斯消去法的基本思想高斯消去法的基本思想 高斯消去法是最古老的求解线性代高斯消去法是最古老的求解线性代数方程组

4、的方法之一，高斯消去法是消数方程组的方法之一，高斯消去法是消去法的一种特殊形式，它包括消元与回去法的一种特殊形式，它包括消元与回代两个过程。代两个过程。第5页，本讲稿共119页设设 ，对行计算乘数，对行计算乘数用用 乘以第乘以第 个方程后加到第个方程后加到第 个方程到第个方程到第 个方程中，消去个方程中，消去 第个第个方程到第方程到第 个方程的未知数个方程的未知数 ，得到，得到第6页，本讲稿共119页 只只要要设设 就就可可以以继继续续进进行行消消元元，直直到到经经过过 次次消消元元后后，将将线线性性方程组方程组(4.1)(4.1)化为化为(4.2)(4.2)所示上三角方程组（以上计算过程称为

5、消元过程）。所示上三角方程组（以上计算过程称为消元过程）。(4.2)(4.2)消消元元过过程程结结束束后后，只只要要设设 ，对对上上三三角角方方程程组组就就可可以以自自下下而而上上逐逐步回代，依次求得步回代，依次求得 ，即，即第7页，本讲稿共119页4.2.2 4.2.2 实现高斯消去法的基本步骤实现高斯消去法的基本步骤（1 1）输入方程组的阶数输入方程组的阶数 n n，系数矩阵，系数矩阵A A和右端常数矩阵和右端常数矩阵b b。（2 2）消元过程消元过程:设设 ，对，对 ，计算：，计算：（3 3）回代过程回代过程（4 4）输出方程组的解。输出方程组的解。（5 5）结束。结束。第8页，本讲稿共

6、119页第9页，本讲稿共119页例例4.1 4.1 用高斯消去法求解线性方程组，式中用高斯消去法求解线性方程组，式中第10页，本讲稿共119页#define N 3/*N#define N 3/*N为方程组系数矩阵的阶数为方程组系数矩阵的阶数*/*/int Gauss(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

7、=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*bi;bj=bj+t*bi;for(k=i;kN;k+)for(k=i;kN;k+)ajk=ajk+t*aik;ajk=ajk+t*aik;return(flag);return(flag);第11页，本讲稿共119页void zg_matric(float aNN,float bN)/*void zg_matric(float aNN,float bN)/*输出增广矩阵输出增广矩阵*/*/in

8、t 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;for(j=i+1;jN;j+)for(j=i+1;jN;j+)xi=xi-aij*xj;xi=xi-aij*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);第13页，本讲稿共119页程序运行结

9、果：程序运行结果：0.100000 0.300000 0.500000 1.000000 0.100000 0.300000 0.500000 1.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.8

10、333344第14页，本讲稿共119页4.3 4.3 列主元高斯消去法列主元高斯消去法4.3.1 4.3.1 列主元高斯消去法的基本思想列主元高斯消去法的基本思想 高斯消去法，一直是在假设高斯消去法，一直是在假设 的情况下进行的，并且是按照方程组中各的情况下进行的，并且是按照方程组中各方程所给定的顺序进行的，所以这种方法又称为顺序消去法。但当方程所给定的顺序进行的，所以这种方法又称为顺序消去法。但当 时，消时，消元过程就无法进行。另外，即使元过程就无法进行。另外，即使 ，但其绝对值很小时，用它作除数时，根，但其绝对值很小时，用它作除数时，根据数值运算中据数值运算中“用绝对值很小的数做除数，舍入

11、误差会增大，而且严重影响计算结果用绝对值很小的数做除数，舍入误差会增大，而且严重影响计算结果的精度的精度”的原则，这种方法在一定程度上也具有局限性。为了克服这一局的原则，这种方法在一定程度上也具有局限性。为了克服这一局限性，产生了列主元高斯消去法。限性，产生了列主元高斯消去法。这种消去法在高斯消去法每次进行消元之前先选取绝对值最大的元素作为主元这种消去法在高斯消去法每次进行消元之前先选取绝对值最大的元素作为主元素（位于主对角线上的在消去过程中用作除数的元素称为主元素，简称主元）进行素（位于主对角线上的在消去过程中用作除数的元素称为主元素，简称主元）进行消去。消去。第15页，本讲稿共119页 列

12、列 主主 元元 高高 斯斯 消消 去去 法法，是是 在在 高高 斯斯 消消 去去 法法 的的 消消 去去 过过 程程 中中，第第 步步()增增 加加 选选 主主 元元 操操 作作，成成 为为：首首 先先 在在 第第 列列 中中，从从 中中选选出出绝绝对对值值最最大大的的元元素素，经经过过行行交交换换(把把绝绝对对值值最最大大者者所所在在的的行行与与第第 行行交交换换)把把绝绝对对值值最最大大的的元元素素交交换换到到 的的单单元元中中；然然后后作作GaussGauss消消去去法法的的消消去去过过程程中中第第 步步，初初等等变变换换产产生生第第 列列的的 个个零零元元素素。回代过程与回代过程与Ga

13、ussGauss消去法的回代过程相同。消去法的回代过程相同。以上就是列主元高斯消去法的操作步骤。由于以上就是列主元高斯消去法的操作步骤。由于 ，可证明可证明 ，中至少有一个元素不为零，从理中至少有一个元素不为零，从理论上保证了列主元高斯消去法的可行性。论上保证了列主元高斯消去法的可行性。第16页，本讲稿共119页4.3.2 4.3.2 实现列主元高斯消去法的基本步骤实现列主元高斯消去法的基本步骤（1 1）输入方程组的阶数）输入方程组的阶数n n，系数矩阵，系数矩阵A A和右端常数矩阵和右端常数矩阵b b 。（2 2）列主元素）列主元素 对对 ，从，从 中选出绝对值最大的元素中选出绝对值最大的元

14、素 ，将，将 行和行和 行交换后，再作第行交换后，再作第 次消元操作。次消元操作。（3 3）消元过程：）消元过程：对对 ,计算：计算：第17页，本讲稿共119页（4）回代过程：（5）输出方程组的解。（6）结束第18页，本讲稿共119页第19页，本讲稿共119页第20页，本讲稿共119页例例 4.2 4.2 用列主元高斯消去法求解线性方程组用列主元高斯消去法求解线性方程组 ，式中，式中第21页，本讲稿共119页void ColGauss(float aNN,float bN)void ColGauss(float aNN,float bN)float t,max_fab;float t,max_

15、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)l=j;max_fab=fabs(aji);l=j;max_fab=fabs(aji);if(il)/*if(il)/*交换交换i i、l l两行两行*/*/t=bi;t=bi;bi=bl;bi=bl;bl=t;bl=t;for(j=i;jN;j+)for(j=i;jN;j+)t

16、=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*bi;bj=bj+t*bi;for(k=i;kN;k+)for(k=i;kN;k+)ajk=ajk+t*aik;ajk=ajk+t*aik;第22页，本讲稿共119页void zg_matric(float aNN,float bN)/*void zg_matric(float aNN,float bN)/*输出增广矩阵输出增广矩阵*/*/int i,j;

17、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);第23页，本讲稿共119页#include#include#define N 4 /*N#define N 4 /*N为方程组系数矩阵的阶数为方程组系数矩阵的阶数*/*/main()main()float aNN=1.003,0.333,1.504,-0.33

18、3,-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,5.113,-4.004,3.332,-1.112 4.329,-1.952,0.006,2.087,5.113,-4.004,3.332,-1.112;float 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

19、 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+1;jN;j+)for(j=i+1;jN;j+)xi=xi-aij*xj;xi=xi-aij*xj;xi=xi/aii;xi=xi/aii;for(i=0;iN;i+)/*for(i=0;iN;i+)/*输出方程组的解输出方程组的解*/*/printf(

20、x%d=%11.7fn,i+1,xi);printf(x%d=%11.7fn,i+1,xi);第24页，本讲稿共119页程序运行结果：程序运行结果：1.003000 0.333000 1.504000 -0.333000 3.005000 1.003000 0.333000 1.504000 -0.333000 3.005000 -2.011000 1.455000 0.506000 2.956000 5.407000 -2.011000 1.455000 0.506000 2.956000 5.407000 4.329000 -1.952000 0.006000 2.087000 0.136

21、000 4.329000 -1.952000 0.006000 2.087000 0.136000 5.113000 -4.004000 3.332000 -1.112000 3.772000 5.113000 -4.004000 3.332000 -1.112000 3.772000 x1=-0.325438 x1=-0.325438 x2=0.327990 x2=0.327990 x3=2.372718 x3=2.372718 x4=1.040163 x4=1.040163第25页，本讲稿共119页 列主元高斯消去法列主元高斯消去法,虽然保证了当前的主元素虽然保证了当前的主元素 是第是第

22、列中列中 以下元素中的绝对值最大者，但还不能保证它是同一行（既第以下元素中的绝对值最大者，但还不能保证它是同一行（既第 行）中的绝对值最大者。因此，经过列选主元后，虽然避免了行）中的绝对值最大者。因此，经过列选主元后，虽然避免了 的情况发生，但其计算过程还是不稳定的，不适于求解较大规模的情况发生，但其计算过程还是不稳定的，不适于求解较大规模的线性方程组，因而这种方法在一定程度上也具有局限性。为了的线性方程组，因而这种方法在一定程度上也具有局限性。为了克服这一局限性，产生了全主元高斯消去法。克服这一局限性，产生了全主元高斯消去法。列主元高斯消去法当变换到第列主元高斯消去法当变换到第 步时，从系数

23、矩阵的右步时，从系数矩阵的右下角下角 阶子矩阵中选取绝对值最大的元素，然后通过阶子矩阵中选取绝对值最大的元素，然后通过行变换与列变换将它交换到主元素行变换与列变换将它交换到主元素 的位置上。回代过程的位置上。回代过程与与GaussGauss消去法的回代过程相同，请读者参考有关书籍。消去法的回代过程相同，请读者参考有关书籍。第26页，本讲稿共119页4.4 4.4 约当（约当（Jordan）消去法）消去法4.4.1 4.4.1 约当消去法的基本思想约当消去法的基本思想 约约当当（JordanJordan）消消去去法法是是一一种种无无回回代代过过程程的的求求解解线线性性方方程程组组的的直直接接解解

24、法法，基基本本思思想想是是，将将线线性性方方程程组组的的系系数数矩矩阵阵通通过过初初等等变变换换，化化为为一一个个等等价价的的对对角角阵阵（只只有有主主对对角角线线上上的的元元素素不不为为0 0，其其余余元元素素均均为为0 0），从而方便地求出方程组的解。从而方便地求出方程组的解。第27页，本讲稿共119页 假如我们将此增广矩阵通过次初等变换，化成下列等价形式（式中每假如我们将此增广矩阵通过次初等变换，化成下列等价形式（式中每个元素右上角圆括号内的数字表示消去计算的次数）个元素右上角圆括号内的数字表示消去计算的次数）由此可立即求得方程组的解：由此可立即求得方程组的解：第28页，本讲稿共119页

25、4.4.2 4.4.2 实现约当消去法的基本步骤实现约当消去法的基本步骤第29页，本讲稿共119页第30页，本讲稿共119页例例 4.3 4.3 用约当消去法求解线性方程组用约当消去法求解线性方程组 ，式中，式中，第31页，本讲稿共119页#define N 3 /*N#define N 3 /*N为方程组的系数矩阵的阶数为方程组的系数矩阵的阶数*/*/int joydan(float aNN+1)int joydan(float aNN+1)float r;float r;int i,j,k,flag=1;int i,j,k,flag=1;for(k=0;kN;k+)for(k=0;kN;k

26、+)if(akk=0)/*if(akk=0)/*当当akkakk为为0 0时，时，约当消去法不能进行约当消去法不能进行*/*/flag=0;flag=0;break;break;else else for(i=0;iN;i+)for(i=0;iN;i+)r=aik/akk;r=aik/akk;if(i!=k)if(i!=k)if(aik!=0)if(aik!=0)for(j=0;j=N;j+)aij=aij-r*akj;for(j=0;j=N;j+)aij=aij-r*akj;return(flag);return(flag);第32页，本讲稿共119页void print_matric(fl

27、oat aNN+1)/*void print_matric(float aNN+1)/*输出增广矩阵输出增广矩阵*/*/int i,j;int i,j;for(i=0;iN;i+)for(i=0;iN;i+)for(j=0;j=N;j+)for(j=0;j=N;j+)printf(%12f,aij);printf(%12f,aij);printf(n);printf(n);printf(n);printf(n);第33页，本讲稿共119页main()main()float aNN+1=1,-1,1,-4,5,-4,3,-12,2,1,1,11;float aNN+1=1,-1,1,-4,5,-

28、4,3,-12,2,1,1,11;float xN;float xN;int i,k;int i,k;clrscr();clrscr();print_matric(a);/*print_matric(a);/*输出增广矩阵的初始数据输出增广矩阵的初始数据*/*/k=joydan(a);k=joydan(a);if(k=1)/*if(k=1)/*当当k=1k=1时时约当消去法正常进行约当消去法正常进行*/*/print_matric(a);/*print_matric(a);/*输出初等变换后的增广矩阵输出初等变换后的增广矩阵*/*/for(i=0;iN;i+)for(i=0;iN;i+)xi=

29、aiN/aii;xi=aiN/aii;printf(x%d=%11.6fn,i+1,xi);/*printf(x%d=%11.6fn,i+1,xi);/*输出方程组的解输出方程组的解*/*/else /*else /*当当k=0k=0时输出时输出约当消去法不能进行的信息约当消去法不能进行的信息*/*/printf(Joydan method does not run!n);printf(Joydan method does not run!n);第34页，本讲稿共119页程序运行结果：程序运行结果：1.000000 -1.000000 1.000000 -4.000000 1.000000 -

30、1.000000 1.000000 -4.000000 5.000000 -4.000000 3.000000 -12.000000 5.000000 -4.000000 3.000000 -12.000000 2.000000 1.000000 1.000000 11.000000 2.000000 1.000000 1.000000 11.000000 1.000000 0.000000 0.000000 3.000000 1.000000 0.000000 0.000000 3.000000 0.000000 1.000000 0.000000 6.000000 0.000000 1.0

31、00000 0.000000 6.000000 0.000000 0.000000 5.000000 -5.000000 0.000000 0.000000 5.000000 -5.000000 x1=3.000000 x1=3.000000 x2=6.000000 x2=6.000000 x3=-1.000000 x3=-1.000000 约当消取法能否进行，取决于每一步能否满足约当消取法能否进行，取决于每一步能否满足 ，否则约当消取法将无，否则约当消取法将无法进行。另外，即使法进行。另外，即使 ，但其绝对值很小时，仍用它作除数时，就会导致，但其绝对值很小时，仍用它作除数时，就会导致舍入误差

32、增大，严重影响计算结果的精度。为解决这一问题，可以采用选主元的技舍入误差增大，严重影响计算结果的精度。为解决这一问题，可以采用选主元的技术进行计算，具体算法可参考有关书籍。术进行计算，具体算法可参考有关书籍。第35页，本讲稿共119页4.5 4.5 三角分解法三角分解法 三三角角分分解解法法是是高高斯斯消消去去法法的的一一种种变变形形解解法法。用用高高斯斯消消去去法法去去解解 阶阶线线性性方方程程组组 ，经经过过 次次消消元元之之后后，得得出出一一个个等等价价的的上上三三角角形形方方程程组组 ，对对上上三三角角形形方方程程组组用用逐逐步步回回代代的方法就可以求出解来。的方法就可以求出解来。上述

33、消元过程也可以通过矩阵分解来实现。上述消元过程也可以通过矩阵分解来实现。第36页，本讲稿共119页4.5.1 4.5.1 三角分解法的基本思想三角分解法的基本思想 设设 阶阶线线性性方方程程组组(4.1)(4.1)的的系系数数矩矩阵阵 非非奇奇异异，并并且且系系数数矩矩阵阵 ，能分解为两个三角形矩阵之积，即能分解为两个三角形矩阵之积，即 (4.3)(4.3)其其中中 为为一一个个下下三三角角阵阵，为为一一个个上上三三角角阵阵，那那么么方方程程组组 就就化为化为 (4.4)(4.4)从而使问题转化为求解两个简单的三角形方程组：从而使问题转化为求解两个简单的三角形方程组：(4.5)(4.5)第37

34、页，本讲稿共119页 矩矩阵阵 分分解解为为两两个个矩矩阵阵 、之之积积，一一般般不不是是惟惟一一的的，但但规规定定了了 、的的形形式式之之后后，由由矩矩阵阵的的乘乘法法运运算算法法则则知知，这这种种分分解解便便是是惟惟一一的的。由由于于 和和 、的的形形式式不不同同，就就得得到到矩矩阵阵的的不不同同解解法，相应地就有解线性代数方程组的几种解法法，相应地就有解线性代数方程组的几种解法.本本节节介介绍绍直直接接三三角角分分解解法法、乔乔里里斯斯基基（CholeskyCholesky）分分解解法、追赶法和对称高斯法（法、追赶法和对称高斯法（分解法）。分解法）。第38页，本讲稿共119页4.5.2

35、4.5.2 直接三角分解法直接三角分解法 设线性方程组设线性方程组(4.1)(4.1)的系数矩阵非奇异，并且系数的系数矩阵非奇异，并且系数矩阵，能分解为同阶的单位下三角形矩阵和同阶的上矩阵，能分解为同阶的单位下三角形矩阵和同阶的上三角矩阵之积，或能分解为同阶的下三角形矩阵和同三角矩阵之积，或能分解为同阶的下三角形矩阵和同阶的单位上三角矩阵之积，即写成矩阵元素的形式如阶的单位上三角矩阵之积，即写成矩阵元素的形式如下：下：第39页，本讲稿共119页或 则称为直接三角分解法。将矩阵 分解为 时，当 是单位下三角形矩阵，是上三角矩阵时，称为Dolittle分解；当 是下三角形矩阵，是单位上三角矩阵时，

36、称为Courant分解。Dolittle分解法和Courant分解法的基本思想相同，本章只介绍Courant分解法。第40页，本讲稿共119页4.5.3 4.5.3 实现实现CourantCourant分解法的基本步骤分解法的基本步骤第41页，本讲稿共119页第42页，本讲稿共119页第43页，本讲稿共119页#define N 3 /*N#define N 3 /*N为方程组系数矩阵的阶数为方程组系数矩阵的阶数*/*/void Lu(float aNN,float LNN,float UNN)/*LUvoid Lu(float aNN,float LNN,float UNN)/*LU分解分解

37、*/*/int i,j,k;int i,j,k;for(k=0;kN;k+)for(k=0;kN;k+)for(i=k;iN;i+)/*for(i=k;iN;i+)/*计算计算L L矩阵的第矩阵的第k k列元素列元素*/*/Lik=aik;Lik=aik;for(j=0;j=k-1;j+)for(j=0;j=k-1;j+)Lik-=(Lij*Ujk);Lik-=(Lij*Ujk);for(j=k+1;jN;j+)/*for(j=k+1;jN;j+)/*计算计算U U矩阵的第矩阵的第k k列元素列元素*/*/Ukj=akj;Ukj=akj;for(i=0;i=k-1;i+)for(i=0;i=k

38、-1;i+)Ukj-=(Lki*Uij);Ukj-=(Lki*Uij);Ukj/=Lkk;Ukj/=Lkk;第44页，本讲稿共119页void solve(float bN,float LNN,float UNN,float xN)void solve(float bN,float LNN,float UNN,float xN)int j,k;int j,k;float yN;float yN;for(k=0;kN;k+)/*for(k=0;kN;k+)/*计算计算Ly=bLy=b中的中的y*/y*/yk=bk;yk=bk;for(j=0;j=k-1;j+)for(j=0;j=0;k-)/*f

39、or(k=N-1;k=0;k-)/*计算计算Ux=yUx=y中的中的x*/x*/xk=yk;xk=yk;for(j=k+1;jN;j+)for(j=k+1;jN;j+)xk=xk-Ukj*xj;xk=xk-Ukj*xj;第45页，本讲稿共119页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;jN;j+)printf(%10f,aij);

40、printf(%10f,aij);printf(%10f,bi);printf(%10f,bi);printf(n);printf(n);printf(n);printf(n);第46页，本讲稿共119页main()main()static float aNN=1,2,-1,1,-1,5,4,1,-2;static float aNN=1,2,-1,1,-1,5,4,1,-2;float bN=3,0,2;float bN=3,0,2;float xN=0,0,0;float xN=0,0,0;float UNN=1,0,0,0,1,0,0,0,1;float UNN=1,0,0,0,1,0,

41、0,0,1;float LNN;float LNN;int i;int i;zg_matric(a,b);zg_matric(a,b);Lu(a,L,U);Lu(a,L,U);solve(b,L,U,x);solve(b,L,U,x);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);第47页，本讲稿共119页程序运行结果：程序运行结果：1.000000 2.000000 -1.000000 3.000000 1.000000 2.000

42、000 -1.000000 3.000000 1.000000 -1.000000 5.000000 0.000000 1.000000 -1.000000 5.000000 0.000000 4.000000 1.000000 -2.000000 2.000000 4.000000 1.000000 -2.000000 2.000000 x1=0.2500000 x1=0.2500000 x2=1.5000000 x2=1.5000000 x3=0.2500000 x3=0.2500000第48页，本讲稿共119页 在在直直接接三三角角分分解解法法中中计计算算时时，如如果果分分母母（甚甚至至

43、可可能能为为零零）绝绝对对值值很很小小，仍仍用用它它作作除除数数时时，就就会会导导致致舍舍入入误误差差增增大大，严严重重影影响响计计算算结结果果的的精精度度。为为解解决决这这一一问问题题，可可以以配配合合按按列列选选主主元元的的技技术术进进行行计计算算，具具体体算算法法可可参参考考有有关关书书籍。籍。第49页，本讲稿共119页4.5.4 4.5.4 乔里斯基（乔里斯基（CholeskyCholesky）分解法）分解法 在在工工程程和和科科学学计计算算中中，有有时时归归结结为为求求一一个个特特殊殊的的线线性性代代数数方方程程组组，例例如如系系数数矩矩阵阵是是一一个个对对称称正正定定矩矩阵阵，为为

44、此此，这这里里介介绍绍系系数数矩矩阵阵为为对对称称正正定定矩矩阵阵的的方方程程组组的的三三角角解解法法乔乔里斯基分解法。里斯基分解法。第50页，本讲稿共119页 对对 阶阶线线性性方方程程组组 ，当当系系数数矩矩阵阵 为为对对称称正正定定时时可可以以惟惟一一地地分分解解为为 ，其其中中 为为上上三三角角矩矩阵。即写成矩阵元素的形式如下：阵。即写成矩阵元素的形式如下：且且 。第51页，本讲稿共119页4.5.5 4.5.5 实现乔里斯基分解法的基本步骤实现乔里斯基分解法的基本步骤第52页，本讲稿共119页第53页，本讲稿共119页实实现现乔乔里里斯斯基基分分解解法法的的N NS S图图第54页，

45、本讲稿共119页第55页，本讲稿共119页例例 4.5 4.5 用乔里斯基分解法求解线性方程组，式中用乔里斯基分解法求解线性方程组，式中第56页，本讲稿共119页#include#include#define N 4 /*N#define N 4 /*N为方程组系数矩阵的阶数为方程组系数矩阵的阶数*/*/void chodec(float aNN,float lNN)void chodec(float aNN,float lNN)int i,j,k;int i,j,k;float s;float s;for(i=0;iN;i+)li0=ai0/sqrt(a00);for(i=0;iN;i+)l

46、i0=ai0/sqrt(a00);for(i=1;iN;i+)for(i=1;iN;i+)s=0.0;s=0.0;for(k=0;k=i-1;k+)for(k=0;k=i-1;k+)s=s+lik*lik;s=s+lik*lik;lii=sqrt(aii-s);lii=sqrt(aii-s);for(j=i+1;jN;j+)for(j=i+1;jN;j+)s=0.0;s=0.0;for(k=0;k=i-1;k+)for(k=0;k=i-1;k+)s=s+ljk*lik;s=s+ljk*lik;lji=(aji-s)/lii;lji=(aji-s)/lii;第57页，本讲稿共119页void s

47、olve(float aNN,float bN,float xN)void solve(float aNN,float bN,float xN)int i,j,k;int i,j,k;float s;float s;for(k=0;kN-1;k+)for(k=0;kN-1;k+)bk/=akk;bk/=akk;for(i=k+1;iN;i+)for(i=k+1;i=0;k-)for(k=N-2;k=0;k-)s=bk;s=bk;for(j=k+1;jN;j+)for(j=k+1;jN;j+)s-=ajk*xj;s-=ajk*xj;xk=s/akk;xk=s/akk;第58页，本讲稿共119页v

48、oid 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;jN;j+)printf(%10f,aij);printf(%10f,aij);printf(%10f,bi);printf(%10f,bi);printf(n);printf(n);printf(n);printf(n);第59页，本讲稿共119页main()main()static flo

49、at aNN=1,2,1,-3,2,5,0,-5,1,0,14,1,static float aNN=1,2,1,-3,2,5,0,-5,1,0,14,1,-3,-5,1,15;-3,-5,1,15;float lNN;float lNN;float bN=1,2,16,8;float bN=1,2,16,8;float xN=0,0,0,0;float xN=0,0,0,0;int i;int i;zg_matric(a,b);zg_matric(a,b);Chodec(a,l);/*Chodec(a,l);/*分解过程分解过程*/*/Solve(l,b,x);/*Solve(l,b,x);

50、/*回代求方程组的解回代求方程组的解*/*/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);第60页，本讲稿共119页程序运行结果：程序运行结果：1.000000 2.000000 1.000000 -3.000000 1.000000 1.000000 2.000000 1.000000 -3.000000 1.000000 2.000000 5.000000 0.000000 -5.000000 2.000000 2.000000