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

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、 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保存進計算機中，再定義2個空矩陣L，U以便保存求出的三角矩陣的值。利用公式，將矩陣A分解為LU，L為單位下三角陣,U為上三角陣。2可知計算方法有三層循環(huán)。先通過公式計算出U矩陣的第一行元素和L矩陣的第一列元素。再根據公式和，和上次的出的值，求出矩陣其余的元素，每次都要三次循環(huán)，求下一個元素需要上一個結果。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、 實驗總結為了求解線性方程組，我們通常需要一定的解法。其中一種解法就是通過矩陣的三角分解來實現的，屬于求解線性方程組的直接法。在不考慮舍入誤差下，直接法可以用有限的運算得到精確解，因此主要適用于求解中小型稠密的線性方程組。1、三角分解法 三角分解法是將A矩陣分解成一個上三角形矩陣U和一個下三角形矩陣L，這樣的分解法又稱為LU分解法。它的用途主要在簡化一個大矩陣的行列式值的計算過程，求反矩陣和求解聯(lián)立方程組。不過要注意這種分解法所得到的上下三角形矩陣并非唯一，還可找到數個不同 的一對上下三角形矩陣，此兩三角形矩陣相乘也會得到原矩陣。 2、 QR分解法 QR分解法是將矩陣分解成一個正規(guī)正交矩陣Q與上三角形矩陣R,所以稱為QR分解法。 在編寫這兩個程序過程中，起初遇到不少麻煩！雖然課上老師反復重復著：“算法不難的，Its so easy!”但是當自己實際操作時，感覺并不是那么容易。畢竟是要把實際的數學問題轉化為計算機能夠識別的編程算法，所以在編寫程序之前我們仔細認真的把所求解的問題逐一進行詳細的分析，最終轉化為程序段。每當遇到問題時，大家或許有些郁悶，但最終還是靜下心來反復仔細的琢磨，一一排除了錯誤，最終完成了本次實驗?；仡^一想原來編個程序其實也沒有想象的那么復雜，只要思路清晰，逐步分析，就可以慢慢搞定了。通過這次實驗，讓我們認知到團隊的作用力度是不容忽視的，以后不管干任何時都要注重團隊合作，遇到不懂得不明白的大家一起討論，越討論越清晰，愈接近最優(yōu)答案。這樣不管干什么都能事半功倍。慶幸自己有這么個團隊，也明白大家一起勞動的果實最珍貴。附 源代碼:LR分解法：#includevoid 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();/顯示矩陣A、bprintf(矩陣A44:n);for(i=0;i4;i+)for(k=0;k4;k+)printf(%-10f,Aik);printf(n);for(i=0;i4;i+)u0i=A0i;for(i=0;i4;i+)li0=Ai0/u00;lii=1;for(r=1;r4;r+)/計算矩陣L和Ufor(i=r;i4;i+)float t=0;for(k=0;kr;k+)t=t+lrk*uki;uri=Ari-t;for(i=r;i4;i+)float t1=0;for(k=0;kr;k+)t1+=lik*ukr;lir=(Air-t1)/urr;y0=b0;for(i=1;i4;i+)float t=0;for(k=0;k=0;i-)float t=0;for(k=i+1;k4;k+)t+=uik*xk;xi=(yi-t)/uii;printf(*n);output_x(x);/輸入矩陣Avoid input_A()float A44;int i,j;printf(input matrix A44:n);for(i=0;i4;i+)for(j=0;j4;j+)scanf(%f,&Aij);/輸入矩陣bvoid input_b()float b4;int i;printf(input matrix b4:n);for(i=0;i4;i+)scanf(%f,&bi);/輸出方程的根xvoid output_x(float x4)int i;printf(方程組的根為:n);for(i=0;i4;i+)printf(x%d=%-13f,i+1,xi);printf(n);QR分解法：#include#include#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(double 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,-7,-7,-3;int i,j;int n=4;printf(待求解的方程組系數矩陣A為：n);for(i=0;in;i+)/顯示矩陣Afor(j=0;jn;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;in;i+)for(j=0;jn;j+)Cij=0;for(k=0;kn;k+)Cij=Cij+Aik*Bkj;void matrix_vec(double AN,double BN,double CN,int n)/矩陣和向量相乘,結果存在向量CNint i,j;for(i=0;in;i+)Ci=0;for(j=0;jn;j+)Ci=Ci+Aij*Bj;double vec_value(double A,int n)/求向量的模double Sum=0;int i;for(i=0;in;i+)Sum=Sum+Ai*Ai;Sum=sqrt(Sum);return Sum;void vec_time(double a,double HN,int n)/兩個向量相乘得一個矩陣；int i,j;for(i=0;in;i+)for(j=0;jn;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;in;i+)for(j=0;jn;j+)if(i=j)Hij=1;elseHij=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的分母部分vec_v=vec_v*vec_v;for(i=m;in;i+)for(j=m;jn;j+)Hij+=-ai*aj*2/vec_v;void matrix_turn(double AN,int n)/求矩陣的轉置double aNN=0;int i,j;for(i=0;in;i+)for(j=0;jn;j+)aij=Aij;for(i=0;in;i+)for(j=0;j=0;i-) for(j=n-1;ji;j-)sum+=Aij*xj;xi=(bi-sum)/Aii;sum=0;printf(矩陣的QR分解求解結果為：n);for(i=0;in;i+)printf(X%d=%-10fn,i+1,xi);void QR(double AN,double b,int n)int i,j,k,t;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;in;i+)/將矩陣A臨時存放在temp中for(j=0;jn;j+)tempij=Aij;for(i=0;in;i+)H2ii=1;/單位陣matrix_turn(temp,n);/矩陣A的轉置 （方便后續(xù)取A的某一列元素）for(i=0;in-1;i+)/Householder的一次變換householder(tempi,H1,n,i);/得到Householder矩陣H1matrix_time(H1,A,Q,n);/矩陣H1和A相乘放在Q中for(k=0;kn;k+)/更新矩陣A，使得第一列元素除第一個元素外，其他全為0for(t=0;tn;t+)Akt=Qkt;for(k=0;kn;k+)/存放臨時矩陣for(t=0;tn;t+)tempkt=Akt;matrix_turn(temp,n);/對矩陣轉置matrix_time(H1,H2,H3,n);for(k=0;kn;k+)for(t=0;tn;t+)H2kt=H3kt;/H2是一次變換與A相乘后的矩陣（H1*A）for(k=0;kn;k+)/R矩陣for(t=0;tn;t+)Rkt=Akt;printf(*n);printf(系數矩陣A進行QR分解后的R矩陣為：n);for(i=0;in;i+)for(j=0;jn;j+)printf(%-10f,