《Gauss消元法解解线性方程组.docx》由会员分享，可在线阅读，更多相关《Gauss消元法解解线性方程组.docx（20页珍藏版）》请在三一办公上搜索。
1、Gauss消元法解解线性方程组摘要 本文叙述了Gauss顺序消元法解线性方程的算法思想以及其求解过程，同时简要叙述了Gauss主元素消元法以及Gauss全主元消元法。紧接着给出了Gauss-Seidel迭代法的算法思想，本文给出了这三个消元方法以及一个迭代法的算法流程图，由于全主元消元法是前两个算法的基础上改进而来，故本文采用第三种方法进行编程计算，前两种方法不再重复编程，然后给出一个实例的计算结果，运行时间，在文章最后分析该实例的计算结果，针对同一实例，又采用Gauss-Seidel方法编程实现，然后对结果进行分析和对比。最后给出了本人在编程时遇到的一些问题和解决办法。 关键词：Gauss顺
2、序消元法 Gauss主元素消元法 Gauss全主元消元法 一、算法的简要描述 1.1Gauss顺序消元法 Gauss消元法在中学里已经学习过，其方法实质，就是运用初等变换，将线性方程组Ax=b转化为同解的上三角矩阵方程组 Ux=L-1b 其中，U为上三角矩阵，L为下三角矩阵。然后对式进行回代求解，即得方程组的解。手算的过程是非常清楚的，现在需回答的是计算机求解，如何实现上述计算过程。 设线性方程组为 a11x1+a12x2+a13x3+a1nxn=b1ax+ax+ax+ax=b2112222332nn2 an1x1+an2x2+an3x3+annxn=bn写成矩阵形式为 a11a12a21a2
3、2a21a22a1mx1b1xba2m2=2 a2mxnbn设线性方程组如上式所示，记A(1)=A，b(1)=b，与是增广矩阵具有形式Ab=A(1)b，此时方程组为A(1)x=b(1)。 (1)第一次消元。设a11为将第二个方程至第n个方程的x1系数ai1(1)消成零，0，(1)构造乘数 ai(1)1 (i=2,3,n, li1=(1)a11用-li1乘以矩阵Ab(1)的第一行加到第i行上(i=2,3,(1)a110b(2)=0(1)a12(2)a22(1)a13(2)a23,n)得 a1(1)n(2)a2nA(2)an2(2)an3(2)annb1(1)(2)b2 (2)bn其中， (2)(
4、1)aij=aij-li1a1(1)j(2)(1)(1)bi=bi-li1b1(i,j=2,3,(i,j=2,3,n),n)假设经过k-1次消元得同解方程组A(k)x=b(k)，此时 (1)a11(k)b=(1)a12(2)a22(1)a13(2)a23a1(1)n(2)a2n(k)akn(k)annA(k)ankb1(1)(2)b2 (1.1. 3(k)bk(k)bn(k)假若akk0，则第k次消元如下：构造乘数 (k)aiklik=(k)akk(k)A用-lik乘以矩阵(i=k+1,n) b(k)第k行加到第i行上(i=k+1,n)得同解方程组A(k+1)x=b(k+1)，其中，A(k+1
5、)b(k+1)中的元素计算公式为 (i,j=2,3,(i,j=2,3,n),n) (k+1)(k)(k)aij=aij-likakj(k+1)=bi(k)-likbk(k)bi如此进行n-1次消元，即k=1,2,n-1，最后得同解方程组 A(n)x=b(n) 其中， (1)a11b(k)=(1)a12(2)a22(1)a13(2)a23a1(1)n(2)a2nAx=b(n)(n)A(n)annb1(1)(2)b2 (n)bn(n)假若ann0，对式进行回代计算得方程组的解 (n)xn=bn/annn(i)(i)(i)xi=(bi-aijxj)/aiij=i+1(i=n-1,2,1)上述全过程称
6、为Gauss顺序消元法。 1.2Gauss主元素消元法 (k) Gauss顺序消元法每一步总有一个要求akk0，否则算法就失败。然而从方程组的理论，只要det(A)0，则方程组解存在且唯一。由此可知，Gauss顺序(k)(k)消元法可执行条件苛刻。不仅如此，从数值计算角度，当akk很小0，但akk(k)时。用akk做分母运算，会引起误差界的放大，数值不稳定现象产生，严重时导致解失真。为克服这一缺点，常采用选主元的消元法。 Gauss列主元消元法主要依据线性方程组任意交换两个方程的次序，方程组的同解性不变，且解的分量次序也不变。于是，第k步在顺序消元法进行之前，(k)(k)从A(k)的第k列元素
7、akk，ak+1k，(k)，ank中选取绝对值最大者，并记录所在行，即 (k)ai(kkk)=maxaik，记l=ik kin(k)A如果lk，则交换矩阵b(k)的第k行与第l行所有对应的元素，然后，再进行第k步顺序消元法。 1.3Gauss全主元消元法 Gauss全主元消元法是在Gauss列主元消元法的基础上进一步改进，目的是更好的提高数值稳定性和解的精度。对那些矩阵性态不太理想的方程组能给出最大可能性的满意解。 (k)全主元算法的第k步是从aij(i=k,k+1,n)中选取绝对值最大者作为主元，并记录主元所在列与所在行，即 (k)ai(kkj)k=maxaij，记l=ik,t=jk ki,
8、jn(k)A如果lk，则交换矩阵(k)则交换矩阵Ab(k)如果tk，的第k行与第l行所有对应的元素；b(k)的第k行与第t列所有对应的元素。值得注意的是，进行行元素的交换，解的分量次序不变，但是进行列元素交换时，解的分量次序将有所改变。 1.4 Gauss-Seidel迭代法 将方程组Ax=b改写成等价方程组 0x1x-a212=a22-a12a110a1nb1aa11x111ab-2nx22a22+a22 (1.4.1) -xnaxn-n1a-an20nnannbnann矩阵形式为 x=GJx+f 其中 0-a12-a1na11a11ab1a1121G=-0-a2nbJa22a222，f=a
9、22 0bn-an1a-an2nnannann任取初始解向量x(0)=(x(0)x(0)1,2,x(0)Tn),带入式(1.4.2)右端得迭代格式 x(k+1)=G(k)Jx+f(k=0,1,2,) 其分量形式为 nxk+1k)i=(bi-aijx(j)/aii(i=1,2,n)(k=0,1,2,) jj=1i此式称为Jacobi迭代法，简称J法。 将上式修改为 xk+1i-1(k+1)i=(bi-aijxj-j=1jnax(k)ijj)/aii(i=1,2,n) =i+1则为Gauss-Seidel迭代法。 二、计算机的基本配置 版本：win7旗舰版 安装内存(RAM):2GB 处理器：In
10、tel(R)Core(TM)i3CPUM350 2.27GHz 系统类型：32位操作系统 三、参数设置及计算结果 (1.4.2) (1.4.3) (1.4. 4 )在本文中选取一个六元方程组 x1+2x2+3x3+4x4+5x5+6x6=442x+3x+4x+5x+7x+8x=732345614x1+5x2+4x3+6x4+9x5+7x6=92 12x+34x+54x+65x+34x+44x=1223456124x1+45x2+89x3+45x4+23x5+43x6=3415x1+23x2+56x3+77x4+43x5+56x6=45这是本人随机写的一个方程组，本文将以之为例，用摘要里说的方法
11、来解方程组。3.1Gauss全主元消元法 3.1.1算法流程图： 1) 输入aij(i,j=1,2,2) z(i)=i(i=1,2,3) 对k=1,2,n),bi(i=1,2,n)； ,n)； ,n-1做 ki,jna) 选主元：aikjk=maxaij，并记l=ik,t=jk，即 amax=akk;l=k;t=k; 对i=k+1,对j=k+1,n做 ,n做 if aijamax then amax=aij,l=i,t=j b) if amax=0，输出奇异标志，停机； c) if lk then akjalj(j=k,k+1,if tk then aikait(i=1,2,d) 对i=k+1
12、,i. ii. iii. ,n);bkbl； ,n);z(k)z(t)； ,n做 aiklik=aik/akk； bi=bi-likbk 对j=k+1,n做aij=aij-likakj； 4) if ann=0 then 输出奇异标志，并停机 else 做 a) bn=bn/ann； b) bi=(bi-j=i+1ab)/a(i=n-1,ijjiin,1); c) x(z(i)=bi(i=1,2,5) 输出解xi(i=1,2,n)； ,n) 3.1.2计算结果及分析： x1=18.5115 x2= 7.2791 x3= -12.0967 x4= -4.6715 x5=-5.9736 x6= 1
13、5.9624 系数矩阵为： 1 2 3 4 5 6442 3 4 5 7 8734 5 4 6 9 792a= b= 12 34 54 65 34 441224 45 89 45 23 433415 23 56 77 43 5645构成增广矩阵为 1 2 3 4 5 62 3 4 5 7 84 5 4 6 9 7A=12 34 54 65 34 4424 45 89 45 23 4315 23 56 77 43 56经过Gauss全主元消元法消元后矩阵a变为 447392 123445 89.00 45.0000 45.0000 23.0000 43.0000 24.000 0 48.6854
14、 -5.3146 28.5281 28.9438 -0.1011 0 0 10.8117 -2.0441 -4.5008 -2.4835a= 0 0 0 6.2806 4.1230 3.7133 0 0 0 0 1.9065 -1.4148 0 0 0 0 0 -0.2568消元后的b变为 34.000 23.6067 -26.9077b= 97.0342 4.2430 -4.7529最终将解得的x带入原方程组，计算出误差矩阵 error=AX-b 下面给出误差矩阵 0.01420.02840.0284error=10-12 0.11370.11370.1137由误差矩阵可知误差已接近于0，可
15、知：这种算法的具有数值稳定，精度高的优点，是比较好的算法，但是缺点是在选主元的过程中，消耗了较多的时间。 该实例的运行时间为0.048152 seconds. 3.2Gauss-Seidel迭代法 3.2.1Gauss-Seidel迭代法算法流程图： 1) 输入：aij(i,j=1,2,n),bi=(i=1,2,n),e; 2) 输入：初始解向量xi(i=1,2,3) 对i=1,2,n)； ,n做 na) xi=(bi-aijxj)/aii； j=1jib) ei=yi-xi； c) xi=yi； 4) fmaxeie then 输出：xi(i=1,2,iin,n)，停机else 返回第三步。
16、 3.2.2 参数的设置： 本题初始解：x(0)=(0,0,0,0,0,0)T; 误差精度e=0.0001; 3.2.3计算结果及分析： 针对这个线性方程组，用迭代法无法求解，因为这个系数矩阵的谱半径不小于1，在本题中谱半径为161.1222远大于1。这个问题的出现从侧面反映了迭代法的适应面会很小，因为现实生活中的很多矩阵都不是会符合谱半径小于1这个条件的，虽然这种方法迭代速度回很快，但是这个谱半径要求小于1弊端是其致命的弱点。 四、计算结果分析及遇到的问题 4.1计算结果的分析 根据误差矩阵来看，此方法解决此类问题具有很好的效果，误差数量级在10-12上面，误差是非常小的，几乎可以说是没有，
17、而且运行时间也很快，而且当矩阵变得很大时，运行时间并没有发生明显的变化，下满给出运行时间与矩阵阶数的变化关系的图表。 运行时间 矩阵阶数 0.042626 2 0.046684 4 0.047888 8 0.044901 16 0.046488 32 0.050995 64 0.084790 128 0.475078 256 6.029018 55.884309 512 1024 通过上表可知:当矩阵阶数太大时，此方法不是一个很好的方法，当矩阵的阶数小于256阶时，采用此方法比较好，因为上述数据是用服从正态分布的随机数的出来的，结果固然具有一定的波动性，但是这种趋势是正确的，我们画出其变化曲线
18、，结果会更加一目了然。 从这幅图像中，我们知道运行时间并不与阶数成正比，而是一种类似二次曲线的效果，所以阶数越大，此种方法越不可行。 4.2遇到的问题 在编这个算法的程序时有一个语句forj=n-1:1 这语句一直没有发现毛病，运行是没有报错，但是因为得出了错误的结果，便开始查程序。查程序时发现没有进入到循环里面去，突然意识到应该写成for犯了一个很低级的错误，害的我找了半天的错误。 j=n-1:-1:1，五、算法的改进 针对算法的改进，本人没有什么好的见解，一般的Gauss消元法，都要采用回代的方法解出x，在这里推荐使用Gauss-Jordan消元法，该算法是不需要回代过程的。Gauss-J
19、ordan消元法在形式上要比Gauss消元法简单一些，没有回代过程，消元完毕就获得解答，11但从工作量角度来看，它大约要O(n3)，比有回代的Gauss消元法O(n3)要231多出O(n3)。 6附录 %Gauss 全主元消元法% tic clc; clear all; a=1 2 3 4 5 6;2 3 4 5 7 8;4 5 4 6 9 7;12 34 54 65 34 44;24 45 89 45 23 43;15 23 56 77 43 56; b=44;73;92;12;34;45; A=a; B=b; n=size(a,1); %A=a b; z=1:n; for i=1:n-1
20、I J=find(abs(a)=max(max(abs(a(i:end,i:end); l=I(1);t=J(1); number=a(l,t); if number(1)=0 disp(method failed); break; end k=i; if i=l temp=a(l,:);temp1=b(l); a(l,:)=a(i,:);b(l)=b(i); a(i,:)=temp;b(i)=temp1; end if i=t temp=a(:,t); a(:,t)=a(:,i); a(:,i)=temp; temp1=z(i); z(i)=z(t); z(t)=temp1; end k=i
21、; ll=a(k+1:end,k)./a(k,k); a(k+1:end,:)=a(k+1:end,:)-ll*a(k,:); b(k+1:end)=b(k+1:end)-b(k).*ll; end if a(n,n)=0 disp(sigular); else b(n)=b(n)./a(n,n); for j=n-1:-1:1 b(j)=(b(j)-a(j,j+1:end)*b(j+1:end)./a(j,j); end end x(z)=b; disp(x); error=A*x-B; disp(error); toc % Gauss-seidel法% clc; clear all; a=
22、1 2 3 4 5 6;2 3 4 5 7 8;4 5 4 6 9 7;12 34 54 65 34 44;24 45 89 45 23 43;15 23 56 77 43 56; b=44;73;92;12;34;45; lamada=0.00001% x=0 0 0 0 0 0;% A=a; n=size(a,1); e=ones(n,1); U V=eig(A); ll=diag(V); G=max(abs(ll); if Glamada for i=1:n s=A(i,:)*x-A(i,i)*x(i); y=(b(i)-s)/A(i,i); e(i)=y-x(i); x(i)=y; end disp(s); end disp(x); else disp(method failed); end
