数值分析实验报告

实验五解线性方程组的直接方法实验5.1（主元的选取与算法的稳定性）问题提出：Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？Gauss消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。实验内容：考虑线性方程组nnnRbRAbAx,,编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss消去过程。实验要求：（1）取矩阵1415157,6816816816bA，则方程有解Tx)1,,1,1(*。取n=10计算矩阵的条件数。让程序自动选取主元，结果如何？（2）现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。思考题一：（Vadermonde矩阵）设niinniiniiniinnnnnnnxxxxbxxxxxxxxxxxxA002010022222121102001111，，其中，nkkxk,,1,0,1.01，（1）对n=2,5,8，计算A的条件数；随n增大，矩阵性态如何变化？（2）对n=5，解方程组Ax=b；设A的最后一个元素有扰动10-4，再求解Ax=b（3）计算（2）扰动相对误差与解的相对偏差，分析它们与条件数的关系。（4）你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因吗？相关MATLAB函数提示：zeros(m,n)生成m行，n列的零矩阵ones(m,n)生成m行，n列的元素全为1的矩阵eye(n)生成n阶单位矩阵rand(m,n)生成m行,n列(0,1)上均匀分布的随机矩阵diag(x)返回由向量x的元素构成的对角矩阵tril(A)提取矩阵A的下三角部分生成下三角矩阵triu(A)提取矩阵A的上三角部分生成上三角矩阵rank(A)返回矩阵A的秩det(A)返回方阵A的行列式inv(A)返回可逆方阵A的逆矩阵[V，D]=eig(A)返回方阵A的特征值和特征向量norm(A,p)矩阵或向量的p范数cond(A,p)矩阵的条件数[L，U，P]=lu(A)选列主元LU分解R=chol(X)平方根分解Hi=hilb(n)生成n阶Hilbert矩阵5.1实验过程：5.1.1程序：functionx=gauss(n,r)n=input('请输入矩阵A的阶数:n=')A=diag(6*ones(1,n))+diag(ones(1,n-1),1)+diag(8*ones(1,n-1),-1)b=A*ones(n,1)fori=1:4p=input('条件数对应的范数是p-范数：p=')pp=cond(A,p)endpause[m,n]=size(A);nb=n+1;Ab=[Ab]r=input('请输入是否为手动，手动输入1，自动输入0：r=')fori=1:n-1ifr==0[pivot,p]=max(abs(Ab(i:n,i)));ip=p+i-1;ifip~=iAb([iip],:)=Ab([ipi],:);disp(Ab);pauseendendifr==1i=iip=input('输入i列所选元素所处的行数：ip=');Ab([iip],:)=Ab([ipi],:);disp(Ab);pauseendpivot=Ab(i,i);fork=i+1:nAb(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb);enddisp(Ab);pauseendx=zeros(n,1);x(n)=Ab(n,nb)/Ab(n,n);fori=n-1:-1:1x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);end5.1.2实验结果如下：1.按照实验要求一：取矩阵A的阶数:n=10且自动选取主元，程序结果运行如下：（2）现选择程序中手动选取主元的功能，观察并记录计算结果。①选取绝对值最大的元素为主元：程序运行开始如第一问的截图也是求范数故这里不在给出。Theansweris：1111111111②选取绝对值最小的元素为主元：Theansweris：1.0e+003*（INF0.0070.0057-0.0410-0.03030.34300.2577-2.7290-2.04632.7308）⑶取矩阵A的阶数:n=20，手动选取主元：选取绝对值最大的元素为主元：Theansweris：11111111111111111111选取绝对值最小的元素为主元：Theansweris：1.0e+007*（-Inf0.00000.0000-0.0000-0.00000.00000.0000-0.0003-0.00020.00220.0016-0.0175-0.01310.13980.1049-1.1185-0.83898.94786.7109-8.9478）⑷修改程序如下：functionx=gaussong(n,r)n=input('请输入矩阵A的阶数:n=')A=hilb(n)b=A*ones(n,1)fori=1:4p=input('条件数对应的范数是p-范数：p=')pp=cond(A,p)endpause[m,n]=size(A);nb=n+1;Ab=[Ab]r=input('请输入是否为手动，手动输入1，自动输入0：r=')fori=1:n-1ifr==0[pivot,p]=max(abs(Ab(i:n,i)));ip=p+i-1;ifip~=iAb([iip],:)=Ab([ipi],:);disp(Ab);pauseendendifr==1i=iip=input('输入i列所选元素所处的行数：ip=');Ab([iip],:)=Ab([ipi],:);disp(Ab);pauseendpivot=Ab(i,i);fork=i+1:nAb(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb);enddisp(Ab);pauseendx=zeros(n,1);x(n)=Ab(n,nb)/Ab(n,n);fori=n-1:-1:1x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);end①所求范数为：自动输入结果为：ans=1.00001.00001.00001.00001.00020.99961.00070.99931.00040.9999②选取绝对值最大的元素为主元结果为：Theansweris：NaNNaNNaNNaNNaNInf-Inf-Inf281.3945-283.3708③选取绝对值最小的元素为主元结果为：Theansweris：NaNNaNNaN-Inf-5.8976-1.9243-2.0291-4.997223.4548-11.10125.1.3对实验结果进行分析：5.1.3.1对实验要求一的结果进行分析：对于Gauss消去法就是用行的初等变换将原线性方程组系数矩阵转化为简单形式，从而进行求解，缺点是迭代次数可能较多，效率不高，且在消去过程中不可以将主元素很小的做除数，否则将导致其他元素数量级的严重增长和舍入误差的扩散，使得计算解不可靠。5.1.3.2对实验要求二的结果进行分析：通过每次选取最大或最小的主元可以发现取绝对值大的元素作为主元比取绝对值小的元素作为主元时产生的结果比较准确，即选取绝对值小的主元时结果产生了较大的误差，条件数越大产生的误差就越大。所以应尽量避免很小的数作为除数。5.1.3.3对实验要求三的结果进行分析：此要求是对要求一和要求二的一个延续，通过实验结果可以看出若采用很小的数作为主元迭代次数越多导致的结果越不可靠，甚至出现错误。5.1.3.4对实验要求四的结果进行分析：对新矩阵进行实验发现依然符合上述规律，可以知道，在进行迭代时主元的选择与算法的稳定性有密切的联系选取绝对值大的元素作为主元比绝对值小的元素作为主元时对结果产生的误差较小。条件数越大对用gauss消去法解线性方程组时，对结果产生的误差就越大。5.1.4实验总结：1.在用gauss消去法解线性方程组时，主元的选取与算法的稳定性有密切的联系，选取适当的主元有利于得出稳定的算法，2.在算法的过程中，选取绝对值较大的主元比选取绝对值较小的主元更有利于算法的稳定，选取绝对值最大的元素作为主元时，得出的结果相对较准确较稳定。3.条件数越小，对用这种方法得出的结果更准确。4.在算除法的过程中要尽量避免使用较小的数做为除数，以免发生结果数量级加大，使大数吃掉小数，产生舍入误差。实验5.2（线性代数方程组的性态与条件数的估计）问题提出：理论上，线性代数方程组bAx的摄动满足bbAAAAAcondxx11)(矩阵的条件数确实是对矩阵病态性的刻画，但在实际应用中直接计算它显然不现实，因为计算1A通常要比求解方程bAx还困难。实验内容：Matlab中提供有函数“condest”可以用来估计矩阵的条件数，它给出的是按1-范数的条件数。首先构造非奇异矩阵A和右端，使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动bA和，使得bA和充分小。实验要求：（1）假设方程Ax=b的解为x，求解方程bbxAAˆ)(，以1-范数，给出xxxxxˆ的计算结果。（2）选择一系列维数递增的矩阵（可以是随机生成的），比较函数“condest”所需机器时间的差别.考虑若干逆是已知的矩阵，借助函数“eig”很容易给出cond2(A)的数值。将它与函数“cond(A,2)”所得到的结果进行比较。（3）利用“condest”给出矩阵A条件数的估计，针对（1）中的结果给出xx的理论估计，并将它与（1）给出的计算结果进行比较，分析所得结果。注意，如果给出了cond(A)和A的估计，马上就可以给出1A的估计。（4）估计著名的Hilbert矩阵的条件数。njijihhHjinnji,,2,1,,11,)(,,5.2实验过程如下：5.2.1.1实验要求一的程序如下：functionn=jisuan(n)a=fix(100*rand(n))+1x=ones(n,1)b=a*xdata=rand(n)*0.00001datb=rand(n,1)*0.00001A=a+dataB=b+datbx0=get(A,B)x1=norm(x0-x,1)/norm(x,1)functionx=get(A,B)[m,n]=size(A);nb=n+1;AB=[AB];fori=1:n-1pivot=AB(i,i);fork=i+1:nAB(k,i:nb)=AB(k,i:nb)-(AB(k,i)/pivot)*AB(i,i:nb);endendx=zeros(n,1);x(n)=AB(n,nb)/AB(n,n);fori=n-1:-1:1x(i)=(AB(i,nb)-AB(i,i+1:n)*x(i+1:n))/AB(i,i);End5.2.1.2实验要求一程序运行结果如下：系数矩阵a为：7018142952536398247346580289074421962620992338538830595879897467437769加扰动后的系数矩阵A为：70.000004618.00