数值分析A实验报告

1数值分析A实验报告清华大学2012年秋季学期2第1章实验3.1（主元的选取与算法的稳定性）1.1问题提出Gauss消去法是我们在线性代数中已经熟悉的，但是由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？Gauss消去法从理论算法到数值算法，其关键是主元的选取。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。1.2实验内容考虑线性方程组𝐀𝒙=𝐛,𝐀∈𝑹𝒏x𝒏,𝒃∈𝑹𝒏编制一个能自动选取主元，又能手动选取主元的求解线性代数方程的Gauss消去过程。1.3实验要求1.3.1取矩阵A=(61861⋱⋱8⋱6816),b=(715⋮1514),则方程有解𝑥∗=(11⋮11)。取n=10计算矩阵的条件数，让程序自动选取主元，结果如何？1.3.2现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验结果。31.3.3取实验阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。1.3.4选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数，重复上述实验，观察记录并分析实验的结果。1.4实验程序1.4.1Matrix31函数函数的主要功能为生成一个与实验3.1要求相符的n维矩阵，输入者仅需要输入matrix31（n）即可获得n维的与实验要求相符的矩阵，其中n为任意自然数。function[A,b]=matrix31(n)A=zeros(n);fori=1:nA(i,i)=6;endfori=1:n-1A(i,i+1)=1;endfori=1:n-1A(i+1,i)=8;endend1.4.2Array31函数函数的主要功能为生成一个与实验3.1要求相符的n维向量，输入者仅需要输入array31（n）即可获得n维的与实验要求相符的向量，其中n为任意自然数。function[B]=array31(n)B=zeros(1,n);B(1,1)=7;B(1,n)=14;4fori=2:n-1B(1,i)=15;endend1.4.3Gaussauto函数函数的主要功能为调用matrix31(n)及array31(n)函数，并使用高斯列主元消去法对方程进行求解。输入者仅需要输入gaussauto(A,b)，便能够获得方程的解。A及b分别为由函数matrix31及array31所生成的方阵与向量。function[x,det,flag]=gaussauto(A,b)[n,m]=size(A);nb=length(b);ifn~=merror('TherowsandcolumnsofmatrixAmustbeequal!');endifm~=nberror('ThecolumnsofAmustbeequalthelengthofb!');endflag='OK';det=1;x=zeros(n,1);fork=1:n-1max1=0;fori=k:nifabs(A(i,k))max1max1=abs(A(i,k));r=i;endendifmax11e-10flag='failure';returnendifrkforj=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det;endfori=k+1:nm=A(i,k)/A(k,k);forj=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k);enddet=det*A(n,n);5ifabs(A(n,n))1e-10flag='failure';return;endfork=n:-1:1forj=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);enddisp('结果X为：')x(k)=b(k)/A(k,k);end1.4.4Gaussmanual函数函数的主要功能为调用matrix31(n)及array31(n)函数，并使用高斯消去法对方程进行求解。输入者仅需要输入gaussauto(A,b)，便能够获得方程的解。在每一次消元过程发生前，系统将会提示输入者选择想要作为主元素的行数，计算机将会按照输入者所输入的行数作为主元素进行结果计算。其中A及b分别为由函数matrix31及array31所生成的方阵与向量。functiongaussmanual(A,b)functionchange(A,i,j)[k,l]=size(A);form=1:lp=A(i,m);A(i,m)=A(j,m);A(j,m)=p;endendfunctionxiaoyuan(A,i,j)[k,l]=size(A);m=A(j,i)/A(i,i);forn=i:lA(j,n)=A(j,n)-m*A(i,n);endend[n,m]=size(A);B=b';C=[AB];X=zeros(n+1,1);M=zeros(n-1,1);Cfori=1:n-1C;disp('这是您第');disp(i);disp('次选取主元，选取主元的行数不能小于');disp(i);j=input('请您输入想要选取的主元所在的行数。');6ifji-1error('请您谨慎选取主元，您所选取的主元所在的行数不能小于选择的次数。')returnendfork=1:m+1p=C(i,k);C(i,k)=C(j,k);C(j,k)=p;endC;ifC(i,i)==0error('请您选择主元不为0的元素，谢谢。')endfora=i:n-1M=C(a+1,i)/C(i,i);foru=i:m+1C(a+1,u)=C(a+1,u)-M*C(i,u);endendCendX(n,1)=C(n,n+1)/C(n,n);fori=n-1:-1:1XX=C(i,:)*X;X(i,1)=(C(i,n+1)-XX)/C(i,i)endend1.4.5Gauss函数Gauss为实验3.1题的主函数，其将matrix31,array31,gaussauto及gaussmanual函数都置放在其下。输入者只需要调用gauss函数便能够对实验3.1题进行全方位求解。Gauss函数调用了matrix31及array31函数来生成了一个方阵及向量，并设置了一个选择环节供实验者进行选择，即自动选择主元或者手动选择主元。functiongaussdisp('此为数值分析A实验3.1专用程序。请您依照步骤依次进行操作。')n=input('请您输入方阵与向量的维数');A=matrix31(n);b=array31(n);option=input('若您想要计算机自动使用列主元消去法请您输入1，若您想要手动选取主元请您输入2。');ifoption==1gaussauto(A,b)endifoption==2gaussmanual(A,b)7endend1.5实验报告1.5.1取矩阵A=(61861⋱⋱8⋱6816),b=(715⋮1514),则方程有解𝑥∗=(11⋮11)。取n=10计算矩阵的条件数，让程序自动选取主元，结果如何？当n=10时，矩阵的2-条件数为1728。我们调用gauss函数，并选择让系统选取最大元素作为主元素，得到的结果为(1,1,1,1,1,1,1,1,1,1,1)𝑇，即满足合理解。1.5.2现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验结果。在这里我们使用∞−范数，并选取最小的∞−范数作为主元。我们调用gauss函数，并选取手动选择主元的选项进行操作。当n=10时，有如下过程：手动选取模最小但是不为0的元素作为主元，可以得到结果依然为(1,1,1,1,1,1,1,1,1,1,1)𝑇。而按照取最大模作为主元的话，依然可以得到结果为(1,1,1,1,1,1,1,1,1,1,1)𝑇。结果分析，由于最大的主元与最小的主元的绝对值没有过于接近于0，都是“比较好的一个数”，然而，如果选择0能够作为模最小的数的话，则计算过程将会崩溃。由此可见，虽然矩阵A的条件数非常大即显得相当“病态”，然而，假设不使用0元素作为主元素的话，不管以最大模或者最小模作为主元8素的话，对于结果都是没有影响的，原因在于不管是哪里一种方法，选取的主元的绝对值与其他元素相较起来都属于“良态”，当然这个是不考虑0元素的存在的情况。在实验要求4的时候，将会举出一个例子来说明如果取非常接近0的一个数作为主元素的话，将会出现问题，而在该例中，其病态就原形毕露了。1.5.3取实验阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。在这里我们使用∞−范数，并选取最小的∞−范数作为主元。我们调用gauss函数，并选取手动选择主元的选项进行操作。当n=20时，有如下过程：手动选取模最小但是不为0的元素作为主元，可以得到结果依然为(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)𝑇。而按照取最大模作为主元的话，依然可以得到结果为(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)𝑇结果分析，由于最大的主元与最小的主元的绝对值没有过于接近于0，都是“比较好的一个数”，然而，如果选择0能够作为模最小的数的话，则计算过程将会崩溃。由此可见，虽然矩阵A的条件数非常大即显得相当“病态”，然而，假设不使用0元素作为主元素的话，不管以最大模或者最小模作为主元素的话，对于结果都是没有影响的，原因在于不管是哪里一种方法，选取的主元的绝对值与其他元素相较起来都属于“良态”，当然这个是不考虑0元素的存在的情况。1.5.4选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数，重复上述实验，观察记录并分析实验的结果。9我们这里直接引用实验2即将出现的Hilbert矩阵作为系数矩阵，并给定其有精确解[1,…,1]’。我们考察六维度的Hilbert系数矩阵，有如下表达，即H为1.00000.50000.33330.25000.20000.16670.50000.33330.25000.20000.16670.14290.33330.25000.20000.16670.14290.12500.25000.20000.16670.14290.12500.11110.20000.16670.14290.12500.11110.10000.16670.14290.12500.11110.10000.0909设有HX=b其中给定X=[111111]’可以计算出b=[2.451.59291.21790.99560.84560.7365]’现在我们已经知道b和H了，便可以使用高斯方法进行求解了！在求解之前，我们可以计算其条件数为1495100，这以之前所计算得到的条件数相比，简直为小巫见大巫了，也足以显示其病态性质。当我们用高斯方法进行求解时，可以解得精确解X=[111111]’。当然，对于维度小于6的Hilbert方程来说，其都是成立的，而且容易获得精确解。现在我们把维度增加到21时，依然使用高斯方法进行求解。我们可以发现21维度的Hilbert矩阵的条件数为4.4080x10^18，可以说是非常病态的方程。对于用高斯法解这样的方程，有如下解：10X=[10.99971.00840.865862.1174-4.269315.159-18.286.252413.83110.008-49.83845.871-27.154.555-