【正文】 實(shí)驗(yàn)報(bào)告 4 通過(guò)重復(fù)求解 N 各線性方程組 ACJ=EJ，其中 J=1， 2， … ， N 來(lái)得到 A1, 則 A[C1 C2 … CN]=[E1 E2 … EN] 而且 A1=[C1 C2 … CN] 保證對(duì) LU 分解只計(jì)算一次 ! 3 設(shè)有如下三角線性方程組，而且系數(shù)矩陣具 有嚴(yán)格對(duì)角優(yōu)勢(shì)： d1x1+c1x2 =b1 a1x1+d2x2+c2x3 =b2 a2x2+d3x3+c3x4 =b3 2 求解線性方程組 AX=B，其中 A=[aij]NN， ,aij=ij1；而且 B=[bij]N1， b11=N，當(dāng) i≥時(shí)， bij=(iN1)/(i1)。 三、 實(shí)驗(yàn)內(nèi)容 1 許多科學(xué)應(yīng)用包含的矩陣帶有很多的零。 研究雅可比迭代法，我們發(fā)現(xiàn)在逐個(gè)求 ）（ 1kX? 的分量時(shí)，當(dāng)計(jì)算到 ）（ 1kiX? 時(shí) ,分量 ）（ 1k1X? ， ......， ）（ 1k1iX? 都已經(jīng)求得，而仍用舊分量 k1X ， ......， ）（ k1iX 計(jì)算 ）（ 1kiX? 。 《數(shù)值方法》實(shí)驗(yàn)報(bào)告 2 三角分解法是將原正方 (square)矩陣分解成一個(gè)上三角形矩陣或是排列 (permuted) 的上三角形矩陣和一個(gè) 下三角形矩陣，這樣的分解法又稱為 LU 分解法。消去過(guò)程實(shí)際上是對(duì)增廣矩陣作行初等變換。 二、 實(shí)驗(yàn)原理 數(shù)學(xué)上，高斯消元法是線性代數(shù)規(guī)劃中的一個(gè)算法，可用來(lái)為線性方程組求解。關(guān)于線性方程組的數(shù)值解法一般有兩類：直接法和迭代法。 《數(shù)值方法》實(shí)驗(yàn)報(bào)告 1 線性方程組 AX=B 的數(shù)值計(jì)算方法實(shí)驗(yàn) 【 摘要 】 在自然科學(xué)與工程技術(shù)中很多問(wèn)題的解決常常歸結(jié)為解線性代數(shù)方程組。 關(guān) 鍵 字 高斯消元法、三角分解法、高斯 賽德 爾迭代、稀疏矩陣 一、 實(shí)驗(yàn)?zāi)康? 、三角分解法、高斯 — 賽德?tīng)柕l(fā)的編程技巧。 高斯（ Gauss）夏鷗按法其實(shí)是將一般的線性方程組變換為三角形（上三角）方程組求解問(wèn)題（消元法），只是步驟規(guī)范，便于編寫(xiě)計(jì)算機(jī)程序。 對(duì)一般的 n 階方程組，消去過(guò)程分 n1 步：第一步消去 11a 下方元素。它的用途主要在簡(jiǎn)化一個(gè)大矩陣的行列式值的計(jì)算過(guò)程，求 反矩陣，和求解聯(lián)立方程組。由于新計(jì)算出的分量比舊分量準(zhǔn)確些，因此設(shè)想一旦新分量 ）（ 1k1X? ， ......， ）（ 1k1iX? 求出，馬上就用新分量 ）（ 1k1X? ， ......， ）（ 1k1iX? 代替雅可比迭代法中 k1X ， ......， ）（ k1iX 來(lái)求 ）（ 1kiX? 這就是高斯 賽德?tīng)?(GaussSeidel)迭代法。在實(shí)際情況中很重要的三角形線性方程組有如下形式： d1x1+c1x2 =b1 a1x1+d2x2+c2x3 =b2 a2x2+d3x3+c3x4 =b3 aN2xN2+dN1xN1+cN1xN =bN1 aN1xN1+dNxN =bN 構(gòu)造一個(gè)程序求解三角形線性方程組。 對(duì) N=3， 7， 11 的 情況分別求解。 算法必須有效地利用系數(shù)矩陣的稀疏性。 a1jxj+ 方程組（ 1） aj1x1+aj2x2+ aNjxj+ N N 算法 流程圖： 《數(shù)值方法》實(shí)驗(yàn)報(bào)告 6 start InputA,B,P,delta,max1 N=length(B)。P))。 1 實(shí)驗(yàn)描述： 本次實(shí)驗(yàn)的解法為使用 LU矩陣求解 X，該解法的內(nèi)容為將系數(shù)矩陣 A分解為一上三角矩陣及一下三角矩陣，且有 A=LU，之后由 LY=B， UX=Y分別求解出 Y，X。 N 的求解方法使用 LU 三角分解法求 《數(shù)值方法》實(shí)驗(yàn)報(bào)告 12 解，用此方法求解出各個(gè) Ek 對(duì)應(yīng)的 Ck，最后以此構(gòu)成 A1； (4)求出 LU 后，應(yīng)判斷矩陣對(duì)角線上是否存在為 0 的元素，若存在，則 A 不存在逆矩陣；若不存在，則可求解逆矩陣 A1； (5)上述方法中的 LU 分解只需要進(jìn)行一次； (6)對(duì)于程序的正確性使用矩陣 1 3 22 6 4251??????及矩陣 1 2 32132 3 1??????進(jìn)行驗(yàn)證，其中，矩陣 1 3 22 6 4251??????不存在逆矩陣，矩陣 1 2 32132 3 1??????的逆矩陣為 66 7 0 .58 33 5 33 3 16 7 0 .25 33 3 0 .08 33 5?????? 實(shí)驗(yàn)結(jié)果： （ 1）輸入為矩陣 1 3 22 6 4251?????? 圖 7 （ 2） 輸入 矩陣為 1 2 32132 3 1?????? 《數(shù)值方法》實(shí)驗(yàn)報(bào)告 13 圖 8 實(shí)驗(yàn)結(jié)論： 如果系數(shù)矩陣能分解為 LU的形式，其中 L為下三角矩陣， U為上三角矩陣，通過(guò)對(duì)系數(shù)矩陣的分解再求解可應(yīng)用簡(jiǎn)單的迭代進(jìn)行求解 x。 算法流程圖： start [N,N]=size(A)。 p=1 N=length(B)。 relerr=err/(norm(X)+eps)。 A(N1,N)*P(N))/A(N1,N1) X(j)=(B(j)A(j,j1)*X(j1)39。 float *x=new float[N1]。n++) //進(jìn)行高斯 消元法計(jì)算 x[k]， k=1， 2，3 A[n][n1]=0。 for(n=N2。 } return x。 int i,k。 //生成動(dòng)態(tài)增廣矩陣 for(i=0。 //輸入增廣矩陣的值 for(i=0。k++) cinA[i][k]。iN。 } 1 includeiostream include using namespace std。 //輸入矩陣的階數(shù)，用于生成動(dòng)態(tài)矩陣 cinN。 float **A=new float*[N]。 cout請(qǐng)輸入矩陣 A的值： endl。k=N。i=N。 for(i=0。 return 0。 //生成二維數(shù)組 U float *x=new float [N]。i=N。 for(k=i。 } } for(i=0。j++) B[i]=B[i]y[j]*U[i][j]。i) //計(jì)算解 X的值 { for(j=N。 } return x。 double **A。 double *lufact(double**,double*,int)。 A=buildA(N)。 for(i=0。 return 0。 //生成二維數(shù)組 U，用于儲(chǔ)存 Lamp。 for(i=0。j++) { c=U[j][i]/U[i][i]。 U[j][i]=c。ji。i=0。 x[i]=y[i]/U[i][i]。 //生成二維動(dòng)態(tài)數(shù)組A[N1][N1] for(i=0。i=N。 return A。 《數(shù)值方法》實(shí)驗(yàn)報(bào)告 25 B[0]=N。 return B。 float **luchange(float**,int*,int)。 N=N1。 //生成用于存放矩陣 A的