高斯消元法_实验报告

华中科技大学数值分析实验报告系、年级研究生院2012级学号姓名类别硕士2013年5月6日华中科技大学数值分析实验报告（实验6.1）2实验6.1实验要求：根据教材实验6.1做出相应改编：分别使用Gauss消元、列选主元。全选主元的方法求解线性方程组，分别比较三种消元方法的结果和算法的区别，并说明主元的选取在Gauss消元的中的作用。问题提出：Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？Gauss消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。一般来说书本上采用的列选主元的办法对其线性方程组进行求解的，那么我们是否可以选择一种行列都选取主元消去的办法来减小相应的误差呢？全主元消元法和列主元消元法一样都是由高斯消元法演变而来。只不过选取主元的范围有所加大。全选主元相对于列选主元的更加复杂化了，因为在运算的过程中导致了元的位置发生了变化，这样我们就不得不追踪每个元的位置。本次实验就几个问题进行了matlab实验分析，比较几种计算方法的优劣性。实验内容：考虑线性方程组nnnRbRAbAx,,编制一个程序：分别能进行Gauss消去、列选主元Gauss消去、全选主元Gauss消去法进行解线性方程组。对三种算法所得到的结果进行比较，分析三种计算方法的准确性。华中科技大学数值分析实验报告（实验6.1）3具体内容：（1）取矩阵1415157,6816816816bA，则方程有解Tx)1,,1,1(*。取n=10、n=20计算矩阵的条件数。分别编写利用matlab编写运算程序，实现Gauss消去、列选主元消去以及全选主元消去的方法。比较三种计算方法的运算结果。在列选主元的过程中分别采用每步消去过程总选取按模最小或按模尽可能小的元素作为主元或每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。1采用普通Gauss消元法进行计算Gauss消去法的基本思想是，通过将一个方程乘或除某个数以及两个方程相加减这两种运算手续，逐步减少方程组中变元的数目，最终使某个方程只含有一个变元，从而得出所求的解。对于𝐴𝑥=𝑏，Gauss消去法的求解思路为：（1）若𝑎11(1)≠0，先让第一个方程组保持不变，利用它消去其余方程组中的𝑥1,使之变成一个关于变元𝑥2,𝑥3……𝑥𝑛的n-1阶方程组。（2）按照（1）中的思路继续运算得到更为低阶的方程组。（3）经过n-1步的消元后，得到一个三角方程。（4）利用求解公式回代得到线性方程组的解。根据这个思路编写matlab程序如下：functionx=gauss(n)disp('请输入构造的矩阵的阶数（10/20）')n=input('');disp('构造矩阵为A=')A=diag(6*ones(1,n))+diag(ones(1,n-1),1)+diag(8*ones(1,n-1),-1)disp('由A构造出矩阵b=')b=A*ones(n,1)[m,n]=size(A);disp('增广矩阵为:')Ab=[Ab]fori=1:n-1yuan=Ab(i,i);华中科技大学数值分析实验报告（实验6.1）4fork=i+1:nAb(k,i:(n+1))=Ab(k,i:(n+1))-(Ab(k,i)/yuan)*Ab(i,i:(n+1));enddisp(Ab);pauseendx=zeros(n,1);x(n)=Ab(n,(n+1))/Ab(n,n);fori=n-1:-1:1x(i)=(Ab(i,(n+1))-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);end取n=10运行的结果为：Ab=61000000007861000000015086100000015008610000015000861000015000086100015000008610015000000861015000000086115000000008614消去结束后的矩阵为：解得线性方程组的解为：ans=1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000取n=20运行的结果为：Ab=华中科技大学数值分析实验报告（实验6.1）5消去结束后的矩阵为：解得线性方程组的解为：ans=1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00002采用列选主元消去法进行计算在采用Gauss消去法的时候，主元绝对值的大小将影响到计算结果，主元的绝对值越大，算法的稳定性越好。列选主元消去法matlab程序的计算思路为：华中科技大学数值分析实验报告（实验6.1）6（1）先构造需要计算的矩阵，得到增广矩阵Ab。（2）将系数矩阵A的每一列的绝对值最大的元素换至对角线上，矩阵b中的元素也随之改变，先判断主元是否为0。然后利用此主元逐行消去此主元所在的列中的元素，矩阵b中的元素也随之改变。（3）经过n-1步运算过后，矩阵A就变换成为一个三角矩阵。。（4）逐次回代，就能计算出方程组的解。分别每步消去过程总选取按模最大的元素作为主元和每步消去过程总选取按模最小的元素作为主元选取列选主元消去法的matlab程序为：functionx=gauss(n)disp('请输入构造的矩阵的阶数（10/20）')n=input('');disp('构造矩阵为A=')A=diag(6*ones(1,n))+diag(ones(1,n-1),1)+diag(8*ones(1,n-1),-1)disp('由A构造出矩阵b=')b=A*ones(n,1)[m,n]=size(A);disp('增广矩阵为:')Ab=[Ab]fori=1:n-1ifway==1i=idisp('输入每一列的主元所在的行')hang=input('');Ab([ihang],:)=Ab([hangi],:);disp(Ab);pauseendzhuyuan=Ab(i,i);fork=i+1:nAb(k,i:(n+1))=Ab(k,i:(n+1))-(Ab(k,i)/zhuyuan)*Ab(i,i:(n+1));enddisp(Ab);pauseendx=zeros(n,1);x(n)=Ab(n,(n+1))/Ab(n,n);fori=n-1:-1:1x(i)=(Ab(i,(n+1))-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);end华中科技大学数值分析实验报告（实验6.1）7（i）：每步消去过程总选取按模最大的元素作为主元:取n=10手动执行该程序：消去结束后的矩阵为：解得线性方程组的解为：ans=1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000取n=20运行的结果为：Ab=解得线性方程组的解为：华中科技大学数值分析实验报告（实验6.1）8ans=1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000（ii）每步消去过程总选取按模最小的元素作为主元:取n=10手动执行该程序：消去结束后的矩阵为：解得线性方程组的解为：ans=1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000取n=20运行的结果为：消去结束后的矩阵为：华中科技大学数值分析实验报告（实验6.1）9解得线性方程组的解为：ans=1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000由上述两种选取主元的方法比较可知，采用每一步消去过程总选取按模最大或按模尽可能大的元素作为主元的选取方法最终得到的方程组的解更加精确，这是因为在计算过程中人为的避免了大数除小数而造成的误差的扩大，而每一步消去过程总选取按模最小或按模尽可能小的元素作为主元的选取方法最终得到的解会有误差，而误差的产生就是由于上诉原因而造成了。所以在实际计算过程中，必须考察计算模型的可行性。3采用全选主元消去进行计算全选主元和列选主元消去法的区别在于全选主元不仅要用到列选主元，而且需要行选绝对值最大的元，这样的麻烦之处就在于不仅需要进行线性方程顺序的变换，而且在解方程的过程中未知数的顺序也发生了相应的改变，这就给程序的编写和运算的时间的上带来了很大的麻烦。全选主元消去法的基本思路为：华中科技大学数值分析实验报告（实验6.1）10（1）先构造需要计算的矩阵，得到增广矩阵Ab（2）将每次更新后系数矩阵的每一列的绝对值最大的元素换至对角线上（绝对值最大的元素必须在对角线的下方），矩阵b中的元素也随之改变。（3）然后将更新后系数矩阵的每行中的绝对值最大元素换至对角线上（该绝对值最大元素必须要在主元的右侧）。先判断主元是否为0，然后利用此主元逐行消去此主元所在的列中的元素，矩阵b中的元素也随之改变。（4）经过n-1步运算过后，矩阵A就变换成为一个三角矩阵。。（5）逐次回代，就能计算出方程组的解。采用全选主元Gauss消去法编写的matlab程序为：functionx=gauss(n)disp('请输入构造的矩阵的阶数（10/20）')n=input('');disp('构造矩阵为A=')A=diag(6*ones(1,n))+diag(ones(1,n-1),1)+diag(8*ones(1,n-1),-1)disp('由A构造出矩阵b=')b=A*ones(n,1)[m,n]=size(A);disp('增广矩阵为:')Ab=[Ab]fori=1:n-1zhuyuan=max(max(abs(A(i:n,i:n))));forr=i:nfort=i:nifzhuyuan==abs(A(r,t))zhuhang=r;zhulie=t;endendendph=A(i,:);A(i,:)=A(zhuhang,:);A(zhuhang,:)=ph;b=b';bb=b(i);b(i)=b(zhuhang);b(zhuhang)=bb;b=b';pl=A(:,i);A(:,i)=A(:,zhulie);A(:,zhulie)=pl;华中科技大学数值分析实验报告（实验6.1）11Ab=[A,b];zhuyuan=Ab(i,i);fork=i+1:nAb(k,i:(n+1))=Ab(k,i:(n+1))-(Ab(k,i)/zhuyuan)*Ab(i,i:(n+1));endA=Ab(:,[1:n]);b=Ab(:,[(n+1)]);disp(Ab);pauseendx=zeros(n,1);x(n)=Ab(n,(n+1))/Ab(n,n);fori=n-1:-1:1x(i)=(Ab(i,(n+1))-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);end取n=10手动执行该程序：消去结束后的矩阵为：解得线性方程组的解为：ans=1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000取n=20运行的结果为：Ab=华中科技大学数值分析实验报告（实验6.1）12解得线性方程组的解为：ans