線性代數(shù)方程組求解參考模板

1、線性代數(shù)方程組求解 一、實(shí)驗(yàn)要求編程求解方程組： 方程組1： 方程組2： 方程組3： 要求：用C/C+語(yǔ)言實(shí)現(xiàn)如下函數(shù)：1. bool lu(double* a, int* pivot, int n);實(shí)現(xiàn)矩陣的LU分解。pivot為輸出參數(shù)，pivot0,n) 中存放主元的位置排列。函數(shù)成功時(shí)返回false，否則返回true。2. bool guass(double const* lu, int const* p, double* b, int n);1 / 17求線代數(shù)方程組的解設(shè)矩陣Lunxn為某個(gè)矩陣anxn的LU分解，在內(nèi)存中按行優(yōu)先次序存放。p0,n)為L(zhǎng)U分解的主元排列。b為方程

2、組Ax=b的右端向量。此函數(shù)計(jì)算方程組Ax=b的解，并將結(jié)果存放在數(shù)組b0,n)中。函數(shù)成功時(shí)返回false，否則返回true。3. void qr(double* a, double* d, int n);矩陣的QR分解假設(shè)數(shù)組anxn在內(nèi)存中按行優(yōu)先次序存放。此函數(shù)使用HouseHolder變換將其就地進(jìn)行QR分解。d為輸出參數(shù)，d 0,n) 中存放QR分解的上三角對(duì)角線元素。4. bool hshld(double const*qr, double const*d, double*b, int n); 求線代數(shù)方程組的解設(shè)矩陣qrnxn為某個(gè)矩陣anxn的QR分解，在內(nèi)存中按行優(yōu)先次序存

3、放。d 0,n) 為QR分解的上三角對(duì)角線元素。b為方程組Ax=b的右端向量。函數(shù)計(jì)算方程組Ax=b的解，并將結(jié)果存放在數(shù)組b0,n)中。函數(shù)成功時(shí)返回false，否則返回true。二、問(wèn)題分析求解線性方程組Ax=b，其實(shí)質(zhì)就是把它的系數(shù)矩陣A通過(guò)各種變換成一個(gè)下三角或上三角矩陣，從而簡(jiǎn)化方程組的求解。因此，在求解線性方程組的過(guò)程中，把系數(shù)矩陣A變換成上三角或下三角矩陣顯得尤為重要，然而矩陣A的變換通常有兩種分解方法：LU分解法和QR分解法。1、 LU分解法：將A分解為一個(gè)下三角矩陣L和一個(gè)上三角矩陣U，即：A=LU,其中 L=， U=2、 QR分解法：將A分解為一個(gè)正交矩陣Q和一個(gè)上三角矩陣

4、R,即：A=QR三、實(shí)驗(yàn)原理解Ax=b 的問(wèn)題就等價(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ì)算公式； 四、實(shí)驗(yàn)步驟1>將矩陣A保存進(jìn)計(jì)算機(jī)中，再定義2個(gè)空矩陣

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

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

7、陣Q與上三角形矩陣R,所以稱為QR分解法。 在編寫這兩個(gè)程序過(guò)程中，起初遇到不少麻煩！雖然課上老師反復(fù)重復(fù)著：“算法不難的，It's so easy!”但是當(dāng)自己實(shí)際操作時(shí)，感覺(jué)并不是那么容易。畢竟是要把實(shí)際的數(shù)學(xué)問(wèn)題轉(zhuǎn)化為計(jì)算機(jī)能夠識(shí)別的編程算法，所以在編寫程序之前我們仔細(xì)認(rèn)真的把所求解的問(wèn)題逐一進(jìn)行詳細(xì)的分析，最終轉(zhuǎn)化為程序段。每當(dāng)遇到問(wèn)題時(shí)，大家或許有些郁悶，但最終還是靜下心來(lái)反復(fù)仔細(xì)的琢磨，一一排除了錯(cuò)誤，最終完成了本次實(shí)驗(yàn)?；仡^一想原來(lái)編個(gè)程序其實(shí)也沒(méi)有想象的那么復(fù)雜，只要思路清晰，逐步分析，就可以慢慢搞定了。附 源代碼:#include <iostream>#i

8、nclude <stdio.h>#include <math.h> #include <stdlib.h>#include <conio.h> using namespace std;bool lu(double* a, int* pivot, int n);/矩陣LU分解bool guass(double const* lu, int const* p, double* b, int n);/求線性代數(shù)方程組的解void qr(double* a, double* d, int n); /矩陣的QR分解bool householder(doub

9、le const*qr, double const*d, double*b, int n); int main() int n=0; int temp=0; bool flag = false; double expct=0;/誤差期望值 double devsq=0;/誤差的方差 int * P= NULL; double * D= NULL; double A= 1, 1/2.0, 1/3.0, 1/4.0, 1/5.0, 1/6.0, 1/2.0, 1/3.0, 1/4.0, 1/5.0, 1/6.0, 1/7.0, 1/3.0, 1/4.0, 1/5.0, 1/6.0, 1/7.0,

10、1/8.0, 1/4.0, 1/5.0, 1/6.0, 1/7.0, 1/8.0, 1/9.0, 1/5.0, 1/6.0, 1/7.0, 1/8.0, 1/9.0, 1/10.0, 1/6.0, 1/7.0, 1/8.0, 1/9.0, 1/10.0,1/11.0 ; double B= 1+ 1/2.0+ 1/3.0+ 1/4.0+ 1/5.0+ 1/6.0, 1/2.0+ 1/3.0+ 1/4.0+ 1/5.0+ 1/6.0+ 1/7.0, 1/3.0+ 1/4.0+ 1/5.0+ 1/6.0+ 1/7.0+ 1/8.0, 1/4.0+ 1/5.0+ 1/6.0+ 1/7.0+ 1/8.

11、0+ 1/9.0, 1/5.0+ 1/6.0+ 1/7.0+ 1/8.0+ 1/9.0+ 1/10.0, 1/6.0+ 1/7.0+ 1/8.0+ 1/9.0+ 1/10.0+1/11.0 ; n = 6 ; P = (int*)malloc(sizeof(int)*n); D = (double*)malloc(sizeof(double)*n); for (int i=0; i<n; i+) Pi=Di=0; cout<<"Please choose the way to solve (0: lu,1:qr):" cin>>flag; if

12、(!flag) cout<<"矩陣LU分解: n" lu(A,P,n);/矩陣LU分解 for (int i=0; i<n; i+) for(int j=0; j<n; j+) printf("%ft",Ai); cout<<endl; cout<<"矩陣LU分解的主元位置: n" for (int i=0; i<n; i+) cout<<Pi<<"t" cout<<"nGuass求線性代數(shù)方程組求解:n"

13、 guass(A,P,B,n);/求線性代數(shù)方程組的解 for (int i=0; i<n; i+) if(i=3) cout<<endl; printf("%.16ft",Bi); cout<<"nGuass解法的誤差為:n" for (int i=0; i<n; i+) if(i=3) cout<<endl; printf("%.16ft",Bi-1); expct=expct + Bi-1; printf("n誤差期望值為%.16ftn",expct/6); c

14、out<<endl; else cout<<"矩陣QR分解: n" qr(A,D,n);/矩陣qr分解 for (int i=0; i<n; i+) for(int j=0; j<n; j+) printf("%dt",Ai); cout<<endl; cout<<"QR分解的上三角矩陣對(duì)角線元素: n" for (int i=0; i<n; i+) printf("%ft",Di); cout<<"nHouseholder線性

15、代數(shù)方程組求解:n" householder(A,D,B,n);/求線性代數(shù)方程組的解 for (int i=0; i<n; i+) if(i=3) cout<<endl; printf("%.16ft",Bi); cout<<"nHouseholder解法的誤差為:n" for (int i=0; i<n; i+) if(i=3) cout<<endl; printf("%.16ft",Bi-1); expct=expct + Bi-1; printf("n誤差期望

16、值為%.16ftn",expct/6); cout<<endl; return 0 ; bool lu(double* a, int* pivot, int n)/矩陣LU分解 int i,j,k; double max,temp; max = 0; temp = 0; for (i=0; i<n-1; i+)/依次對(duì)第i列進(jìn)行處理 / 選出i列的主元,記錄主元位置 max = fabs(an*i + i); pivoti=i; for(j=i+1; j<n; j+)/j表示第j行 if( fabs(an*j + i)>max) max= fabs(an

17、*j + i) ; pivoti=j; / 對(duì)第i列進(jìn)行行變換，使得主元在對(duì)角線上 if(pivoti!=i) for(j=i; j<n; j+)/ij與pivotij換 只用對(duì)上三角進(jìn)行處理 temp=an*i + j; an*i + j=an*pivoti+ j; an*pivoti+ j=temp; for(j=i+1; j<n; j+)/Pi 部分下三角L an*j + i=an*j+i/an*i+i; for(j=i+1; j<n; j+)/計(jì)算上三角U for(k=i+1; k<n; k+) an*j + k=an*j+k- an*j+i*an*i+k; /

18、計(jì)算下三角 L for(i=0; i<n-2; i+)/i行k列 for(k=i+1; k<n-1;k+) temp=an*pivotk + i; an*pivotk + i=ak*n + i; ak*n + i=temp; return false ; bool guass(double const* lu, int const* p, double* b, int n)/求線性代數(shù)方程組的解 int i,j,k; double temp; /按qivot對(duì)b行變換，與LU匹配 for(i=0; i<n-1; i+) /貌似錯(cuò)誤在這里哦 temp = bpi; bpi =

19、bi; bi=temp; /Ly=b,將y的內(nèi)容放入b for(i=0; i<n; i+) for(j=0; j<i; j+) bi=bi-lun*i+j*bj; /Uy=x,將x的內(nèi)容放入b for(i=n-1; i>=0; i-) for(j=n-1; j>i; j-) bi=bi-lun*i+j*bj; bi=bi/lun*i+i; return false; void qr(double* a, double* d, int n) /矩陣的QR分解 int i,j,l,k; double tem,m; double *temp; temp = (double *

20、)malloc(sizeof(double)*n); for (i=0; i<n-1; i+)/依次對(duì)第i列進(jìn)行處理 /獲得tem值 m = 0 ; for(j=i; j<n; j+)/j表示第j行 m = m +an*j + i*an*j + i ; if(an*i + i >0) m = -sqrt(m); else m = sqrt(m); /獲得temp放入矩陣，并存主元d tem = 0 ; di = m ; an*i +i = an*i +i - m; for(j=i; j<=n-1; j+) tem=tem + an*j +i*an*j +i; tem= s