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

《LU和QR分解法解线性方程组》由会员分享，可在线阅读，更多相关《LU和QR分解法解线性方程组（19页珍藏版）》请在装配图网上搜索。

1、LU和QR法解线性方程组 一、 问题描述求解方程组=,规定:1、 编写用三角（LU）分解法求解线性方程组；2、 编写用正交三角（R)分解法求解线性方程组。二、 问题分析求解线性方程组A=，其实质就是把它的系数矩阵A通过多种变换成一种下三角或上三角矩阵，从而简化方程组的求解。因此，在求解线性方程组的过程中,把系数矩阵变换成上三角或下三角矩阵显得尤为重要，然而矩阵A的变换一般有两种分解措施：LU分解法和QR分解法。1、 LU分解法:将A分解为一种下三角矩阵L和一种上三角矩阵，即：ALU，其中 L=, 2、 QR分解法:将A分解为一种正交矩阵Q和一种上三角矩阵R，即:A=Q三、实验原理1、 L分解法

2、解Ax=b 的问题就等价于规定解两个三角形方程组：y=b,求y; x=y,求x.设A为非奇异矩阵,且有分解式A=U,为单位下三角阵，为上三角阵。L，的元素可以有n步直接计算定出。用直接三角分解法解Ax=b（规定的所有顺序主子式都不为零)的计算公式: , ，i,，n.计算U的第r行,L的第r列元素（i2,3,n)： ,i=r,r+1，,n； , ir+1，,且rn.求解Ly=b，U=y的计算公式; 2、 Q分解法四、 实验环节1、LU分解法1将矩阵保存进计算机中,再定义2个空矩阵L，U以便保存求出的三角矩阵的值。运用公式，将矩阵A分解为LU,L为单位下三角阵,U为上三角阵。2可知计算措施有三层循

3、环。先通过公式计算出U矩阵的第一行元素和L矩阵的第一列元素。再根据公式和,和上次的出的值，求出矩阵其他的元素,每次都要三次循环,求下一种元素需要上一种成果。3先由公式,Ly 求出y，由于L为下三角矩阵,因此由第一行开始求.4再由公式，Ux求出x, 由于U为上三角矩阵,因此由最后一行开始求x.2、QR分解法五、 程序流程图1、L分解法2、Q分解法六、 实验成果1、 L分解法2、QR分解法七、 实验总结为了求解线性方程组,我们一般需要一定的解法。其中一种解法就是通过矩阵的三角分解来实现的，属于求解线性方程组的直接法。在不考虑舍入误差下，直接法可以用有限的运算得到精确解，因此重要合用于求解中小型稠密

4、的线性方程组。1、三角分解法 三角分解法是将A矩阵分解成一种上三角形矩阵U和一种下三角形矩阵L，这样的分解法又称为分解法。它的用途重要在简化一种大矩阵的行列式值的计算过程，求反矩阵和求解联立方程组。但是要注意这种分解法所得到的上下三角形矩阵并非唯一,还可找到数个不同 的一对上下三角形矩阵,此两三角形矩阵相乘也会得到原矩阵。 2、 分解法 R分解法是将矩阵分解成一种正规正交矩阵Q与上三角形矩阵R,因此称为QR分解法。 在编写这两个程序过程中,起初遇到不少麻烦!虽然课上教师反复反复着:“算法不难的，Its esy!”但是当自己实际操作时,感觉并不是那么容易。毕竟是要把实际的数学问题转化为计算机可以

5、辨认的编程算法,因此在编写程序之前我们仔细认真的把所求解的问题逐个进行具体的分析，最后转化为程序段。每当遇到问题时,人们或许有些郁闷,但最后还是静下心来反复仔细的揣摩,一一排除了错误，最后完毕了本次实验。 回头一想本来编个程序其实也没有想象的那么复杂,只要思路清晰，逐渐分析，就可以慢慢搞定了。通过这次实验,让我们认知到团队的作用力度是不容忽视的,后来不管干任何时都要注重团队合伙,遇到不懂得不明白的人们一起讨论，越讨论越清晰，愈接近最优答案。这样不管干什么都能事半功倍。庆幸自己有这样个团队，也明白人们一起劳动的果实最贵重。附 源代码:LR分解法：inludevod input_A();/输入矩阵

6、Avoid inp_b（）；/输入矩阵bvdoutputx(float x4);/输出方程组的根oid i()it i,k,r;flotA44=4,，1,5,8,7,2,1，4,8,3，6,2，6，1,20,b4=-2,-7,-7，-3,44，44，y4;/给定的方程组/inputA()；/npt_b();显示矩阵A、britf（矩阵A44:n）;for(i=0;4;i+）or(k=0；k4;k+)prif（%10f，Ai);rint();for(i0；4;+)u0i=A0i;for(i=0;i4；i+)i0Ai0/u00;l=1；or（r=1;r；r+)/计算矩阵L和Ufor(r;i4；i+

7、)l =；for(k=；kr;+)t=t+lkki;uri=A-t;fr（i=r;；+)fa t10;for(k=0;r;k)t+=li*ur;lr=(Air-t1)/ur；y0=b0；for(i=1；i;i+)fltt0;or(k0;k=0;i-)float t=;for(ki+1;k4;k+)t+=uik*x;i(y-t)ui;printf(*n)；outt_x(x）;/输入矩阵Avoid iut_A()floa A4;it i,j；prinf(input mtriA44:）；r(i=;i4;i+)for（j;j+)scn(%,&Aij)；/输入矩阵bod inut_b（）flo b;in

8、t ;it(inut mar4：n);for（i=0；i4;i+)scnf(%f，&bi）;/输出方程的根xvoid out_x(lot x4)inti；rinf(方程组的根为：n);or(i;+)printf(x%d=%-13f，+1,xi);pritf();QR分解法：#iclde#includedee N 4oiatrix_tme（doubeA,double B,double CN，int n）;/两个矩阵相乘，成果存在矩阵Cvoid atrixvec(oubl AN,oule BN,oubleCN，int n)；/矩阵和向量相乘,成果存在向量Cduble vec_va(dole A,i

9、n n）;/求向量的模vo vetie(oublea,doubleH,int n)；/两个向量相乘得一种矩阵;void oueolder(double *，duble N,int n， int );/求解Huehlder矩阵函数oidmatri_turn(double AN,in n);/求矩阵装置void soluio(ouble ，dule ,i n);求解上三角方程的解void QR(dole AN，ube b,nt n);voiin（)oubl A44=4,2,5，8,7，2，10，4,3,2,，11,20;obeb4=-2，-，-7,3;int i,j；int n=4;ritf(待求

10、解的方程组系数矩阵为:n);r(i;in；+）/显示矩阵Afo(j=0；jn;)pif(%-10f,ij)；printf（n）；QR(A,b，)；void mri_tim(doube A,dolBN，dole CN,in n)/两个矩阵相乘，成果存在矩阵CNin ,；or(i0;n;+)for(j=0；n;j+)Cij0；for(=0；k;+)Cijij+Ak*Bkj;vid matrx_ec(obe A,double B，oube N,int n)/矩阵和向量相乘,成果存在向量CNint ,;or(i=0;in;+)i=0;or（；n;j+）C=C+ij*Bj;duble v_value（d

11、uleA，nt n)/求向量的模oub u=0;it i;for(=0;in;i+)SumSum+Ai*Ai；Sum=qt(u)；tn ;vodv_time(oula,oubl HN,int )两个向量相乘得一种矩阵;i,j;fo(i=0;in;i+）for(j=0；jn；j+)Hij=aiaj;vid househlder(uble *a,dbl ,intn,inm)/计算Householer矩阵int i,j;double fist;/寄存首元素doulee_v;/寄存向量的模doule a1N=0；or（i=0;in；i+)for(j=;jn;+)if（i=j）H=1;elseHij=;

12、fist=am；/取矩阵A转置的第m行(取矩阵A的第m列数)ec_v=ve_value(am,n-)；计算m列元素所构成向量的模m=amec_v;/的分子部分c_v=vcvale(m，n-)；/的分母部分ec_=ve_v*c_;for（i=m；n;+)for(;jn；+)Hij=ai*aj*2vcv;vod rixun(doubeA，int n)/求矩阵的转置dube aNN=0;int i,j;fo(i=0;in；i+）fo(j=0;n;j+)aij=Aij;for（i=；in;i+)or(j0;j=0;i-)fo(jn-1;ji；j-)sum+ij*xj;=（bi-sm)/Aii；su0;

13、rntf(矩阵的R分解求解成果为:);fo（；in;+)printf(X%d=%10f,i1,xi）;od QR(duleAN,double b,t )int i，j,k,t;doube HNN=0,2N0,H3NN=0;/H1,H2寄存相邻两次的Hoeler矩阵,H3为中间变量矩阵doul tempN=0；doube QNN=，RNN=,Q1N=0;ouleb1N=；fo（i0;n；i+)/将矩阵A临时寄存在emp中or（=0；jn；j+)temp=Aij;for(i=0；in;i+)H2ii=；/单位阵matrx_urn(tep,n）；/矩阵A的转置 (以便后续取A的某一列元素)fr(i;

14、n-1；i+）/Houseoder的一次变换ouseholder(temi，,n，i）;/得到Huseolr矩阵1arix_tim(H,Q,n);/矩阵H和A相乘放在Q中for(k=；kn;k+)/更新矩阵A,使得第一列元素除第一种元素外，其她全为0for(t=0;tn;t+）Akkt;o(k=0;kn;+)/寄存临时矩阵for(t=;t+)emtAkt;mrixtrn（emp，n);/对矩阵转置attime(H1,2,H3，n）;for(k0;kn；k+)fo(t=;tn;t+)H2kt=H3kt;/2是一次变换与A相乘后的矩阵(H*A）fr(0;kn;k+)/R矩阵for（=0;tn;+)

15、Rkt=Akt；printf(*n);rintf(系数矩阵A进行QR分解后的矩阵为：n);fo(=0；in;+)for(j0;n;j+)printf(%-0,Rj);printf();for(k;+）for(t=0;tn;t+)Qkt=kt;/此时Q为最后正交矩阵Q的转置matrix_vec（Q，b,b1，n)；/Q矩阵为正交阵，故的逆就等于Q装置,求解的逆matixurn(,n）;/转置,求解pritf（*n）;prntf(系数矩阵A进行QR分解后的Q矩阵为：);for(i;in；i)or(=0;jn;j+）pintf(%-10,Qi);pritf（n);slut(R,b1,n);/求解上三角方程