LU解线性方程组

1/12LU分解法解线性方程组一．实验目的1.熟练运用已学计算方法求解方程组2.加深对计算方法技巧，选择正确的计算方法来求解各种方程组3.培养使用电子计算机进行科学计算和解决问题的能力二．实验环境VC++6.0语言：C++三．实验内容编写一个解线性方程组的LU算法。算法包括：1.矩阵的输入输出2.判断矩阵是否奇异3.进行分解4.输出L、U5.回带解出y、x四．实验原理1.若一个线性方程组系数矩阵为n阶方阵A且各阶顺序主子式均不为零，则A的LU分解存在且唯一。2.在满足1的条件下可推导得出以下公式(1):11()jijijikijijklaluu11jjijijkkikualu(2):11iiiikkkybly1()niiikkiikixyuxu2/123.公式(1)用于求解矩阵L、U，公式(2)用于回带求解y、x。从公式中可以看出：L对角线上元素为1，U第一行与A第一行相同。4.LU分解的具体过程和计算顺序如下：(1)第一步分解：1111ua(2)第二步分解：212111lau1212ua22222112ualu(3)第三步分解：313111lau3232311222()laluu1313ua23232113ualu333331133223ualulu(n).第n步分解：依次计算：1.1..nnninnnlluu1112121222121.1.......1nnnnnnuuuluullu五.实验流程图3/121.程序总的流程图初始化系数矩阵与右端项开始程序判断矩阵是否奇异做l-u分解输出无解信息若非奇异回带解出y、x并输出程序结束2.考虑主要目的这里只给出LU分解的流程图：4/12i=1i+n?回带求值语句J=1J=i-1?Temp=0K=1K=J-1Temp=Temp+Aik*AkjK=++是是是Lij=（Aij-Temp）/UjjJ++否M=0P《=j-1M=M+Ajp*ApiP++是P=1i++否，转到其他语句Uji=Aji-M5/12六．VC++源程序#includeiostreamusingnamespacestd;doubleFun(intn1,doublea1[11][11]);//声明求矩阵行列式的递归函数函数doubleFun(intn1,doublea1[11][11]){inti_1,j_1,c;//c为数组b的行doubleb[11][11];//用于存放余子式intp=0,q=0;intsum=0;if(n1==1)returna1[1][1];for(i_1=1;i_1=n1;i_1++){for(c=1;c=n1-1;c++){if(ci_1)p=0;elsep=1;for(j_1=1;j_1=n1-1;j_1++){b[c][j_1]=a1[c+p][j_1+1];}}if((i_1-1)%2==0)q=1;elseq=(-1);sum=sum+a1[i_1][1]*q*Fun(n1-1,b);}returnsum;}//主程序，用于矩阵输入输出，判断，分解，回带voidmain(){intn,i,j,k;doubletemp=1.0;doublea[11][11]={0};doubleL[11][11]={0};6/12doubleU[11][11]={0};doubleb[11]={1};doublex[11]={1};doublem=0.0;doubley[11]={1};doubletemp1=0.0;doubletemp2=0.0;//完成系数矩阵与右端项的输入与输出cout!!!请输入待解方程组未知数的个数:nendl;cinn;cout!!!请输入待解方程组的系数endl;for(i=1;i=n;i++){for(intj=1;j=n;j++){cina[i][j];}}cout!!!请输入待解方程组的系数矩阵如下:endl;for(i=1;i=n;i++){for(j=1;j=n;j++){couta[i][j];}coutendl;}cout!!!请输入待解方程组的右端项endl;//初始化LUfor(intp=1;p=n;p++){L[p][p]=1;}for(p=1;p=n;p++){7/12U[1][p]=a[1][p];}coutendl;for(i=1;i=n;i++){cinb[i];}coutb=endl;for(i=1;i=n;i++){coutb[i]endl;}for(i=2;i=n;i++){temp=Fun(i,a);if(temp==0){coutendlendl;coutendloutput：temp=tempendl;cout矩阵奇异;break;}}coutendl;//LU分解if(temp!=0){for(i=1;i=n;i++){for(j=1;j=i-1;j++){if(a[j][j]==0)cout请重新输入系数矩阵与右端项：endl;elsem=0.0;8/12for(k=1;k=j-1;k++){m=m+L[i][k]*U[k][j];}L[i][j]=(a[i][j]-m)/U[j][j];}for(j=1;j=i;j++){m=0.0;for(k=1;k=j-1;k++){m=m+L[j][k]*U[k][i];}U[j][i]=a[j][i]-m;}}}if(temp!=0){coutL=endl;for(i=1;i=n;i++){for(j=1;j=n;j++){coutL[i][j];}coutendlendl;}coutendlendlU=endl;for(i=1;i=n;i++){for(j=1;j=n;j++){coutU[i][j];}coutendlendl;}9/12//回带求解y[]couty=endl;for(ints=1;s=n;s++){for(intt=1;t=s-1;t++){temp1=temp1+L[s][t]*y[t];}y[s]=b[s]-temp1;couty[s]endl;temp1=0.0;}coutendl;//回带求解x[]x[1]=y[1];for(s=n;s=1;s--){for(intt=s+1;t=n;t++){temp2=temp2+U[s][t]*x[t];}x[s]=(y[s]-temp2)/U[s][s];temp2=0.0;}coutx=endl;for(i=1;i=n;i++)coutx[i]endl;}}10/12七．程序执行结果1.系数矩阵主子式有为0时，随便输入一个矩阵结果如下：2.系数矩阵各阶主子式不为0时，输入课本70页习题12所示矩阵系数5791068109710895765A与右端项2618229b，顺序执行程序，结果经matlab验证无误，结果如下：11/12验证：12/12七．实验优缺点此算法可以一定程度上判断是否可解并解出其解，但由于未添加主元选取的步骤，导致对于有的方程组会判断失误并认为方程无解。八．参考文献1.张军等编著.数值计算.清华大学出版社2008年7月第一版2.罗建军等编著.C++程序设计教程（第二版）.高等教育出版社2007年8月第二版3.施吉林刘淑珍陈桂芝.计算机数值方法[M].北京:高等教育出版社,1999