LU和QR分解法解線性方程組

1、LU和QR法解線性方程組 一、 問題描述求解方程組=,要求：1、 編寫用三角LU分解法求解線性方程組；2、 編寫用正交三角QR分解法求解線性方程組。2、 問題分析求解線性方程組Ax=b，其實(shí)質(zhì)就是把它的系數(shù)矩陣A通過各種變換成一個(gè)下三角或上三角矩陣，從而簡化方程組的求解。因此，在求解線性方程組的過程中，把系數(shù)矩陣A變換成上三角或下三角矩陣顯得尤為重要，然而矩陣A的變換通常有兩種分解方法：LU分解法和QR分解法。1、 LU分解法：將A分解為一個(gè)下三角矩陣L和一個(gè)上三角矩陣U，即：A=LU,其中 L=， U=2、 QR分解法：將A分解為一個(gè)正交矩陣Q和一個(gè)上三角矩陣R,即：A=QR三、實(shí)驗(yàn)原理1、

2、 LU分解法解Ax=b 的問題就等價(jià)于要求解兩個(gè)三角形方程組： Ly=b,求y;    Ux=y,求x.設(shè)A為非奇異矩陣,且有分解式A=LU, L為單位下三角陣,U為上三角陣。L,U的元素可以有n步直接計(jì)算定出。用直接三角分解法解Ax=b要求A的所有順序主子式都不為零的計(jì)算公式： , ,i=2,3,，n.計(jì)算U的第r行，L的第r列元素(i=2,3,n)：     ， i=r,r+1,n; , i=r+1,n,且rn.求解Ly=b，Ux=y的計(jì)算公式； 2、 QR分解法4、 實(shí)驗(yàn)步驟1、LU分解法1>將矩陣A保存進(jìn)計(jì)算機(jī)中，再定義

3、2個(gè)空矩陣L，U以便保存求出的三角矩陣的值。利用公式，將矩陣A分解為LU，L為單位下三角陣,U為上三角陣。2>可知計(jì)算方法有三層循環(huán)。先通過公式計(jì)算出U矩陣的第一行元素 和L矩陣的第一列元素。再根據(jù)公式和，和上次的出的值，求出矩陣其余的元素，每次都要三次循環(huán)，求下一個(gè)元素需要上一個(gè)結(jié)果。3>先由公式 ，Ly=b 求出y，因?yàn)長為下三角矩陣，所以由第一行開始求y.4>再由公式，Ux=y求出x, 因?yàn)閁為上三角矩陣，所以由最后一行開始求x.2、QR分解法5、 程序流程圖1、LU分解法2、QR分解法6、 實(shí)驗(yàn)結(jié)果1、 LU分解法2、QR分解法7、 實(shí)驗(yàn)總結(jié)為了求解線性方程

4、組，我們通常需要一定的解法。其中一種解法就是通過矩陣的三角分解來實(shí)現(xiàn)的，屬于求解線性方程組的直接法。在不考慮舍入誤差下，直接法可以用有限的運(yùn)算得到精確解，因此主要適用于求解中小型稠密的線性方程組。1、三角分解法 三角分解法是將A矩陣分解成一個(gè)上三角形矩陣U和一個(gè)下三角形矩陣L，這樣的分解法又稱為LU分解法。它的用途主要在簡化一個(gè)大矩陣的行列式值的計(jì)算過程，求反矩陣和求解聯(lián)立方程組。不過要注意這種分解法所得到的上下三角形矩陣并非唯一，還可找到數(shù)個(gè)不同 的一對上下三角形矩陣，此兩三角形矩陣相乘也會得到原矩陣。 2、 QR分解法 QR分解法是將矩陣分解成一個(gè)正規(guī)正交矩陣Q與上三角形矩陣R,所以稱為Q

5、R分解法。 在編寫這兩個(gè)程序過程中，起初遇到不少麻煩！雖然課上老師反復(fù)重復(fù)著：“算法不難的，It's so easy!但是當(dāng)自己實(shí)際操作時(shí)，感覺并不是那么容易。畢竟是要把實(shí)際的數(shù)學(xué)問題轉(zhuǎn)化為計(jì)算機(jī)能夠識別的編程算法，所以在編寫程序之前我們仔細(xì)認(rèn)真的把所求解的問題逐一進(jìn)行詳細(xì)的分析，最終轉(zhuǎn)化為程序段。每當(dāng)遇到問題時(shí)，大家或許有些郁悶，但最終還是靜下心來反復(fù)仔細(xì)的琢磨，一一排除了錯(cuò)誤，最終完成了本次實(shí)驗(yàn)?；仡^一想原來編個(gè)程序其實(shí)也沒有想象的那么復(fù)雜，只要思路清晰，逐步分析，就可以慢慢搞定了。通過這次實(shí)驗(yàn)，讓我們認(rèn)知到團(tuán)隊(duì)的作用力度是不容無視的，以后不管干任何時(shí)都要注重團(tuán)隊(duì)合作，遇到不懂得不

6、明白的大家一起討論，越討論越清晰，愈接近最優(yōu)答案。這樣不管干什么都能事半功倍。慶幸自己有這么個(gè)團(tuán)隊(duì)，也明白大家一起勞動(dòng)的果實(shí)最珍貴。附 源代碼: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+)/計(jì)算矩陣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);/兩個(gè)矩陣相乘,結(jié)果存在矩陣CNvoid matrix_vec(double AN,double BN,double CN,int n);/矩陣和向量相乘,結(jié)果存在向量CNdouble vec_value(double A,int n);/求向量的模void vec_time(dou

11、ble a,double HN,int n);/兩個(gè)向量相乘得一個(gè)矩陣；void householder(double *a,double HN,int n, int m);/求解Householder矩陣函數(shù)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("待求解的方程組系數(shù)矩陣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)/兩個(gè)矩陣相乘，結(jié)果存在矩陣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)/矩陣和向量相乘,結(jié)果存在向量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)/兩個(gè)向量相乘得一個(gè)矩陣；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)/計(jì)算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轉(zhuǎn)置的第m行取矩陣A的第m列數(shù)vec_v=vec_value

15、(&am,n-m);/計(jì)算m列元素所構(gòu)成向量的模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)/求矩陣的轉(zhuǎ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分解求解結(jié)果為：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臨時(shí)存放在temp中for(j=0;j<n;j+)tempij=Aij;for(i=0;i<n;i+)H2ii=1;/單位陣matrix_turn(temp,n);/矩陣A的轉(zhuǎn)置 方便后續(xù)取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，使得第一列元素除第一個(gè)元素外，其他全為0for(t=0;t<n;t+)Akt=Qkt;for(k=0;k<n;k+)/存放臨時(shí)矩陣for(t=0;t<n;t+)tempkt=Akt;matrix_turn(temp,n);/對矩陣轉(zhuǎ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("系數(shù)矩陣A進(jìn)行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