科学仿真实验五正文NO1

消去法等解线性方程组（二）第1页共11页一引言在古代数学和近代数学数值计算和工程应用中，求解线性方程组都是重要的课题。西方的Gauss消元法在求解未知数多的大型线形方程组则是在20世纪中叶电子计算机问世后才成为现实。这次实验将充分运用各种数值计算方法，及计算机的强大数值计算功能来解决所遇到的问题。因而使数值计算方法和相应的数学软件包都产生了变化，。本文讨论使用不列主元、列主元和Guass消去法解线性方程组的一般子程序，并用C语言算法描述，用Visualc++6.0编译，用消去法求解具体的线性方程组的解并成功运行程序求解方程。二模型建立与数学原理：（1）Gauss消去法Gauss消去法即是把一个给定的矩阵分解为一个下三角L和一个上三角U的乘积。为此我们找到了一个最简单的下三角矩阵其中即这种类型的下三角称作Gauss变换。而称向量kl为Gauss向量。对于一个给定的向量，我们有只要取，便有。于是将Gauss变换作用于一个矩阵，就将对应的列向量相应项变为0。一般的阶矩阵A,在一定条件下，我们也可以计算n-1个Gauss变换使为上三角矩阵。ALU，这种计算主元分解的方法称做Gauss消去法。LU分解的可行性分析：消去法等解线性方程组（二）第2页共11页若的顺序主子阵L均非奇异，则存在唯一的单位下三角和上三角阵使得ALU对于题1中，我们可以通过高等代数的方法计算出对于,也即存在唯一性。(2)列主元Gauss消去法由于可能存在小主元现象，在Gauss消去的过程除数的绝对值过小，则导致其他变量量级的巨大增长和舍入误差扩散。因而要保持数值解得精度应避免选取绝对值过小的元素作主元。(3)平方根法Cholesky分解定理：若对称正定，则存在一个对角元均为正数的下三角阵，使得。Cholesky分解可用不选主元的Gauss消去法来实现，更简单的方法是比较两边对应元素来计算L的，设比较两边对应元素,得关系式。经计算得这样便求出了L的第k列元素。这种方法称为平方根法。三实验环境：VisualC++6.0Matlab7.0四实验过程及实验分析(1)高斯消去法算法原理在1中对题目用高斯消去法C++编程后运行，结果用matlab作图后显示可以消去法等解线性方程组（二）第3页共11页看出，在初始分量，效果还是比较理想的逼近于x=1,到后面出现了震荡。可见得出的解数值不太稳定。这是由于系数矩阵的条件数非常大，用matlab计算可得，已经非常大，也即该线形方程组得求解问题是病态的，我们得到的数值解是不稳定的。而应用列主元高斯消去法，当矩阵是84阶时，不能得出结果，这是因为条件数很大时矩阵接近于一个奇异矩阵，高斯消去法不一定能够进行下去。当为30阶时，可以看出也得到了精确解（附录中只提供了84阶矩阵，相应30阶读者可以自行改动）。（2）Gauss-seidel迭代算法原理线性方程组系数矩阵是严格对角占优的，条件数,比较小，可知此问题数值比较稳定。由于细数矩阵严格对角的方程组Gauss-seidel迭代收敛，所以可以应用此算法来解决（具体程序见附录）。从计算结果的图形表示我们可以看出，此种方法达到了数值要求，比较好的收敛于线性方程的解。（3）平方根分解法算法原理应用平方根法时的程序见附录四。由于b是计算机随机选取的，所以得到的解是人为无法确定的，我们在这里只画出了一种解的情况下的图象用改进的平方根法解Hilbert矩阵，有趣的是，当我们用平方根法解十并不会得到结果。应着重指出，正定矩阵的分解具有重要的性质，即，这说明L的元素的绝对值不会变得很大，这个性质对计算方法的稳定性是非常重要的。平方根法可以免去高斯消元过程。但是只限于对称方程组。高斯消元法是比较经典的求解稠密线性方程组的方法。为避免小主元，应选用高斯列主元法。（4）列主元消去算法原理设有线性方程组Ax=b，其中设A为非奇异矩阵。方程组的增广矩阵为消去法等解线性方程组（二）第4页共11页111211212222112[,]nnnnlnnnaaabaaabaaaAabb第1步（k=1）：首先在A的第一列中选取绝对值最大的元素1ia，作为第一步的主元素:1,111max0iiinaa，然后交换（A，b）的第1行与第i1行元素，再进行消元计算。设列主元素消去法已经完成第1步到第k-1步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组()()kkAxb(1)(1)(1)(1)(1)1112111(2)(2)(2)(2)22222()()()()()()()1,1()1,()()()[,][,]kkkknknkkkkkkkkknkkkknkkknnnnkaaaabaaabaabaabaabAbAb第k步计算如下：对于k=1，2，…，n-1（1）按列选主元：即确定ik使,max0ikkikkinaa（2）如果,0ikka，则A为非奇异矩阵，停止计算。（3）如果ik≠k，则交换[A，b]第ik行与第k行元素。（4）消元计算,(1,,)ikikikkkaamikna,(,1,,),(1,,)ijijikkjiiikkaamaijknbbmbikn消去法等解线性方程组（二）第5页共11页消元乘数ikm满足:1,(1,,)ikikkkamikna（5）回代求解1/(0)()/,1,,2,1n当nnnnnniiijjiijibbaabbabain计算解1(,,)nxx在常数项b(n)内得到。五总结实验的数学原理很容易理解，也容易上手。把运算的结果带入原方程组，可以发现符合的还是比较好。这说明列主元消去法计算这类方程的有效性。当A可逆时，能够将计算进行到底，列主元法就能确保算法的稳定，而且计算量不大。直接三角消去过程，实质上是将A分解为两个三角矩阵的乘积A=LU，并求解Ly=b的过程。回带过程就是求解上三角方程组Ux=y。所以在实际的运算中，矩阵LU可以直接计算出，而不需要任何中间步骤，从而在计算过程中将高斯消去法的步骤进行了进一步的简略，大大提高了运算速度，这就是三角分解法。通过以上的计算比较，方程组具有严重的病态性。当系数矩阵有微小的变化时，所得的解与原方程组的解有很大的相对误差。当系数矩阵A和b有微小变化时，所得的解与方程组的解没有相对误差。用MATLAB内部函数inv通过求逆矩阵，然后通过x=inv（A）*b也可以求出方程组的解，但是没有列主元高斯消去法具有良好的稳定性。det函数求方程组系数矩阵的行列式时所得结果和高斯消去法和所得结果相同，具有方便快捷的优点。消去法等解线性方程组（二）第6页共11页参考文献[1]严蔚敏吴伟民．数据结构(C语言版)[M].北京：清华大学出版社，1997.4[2]李庆扬，王能超，易大义.数值分析[M].武汉.华中科技大学出版社，2006.7.[3]清华大学、北京大学计算方法编写组。计算方法[M]。北京。科学出版社，1980消去法等解线性方程组（二）第7页共11页附录程序源代码附录一：Gauss消去法#includeiostream.h#includestdlib.h#includemath.h#ifndefEpsilon#defineEpsilon0.00000001#endif#defineDIM84intgcplim(intprocess,doubleA[DIM][DIM],doublexx[DIM]){intk,i,j;//doublepelement;if(process==1)coutTheprocessofelimination:\n;for(k=0kDIMk++){/*pelement=fabs(A[k][k])i0=kfor(i=k;iDIM;i++)if(fabs(A[i][k])pelement){pelement=fabs(A[i][k]);i0=i;}if(i0!=k){for(j=0;jDIM;消去法等解线性方程组（二）第8页共11页j++){pelement=A[k][j];A[k][j]=A[i0][j];A[i0][j]=pelement;}pelement=xx[k];xx[k]=xx[i0];xx[i0]=pelement;}*/if(fabs(A[k][k])Epsilon)return(1);for(i=k+1iDIMi++){A[i][k]=A[i][k]/A[k][k];for(j=k+1;jDIM;j++)A[i][j]=A[i][j]-A[i][k]*A[k][j];xx[i]=xx[i]-A[i][k]*xx[k];}if(process==1)for(i=0;iDIM;i++){for(j=0;jDIM;j++)coutA[i][j]'\t';cout|xx[i];coutendl;}coutendl;消去法等解线性方程组（二）第9页共11页}}//for(i=DIM-1;i=0i--)//{xx[DIM-1]=xx[DIM-1]/A[DIM-1][DIM-1]for(i=DIM-2;i=0;i--){for(j=i+1;jDIM;j++)xx[i]=xx[i]-A[i][j]*xx[j];xx[i]=xx[i]/A[i][i];}return(0);}voidmain(){inti;staticdoubleA[DIM][DIM];for(i=1;iDIM;i++){A[i][i]=6;A[i][i+1]=1;A[i][i-1]=8;}A[0][0]=6;A[0][1]=1;A[DIM-1][DIM-1]=6;A[DIM-1][DIM-2]=8;staticdoubleb[DIM];for(i=1;iDIM-1;i++)b[i]=15;b[0]=7;b[DIM-1]=14;coutColumeendl;消去法等解线性方程组（二）第10页共11页if(gcplim(1,A,b)==1){coutThelinearsystemhasn'tsolution!endl;coutStrikeanykeytoexit!endl;exit(1);}coutColumnprincipleeliminationforthesolution:\n;for(i=0;iDIM;i++)coutb[i]endl;}解的图象如下所示：附录二：Gauss列主元消去法#includeiostream.h#includestdlib.h#includemath.h#ifndefEpsilon#defineEpsilon0.00000001#endif#defineDIM84intgcplim(intprocess,doubleA[DIM][DIM],doublexx[DIM]){intk,i,j,i0;doublepelement;if(process==1)coutTheprocessofelimination:\n;for(k=0;kDIM;k++){pelement=fabs(A[k][k]);消去法等解线性方程组（二）第11页共11页i0=k;for(i=k;iDIM;i++)if(fabs(A[i][k])pelement){pelement=fabs(A[i][k]);i0=i;}if(i0!=k){for(j=0;jDIM;j++){pelement=A[k][j];A[k][j]=A[i0][j];A[i0][j]=pelement;}pelement=xx[k];xx[k]=xx[i0];xx[i0]=pelement;}if(fabs(A[k][k])Epsilon)return(1);for(i=k+1;iDIM;i++){A[i][k]=A[i][k]/A[k][k];for(j=k+1;jDIM;j++)A[i]