三种消元法程序——数学软件实验报告

北京大学数学软件实验任务书课程名称数学软件实验班级数学0901 学号0912020112 姓名****实验课题线性方程组高斯消去法，高斯列主元消去法，高斯全主元消去法实验目的熟悉线性代数方程组高斯消去法，高斯列主元消去法，高斯全主元消去法实验要求运用Matlab/C/C++/Java/Maple/Mathematica等其中一种语言完成实验内容线性方程组高斯消去法线性方程组高斯列主元消去法线性方程组高斯全主元消去法成绩100 教师****实验一 高斯消去法求解线性方程组【实验名称】高斯消去法求解线性方程组【实验目的】进一步熟悉理解Guass 消元法解法思路，提高matlab 编程能力。

【实验要求】已知线性方程矩阵，利用软件求解线性方程组的解。

【实验原理】消元过程：设0)0(11≠a ，令乘数)0(11)0(11/a a m i i -=，完成以下步骤： 首先，进行消元操作（消去第i 个方程组的i x ）1i m ×第1个方程+第i 个方程（i=2，3，.....n ），则第i 个方程变为1)1(2)1(2...i n in i b x a x a =++；这样消去第2，3，... ，n 个方程的变元i x 后。

原线性方程组变为：⎪⎪⎪⎩⎪⎪⎪⎨⎧=++=++=++)1()1(2)1(2)1(2)1(22)1(22)0(1)0(11)0(11... . .... ...n n nn n n n n n b x a x a b x a x a b x a x a这样就完成了第一步消元。

回代过程：在最后的一方程中解出n x ，得：)1()1(/--=n nnn n n a b x ；再将n x 的值代入倒数第二个方程，解出1-n x ，依次往上反推，即可求出方程组的解，其通项为：3,...,1-n 2,-n k /)()1(1)1()1(=-=-+=--∑k kk nk j j k kjk kk a x abx【程序设计】function maintest1 clc clear allnum=input('please input the order n=') A=zeros(num,num); for i=1:num for j=1:numA(i,j)=input(''); endend Aif det(A)~=0for i=1:numb(i)=input('');endb=b'num=length(b)for k=1:num-1for i=k+1:numif A(k,k)~=0l=A(i,k)/A(k,k);A(i,:)=A(i,:)-A(k,:).*l;b(i)=b(i)-b(k)*l;endendendAB%回代求xx(num)=b(num)/A(num,num);for i=num-1:-1:1sum=0;for j=i+1:numsum=sum+A(i,j)*x(j);endx(i)=(b(i)-sum)/A(i,i);endxend【实验结论】高斯消元能很好的求解线性方程组，和用克莱姆法则求解方程组该算法简单而且求解次数少。

但是高斯顺序消元是有缺陷的，它并没有考虑主对角元素相对其他同列元素绝对值很小，有时会导致用绝对值较小的数字做分母而出现较大的数字与较小数字进行相加减运算而被“吞掉”的现象，以至于产出某一未知元较小的误差却引起其他元更大的误差。

实验二列主元消去法求解线性方程组【实验名称】列主元消去法求解线性方程组【实验目的】进一步熟悉理解列主元消去法求解线性方程组解法思路，提高matlab编程能力。

【实验要求】已知线性方程矩阵，利用软件求解线性方程组的解。

【实验原理】列主元消去法基本思路:设有线性方程组x b A =，设A 是可逆矩阵。

列主元消去法的基本思想就是通过列主元的选取将初等行变换作用于方程组的增广矩阵[]A |b =A ，将其中的A 变换成一个上三角矩阵，然后求解这个三角形方程组。

将方程组用增广矩阵[]()(1)|ij n n b a ⨯+A =A =表示。

1 消元过程Step1:选主元，找{},1,,k i k k n ∈+ 使得,max k i k ik k i na a ≤≤=,1,2,,1k n =-Step2:如果,0k i k a =，则矩阵A 奇异，程序结束；否则执行step3；Step3:如果k i k ≠，则交换第k 行与第k i 行对应元素位置，k kj i j a a ↔， ,,1j k n =+ ； Step4:消元，对,,i k n = ，求解/,ik ik kk l a a =对1,,1j k n =++ ，求解.ij ij ik kj a a l a =-2 回代过程：Step1:若0,nn a =则矩阵奇异，程序结束；否则执行step2 Step2:,1/;n n n nn x a a +=对1,,2,1i n =- ，求解：,11/n i i n ij j ii j i x a a x a +=+⎛⎫=- ⎪⎝⎭∑【程序设计】function [x,XA]=GaussLiezhuYuan(A,b) N=size(A) n=N(1); index = 0; for i=1:(n-1)me=max(abs(A(1:n,i))); %确定每一列绝对值最大的数值 for k=i:nif(abs(A(k,i))==me) index=k; break; end endtemp=A(i,1:n);A(i,1:n)=A(index,1:n);A(index,1:n)=temp; %与主元所在的行进行交换 bb=b(index); b(index)=b(i);b(i)=bb;for j=(i+1):nif(A(i,i)==0)disp('对角元素为0！');return;endl=A(j,i);m=A(i,i);A(j,1:n)=A(j,1:n)-l*A(i,1:n)/m;b(j)=b(j)-l*b(i)/m;endendx=SolveUpTriangle(A,b);XA=A;end【实验结论】将向量求解得到的解向量代入原方程，可以发现符合的还是比较好的。

这说明了列主元消去法求解这一类方程的有效性。

事实上，对于普通Gauss消去法的修正，还可以有行主元消去法，全主元消去法。

但是全主元消去法相对于行、列主元消去法的工作量要大的多。

所以列主元消去法是解线性方程组实用的方法之一。

实验三全主元消去法求解线性方程组【实验名称】全主元消去法求解线性方程组【实验目的】进一步熟悉理解全主元消去法求解线性方程组解法思路，提高matlab编程能力。

【实验要求】已知线性方程矩阵，利用软件求解线性方程组的解。

【实验原理】主元就是在每一步里面都要在被处理的矩阵里选取绝对值最大的元素作为主元，从而使高斯消去法具有更好的数据稳定性的方法。

虽然全主元消去法的求解结果更加可靠，但由于全主元每步耗时多，而且进行列变换，就会打乱未知变元的顺序，因此，比较繁琐。

（参考程序4）【程序设计】%利用高斯全主元消去法求AX=b方程组的解A=[3 2 0 0;1 4 2 0;0 1 2 1;0 0 -1 1]b=[5 7 2 -2]'N=size(A);n = N(1);index_l = 0;index_r = 0;order = 1:n; %记录未知数顺序的向量for i=1:(n-1)me = max(max(abs(A(i:n,i:n)))) %选取全主元for k=i:nfor r=i:nif(abs(A(k,r))==me)index_l = k;index_r = r; %保存主元所在的行和列k=n;break;endendendtemp = A(i,1:n);A(i,1:n) = A(index_l,1:n);A(index_l,1:n) = temp;bb = b(index_l);b(index_l)=b(i);b(i) = bb; %交换主行temp = A(1:n,i);A(1:n,i) = A(1:n,index_r);A(1:n,index_r) = temp; %交换主列pos = order(i);order(i) = order(index_r);order(index_r) = pos; %主列的交换会造成未知数顺序的变化for j=(i+1):nif(A(i,i)==0)disp('对角元素为0！');return;endl = A(j,i);m = A(i,i);A(j,1:n)=A(j,1:n)-l*A(i,1:n)/mb(j)=b(j)-l*b(i)/mendendx=SolveUpTriangle(A,b); %调用上三角系数矩阵求线性方程组的函数y=zeros(n,1);for i=1:nfor j=1:nif(order(j)==i)y(i)=x(j);endendend %恢复未知数原来的顺序x=y;XA = A;x【实验结论】虽然全主元消去法的求解结果更加可靠，高斯消去具有更好的数值稳定性。

但是由于全选主元每步消耗的时间更多，而且要进行列交换，那么所求未知量x1,x2,……xn的顺序就会被打乱，因此，实际应用中一般使用列主元消去法。