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

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

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

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

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

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

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

6、）结束。结束。本讲稿第八页，共一百一十九页本讲稿第九页，共一百一十九页例例4.1 4.1 用高斯消去法求解线性方程组，式中用高斯消去法求解线性方程组，式中本讲稿第十页，共一百一十九页#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

7、=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*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);本讲稿第十一页，共一百一十九页void zg_matric(float aNN,float bN)/*void zg_matric(float aNN,float b

8、N)/*输出增广矩阵输出增广矩阵*/*/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;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,

9、xi);本讲稿第十三页，共一百一十九页程序运行结果：程序运行结果：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.541

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

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

12、，简称主元）进行消去。本讲稿第十五页，共一百一十九页 列列 主主 元元 高高 斯斯 消消 去去 法法，是是 在在 高高 斯斯 消消 去去 法法 的的 消消 去去 过过 程程 中中，第第 步步()增增 加加 选选 主主 元元 操操 作作，成成 为为：首首 先先 在在 第第 列列 中中，从从 中中选选出出绝绝对对值值最最大大的的元元素素，经经过过行行交交换换(把把绝绝对对值值最最大大者者所所在在的的行行与与第第 行行交交换换)把把绝绝对对值值最最大大的的元元素素交交换换到到 的的单单元元中中；然然后后作作GaussGauss消消去去法法的的消消去去过过程程中中第第 步步，初初等等变变换换产产生生第

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

14、素 对对 ，从，从 中选出绝对值最大的元素中选出绝对值最大的元素 ，将将 行和行和 行交换后，再作第行交换后，再作第 次消元操作。次消元操作。（3 3）消元过程：）消元过程：对对 ,计算：计算：本讲稿第十七页，共一百一十九页（4）回代过程：（5）输出方程组的解。（6）结束本讲稿第十八页，共一百一十九页本讲稿第十九页，共一百一十九页本讲稿第二十页，共一百一十九页例例 4.2 4.2 用列主元高斯消去法求解线性方程组用列主元高斯消去法求解线性方程组 ，式中，式中本讲稿第二十一页，共一百一十九页void ColGauss(float aNN,float bN)void ColGauss(float

15、aNN,float 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)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=b

16、l;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*bi;bj=bj+t*bi;for(k=i;kN;k+)for(k=i;kN;k+)ajk=ajk+t*aik;ajk=ajk+t*aik;本讲稿第二十二页，共一百一十九页void zg_matric(float aNN,float bN)/*void zg_matric(f

17、loat 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#define N 4 /*N#define N 4 /*N为方程组系数矩阵的阶数为方程组系数矩阵的阶数*/

18、*/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 bN=3.005,5.407,0.136,3.772;float bN=3.005,5.407,0.136,3.772

19、;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+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;

20、i+)/*for(i=0;iN;i+)/*输出方程组的解输出方程组的解*/*/printf(x%d=%11.7fn,i+1,xi);printf(x%d=%11.7fn,i+1,xi);本讲稿第二十四页，共一百一十九页程序运行结果：程序运行结果：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.

21、407000 4.329000 -1.952000 0.006000 2.087000 0.136000 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本讲稿第二十五页，共

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

23、局限性，产生了全主元高斯消去法。去法。列主元高斯消去法当变换到第列主元高斯消去法当变换到第 步时，从系数矩阵的右下角步时，从系数矩阵的右下角 阶子矩阵中选取绝对值最大的元素，然后通过行变换与列变换将阶子矩阵中选取绝对值最大的元素，然后通过行变换与列变换将它交换到主元素它交换到主元素 的位置上。回代过程与的位置上。回代过程与GaussGauss消去法的回代过消去法的回代过程相同，请读者参考有关书籍。程相同，请读者参考有关书籍。本讲稿第二十六页，共一百一十九页4.4 4.4 约当（约当（Jordan）消去法）消去法4.4.1 4.4.1 约当消去法的基本思想约当消去法的基本思想 约约当当（Jord

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

25、的次数）圆括号内的数字表示消去计算的次数）由此可立即求得方程组的解：由此可立即求得方程组的解：本讲稿第二十八页，共一百一十九页4.4.2 4.4.2 实现约当消去法的基本步骤实现约当消去法的基本步骤本讲稿第二十九页，共一百一十九页本讲稿第三十页，共一百一十九页例例 4.3 4.3 用约当消去法求解线性方程组用约当消去法求解线性方程组 ，式中，式中，本讲稿第三十一页，共一百一十九页#define N 3 /*N#define N 3 /*N为方程组的系数矩阵的阶数为方程组的系数矩阵的阶数*/*/int joydan(float aNN+1)int joydan(float aNN+1)float

26、 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+)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+)ai

27、j=aij-r*akj;return(flag);return(flag);本讲稿第三十二页，共一百一十九页void print_matric(float 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);本讲稿第三十三页，共一百一十九页ma

28、in()main()float aNN+1=1,-1,1,-4,5,-4,3,-12,2,1,1,11;float aNN+1=1,-1,1,-4,5,-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);/*

29、print_matric(a);/*输出初等变换后的增广矩阵输出初等变换后的增广矩阵*/*/for(i=0;iN;i+)for(i=0;iN;i+)xi=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

30、);本讲稿第三十四页，共一百一十九页程序运行结果：程序运行结果：1.000000 -1.000000 1.000000 -4.000000 1.000000 -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.00

31、0000 0.000000 0.000000 3.000000 0.000000 1.000000 0.000000 6.000000 0.000000 1.000000 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、法将无法进，否则约当消取法将无法进行。另外，即使行。另外，即使 ，但其绝对值很小时，仍用它作除数时，就会导致舍入误差，但其绝对值很小时，仍用它作除数时，就会导致舍入误差增大，严重影响计算结果的精度。为解决这一问题，可以采用选主元的技术进行计算，增大，严重影响计算结果的精度。为解决这一问题，可以采用选主元的技术进行计算，具体算法可参考有关书籍。具体算法可参考有关书籍。本讲稿第三十五页，共一百一十九页4.5 4.5 三角分解法三角分解法 三三角角分分解解法法是是高高斯斯消消去去法法的的一一种种变变形形解解法法。用用高高斯斯消消去去法法去去解解 阶阶线线性性方方程程组组 ，经经过过 次次消消元元之之

33、后后，得得出出一一个个等等价价的的上上三三角角形形方方程程组组 ，对对上上三三角角形形方方程程组组用用逐逐步步回代的方法就可以求出解来。回代的方法就可以求出解来。上述消元过程也可以通过矩阵分解来实现。上述消元过程也可以通过矩阵分解来实现。本讲稿第三十六页，共一百一十九页4.5.1 4.5.1 三角分解法的基本思想三角分解法的基本思想 设设 阶阶线线性性方方程程组组(4.1)(4.1)的的系系数数矩矩阵阵 非非奇奇异异，并并且且系系数数矩阵矩阵 ，能分解为两个三角形矩阵之积，即，能分解为两个三角形矩阵之积，即 (4.3)(4.3)其其中中 为为一一个个下下三三角角阵阵，为为一一个个上上三三角角阵

34、阵，那那么么方方程程组组 就化为就化为 (4.4)(4.4)从而使问题转化为求解两个简单的三角形方程组：从而使问题转化为求解两个简单的三角形方程组：(4.5)(4.5)本讲稿第三十七页，共一百一十九页 矩矩阵阵 分分解解为为两两个个矩矩阵阵 、之之积积，一一般般不不是是惟惟一一的的，但但规规定定了了 、的的形形式式之之后后，由由矩矩阵阵的的乘乘法法运运算算法法则则知知，这这种种分分解解便便是是惟惟一一的的。由由于于 和和 、的的形形式式不不同同，就就得得到到矩矩阵阵的的不不同同解法，相应地就有解线性代数方程组的几种解法解法，相应地就有解线性代数方程组的几种解法.本本节节介介绍绍直直接接三三角角

35、分分解解法法、乔乔里里斯斯基基（CholeskyCholesky）分分解解法法、追赶法和对称高斯法（追赶法和对称高斯法（分解法）。分解法）。本讲稿第三十八页，共一百一十九页4.5.2 4.5.2 直接三角分解法直接三角分解法 设线性方程组设线性方程组(4.1)(4.1)的系数矩阵非奇异，并且系数矩的系数矩阵非奇异，并且系数矩阵，能分解为同阶的单位下三角形矩阵和同阶的上三角阵，能分解为同阶的单位下三角形矩阵和同阶的上三角矩阵之积，或能分解为同阶的下三角形矩阵和同阶的单矩阵之积，或能分解为同阶的下三角形矩阵和同阶的单位上三角矩阵之积，即写成矩阵元素的形式如下：位上三角矩阵之积，即写成矩阵元素的形式

36、如下：本讲稿第三十九页，共一百一十九页或 则称为直接三角分解法。将矩阵 分解为 时，当 是单位下三角形矩阵，是上三角矩阵时，称为Dolittle分解；当 是下三角形矩阵，是单位上三角矩阵时，称为Courant分解。Dolittle分解法和Courant分解法的基本思想相同，本章只介绍Courant分解法。本讲稿第四十页，共一百一十九页4.5.3 4.5.3 实现实现CourantCourant分解法的基本步骤分解法的基本步骤本讲稿第四十一页，共一百一十九页本讲稿第四十二页，共一百一十九页本讲稿第四十三页，共一百一十九页#define N 3 /*N#define N 3 /*N为方程组系数矩阵

37、的阶数为方程组系数矩阵的阶数*/*/void Lu(float aNN,float LNN,float UNN)/*LUvoid Lu(float aNN,float LNN,float UNN)/*LU分解分解*/*/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);fo

38、r(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-1;i+)Ukj-=(Lki*Uij);Ukj-=(Lki*Uij);Ukj/=Lkk;Ukj/=Lkk;本讲稿第四十四页，共一百一十九页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;f

39、loat 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-)/*for(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;本讲稿第四十五页，共一百一十九页void zg_matric(float aNN,float bN)/*void zg_matric(float aNN,floa

40、t 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);本讲稿第四十六页，共一百一十九页main()main()static float aNN=1,2,-1,1,-1,5,4,1,-2;static float aNN=1,2,-1,1,-1,5,4,

41、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,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.7f

42、n,i+1,xi);printf(x%d=%11.7fn,i+1,xi);本讲稿第四十七页，共一百一十九页程序运行结果：程序运行结果：1.000000 2.000000 -1.000000 3.000000 1.000000 2.000000 -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.250000

43、0 x1=0.2500000 x2=1.5000000 x2=1.5000000 x3=0.2500000 x3=0.2500000本讲稿第四十八页，共一百一十九页 在在直直接接三三角角分分解解法法中中计计算算时时，如如果果分分母母（甚甚至至可可能能为为零零）绝绝对对值值很很小小，仍仍用用它它作作除除数数时时，就就会会导导致致舍舍入入误误差差增增大大，严严重重影影响响计计算算结结果果的的精精度度。为为解解决决这这一一问问题题，可可以以配配合合按按列列选选主主元元的的技技术术进进行行计计算算，具具体体算算法法可可参参考考有有关关书书籍。籍。本讲稿第四十九页，共一百一十九页4.5.4 4.5.4

44、乔里斯基（乔里斯基（CholeskyCholesky）分解法）分解法 在在工工程程和和科科学学计计算算中中，有有时时归归结结为为求求一一个个特特殊殊的的线线性性代代数数方方程程组组，例例如如系系数数矩矩阵阵是是一一个个对对称称正正定定矩矩阵阵，为为此此，这这里里介介绍绍系系数数矩矩阵阵为为对对称称正正定定矩矩阵阵的的方方程程组组的的三三角角解解法法乔乔里斯基分解法。里斯基分解法。本讲稿第五十页，共一百一十九页 对对 阶阶线线性性方方程程组组 ，当当系系数数矩矩阵阵 为为对对称称正正定定时时可可以以惟惟一一地地分分解解为为 ，其其中中 为为上上三角矩阵。即写成矩阵元素的形式如下：三角矩阵。即写成

45、矩阵元素的形式如下：且且 。本讲稿第五十一页，共一百一十九页4.5.5 4.5.5 实现乔里斯基分解法的基本步骤实现乔里斯基分解法的基本步骤本讲稿第五十二页，共一百一十九页本讲稿第五十三页，共一百一十九页实实现现乔乔里里斯斯基基分分解解法法的的N NS S图图本讲稿第五十四页，共一百一十九页本讲稿第五十五页，共一百一十九页例例 4.5 4.5 用乔里斯基分解法求解线性方程组，式中用乔里斯基分解法求解线性方程组，式中本讲稿第五十六页，共一百一十九页#include#include#define N 4 /*N#define N 4 /*N为方程组系数矩阵的阶数为方程组系数矩阵的阶数*/*/voi

46、d 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+)li0=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

47、;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;本讲稿第五十七页，共一百一十九页void solve(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/=a

48、kk;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;本讲稿第五十八页，共一百一十九页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

49、=0;jN;j+)printf(%10f,aij);printf(%10f,aij);printf(%10f,bi);printf(%10f,bi);printf(n);printf(n);printf(n);printf(n);本讲稿第五十九页，共一百一十九页main()main()static float 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=

50、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);/*回代求方程组的解回代求方程组的解*/*/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);本讲稿第六十页，共一百一十九页程序运行结果