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

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

2、 LU分解法解Ax=b 的问题就等价于要求解两个三角形方程组： Ly=b,求y;    Ux=y,求x.设A为非奇异矩阵,且有分解式A=LU, L为单位下三角阵,U为上三角阵。L,U的元素可以有n步直接计算定出。用直接三角分解法解Ax=b要求A的所有顺序主子式都不为零的计算公式： , ,i=2,3,，n.计算U的第r行，L的第r列元素(i=2,3,n)：     ， i=r,r+1,n; , i=r+1,n,且rn.求解Ly=b，Ux=y的计算公式； 2、 QR分解法4、 实验步骤1、LU分解法1>将矩阵A保存进计算机中，再定义

3、2个空矩阵L，U以便保存求出的三角矩阵的值。利用公式，将矩阵A分解为LU，L为单位下三角阵,U为上三角阵。2>可知计算方法有三层循环。先通过公式计算出U矩阵的第一行元素 和L矩阵的第一列元素。再根据公式和，和上次的出的值，求出矩阵其余的元素，每次都要三次循环，求下一个元素需要上一个结果。3>先由公式 ，Ly=b 求出y，因为L为下三角矩阵，所以由第一行开始求y.4>再由公式，Ux=y求出x, 因为U为上三角矩阵，所以由最后一行开始求x.2、QR分解法5、 程序流程图1、LU分解法2、QR分解法6、 实验结果1、 LU分解法2、QR分解法7、 实验总结为了求解线性方程

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

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

6、明白的大家一起讨论，越讨论越清晰，愈接近最优答案。这样不管干什么都能事半功倍。庆幸自己有这么个团队，也明白大家一起劳动的果实最珍贵。附 源代码:LR分解法：#include<stdio.h>void input_A();/输入矩阵Avoid input_b();/输入矩阵bvoid output_x(float x4);/输出方程组的根void main()int i,k,r;float A44=4,2,1,5,8,7,2,10,4,8,3,6,12,6,11,20,b4=-2,-7,-7,-3,x4,u44,l44,y4;/给定的方程组/input_A();/input_b();

7、/显示矩阵A、bprintf("矩阵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;lii=1;for(r=1;r<4;r+)/计算矩阵L和Ufor(i=r;i<4;i+)float t=0;for(k=0;k<r;k+)t=t+lrk*uki;uri=Ari-t;for(i=r;i<4;i+)fl

8、oat t1=0;for(k=0;k<r;k+)t1+=lik*ukr;lir=(Air-t1)/urr;y0=b0;for(i=1;i<4;i+)float t=0;for(k=0;k<i;k+)t+=lik*yk;yi=bi-t;for(i=3;i>=0;i-)float t=0;for(k=i+1;k<4;k+)t+=uik*xk;xi=(yi-t)/uii;printf("*n");output_x(x);/输入矩阵Avoid input_A()float A44;int i,j;printf("input matrix A4

9、4:n");for(i=0;i<4;i+)for(j=0;j<4;j+)scanf("%f",&Aij);/输入矩阵bvoid input_b()float b4;int i;printf("input matrix b4:n");for(i=0;i<4;i+)scanf("%f",&bi);/输出方程的根xvoid output_x(float x4)int i;printf("方程组的根为:n");for(i=0;i<4;i+)printf("x%d=

10、%-13f",i+1,xi);printf("n");QR分解法：#include<stdio.h>#include<math.h>#define N 4void matrix_time(double AN,double BN,double CN,int n);/两个矩阵相乘,结果存在矩阵CNvoid matrix_vec(double AN,double BN,double CN,int n);/矩阵和向量相乘,结果存在向量CNdouble vec_value(double A,int n);/求向量的模void vec_time(dou

11、ble a,double HN,int n);/两个向量相乘得一个矩阵；void householder(double *a,double HN,int n, int m);/求解Householder矩阵函数void matrix_turn(double AN,int n);/求矩阵装置void solution(double AN,double b,int n);/求解上三角方程的解void QR(double AN,double b,int n);void main()double A44=4,2,1,5,8,7,2,10,4,8,3,6,12,6,11,20;double b4=-2,

12、-7,-7,-3;int i,j;int n=4;printf("待求解的方程组系数矩阵A为：n");for(i=0;i<n;i+)/显示矩阵Afor(j=0;j<n;j+)printf("%-10f",Aij);printf("n");QR(A,b,n);void matrix_time(double AN,double BN,double CN,int n)/两个矩阵相乘，结果存在矩阵CNint i,j,k;for(i=0;i<n;i+)for(j=0;j<n;j+)Cij=0;for(k=0;k<n

13、;k+)Cij=Cij+Aik*Bkj;void matrix_vec(double AN,double BN,double CN,int n)/矩阵和向量相乘,结果存在向量CNint i,j;for(i=0;i<n;i+)Ci=0;for(j=0;j<n;j+)Ci=Ci+Aij*Bj;double vec_value(double A,int n)/求向量的模double Sum=0;int i;for(i=0;i<n;i+)Sum=Sum+Ai*Ai;Sum=sqrt(Sum);return Sum;void vec_time(double a,double HN,in

14、t n)/两个向量相乘得一个矩阵；int i,j;for(i=0;i<n;i+)for(j=0;j<n;j+)Hij=ai*aj;void householder(double *a,double HN,int n, int m)/计算Householder矩阵int i,j;double first;/存放首元素double vec_v;/存放向量的模double a1N=0;for(i=0;i<n;i+)for(j=0;j<n;j+)if(i=j)Hij=1;elseHij=0;first=am;/取矩阵A转置的第m行取矩阵A的第m列数vec_v=vec_value

15、(&am,n-m);/计算m列元素所构成向量的模am=am-vec_v;/w的分子局部vec_v=vec_value(&am,n-m);/w的分母局部vec_v=vec_v*vec_v;for(i=m;i<n;i+)for(j=m;j<n;j+)Hij+=-ai*aj*2/vec_v;void matrix_turn(double AN,int n)/求矩阵的转置double aNN=0;int i,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+)Ai

16、j=aji;void solution(double AN,double b,int n)/求解上三角方程的解int i,j;double xN=0;double sum=0;for(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);void QR(double AN,double b,int n)int i,j,k,t;

17、double H1NN=0,H2NN=0,H3NN=0;/H1,H2存放相邻两次的Householder矩阵，H3为中间变量矩阵double tempNN=0;double QNN=0,RNN=0,Q1NN=0;double b1N=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(t

18、empi,H1,n,i);/得到Householder矩阵H1matrix_time(H1,A,Q,n);/矩阵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+)tempkt=Akt;matrix_turn(temp,n);/对矩阵转置matrix_time(H1,H2,H3,n);for(k=0;k<n;k+)for(t=0;t<n;t+)H2kt=H3kt;/H2是一次变换与A相乘后的矩阵H1*Afor(k=0;k<n;k+)/R矩阵for(t=0;t<n;t+)Rkt=Akt;printf("*n");printf("系数矩阵A进行QR分解后的R矩阵为：n");for(i=0;i<n;i+)for(j=0;j<n;j+)printf("%-10f",Rij);printf("n");for(k=0;k&l