LU和QR分解法解线性方程组

LU和QR法解线性方程组一、问题描述求解方程组201161263841027851244321xxxx==3772,要求：1、编写用三角（LU）分解法求解线性方程组；2、编写用正交三角（QR）分解法求解线性方程组。二、问题分析求解线性方程组Ax=b，其实质就是把它的系数矩阵A通过各种变换成一个下三角或上三角矩阵，从而简化方程组的求解。因此，在求解线性方程组的过程中，把系数矩阵A变换成上三角或下三角矩阵显得尤为重要，然而矩阵A的变换通常有两种分解方法：LU分解法和QR分解法。1、LU分解法：将A分解为一个下三角矩阵L和一个上三角矩阵U，即：A=LU,其中L=10010012121nnlll，U=nnnnuuuuuu00000222112112、QR分解法：将A分解为一个正交矩阵Q和一个上三角矩阵R,即：A=QR三、实验原理1、LU分解法解Ax=b的问题就等价于要求解两个三角形方程组：⑴Ly=b,求y;⑵Ux=y,求x.设A为非奇异矩阵,且有分解式A=LU,L为单位下三角阵,U为上三角阵。L,U的元素可以有n步直接计算定出。用直接三角分解法解Ax=b（要求A的所有顺序主子式都不为零）的计算公式：①),,2,1(niaulili,11/ualilil,i=2,3,…，n.计算U的第r行，L的第r列元素(i=2,3,…,n)：②11rkkirkririulau，i=r,r+1,…,n;③rrrkkrikiriruulal/)(11,i=r+1,…,n,且r≠n.求解Ly=b，Ux=y的计算公式；④:,3,2,,1111niylbybyikkikii⑤.1,,2,1,/)(,/1nniuxuyxuyxiinikkikiinnnn2、QR分解法四、实验步骤1、LU分解法1将矩阵A保存进计算机中，再定义2个空矩阵L，U以便保存求出的三角矩阵的值。利用公式①，②，③将矩阵A分解为LU，L为单位下三角阵,U为上三角阵。2可知计算方法有三层循环。先通过公式①计算出U矩阵的第一行元素liu和L矩阵的第一列元素ill。再根据公式②和③，和上次的出的值，求出矩阵其余的元素，每次都要三次循环，求下一个元素需要上一个结果。3先由公式④，Ly=b:,3,2,,1111niylbybyikkikii求出y，因为L为下三角矩阵，所以由第一行开始求y.4再由公式⑤，Ux=y.1,,2,1,/)(,/1nniuxuyxuyxiinikkikiinnnn求出x,因为U为上三角矩阵，所以由最后一行开始求x.2、QR分解法五、程序流程图1、LU分解法开始输入系数矩阵A，常数项b及ndet←1K=1,n-1,1调选列主元子程序i=k+1,n,1aik←aik/akk(即求乘数mik)j=k+1,n,1aik←aij-aik*akjjbi←bi←aik*bkidet←akkdetkbn←bn/am(即求出xn)i=n-1,1,-1s←0j=i+1,n,1s←s+aij*bjjbi←(bi-s)/aij输出b及det结束idet←amndet回代过程消元过程由主程序转来d←akk,p←ki=k+1,n,1|aik||d|d←aik,p←iid=0P=kj=k,n,1t←apj,apj←akj,akj←tjt←bp,bp←bk,bk←tdet←det返回主程序返回主程序输出失败标志，|A|=0det←0结束选列主元≤=≠=2、QR分解法开始输入数据A，b矩阵A的转置Householder变换转置矩阵相乘输出上三角阵R=H3*H2*H1*A输入正交阵Q=-（H3*H2*H1）T解上三角阵输出结果结束六、实验结果1、LU分解法2、QR分解法七、实验总结为了求解线性方程组，我们通常需要一定的解法。其中一种解法就是通过矩阵的三角分解来实现的，属于求解线性方程组的直接法。在不考虑舍入误差下，直接法可以用有限的运算得到精确解，因此主要适用于求解中小型稠密的线性方程组。1、三角分解法三角分解法是将A矩阵分解成一个上三角形矩阵U和一个下三角形矩阵L，这样的分解法又称为LU分解法。它的用途主要在简化一个大矩阵的行列式值的计算过程，求反矩阵和求解联立方程组。不过要注意这种分解法所得到的上下三角形矩阵并非唯一，还可找到数个不同的一对上下三角形矩阵，此两三角形矩阵相乘也会得到原矩阵。2、QR分解法QR分解法是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R,所以称为QR分解法。在编写这两个程序过程中，起初遇到不少麻烦！虽然课上老师反复重复着：“算法不难的，It'ssoeasy!”但是当自己实际操作时，感觉并不是那么容易。毕竟是要把实际的数学问题转化为计算机能够识别的编程算法，所以在编写程序之前我们仔细认真的把所求解的问题逐一进行详细的分析，最终转化为程序段。每当遇到问题时，大家或许有些郁闷，但最终还是静下心来反复仔细的琢磨，一一排除了错误，最终完成了本次实验。回头一想原来编个程序其实也没有想象的那么复杂，只要思路清晰，逐步分析，就可以慢慢搞定了。通过这次实验，让我们认知到团队的作用力度是不容忽视的，以后不管干任何时都要注重团队合作，遇到不懂得不明白的大家一起讨论，越讨论越清晰，愈接近最优答案。这样不管干什么都能事半功倍。庆幸自己有这么个团队，也明白大家一起劳动的果实最珍贵。附源代码:LR分解法：#includestdio.hvoidinput_A();//输入矩阵Avoidinput_b();//输入矩阵bvoidoutput_x(floatx[4]);//输出方程组的根voidmain(){inti,k,r;floatA[4][4]={{4,2,1,5},{8,7,2,10},{4,8,3,6},{12,6,11,20}},b[4]={-2,-7,-7,-3},x[4],u[4][4],l[4][4],y[4];//给定的方程组//input_A();//input_b();//显示矩阵A、bprintf(矩阵A[4][4]:\n);for(i=0;i4;i++){for(k=0;k4;k++)printf(%-10f,A[i][k]);printf(\n);}for(i=0;i4;i++)u[0][i]=A[0][i];for(i=0;i4;i++){l[i][0]=A[i][0]/u[0][0];l[i][i]=1;}for(r=1;r4;r++)//计算矩阵L和U{for(i=r;i4;i++){floatt=0;for(k=0;kr;k++)t=t+l[r][k]*u[k][i];u[r][i]=A[r][i]-t;}for(i=r;i4;i++){floatt1=0;for(k=0;kr;k++)t1+=l[i][k]*u[k][r];l[i][r]=(A[i][r]-t1)/u[r][r];}}y[0]=b[0];for(i=1;i4;i++){floatt=0;for(k=0;ki;k++)t+=l[i][k]*y[k];y[i]=b[i]-t;}for(i=3;i=0;i--){floatt=0;for(k=i+1;k4;k++)t+=u[i][k]*x[k];x[i]=(y[i]-t)/u[i][i];}printf(*****************************\n);output_x(x);}//输入矩阵Avoidinput_A(){floatA[4][4];inti,j;printf(inputmatrixA[4][4]:\n);for(i=0;i4;i++)for(j=0;j4;j++)scanf(%f,&A[i][j]);}//输入矩阵bvoidinput_b(){floatb[4];inti;printf(inputmatrixb[4]:\n);for(i=0;i4;i++)scanf(%f,&b[i]);}//输出方程的根xvoidoutput_x(floatx[4]){inti;printf(方程组的根为:\n);for(i=0;i4;i++)printf(x[%d]=%-13f,i+1,x[i]);printf(\n);}QR分解法：#includestdio.h#includemath.h#defineN4voidmatrix_time(doubleA[][N],doubleB[][N],doubleC[][N],intn);//两个矩阵相乘,结果存在矩阵C[][N]voidmatrix_vec(doubleA[][N],doubleB[N],doubleC[N],intn);//矩阵和向量相乘,结果存在向量C[N]doublevec_value(doubleA[],intn);//求向量的模voidvec_time(doublea[],doubleH[][N],intn);//两个向量相乘得一个矩阵；voidhouseholder(double*a,doubleH[][N],intn,intm);//求解Householder矩阵函数voidmatrix_turn(doubleA[][N],intn);//求矩阵装置voidsolution(doubleA[][N],doubleb[],intn);//求解上三角方程的解voidQR(doubleA[][N],doubleb[],intn);voidmain(){doubleA[4][4]={{4,2,1,5},{8,7,2,10},{4,8,3,6},{12,6,11,20}};doubleb[4]={-2,-7,-7,-3};inti,j;intn=4;printf(待求解的方程组系数矩阵A为：\n);for(i=0;in;i++)//显示矩阵A{for(j=0;jn;j++)printf(%-10f,A[i][j]);printf(\n);}QR(A,b,n);}voidmatrix_time(doubleA[][N],doubleB[][N],doubleC[][N],intn)//两个矩阵相乘，结果存在矩阵C[][N]{inti,j,k;for(i=0;in;i++)for(j=0;jn;j++){C[i][j]=0;for(k=0;kn;k++)C[i][j]=C[i][j]+A[i][k]*B[k][j];}}voidmatrix_vec(doubleA[][N],doubleB[N],doubleC[N],intn)//矩阵和向量相乘,结果存在向量C[N]{inti,j;for(i=0;in;i++){C[i]=0;for(j=0;jn;j++)C[i]=C[i]+A[i][j]*B[j];}}doublevec_value(doubleA[],intn)//求向量的模{doubleSum=0;inti;for(i=0;in;i++)Sum=Sum+A[i]*A[i];Sum=sqrt(Sum);returnSum;}voidvec_time(doublea[],doubleH[][N],intn)//两个向量相乘得一个矩阵；{inti,j;for(i=0;in;i++)for(j=0;jn;j++)H[i][j]=a[i]*a[j];}voidhouseholder(double*a,doubleH[][N],intn,intm)//计算Householder矩阵{inti,j;doublefirst;//存放首元素doublevec_v;//存放向量的模doublea1[N]={0};for(i=0;in;i++)for(j=0;jn;j++){if(i==j)H[i][j]=1;elseH[i][j]=0;}first=a[m];//取矩阵