消去法等解线性方程组(一)

1、精选优质文档-倾情为你奉上消去法等解线性方程组<一）中文摘要Gauss消去法是目前求解中小规模线性方程组<即阶数不要太高，例如不超过1000）常用的方法，它一般用于系数矩阵稠密<即矩阵的绝大多数元素都是非零的）而没有任何特殊结构的线性方程组。本文对GAUSS、列主元及完全主元消去法这几种方法进行了比较。求解n阶线性方程组Ax=b的基本思想是在逐步消元的过程中，把方程组的系数矩阵化为上三角矩阵，从而将原方程组约化为容易求解的等价三角方程组，然后进行回代求解。b5E2RGbCAP关键词：Gauss消去法； C+编程ELIMINATION METHOD FOR SOLVING LI

2、NEAR EQUATIONS (一>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 matrix dense (ie, most of the matrix elements are

3、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 equations of the coefficient matrix into upper

4、 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目录消去法等解线性方程组<一）1目录3一、问题背景4二、问题提出5三、数学原理与算法分析53.1不选主元Gauss消去法53.1.1不选主元Gauss消去

5、法数学原理53.1.2不选主元Gauss消去法算法步骤63.2列主元Gauss消去法73.2.1列主元Gauss消去法数学原理73.2.1列主元Gauss消去法算法描述83.3完全主元Gauss消去法83.3完全主元Gauss消去法数学原理83.3.1算法描述如下103.4计算结果及敏度分析11四、其他方法13五、参考文献14六、附录15一、问题背景在古代数学和近代数学数值计算和工程应用中，求解线性方程组都是重要的课题。西方的Gauss消元法在求解未知数多的大型线形方程组则是在20世纪中叶电子计算机问世后才成为现实。这次实验将充分运用各种数值计算方法，及计算机的强大数值计算功能来解决所遇到的问

6、题。5PCzVD7HxA数学上，高斯消元法，是线性代数中的一个算法，可用来为线性方程组求解，求出矩阵的秩，以及求出可逆方阵的逆矩阵。其基本思路如下：jLBHrnAILg设有方程组Ax= b，设A是可逆矩阵。将矩阵的初等行变换作用于方程组的增广矩阵，将其中的A变换成一个上三角矩阵，然后求解这个三角形方程组xHAQX74J0X高斯消元法的算法复杂度是O(n3>，这就是说，如果系数矩阵的是n × n，那么高斯消元法所需要的计算量大约与n3成比例。LDAYtRyKfE高斯消元法对于一些矩阵来说是稳定的。对于普遍的矩阵来说，高斯消元法在应用上通常也是稳定的，所以本次课程设计采用GAUSS

7、、列主元及完全主元消去法求解线性方程组。Zzz6ZB2Ltk二、问题提出先用你所熟悉的计算机语言将不选主元、列主元和完全主元Gauss消去法编写成通用的子程序，然后用你编写的程序求解下面的方程组<考虑从70到80）dvzfvkwMI1对上述方程组还可以采用哪些方法求解？选择其中一些方法编程上机求解上述方程组，说明最适合的是什么方法；将计算结果进行比较分析，谈谈你对这些方法的看法。rqyn14ZNXI三、数学原理与算法分析3.1不选主元Gauss消去法3.1.1不选主元Gauss消去法数学原理Gauss消去法即是把一个给定的矩阵分解为一个下三角L和一个上三角U的乘积。为此我们找到了一个最简

8、单的下三角矩阵.其中EmxvxOtOco即<1）这种类型的下三角称作Gauss变换。而称向量为Gauss向量。对于一个给定的向量，我们有只要取，便有。于是将Gauss变换作用于一个矩阵，就将对应的列向量相应项变为0。一般的阶矩阵A,在一定条件下，我们也可以计算n-1个Gauss变换使SixE2yXPq5为上三角矩阵。，这种计算主元分解的方法称做Gauss消去法。LU分解的可行性分析：若的顺序主子阵均非奇异，则存在唯一的单位下三角和上三角阵使得对于题1中，我们可以通过高等代数的方法计算出对于*.0,也即存在唯一性。6ewMyirQFL3.1.2不选主元Gauss消去法算法步骤步骤1，首先将

9、方程组用增广矩阵表示。步骤2: <1）输入数据A和b，置det=1；<2）消元过程，对a、选主元，找使得；b、如果，则矩阵奇异，程序结束；否则执行<3）。c、如果，则交换第行与第行对应元素位置，。d、消元，对，计算对，计算：<2）<3）e、置det=*det步骤3：回代过程：<1）若则矩阵奇异，程序结束；否则执行<2）。<2）对，计算3.2列主元Gauss消去法3.2.1列主元Gauss消去法数学原理每次消去中，选取每列的绝对值最大的元素作为消去对象并作为主元素。然后换行使之变到主元位子上，在进行消元计算。设，确定第k列主元所在位置，在交换行和k

10、行后，在进行消元。根据矩阵理论，交换和k两个方程的位置，列主元素的消去过程相当于对交换后的新矩阵进行消元，即kavU42VRUs <4）同时，右端向量变化为 <5）3.2.1列主元Gauss消去法算法描述如果在高斯顺序消去法消去过程进行到第i步时，现选取中绝对值最大的元素，设为第j行的元素，把矩阵的第i行和第j行互换，这时变为，然后将第i+1行至第n行中的每一行减去第i行乘以<k代表行号），依次进行消元。y6v3ALoS89Gauss列主元消去法的算法步骤如下：步骤1：将方程组写成以下的增广矩阵的形式： <6）步骤2：对k=1，2，3,.,n-1,令；交换增广矩阵的第k

11、行与第p行；步骤3：对j=k+1，k+2，.,n,计算<m=看，k+1，.，n）；步骤4：算法结束。3.3完全主元Gauss消去法3.3完全主元Gauss消去法数学原理设方程组的增广矩阵为<7）首先在中选取绝对值最大的元素作为主元素，例如，然后交换的第行与第行，经第一次消元计算得<8）重复上述过程，已完成第步的选主元素，交换两行及交换两列，消元计算，约化为<9）其中元素仍记作，元素仍记作。第步选主元素<在右下角方框内选），即确定，使<10）交换第行与行元素交换第列与列元素，将调到位置，再进行消元计算，最后将原方程化为<11）其中，的次序为未知数，调换后

12、的次序。回代求解得<12）3.3.1算法描述如下设。本算法用的带有行、列交换的Gauss消去法，消元结果冲掉，乘数冲掉，计算解冲掉常数项，用表示对的消元次数。用一整型数组开始记录未知数，的次序<即下标，2，），最后记录调换后未知数的下标。M2ub6vSTnP步骤1 对于=1,2，有；对于，2，做到步6；步骤2 选主元素；步骤3 如果，则计算停止<这时）；步骤4 <1）如果，则转<2），否则换行：，；<2）如果，则转步5，否则换行：，；步骤5 计算乘数；步骤6 消元计算；。步骤7 回代求解<1）；<2）对于，。步骤8 调整未知数的次序<1）对

13、于；，；<2）对于；。3.4计算结果及敏度分析在对题目用高斯消去法C+编程<见附录一）后运行，结果用matlab作图后显示可以看出，在初始分量，效果还是比较理想的逼近于x=1,到后面出现了震荡。可见得出的解数值不太稳定。这是由于系数矩阵的条件0YujCfmUCw图1 N=84 Gauss消去法<不选主元）数k(A)非常大，用matlab计算<程序见附录六）可得k(A)=cond(A>= 4.1779e+016，已经非常大，也即该线形方程组得求解问题是病态的，我们得到的数值解是不稳定的。而应用列主元高斯消去法，当矩阵是84阶时，不能得出结果，这是因为条件数很大时矩阵

14、接近于一个奇异矩阵，高斯消去法不一定能eUts8ZQVRd图2 N=84 Gauss全主元消去法够进行下去。当为30阶时，可以看出也得到了精确解。当DIM=84时程序最后输出结果为：The linear system hasn't solution!Strike any key to exit!Press any key to continue而且，对于方程系数矩阵是严格对角占优的，条件数k(A)=Cond(A>=1.4997,比较小，可知此问题数值比较稳定。由于系数矩阵严格对角的方程组Gauss-seidel迭代收敛，所以可以应用此算法来解决<具体程序见附录）。从计算结果

15、的图形表示我们可以看出，此种方法达到了数值要求，比较好的收敛于线性方程的解。sQsAEJkW5T四、其他方法平方根法就是利用对称正定矩阵的三角分解而得到的求解对称正定方程组的一种有效方法。平方根法递推公式可以证明对于对称正定矩阵A,可以唯一地分解成A=LLT,其中L是非奇异下三角形矩阵。GMsIasNXkACholesky分解定理：若对称正定，则存在一个对角元均为正数的下三角阵，使得。Cholesky分解可用不选主元的Gauss消去法来实现，更简单的方法是比较A=LLT两边对应元素来计算L的，设TIrRGchYzg<13）比较两边对应元素,得关系式。经计算得，这样便求出了L的第k列元素。

16、这种方法称为平方根法。图3 N=84 Gauss平方根法应用平方根法时的程序见附录四。由于b是计算机随机选取的，所以得到的解是人为无法确定的，我们在这里只画出了一种解的情况下的图象<见图2）7EqZcWLZNX用改进的平方根法解Hilbert矩阵的程序见附录五，相应的图形见上图。有趣的是，当应用平方根解此方程时并不会得到结果。lzq7IGf02E应着重指出，正定矩阵的分解具有重要的性质，即，这说明的元素的绝对值不会变得很大，这个性质对计算方法的稳定性是非常重要的。zvpgeqJ1hk平方根法可以免去高斯消元过程。但是只限于对称方程组。高斯消元法是比较经典的求解稠密线性方程组的方法。为避免

17、小主元，应选用高斯列主元法。NrpoJac3v1五、参考文献1李庆扬，王能超，易大义.数值分析M.武汉：华中科技大学出版社.2006.72萧树铁.数学实验M.北京：高等教育出版社.19903莱<Lay D.C.），刘深泉.数值线性代数及其应用M.机械工业出版社.2005.84 黄友谦，程诗杰，陈浙鹏，数值实验，北京：高等教育出版社，19895 蔡大用，数值分析与实验学习指导，北京：清华大学出版社与施普林格出版社，20016 肖筱南，数值计算方法与上机实习指导，北京：北京大学出版社，2004六、附录专心-专注-专业附录一：Gauss消去法<不选主元）#include<iostr

18、eam.h>#include<stdlib.h>#include<math.h>#ifndef Epsilon#define Epsilon 0.#endif#define DIM 84intgcplim(intprocess,double ADIMDIM,double xxDIM>1nowfTG4KIintk,i,j。/double pelement。if(process=1>cout<<"The process of elimination:n"。for(k=0。k<DIM。k+>/*pelement=f

19、abs(Akk>。i0=k。for(i=k。i<DIM。i+>if(fabs(Aik>>pelement>pelement=fabs(Aik>。i0=i。fjnFLDa5Zoif(i0!=k>for(j=0。j<DIM。j+>pelement=Akj。Akj=Ai0j。Ai0j=pelement。tfnNhnE6e5pelement=xxk。xxk=xxi0。xxi0=pelement。*/if(fabs(Akk><Epsilon> return(1>。for(i=k+1。i<DIM。i+>Aik=

20、Aik/Akk。for(j=k+1。j<DIM。j+>Aij=Aij-Aik*Akj。xxi=xxi-Aik*xxk。if(process=1>for(i=0。i<DIM。i+>for(j=0。j<DIM。j+>cout<<Aij<<'t'。cout<<" |"<<xxi。cout<<endl。cout<<endl。/for(i=DIM-1。i>=0。i->/xxDIM-1=xxDIM-1/ADIM-1DIM-1。for(i=DIM-

21、2。i>=0。i->for(j=i+1。j<DIM。j+>xxi=xxi-Aij*xxj。xxi=xxi/Aii。return(0>。void main(>int i。static double ADIMDIM。for(i=1。i<DIM。i+>Aii=6。Aii+1=1。Aii-1=8。A00=6。A01=1。ADIM-1DIM-1=6。ADIM-1DIM-2=8。HbmVN777sLstatic double bDIM。for(i=1。i<DIM-1。i+>bi=15。b0=7。bDIM-1=14。cout<<&quo

22、t;Colume"<<endl。if(gcplim(1,A,b>=1>cout<<"The linear system hasn't solution!"<<endl。cout<<"Strike any key to exit!"<<endl。exit(1>。cout<<"Column principle elimination for the solution:n"。V7l4jRB8Hsfor(i=0。i<DIM。i+&

23、gt;cout<<bi<<endl。附录二：Gauss列主元消去法#include<iostream.h>#include<stdlib.h>#include<math.h>#ifndef Epsilon#define Epsilon 0.#endif#define DIM 84intgcplim(intprocess,double ADIMDIM,double xxDIM>83lcPA59W9int k,i,j,i0。doublepelement。if(process=1>cout<<"The pr

24、ocess of elimination:n"。for(k=0。k<DIM。k+>pelement=fabs(Akk>。i0=k。for(i=k。i<DIM。i+>if(fabs(Aik>>pelement>pelement=fabs(Aik>。i0=i。mZkklkzaaPif(i0!=k>for(j=0。j<DIM。j+>pelement=Akj。Akj=Ai0j。Ai0j=pelement。AVktR43bpwpelement=xxk。xxk=xxi0。xxi0=pelement。if(fabs(Akk&g

25、t;<Epsilon> return(1>。for(i=k+1。i<DIM。i+>Aik=Aik/Akk。for(j=k+1。j<DIM。j+>Aij=Aij-Aik*Akj。xxi=xxi-Aik*xxk。if(process=1>for(i=0。i<DIM。i+>for(j=0。j<DIM。j+>cout<<Aij<<'t'。cout<<" |"<<xxi。cout<<endl。cout<<endl。for(i=

26、DIM-1。i>=0。i->/xxDIM-1=xxDIM-1/ADIM-1DIM-1。/for(i=DIM-2。i>=0。i->/for(j=i+1。j<DIM。j+>xxi=xxi-Aij*xxj。xxi=xxi/Aii。return(0>。void main(>int i。static double ADIMDIM。for(i=1。i<DIM。i+>Aii=6。Aii+1=1。Aii-1=8。A00=6。A01=1。ADIM-1DIM-1=6。ADIM-1DIM-2=8。ORjBnOwcEdstatic double bDIM。f

27、or(i=1。i<DIM-1。i+>bi=15。b0=7。bDIM-1=14。cout<<"Colume"<<endl。if(gcplim(1,A,b>=1>cout<<"The linear system hasn't solution!"<<endl。cout<<"Strike any key to exit!"<<endl。exit(1>。cout<<"Column principle elimin

28、ation for the solution:n"。2MiJTy0dTTfor(i=0。i<DIM。i+>cout<<bi<<endl。当DIM=84时程序最后输出结果为：The linear system hasn't solution!Strike any key to exit!Press any key to continue附录三：全主元Gauss消去法#include <stdio.h>#include <math.h>#define N 4 /*方程阶数*/void printf_ab(double a

29、N,double *b>。/*打印运算矩阵a,b*/gIiSpiue7Avoid get_CEMEresult(double aN,double *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/-/-测试数组 /*double

30、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 。/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

31、,2,-1, 0,0,0,-1,2。double bN = 1,0,0,0,0 。 /yes 5double aNN= 0.4096, 0.1234, 0.3678, 0.2943, 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

32、 aNN=2,1,1,0, 4,3,3,1, 8,7,9,5, 6,7,9,8。int LabelN。 /用于保存未知数x的次序printf("变换前，增广矩阵A|B:n">。printf_ab(a,b>。 CEME(a,b,Label>。/完全选主元素的消去法printf("变换后，增广矩阵A|B:n">。printf_ab(a,b>。get_CEMEresult(a,b,Label>。/计算线性方程组的解，并顺序打印/打印运算增广矩阵a|b/voidprintf_ab(double aN,double *b>

33、inti,j。for(i=0。i<N。i+> for(j=0。j<N。j+> printf("%4.3f ",aij>。 printf("%4.3f ",bi>。printf("n">。 /计算线性方程组的解，并顺序打印/ /voidget_CEMEresult(double aN,double *b,int *Label>ooeyYZTjj1inti,j。double xN, sum 。xN-1=bN-1/aN-1N-1。for(i=N-2。i>=0。i-> sum=0。f

34、or(j=i+1。j<N。j+>sum+=aij*xj。xi=(bi-sum>/aii。 printf("计算并调整x的次序后，结果为:n">。for(i=0。i<N。i+> for(j=0。j<N。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 *Label>int

35、 i=0,j=0,k=0。double l=0.0。/-for(k=0。k<N。k+>Labelk = k。 /初始化/-for(k=0。k<N-1。k+> choose_ME(a,b,Label,k>。for(i=k+1。i<N。i+> l=aik/akk。bi-=l*bk。for(j=k。j<N。j+> aij -= l*akj。 /选取主元素，并进行行、列的交换/voidchoose_ME(double aN,double *b,int *Label,int k>BkeGuInkxIinti,j,Di,Dj。doubletemp

36、a。int temp 。 Di=Dj=k。tempa=akk。/akkfor(i=k。i<N。i+> for(j=k。j<N。j+> if(fabs(tempa><fabs(aij>> tempa=aij。 Di=i。/获取列主元素Dj=j。 / printf("nk=%d Di=%d,Dj=%d ",k,Di,Dj>。if(tempa>-0.0001 && tempa <0.0001> printf("输入数组有问题，因为detA=0n">。return 。i

37、f(Di!=k> for(j=k。j<N。j+> /at行 <-> ak , aki =0 (i=0,1,.,k-1>PgdO0sRlMo tempa=akj。akj=aDij。aDij=tempa。 /bk行 <-> attempa=bk。bk=bDi。bDi=tempa。 if(Dj!=k>/need change Label for(j=0。j<N。j+> /a*k列 <-> a*Dj列 tempa=ajk。ajk=ajDj。ajDj=tempa。 / 记录未知数x的次序temp = Labelk 。Labelk = LabelDj。LabelDj = temp 。 return 。附录四：平方根法程序：#include<iostream.h>#include<stdlib.h>#include<math.h>#ifndef EPSILON#defi