LU分解程序部分全部.doc

实验报告2：线性方程组的求解姓名：祁晓 学号：08012325 班级：勘查三班姓名：惠婉婷 学号：08014122 班级：勘查三班一上机题目 已知方程组为：x1-2*x2+2*x3=-2;2*x1-3*x2-3*x3=4;4*x1+x2+6*x3=3；分别用矩阵的三角分解法和高斯列主元消去法求方程组的解二目的要求 掌握用矩阵的三角分解法和高斯列主元消去法设计程序，从而实现对线性方程组的求解。三方法原理1高斯消去法是通过逐步消元的方法把原方程组化为等价的上三角形方程组，然后回代的求解过程。高斯列主元消去法是在高斯消去法第k步时，不取akk作为主元，而是取满足|ark|=max|aik|(k=i=n)的ark作为主元，若有多个r满足则取最小的，然后再去交换第r行与第k行；最后利用高斯消元法求解。2. 矩阵的三角分解法是方程组的系数矩阵A可以分解为一个下三角阵L和一个上三角阵U的乘积，即A=LU，则AX=b为LUX=b,根据Ly=b和UX=y求出方程组的解。其中求L,U的过程为：先计算U的第一行和L的第一列U1j=a1j(1=j=n),Li1=ai1/U11;然后根据 Ukj=akj-Lkr*Urj(1=r=k-1;j=k,k+1,n)和 Lik=(aik-Lir*Urk)/Ukk(1=rk-1;i=k+1,n)计算出U的第k行和L的第k列，从而求出L和U矩阵四算法步骤（N-S流程图）1矩阵的三角分解法算法步骤如下：Step1:计算U的第一行和L的第一列 U的第一行：i=1时a1j=L11*U1j则U1j=a1j L11=1 (j=1,2,n ); L的第一列：j=1时 ai1=Li1*U11则Li1=ai1/a1j (i=2,3,n);Step2:计算U的第二行和L的第二列 U的第二行：i=2时a2j=L21*U1j+L22*U2j L22=1则U2j=a2j L21*U1j (j=2,3,n); L的第二列：j=2时 ai2=Li1*U12+Li2*U22则Li2=(ai2-Li1*U12)/U22 (i=3,4,n);Step3:假设已经进行了（k-1）步，得到了U的前k-1行和L的前k-1列，则计算U的第k行和L的第k列U的第k行：i=k时 Ukj=akj-Lkr*Urj (1=r=k-1;j=k,k+1,n)L的第k列：j=k时 Lik=(aik-Lir*Urk)/Ukk (1=rk-1;i=k+1,n)Step4:利用L，U求解方程组的解 根据A=L*U,A*X=y,从而得到 Ly=b和UX=y，然后根据高斯消元法的回代部分即可分别求出y和x2. 高斯列主元消去法算法步骤如下：Step1:选出列主元ark=max|aik| (k=i=n);Step2:判断方程组是否适合使用该方法求解若arke,则方程组不适合使用该方法求解，程序结束，否则进行下一步Step3:若r不等于k，交换r,k行Step4:高斯消元求解过程 消元过程：当k=0,i=1时，令a10=a10/a00; b1=b1-a10*b1； j=1时a11=a11-a10*a01当消元进行到第k次时，令akk-1=akk-1/ak-1k-1; bk=bk-akk-1*bk； j=k时akk=akk-akk-1*akk-1 回代求解过程：bn-1=bn-1/an-1n-1;当k=n-1时，bk=(bk-akj*bj)/akk(k+1=j=n;k=n-1,n-2,.1);返回值bk即为所求解。五C语言程序及变量说明利用矩阵的三角分解法即LU分解法求解线性方程组的解程序部分：/*矩阵的三角分解法即Lu分解法*/#include#includeint i,j,k,r;double m=0,p=0;double a33;void lu(double a33) for(i=1;i=2;i+) if(a00!=0) ai0=ai0/a00;/*先求出第一行和第一列的仍就放在a数组中*/for(k=1;k=2;k+)for(j=k;j=2;j+) for(r=0;r=k-1;r+)m=m+akr*arj; akj=akj-m;m=0; for(i=k+1;i=2;i+) for(r=0;r=k-1;r+) p=p+air*ark; aik=(aik-p)/akk; p=0;/*根据求行与列计算公式,做五个FOR循环,依次求出各行各列的数据放在a数组中*/void main()static double a33=1,-2,2,2,-3,-3,4,1,6;static double b3=-2,4,3;double c3;double d3;double f33;double m=0;double n=0;int r;int i,j;lu(a);printf(输出U的矩阵为n);for(i=0;i=2;i+)for(j=i;j=2;j+)fij=aij;printf( %f,fij);printf(n);/*为防止输出L矩阵时U被覆盖掉,先输出U矩阵,取其上三角部分,即为所求*/printf(输出L的矩阵为n);for(i=0;i=2;i+)for(j=0;j=i;j+)if(i=j)aij=1;printf( %f,aij);elseprintf( %f,aij);printf(n);c0=b0;for(i=1;i=2;i+) for(r=0;r=0;i=i-1)for(r=2;ri;r=r-1)n=n+fir*dr;di=(ci-n)/fii;n=0;printf(所求方程组解为x1=%f, x2=%f, x3=%f,d0,d1,d2);/*根据LU分解所得两个矩阵及求解步骤计算所求X一组解*/程序结束.运行结果如下:利用高斯列主元消去法求解线性方程组的解程序部分：/*高斯列主元消去法*/ #include#includedouble f(double a33,double b3,int n,int c,double e) double t,d,s;int i,r,j,k;/*通过对每列元素绝对值大小的比较选出最大的，从而选出列主元*/for(k=0;k=n-1;k+)d=fabs(akk); r=k;for(i=k+1;id) d=fabs(aik); r=i;if(de)/*判断所给方程组是否为病态的*/c=-1;break;if(r!=k)/*若r不等于k,交换k,r行*/for(j=k;j=n-1;j+) t=arj; arj=akj; akj=t; t=br;br=bk;bk=t;/*高斯消元过程*/ for(k=0;k=n-2;k+) for(i=k+1;i=n-1;i+) aik=aik/akk; bi=bi-aik*bk; for(j=k+1;j=0;k-) s=0;for(j=k+1;j=n-1;j+)s=s+akj*bj;bk=(bk-s)/akk; return (b0,b1,b2);void main()static double a33=1,-2,2,2,-3,-3,4,1,6;static double b3=-2,4,3; double e=0.01;int c;f(a,b,3,c,e);if(c=-1)printf(This method is failur);else printf(the answer isnX1=%fnX2=%fnX3=%fn,b0,b1,b2);程序结束.运行结果如下:六.运行结果及分析对使用矩阵的三角分解法分析运行结果：由矩阵的三角分解法所编程序运行结果可以看出：利用LU分解法求解方程组的解，将A分解为L和U之后，很容易就可以得方程组的解。 对使用高斯列主元消去法分析运行结果：由高斯列主元消去法所编程序运行结果可以看出：高斯列主元消去法比一般消去法多了一些比较运算，但是结果更稳定一些，该方法是目前首选的直接法。七.对程序评价 通过运用两种方法求解方程组的解，以及所得的运行结果和对其进