消去法等解線性方程組(一)

1、消去法等解線性方程組一）中文摘要Gauss消去法是目前求解中小規(guī)模線性方程組即階數(shù)不要太高，例如不超過1000）常用的方法，它一般用于系數(shù)矩陣稠密p1EanqFDPwENGLISH ABSTRACTGauss elimination method is the small size of solving linear equations (i.e. the order is not too high, for example, not more than 1000 commonly used method, it is generally used for the coefficient ma

2、trix dense (ie, most of the matrix elements are nonzero without any special structure of linear equations. In this paper, GAUSS, the main element and completely Elimination PCA these methods were compared. Solving n linear equations Ax = b is the basic idea of the gradual elimination process, the eq

3、uations of the coefficient matrix into upper triangular matrix, thus reduced to the original equations easy to solve trigonometric equations equivalence group then back substitution to solve.DXDiTa9E3dKeywords: Gauss elimination method。 C + + ProgrammingRTCrpUDGiT目錄消去法等解線性方程組，這就是說，如果系數(shù)矩陣的是n n，那么高斯消元

4、法所需要的計(jì)算量大約與n3成比例。LDAYtRyKfE高斯消元法對(duì)于一些矩陣來說是穩(wěn)定的。對(duì)于普遍的矩陣來說，高斯消元法在應(yīng)用上通常也是穩(wěn)定的，所以本次課程設(shè)計(jì)采用GAUSS、列主元及完全主元消去法求解線性方程組。Zzz6ZB2Ltk二、問題提出先用你所熟悉的計(jì)算機(jī)語言將不選主元、列主元和完全主元Gauss消去法編寫成通用的子程序，然后用你編寫的程序求解下面的方程組考慮從70到80）dvzfvkwMI1對(duì)上述方程組還可以采用哪些方法求解？選擇其中一些方法編程上機(jī)求解上述方程組，說明最適合的是什么方法；將計(jì)算結(jié)果進(jìn)行比較分析，談?wù)勀銓?duì)這些方法的看法。rqyn14ZNXI三、數(shù)學(xué)原理與算法分析3.

5、1不選主元Gauss消去法3.1.1不選主元Gauss消去法數(shù)學(xué)原理Gauss消去法即是把一個(gè)給定的矩陣分解為一個(gè)下三角L和一個(gè)上三角U的乘積。為此我們找到了一個(gè)最簡(jiǎn)單的下三角矩陣.其中EmxvxOtOco即1）這種類型的下三角稱作Gauss變換。而稱向量為Gauss向量。對(duì)于一個(gè)給定的向量，我們有只要取，便有。于是將Gauss變換作用于一個(gè)矩陣，就將對(duì)應(yīng)的列向量相應(yīng)項(xiàng)變?yōu)?。一般的階矩陣A,在一定條件下，我們也可以計(jì)算n-1個(gè)Gauss變換使SixE2yXPq5為上三角矩陣。，這種計(jì)算主元分解的方法稱做Gauss消去法。LU分解的可行性分析：若的順序主子陣均非奇異，則存在唯一的單位下三角和上

6、三角陣使得對(duì)于題1中，我們可以通過高等代數(shù)的方法計(jì)算出對(duì)于*.0,也即存在唯一性。6ewMyirQFL3.1.2不選主元Gauss消去法算法步驟步驟1，首先將方程組用增廣矩陣表示。步驟2: 1）輸入數(shù)據(jù)A和b，置det=1；2）消元過程，對(duì)a、選主元，找使得；b、如果，則矩陣奇異，程序結(jié)束；否則執(zhí)行3）。c、如果，則交換第行與第行對(duì)應(yīng)元素位置，。d、消元，對(duì)，計(jì)算對(duì)，計(jì)算：2）3）e、置det=*det步驟3：回代過程：1）若則矩陣奇異，程序結(jié)束；否則執(zhí)行2）。2）對(duì)，計(jì)算3.2列主元Gauss消去法3.2.1列主元Gauss消去法數(shù)學(xué)原理每次消去中，選取每列的絕對(duì)值最大的元素作為消去對(duì)象并作

7、為主元素。然后換行使之變到主元位子上，在進(jìn)行消元計(jì)算。設(shè)，確定第k列主元所在位置，在交換行和k行后，在進(jìn)行消元。根據(jù)矩陣?yán)碚?，交換和k兩個(gè)方程的位置，列主元素的消去過程相當(dāng)于對(duì)交換后的新矩陣進(jìn)行消元，即kavU42VRUs 4）同時(shí)，右端向量變化為 5）3.2.1列主元Gauss消去法算法描述如果在高斯順序消去法消去過程進(jìn)行到第i步時(shí)，現(xiàn)選取中絕對(duì)值最大的元素，設(shè)為第j行的元素，把矩陣的第i行和第j行互換，這時(shí)變?yōu)?，然后將第i+1行至第n行中的每一行減去第i行乘以k代表行號(hào)），依次進(jìn)行消元。y6v3ALoS89Gauss列主元消去法的算法步驟如下：步驟1：將方程組寫成以下的增廣矩陣的形式： 6

8、）步驟2：對(duì)k=1，2，3,.,n-1,令；交換增廣矩陣的第k行與第p行；步驟3：對(duì)j=k+1，k+2，.,n,計(jì)算m=看，k+1，.，n）；步驟4：算法結(jié)束。3.3完全主元Gauss消去法3.3完全主元Gauss消去法數(shù)學(xué)原理設(shè)方程組的增廣矩陣為7）首先在中選取絕對(duì)值最大的元素作為主元素，例如，然后交換的第行與第行，經(jīng)第一次消元計(jì)算得8）重復(fù)上述過程，已完成第步的選主元素，交換兩行及交換兩列，消元計(jì)算，約化為9）其中元素仍記作，元素仍記作。第步選主元素在右下角方框內(nèi)選），即確定，使10）交換第行與行元素交換第列與列元素，將調(diào)到位置，再進(jìn)行消元計(jì)算，最后將原方程化為11）其中，的次序?yàn)槲粗獢?shù)，

9、調(diào)換后的次序。回代求解得12）3.3.1算法描述如下設(shè)。本算法用的帶有行、列交換的Gauss消去法，消元結(jié)果沖掉，乘數(shù)沖掉，計(jì)算解沖掉常數(shù)項(xiàng)，用表示對(duì)的消元次數(shù)。用一整型數(shù)組開始記錄未知數(shù)，的次序即下標(biāo)，2，），最后記錄調(diào)換后未知數(shù)的下標(biāo)。M2ub6vSTnP步驟1 對(duì)于=1,2，有；對(duì)于，2，做到步6；步驟2 選主元素；步驟3 如果，則計(jì)算停止這時(shí)）；步驟4 1）如果，則轉(zhuǎn)2），否則換行：，；2）如果，則轉(zhuǎn)步5，否則換行：，；步驟5 計(jì)算乘數(shù)；步驟6 消元計(jì)算；。步驟7 回代求解1）；2）對(duì)于，。步驟8 調(diào)整未知數(shù)的次序1）對(duì)于；，；2）對(duì)于；。3.4計(jì)算結(jié)果及敏度分析在對(duì)題目用高斯消去法C

10、+編程見附錄一）后運(yùn)行，結(jié)果用matlab作圖后顯示可以看出，在初始分量，效果還是比較理想的逼近于x=1,到后面出現(xiàn)了震蕩?？梢姷贸龅慕鈹?shù)值不太穩(wěn)定。這是由于系數(shù)矩陣的條件0YujCfmUCw圖1 N=84 Gauss消去法不選主元）數(shù)k(A)非常大，用matlab計(jì)算= 4.1779e+016，已經(jīng)非常大，也即該線形方程組得求解問題是病態(tài)的，我們得到的數(shù)值解是不穩(wěn)定的。而應(yīng)用列主元高斯消去法，當(dāng)矩陣是84階時(shí)，不能得出結(jié)果，這是因?yàn)闂l件數(shù)很大時(shí)矩陣接近于一個(gè)奇異矩陣，高斯消去法不一定能eUts8ZQVRd圖2 N=84 Gauss全主元消去法夠進(jìn)行下去。當(dāng)為30階時(shí)，可以看出也得到了精確解。

11、當(dāng)DIM=84時(shí)程序最后輸出結(jié)果為：The linear system hasnt solution!Strike any key to exit!Press any key to continue而且，對(duì)于方程系數(shù)矩陣是嚴(yán)格對(duì)角占優(yōu)的，條件數(shù)k(A)=Cond(A=1.4997,比較小，可知此問題數(shù)值比較穩(wěn)定。由于系數(shù)矩陣嚴(yán)格對(duì)角的方程組Gauss-seidel迭代收斂，所以可以應(yīng)用此算法來解決具體程序見附錄）。從計(jì)算結(jié)果的圖形表示我們可以看出，此種方法達(dá)到了數(shù)值要求，比較好的收斂于線性方程的解。sQsAEJkW5T四、其他方法平方根法就是利用對(duì)稱正定矩陣的三角分解而得到的求解對(duì)稱正定方程組

12、的一種有效方法。平方根法遞推公式可以證明對(duì)于對(duì)稱正定矩陣A,可以唯一地分解成A=LLT,其中L是非奇異下三角形矩陣。GMsIasNXkACholesky分解定理：若對(duì)稱正定，則存在一個(gè)對(duì)角元均為正數(shù)的下三角陣，使得。Cholesky分解可用不選主元的Gauss消去法來實(shí)現(xiàn)，更簡(jiǎn)單的方法是比較A=LLT兩邊對(duì)應(yīng)元素來計(jì)算L的，設(shè)TIrRGchYzg13）比較兩邊對(duì)應(yīng)元素,得關(guān)系式。經(jīng)計(jì)算得，這樣便求出了L的第k列元素。這種方法稱為平方根法。圖3 N=84 Gauss平方根法應(yīng)用平方根法時(shí)的程序見附錄四。由于b是計(jì)算機(jī)隨機(jī)選取的，所以得到的解是人為無法確定的，我們?cè)谶@里只畫出了一種解的情況下的圖象

13、見圖2）7EqZcWLZNX用改進(jìn)的平方根法解Hilbert矩陣的程序見附錄五，相應(yīng)的圖形見上圖。有趣的是，當(dāng)應(yīng)用平方根解此方程時(shí)并不會(huì)得到結(jié)果。lzq7IGf02E應(yīng)著重指出，正定矩陣的分解具有重要的性質(zhì)，即，這說明的元素的絕對(duì)值不會(huì)變得很大，這個(gè)性質(zhì)對(duì)計(jì)算方法的穩(wěn)定性是非常重要的。zvpgeqJ1hk平方根法可以免去高斯消元過程。但是只限于對(duì)稱方程組。高斯消元法是比較經(jīng)典的求解稠密線性方程組的方法。為避免小主元，應(yīng)選用高斯列主元法。NrpoJac3v1五、參考文獻(xiàn)1李慶揚(yáng)，王能超，易大義.數(shù)值分析M.武漢：華中科技大學(xué)出版社.2006.72蕭樹鐵.數(shù)學(xué)實(shí)驗(yàn)M.北京：高等教育出版社.1990

14、3萊Lay D.C.），劉深泉.數(shù)值線性代數(shù)及其應(yīng)用M.機(jī)械工業(yè)出版社.2005.84 黃友謙，程詩(shī)杰，陳浙鵬，數(shù)值實(shí)驗(yàn)，北京：高等教育出版社，19895 蔡大用，數(shù)值分析與實(shí)驗(yàn)學(xué)習(xí)指導(dǎo)，北京：清華大學(xué)出版社與施普林格出版社，20016 肖筱南，數(shù)值計(jì)算方法與上機(jī)實(shí)習(xí)指導(dǎo)，北京：北京大學(xué)出版社，2004六、附錄附錄一：Gauss消去法不選主元）#include#include#include#ifndef Epsilon#define Epsilon 0.#endif#define DIM 84intgcplim(intprocess,double ADIMDIM,double xxDIM1n

15、owfTG4KIintk,i,j。/double pelement。if(process=1coutThe process of elimination:n。for(k=0。k/*pelement=fabs(Akk。i0=k。for(i=k。iif(fabs(Aikpelementpelement=fabs(Aik。i0=i。fjnFLDa5Zoif(i0!=kfor(j=0。jpelement=Akj。Akj=Ai0j。Ai0j=pelement。tfnNhnE6e5pelement=xxk。xxk=xxi0。xxi0=pelement。*/if(fabs(Akk return(1。for(

16、i=k+1。iAik=Aik/Akk。for(j=k+1。jAij=Aij-Aik*Akj。xxi=xxi-Aik*xxk。if(process=1for(i=0。ifor(j=0。jcoutAijt。cout |xxi。coutendl。cout=0。i-/xxDIM-1=xxDIM-1/ADIM-1DIM-1。for(i=DIM-2。i=0。i-for(j=i+1。jxxi=xxi-Aij*xxj。xxi=xxi/Aii。return(0。void main(int i。static double ADIMDIM。for(i=1。iAii=6。Aii+1=1。Aii-1=8。A00=6。A

17、01=1。ADIM-1DIM-1=6。ADIM-1DIM-2=8。HbmVN777sLstatic double bDIM。for(i=1。ibi=15。b0=7。bDIM-1=14。coutColume=1coutThe linear system hasnt solution!endl。coutStrike any key to exit!。coutColumn principle elimination for the solution:n。V7l4jRB8Hsfor(i=0。icoutbiendl。附錄二：Gauss列主元消去法#include#include#include#ifnd

18、ef Epsilon#define Epsilon 0.#endif#define DIM 84intgcplim(intprocess,double ADIMDIM,double xxDIM83lcPA59W9int k,i,j,i0。doublepelement。if(process=1coutThe process of elimination:n。for(k=0。kpelement=fabs(Akk。i0=k。for(i=k。iif(fabs(Aikpelementpelement=fabs(Aik。i0=i。mZkklkzaaPif(i0!=kfor(j=0。jpelement=Ak

19、j。Akj=Ai0j。Ai0j=pelement。AVktR43bpwpelement=xxk。xxk=xxi0。xxi0=pelement。if(fabs(Akk return(1。for(i=k+1。iAik=Aik/Akk。for(j=k+1。jAij=Aij-Aik*Akj。xxi=xxi-Aik*xxk。if(process=1for(i=0。ifor(j=0。jcoutAijt。cout |xxi。coutendl。cout=0。i-/xxDIM-1=xxDIM-1/ADIM-1DIM-1。/for(i=DIM-2。i=0。i-/for(j=i+1。jxxi=xxi-Aij*xxj

20、。xxi=xxi/Aii。return(0。void main(int i。static double ADIMDIM。for(i=1。iAii=6。Aii+1=1。Aii-1=8。A00=6。A01=1。ADIM-1DIM-1=6。ADIM-1DIM-2=8。ORjBnOwcEdstatic double bDIM。for(i=1。ibi=15。b0=7。bDIM-1=14。coutColume=1coutThe linear system hasnt solution!endl。coutStrike any key to exit!。coutColumn principle elimina

21、tion for the solution:n。2MiJTy0dTTfor(i=0。icoutbiendl。當(dāng)DIM=84時(shí)程序最后輸出結(jié)果為：The linear system hasnt solution!Strike any key to exit!Press any key to continue附錄三：全主元Gauss消去法#include #include #define N 4 /*方程階數(shù)*/void printf_ab(double aN,double *b。/*打印運(yùn)算矩陣a,b*/gIiSpiue7Avoid get_CEMEresult(double aN,double

22、 *b,int *Label。/*回代求解*/uEh0U1Yfmhvoid CEME(double aN,double *b,int* Label。/完全選主元素的消去法 Complete elimination of major elementsIAg9qLsgBXvoidchoose_ME(double aN,double *b,int *Label,int k。WwghWvVhPE/-/-測(cè)試數(shù)組 /*double aNN= 2,1, 2,2 。double bN = 3,4。/yes 2double aNN=3,2,1, 5,6,4, 7,8,9 。double bN = 5,6,7

23、。/yes 3double aNN=4,-1,1, -1,4.25,2.75, 1,2.75,3.5 。double bN = 6,-0.5,1.25 。/yes 3double aNN=2,1,1,0, 4,3,3,1, 8,7,9,5, 6,7,9,8。double bN = 1,2,2,-1 。/yes 4double aNN= 2,-1,0,0,0, -1,2,-1,0,0, 0,-1,2,-1,0, 0,0,-1,2,-1, 0,0,0,-1,2。double bN = 1,0,0,0,0 。 /yes 5double aNN= 0.4096, 0.1234, 0.3678, 0.2

24、943, 0.2246, 0.3872, 0.4015, 0.1129, 0.3645, 0.1920, 0.3781, 0.0643, 0.1784, 0.4002, 0.2786, 0.3927。double bN = 0.4043, 0.1550, 0.4240, -0.2557。 /yes 4asfpsfpi4k */void main( double bN = 1,2,2,-1 。/yes 4double aNN=2,1,1,0, 4,3,3,1, 8,7,9,5, 6,7,9,8。int LabelN。 /用于保存未知數(shù)x的次序printf(變換前，增廣矩陣A|B:n。printf

25、_ab(a,b。 CEME(a,b,Label。/完全選主元素的消去法printf(變換后，增廣矩陣A|B:n。printf_ab(a,b。get_CEMEresult(a,b,Label。/計(jì)算線性方程組的解，并順序打印/打印運(yùn)算增廣矩陣a|b/voidprintf_ab(double aN,double *binti,j。for(i=0。i for(j=0。j printf(%4.3f ,aij。 printf(%4.3f ,bi。printf(n。 /計(jì)算線性方程組的解，并順序打印/ /voidget_CEMEresult(double aN,double *b,int *Labeloo

26、eyYZTjj1inti,j。double xN, sum 。xN-1=bN-1/aN-1N-1。for(i=N-2。i=0。i- sum=0。for(j=i+1。jsum+=aij*xj。xi=(bi-sum/aii。 printf(計(jì)算并調(diào)整x的次序后，結(jié)果為:n。for(i=0。i for(j=0。j if(i=Labelj printf(x%d=%6.5fn,i+1,xj。break。 /Complete elimination of major elements/完全主元素消去法/void CEME(double aN,double *b , int *Labelint i=0,j=

27、0,k=0。double l=0.0。/-for(k=0。kLabelk = k。 /初始化/-for(k=0。k choose_ME(a,b,Label,k。for(i=k+1。i l=aik/akk。bi-=l*bk。for(j=k。j aij -= l*akj。 /選取主元素，并進(jìn)行行、列的交換/voidchoose_ME(double aN,double *b,int *Label,int kBkeGuInkxIinti,j,Di,Dj。doubletempa。int temp 。 Di=Dj=k。tempa=akk。/akkfor(i=k。i for(j=k。j if(fabs(tempa tempa=aij。 Di=i。/獲取列主元素Dj=j。 / printf(nk=%d Di=%d,Dj=%d ,k,Di,Dj。if(tempa-0.0001 & tempa printf(輸入數(shù)組有問題，因?yàn)閐etA=0n。return 。if(Di!=k for(j=k。j /at行 ak , aki =0 (i=0,1,.,k-1PgdO0sRlMo tempa=akj。akj