实验六报告6LU分解

实验六解线性方程组的直接法——LU分解实验目的:1、了解LU分解的数学依据与基本原理；2、能熟练地应用、构造矩阵的LU分解程序求解方程组。实验内容:根据LU分解的原理与公式自己填写实验题目：（1）编写程序写出下面系数矩阵的LU分解223457245（2）编程序，分别用LU分解法与主元素消去法解方程组第一题：1231231234212223523144xxxxxxxxx#includestdio.h#includestdlib.hvoidmain(){floata[3][3]={{2,2,3},{4,5,7},{-2,4,5}};floatL[3][3],U[3][3];intn=3;intk,i,j;floats,t;for(j=0;jn;j++)a[0][j]=a[0][j];for(i=1;in;i++)a[i][0]=a[i][0]/a[0][0];for(k=1;kn;k++){for(j=k;jn;j++){s=0;for(i=0;ik;i++)s=s+a[k][i]*a[i][j];a[k][j]=a[k][j]-s;}for(i=k+1;in;i++){t=0;for(j=0;jk;j++)t=t+a[i][j]*a[j][k];a[i][k]=(a[i][k]-t)/a[k][k];}}for(i=0;in;i++)for(j=0;jn;j++){if(ij){L[i][j]=a[i][j];U[i][j]=0;}else{U[i][j]=a[i][j];if(i==j)L[i][j]=1;elseL[i][j]=0;}}printf(L[3][3]=\n);for(i=0;in;i++){for(j=0;jn;j++)printf(%0.0f,L[i][j]);printf(\n);}printf(U[3][3]=\n);for(i=0;in;i++){for(j=0;jn;j++)printf(%0.0f,U[i][j]);printf(\n);}}第二题：#includestdio.h#includestdlib.h#defineN10//矩阵大小范围//求x回调函数floatgetmx(floata[N][N],floatx[N],inti,intn){floatmx=0;intr;for(r=i+1;rn;r++){mx+=a[i][r]*x[r];}returnmx;}//求y回调函数floatgetmy(floata[N][N],floaty[N],inti,intn){floatmy=0;intr;for(r=0;rn;r++){if(i!=r)my+=a[i][r]*y[r];}returnmy;}//计算xfloatgetx(floata[N][N],floatb[N],floatx[N],inti,intn){floatresult;if(i==n-1)//计算最后一个x的值result=(float)(b[i]/a[n-1][n-1]);elseresult=(float)((b[i]-getmx(a,x,i,n))/a[i][i]);returnresult;}//计算yfloatgety(floata[N][N],floatb[N],floaty[N],inti,intn){floatresult;if(i==0)//计算第一个y的值result=float(b[i]/a[i][i]);elseresult=float((b[i]-getmy(a,y,i,n))/a[i][i]);returnresult;}voidmain(){floatl[N][N]={0};//定义L矩阵floatu[N][N]={0};//定义U矩阵floaty[N]={0};//定义数组Yfloatx[N]={0};//定义数组Xfloata[N][N];//定义系数矩阵floatb[N];//定义右端项floatsum=0;inti,j,k;intn;intflag=1;while(flag){printf(请输入系数矩阵的大小：);scanf(%d,&n);if(nN){printf(矩阵过大！\n);continue;}flag=0;}printf(请输入系数矩阵值：\n);for(i=1;in+1;i++){for(j=1;jn+1;j++){printf(a[%d][%d]:,i,j);scanf(%f,&a[i][j]);}}printf(请输入右端项数组：\n);for(i=1;in+1;i++){printf(b[%d]:,i);scanf(%f,&b[i]);}//初始化矩阵lfor(i=0;in;i++){for(j=0;jn;j++){if(i==j)l[i][j]=1;}}//LU分解for(i=0;in;i++){u[0][i]=(float)(a[0][i]/l[0][0]);}for(i=0;in-1;i++){for(j=i+1;jn;j++){for(k=0,sum=0;kn;k++){if(k!=i)sum+=l[j][k]*u[k][i];}l[j][i]=(float)((a[j][i]-sum)/u[i][i]);}for(j=i+1;jn;j++){for(k=0,sum=0;kn;k++){if(k!=i+1)sum+=l[i+1][k]*u[k][j];}u[i+1][j]=(float)((a[i+1][j]-sum));}}//回代方式计算数组for(i=0;in;i++){y[i]=gety(l,b,y,i,n);}//回代方式计算数组Xfor(i=n-1;i=0;i--){x[i]=getx(u,y,x,i,n);}//显示数组X即结果printf(\n\n方程组解X：\n);for(i=0;in;i++){printf(x%d=%0.3f\n,i+1,x[i]);}}