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

1、LU和QR法解线性方程组问题描述求解方程组一4810X2-21-7121120X3-734一二3一要求:1、编写用三角（LU）分解法求解线性方程组;2、编写用正交三角（QR）分解法求解线性方程组。二、问题分析求解线性方程组Ax=b,其实质就是把它的系数矩阵A通过各种变换成一个下三角或上三角矩阵，从而简化方程组的求解。因此，在求解线性方程组的过程中，把系数矩阵变换成上三角或下三角矩阵显得尤为重要，然而矩阵A的变换通常有两种分解方法:LU分解法和QR分解法。将A分解为一个下三角矩阵即：A=LU,1、LU分解法:一1001U11U12u1nl21100U22U2n99+,U=+a-000-Jn1ln

2、21_J0UnnL和一个上三角矩阵U,其中L=2、QR分解法:将A分解为一个正交矩阵Q和一个上三角矩阵R,即:A=QR三、实验原理1、LU分解法解Ax=b的问题就等价于要求解两个三角形方程组:Ly=b,求y;Ux=y,求x.设A为非奇异矩阵，且有分解式A=LU,1为单位下三角阵，U为上三角阵。L,U的元素可以有n步直接计算定出。用直接三角分解法解Ax=b（要求A的所有顺序主子式都不为零）的计算公式： Uli=aii（i=1,2,，n）,li|=a/5，i=2,3,，n.计算U的第r行，L的第r列元素（i=2,3，,n）：r-1 Uri=ai-ZlrkUki,i=r,r+1,n;k4r1 lir

3、=（air-£likUkr）/Urr,i=r+1,，n,且rwn.k4求解Ly=b,Ux=y的计算公式；y1=bi, iX-likyj=2,3,n:kWXn=yn/Unn,nX=（yi-%UkXk）/Uii,i=n-1,n-2,1.k文;12、QR分解法四、实验步骤1、LU分解法1将矩阵A保存进计算机中，再定义2个空矩阵L,U以便保存求出的三角矩阵的值。利用公式，将矩阵A分解为LU,L为单位下三角阵,U为上三角阵。2可知计算方法有三层循环。先通过公式计算出U矩阵的第一行元素U"和L矩阵的第一列元素八。再根据公式和，和上次的出的值，求出矩阵其余的元素，每次都要三次循环，求下一

4、个元素需要上一个结果。3先由公式，Ly=byi=h,iJV"bilikyk,i=2,3,n:k1求出y,因为L为下三角矩阵，所以由第一行开始求y.4再由公式，Ux=yxn=yn/unn,nx=(yi-'Uikxk)/Uii,i=n-1,n-2,1.k二i求出x,因为U为上三角矩阵，所以由最后一行开始求x.2、QR分解法五、程序流程图1、LU分解法输入系数矩阵A,常数项b及n/det-1K=1,n-1,1调选列主元子程序i=k+1,n,1aik-aik/akk（即求乘数mik）j=k+1,n,1aikaij-aik*akjbibiaik*bkdetakkdet消元过程回代过程返

5、回主程序2、QR分解法dakk,kk=k+1,n,1A<|aik|>|d|二选列主元输出上三角阵R=H3*H2*H1*A输入正交阵Q=-(H3*H2*H1)T解上三角阵TK、实验结果1、LU分解法矩阵鼠414；4.g0呢的2.0000098.000000?.0000004.0000008.00000U12.0000006.0000001.0000002.M00003.09000U11.0000005.000000io.mum&,00003020,000000方程组的根为：X1=1.000000x21=-la000000x3=1.000000x41=-1.000000Pres

6、sanykeytocontinue2、QR分解法七、实验总结为了求解线性方程组，我们通常需要一定的解法。其中一种解法就是通过矩阵的三角分解来实现的，属于求解线性方程组的直接法。在不考虑舍入误差下，直接法可以用有限的运算得到精确解，因此主要适用于求解中小型稠密的线性方程组。1、三角分解法三角分解法是将A矩阵分解成一个上三角形矩阵U和一个下三角形矩阵L,这样的分解法又称为LU分解法。它的用途主要在简化一个大矩阵的行列式值的计算过程，求反矩阵和求解联立方程组。不过要注意这种分解法所得到的上下三角形矩阵并非唯一，还可找到数个不同的一对上下三角形矩阵，此两三角形矩阵相乘也会得到原矩阵。2、QR分解法QR

7、分解法是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R,所以称为QR分解法。在编写这两个程序过程中，起初遇到不少麻烦！虽然课上老师反复重复着：“算法不难的，It'ssoeasy!”但是当自己实际操作时，感觉并不是那么容易。毕竟是要把实际的数学问题转化为计算机能够识别的编程算法，所以在编写程序之前我们仔细认真的把所求解的问题逐一进行详细的分析，最终转化为程序段。每当遇到问题时，大家或许有些郁闷，但最终还是静下心来反复仔细的琢磨，一一排除了错误，最终完成了本次实验。回头一想原来编个程序其实也没有想象的那么复杂，只要思路清晰，逐步分析，就可以慢慢搞定了。通过这次实验，让我们认知到团队的作用力

8、度是不容忽视的，以后不管干任何时都要注重团队合作，遇到不懂得不明白的大家一起讨论，越讨论越清晰，愈接近最优答案。这样不管干什么都能事半功倍。庆幸自己有这么个团队，也明白大家一起劳动的果实最珍贵。附源代码：LR分解法：#include<stdio.h>voidinput_A();voidinput_b();voidoutput_x(floatx4);输入矩阵A输入矩阵b/输出方程组的根voidmain()inti,k,r;floatA44=421,5,8,7210,4,836,12,6,11,20,b4=-2,-7,-7,-3,x4,u44,l44,y4;给定的方程组/input_A

9、();/input_b();/显示矢I阵A、bprintf("矩B$A44:n");for(i=0;i<4;i+)for(k=0;k<4;k+)printf("%-10f",Aik);printf("n");for(i=0;i<4;i+)u0i=A0i;for(i=0;i<4;i+)li0=Ai0/u00;l皿i=1;for(r=1;r<4;r+)计算矩阵L和U(for(i=r;i<4;i+)(floatt=0;for(k=0;k<r;k+)t=t+lrk*uki;uri=Ari-t;)for

10、(i=r;i<4;i+)(floatt1=0;for(k=0;k<r;k+)t1+=lik*ukr;lir=(Air-t1)/urr;y0=b0;for(i=1;i<4;i+)(floatt=0;for(k=0;k<i;k+)t+=lik*yk;yi=bi-t;for(i=3;i>=0;i-)(floatt=0;for(k=i+1;k<4;k+)t+=uik*xk;xi=(yi-t)/uii;printf("*n");output_x(x);输入矩阵Avoidinput_A()(floatA44;inti,j;printf("i

11、nputmatrixA44:n");for(i=0;i<4;i+)for(j=0;j<4;j+)scanf("%f",&Aij);输入矩阵bvoidinput_b()floatb4;inti;printf("inputmatrixb4:n");for(i=0;i<4;i+)scanf("%f",&bi);输出方程的根xvoidoutput_x(floatx4)inti;printf("方程组的根为:n");for(i=0;i<4;i+)printf("x%

12、d=%-13f",i+1,xi);printf("n");QR分解法：#include<stdio.h>#include<math.h>#defineN4voidmatrix_time(doubleAN,doubleBN,doubleCN,intn);两个矩阵相乘结果存在矩阵CNvoidmatrixvec(doubleAN,doubleBN,doubleCN,intn);矩阵和向量相乘,结果存在向量CNdoublevec_value(doubleA,intn);voidvec_time(doublea,doubleHN,intn);void

13、householder(double*a,doubleHN,intvoidmatrix_turn(doubleAN,intn);求向量的模两个向量相乘得一个矩阵；n,intm);求解Householder矩阵函数求矩阵装置求解上三角方程的解voidsolution(doubleAN,doubleb,intn);voidQR(doubleAN,doubleb,intn);voidmain()(doubleA44=4,2,1,5,8,7,2,10,4,8,3,6,12,6,11,20);doubleb4=-2,-7,-7,-3;inti,j;intn=4;printf("待求解的方程组系

14、数矩阵A为：n");for(i=0;i<n;i+)/显示矢I阵Afor(j=0;j<n;j+)printf("%-10f",Aij);printf("n");QR(A,b,n);两个矩阵相乘，结矩阵和向量相乘，结果存voidmatrix_time(doubleAN,doubleBN,doubleCN,intn)果存在矩阵CNinti,j,k;for(i=0;i<n;i+)for(j=0;j<n;j+)Cij=0;for(k=0;k<n;k+)Cij=Cij+Aik*Bkj;voidmatrix_vec(double

15、AN,doubleBN,doubleCN,intn)在向量CNinti,j;for(i=0;i<n;i+)Ci=0;for(j=0;j<n;j+)Ci=Ci+Aij*Bj;doublevec_value(doubleA,intn)求向量的模(doubleSum=0;inti;for(i=0;i<n;i+)Sum=Sum+Ai*Ai;Sum=sqrt(Sum);returnSum;)voidvec_time(doublea,doubleHN,intn)/两个向量相乘得一个矩阵；(inti,j;for(i=0;i<n;i+)for(j=0;j<n;j+)Hij=ai*

16、aj;)voidhouseholder(double*a,doubleHN,intn,intm)计算Householder矩阵(inti,j;doublefirst;/存放首元素doublevec_v;存放向量的模doublea1N=0;for(i=0;i<n;i+)for(j=0;j<n;j+)if(i=j)H皿=1;elseH皿=0;first=am;/取矩阵A转置的第m行(取矩阵A的第m列数)vec_v=vec_value(&am,n-m);/计算m列元素所构成向量的模am=am-vec_v;/w的分子部分vec_v=vec_value(&am,n-m);/w

17、的分母部分vec_v=vec_v*vec_v;for(i=m;i<n;i+)for(j=m;j<n;j+)Hij+=-ai*aj*2/vec_v;)voidmatrix_turn(doubleAN,intn)/求矩阵的转置doubleaNN=0;inti,j;for(i=0;i<n;i+)for(j=0;j<n;j+)aij=Aij;for(i=0;i<n;i+)for(j=0;j<n;j+)Aij=aji;voidsolution(doubleAN,doubleb,intn)求解上三角方程的解inti,j;doublexN=0;doublesum=0;fo

18、r(i=n-1;i>=0;i-)for(j=n-1;j>i;j-)sum+=Aij*xj;xi=(bi卜sum)/Aii;sum=0;printf("矩阵的QR分解求解结果为：n");for(i=0;i<n;i+)printf("X%d=%-10fn”,i+1,xi);voidQR(doubleAN,doubleb口,intn)inti,j,k,t;doubleH1NN=0,H2NN=0,H3NN=0;/H1,H2存放相邻两次的Householder矩阵，H3为中间变量矩阵doubletempNN=0;doubleQNN=0,RNN=0,Q1NN

19、=0;doubleb1N=0;for(i=0;i<n;i+)将矩阵A临时存放在temp中for(j=0;j<n;j+)tempij=Aij;/单位阵for(i=0;i<n;i+)H2ii=1;matrix_turn(temp,n);/矩阵A的转置（方便后续取A的某一列元素）for(i=0;i<n-1;i+)/Householder的一次变换householder(tempi,H1,n,i);matrix_time(H1,A,Q,n);得至UHouseholder矩阵H1矩B$H1和A相乘放在Q中for(k=0;k<n;k+)更新矩阵A,使得第一列元素除第一个元素外,其他全为0for(t=0;t<n;t+)Akt=Qkt;for(k=0;k<n;k+)for(t=0;t<n;t+)te