線性代數(shù)方程組求解

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

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

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

4、QR三、實驗原理解Ax=b 的問題就等價于要求解兩個三角形方程組： Ly=b,求y; Ux=y,求x.設(shè)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的計算公式； 四、實驗步驟1將矩陣A保存進計算機中，再定義2個空矩陣L，U以便保存求出的三角矩陣的值。利用公式，將矩陣A分解為LU，L為單位下三角陣,U為上三角陣。

5、2可知計算方法有三層循環(huán)。先通過公式計算出U矩陣的第一行元素和L矩陣的第一列元素。再根據(jù)公式和，和上次的出的值，求出矩陣其余的元素，每次都要三次循環(huán)，求下一個元素需要上一個結(jié)果。3先由公式 ，Ly=b 求出y，因為L為下三角矩陣，所以由第一行開始求y.4再由公式，Ux=y求出x, 因為U為上三角矩陣，所以由最后一行開始求x.五、程序流程圖1、LU分解法2、QR分解法六、實驗結(jié)果1、 LU分解法方程組1 ：方程組2：方程組3：2、 QR分解法方程組1：方程組2：方程組3：七、實驗總結(jié)為了求解線性方程組，我們通常需要一定的解法。其中一種解法就是通過矩陣的三角分解來實現(xiàn)的，屬于求解線性方程組的直接法

6、。在不考慮舍入誤差下，直接法可以用有限的運算得到精確解，因此主要適用于求解中小型稠密的線性方程組。1、三角分解法 三角分解法是將A矩陣分解成一個上三角形矩陣U和一個下三角形矩陣L，這樣的分解法又稱為LU分解法。它的用途主要在簡化一個大矩陣的行列式值的計算過程，求反矩陣和求解聯(lián)立方程組。不過要注意這種分解法所得到的上下三角形矩陣并非唯一，還可找到數(shù)個不同 的一對上下三角形矩陣，此兩三角形矩陣相乘也會得到原矩陣。 2、 QR分解法QR分解法是將矩陣分解成一個正規(guī)正交矩陣Q與上三角形矩陣R,所以稱為QR分解法。 在編寫這兩個程序過程中，起初遇到不少麻煩！雖然課上老師反復(fù)重復(fù)著：“算法不難的，Its

7、so easy!”但是當(dāng)自己實際操作時，感覺并不是那么容易。畢竟是要把實際的數(shù)學(xué)問題轉(zhuǎn)化為計算機能夠識別的編程算法，所以在編寫程序之前我們仔細(xì)認(rèn)真的把所求解的問題逐一進行詳細(xì)的分析，最終轉(zhuǎn)化為程序段。每當(dāng)遇到問題時，大家或許有些郁悶，但最終還是靜下心來反復(fù)仔細(xì)的琢磨，一一排除了錯誤，最終完成了本次實驗?；仡^一想原來編個程序其實也沒有想象的那么復(fù)雜，只要思路清晰，逐步分析，就可以慢慢搞定了。附 源代碼:#include #include #include #include #include using namespace std;bool lu(double* a, int* pivot, in

8、t 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(double 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

9、= 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, 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+

10、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.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

11、*)malloc(sizeof(double)*n); for (int i=0; in; i+) Pi=Di=0; coutflag; if(!flag) cout矩陣LU分解: n; lu(A,P,n);/矩陣LU分解 for (int i=0; in; i+) for(int j=0; jn; j+) printf(%ft,Ai); coutendl; cout矩陣LU分解的主元位置: n; for (int i=0; in; i+) coutPit; coutnGuass求線性代數(shù)方程組求解:n; guass(A,P,B,n);/求線性代數(shù)方程組的解 for (int i=0; in;

12、 i+) if(i=3) coutendl; printf(%.16ft,Bi); coutnGuass解法的誤差為:n; for (int i=0; in; i+) if(i=3) coutendl; printf(%.16ft,Bi-1); expct=expct + Bi-1; printf(n誤差期望值為%.16ftn,expct/6); coutendl; else cout矩陣QR分解: n; qr(A,D,n);/矩陣qr分解 for (int i=0; in; i+) for(int j=0; jn; j+) printf(%dt,Ai); coutendl; coutQR分解

13、的上三角矩陣對角線元素: n; for (int i=0; in; i+) printf(%ft,Di); coutnHouseholder線性代數(shù)方程組求解:n; householder(A,D,B,n);/求線性代數(shù)方程組的解 for (int i=0; in; i+) if(i=3) coutendl; printf(%.16ft,Bi); coutnHouseholder解法的誤差為:n; for (int i=0; in; i+) if(i=3) coutendl; printf(%.16ft,Bi-1); expct=expct + Bi-1; printf(n誤差期望值為%.16

14、ftn,expct/6); coutendl; 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; in-1; i+)/依次對第i列進行處理 / 選出i列的主元,記錄主元位置 max = fabs(an*i + i); pivoti=i; for(j=i+1; jmax) max= fabs(an*j + i) ; pivoti=j; / 對第i列進行行變換，使得主元在對角線上 if(pivoti!=i) for(j=i;

15、 jn; j+)/ij與pivotij換 只用對上三角進行處理 temp=an*i + j; an*i + j=an*pivoti+ j; an*pivoti+ j=temp; for(j=i+1; jn; j+)/Pi 部分下三角L an*j + i=an*j+i/an*i+i; for(j=i+1; jn; j+)/計算上三角U for(k=i+1; kn; k+) an*j + k=an*j+k- an*j+i*an*i+k; /計算下三角 L for(i=0; in-2; i+)/i行k列 for(k=i+1; kn-1;k+) temp=an*pivotk + i; an*pivot

16、k + 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對b行變換，與LU匹配 for(i=0; in-1; i+) /貌似錯誤在這里哦 temp = bpi; bpi = bi; bi=temp; /Ly=b,將y的內(nèi)容放入b for(i=0; in; i+) for(j=0; j=0; i-) for(j=n-1; ji; j-) bi=bi-lun*

17、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 *)malloc(sizeof(double)*n); for (i=0; in-1; i+)/依次對第i列進行處理 /獲得tem值 m = 0 ; for(j=i; j0) m = -sqrt(m); else m = sqrt(m); /獲得temp放入矩陣，并存主元d tem = 0 ; di = m ; an*i

18、 +i = an*i +i - m; for(j=i; j=n-1; j+) tem=tem + an*j +i*an*j +i; tem= sqrt(tem); for(j=i; j=n-1; j+) an*j +i = an*j +i/tem ; / 調(diào)整矩陣 for(k=i+1;kn;k+) for(j=i; jn; j+) tem = 0 ; for(l=i; ln; l+) tem =tem + an*j + i*an*l + i*an*l + k; tempj = aj*n+k - 2*tem; for(j=i; jn; j+) aj*n+k = tempj; dn-1 = a(n-1)*n+n-1; bool householder(double const*qr, double const*d, double*b, int n)/求