实验五 解线性方程组的直接方法实验5.1 （主元的选取与算法的稳定性） 问题提出：Gauss 消去法是我们在线性代数中已经熟悉的。

但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss 消去法作为数值算法的稳定性呢？Gauss 消去法从理论算法到数值算法，其关键是主元的选择。

主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。

实验内容：考虑线性方程组编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss 消去过程。

实验要求：（1）取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816 b A ，则方程有解T x )1,,1,1(* =。

取n=10计算矩阵的条件数。

让程序自动选取主元，结果如何？（2）现选择程序中手动选取主元的功能。

每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。

若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。

（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。

（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。

重复上述实验，观察记录并分析实验结果。

思考题一：（Vadermonde 矩阵）设⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=∑∑∑∑====n i i n n i i ni i n i i n n n n nn nx x x x b x x x x x x x x x x x x A 002010022222121102001111 ，，其中，n k k x k ,,1,0,1.01 =+=，（1）对n=2,5,8，计算A 的条件数；随n 增大，矩阵性态如何变化？（2）对n=5，解方程组Ax=b ；设A 的最后一个元素有扰动10-4，再求解Ax=b （3）计算（2）扰动相对误差与解的相对偏差，分析它们与条件数的关系。

（4）你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因吗？5.1实验过程：5.1.1程序：function x=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)for i=1:4p=input('条件数对应的范数是p-范数：p=')pp=cond(A,p)endpause[m,n]=size(A);nb=n+1;Ab=[A b]r=input('请输入是否为手动，手动输入1，自动输入0：r=')for i=1:n-1if r==0[pivot,p]=max(abs(Ab(i:n,i)));ip=p+i-1;if ip~=iAb([i ip],:)=Ab([ip i],:);disp(Ab); pauseendendif r==1i=iip=input('输入i列所选元素所处的行数：ip=');Ab([i ip],:)=Ab([ip i],:);disp(Ab); pauseendpivot=Ab(i,i);for k=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);for i=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）现选择程序中手动选取主元的功能，观察并记录计算结果。

①选取绝对值最大的元素为主元：程序运行开始如第一问的截图也是求范数故这里不在给出。

The answer is ：1 1 1 1 1 1 1 1 1 1②选取绝对值最小的元素为主元：The answer is：1.0e+003*（INF 0.007 0.0057 -0.0410 -0.0303 0.3430 0.2577 -2.7290 -2.0463 2.7308）⑶取矩阵A的阶数:n=20，手动选取主元：①选取绝对值最大的元素为主元：The answer is ：1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1②选取绝对值最小的元素为主元：The answer is：1.0e+007*（-Inf 0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0003 -0.0002 0.00220.0016 -0.0175 -0.0131 0.1398 0.1049 -1.1185 -0.8389 8.9478 6.7109 -8.9478）⑷修改程序如下：function x=gaussong(n,r)n=input('请输入矩阵A的阶数:n=')A=hilb(n)b=A*ones(n,1)for i=1:4p=input('条件数对应的范数是p-范数：p=')pp=cond(A,p)endpause[m,n]=size(A);nb=n+1;Ab=[A b]r=input('请输入是否为手动，手动输入1，自动输入0：r=')for i=1:n-1if r==0[pivot,p]=max(abs(Ab(i:n,i)));ip=p+i-1;if ip~=iAb([i ip],:)=Ab([ip i],:);disp(Ab); pauseendendif r==1i=iip=input('输入i列所选元素所处的行数：ip=');Ab([i ip],:)=Ab([ip i],:);disp(Ab); pauseendpivot=Ab(i,i);for k=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);for i=n-1:-1:1x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);end①所求范数为：自动输入结果为：ans =1.0000 1.00001 .0000 1.0000 1.0002 0.9996 1.0007 0.9993 1.0004 0.9999②选取绝对值最大的元素为主元结果为：The answer is ：NaN NaN NaN NaN NaN Inf -Inf -Inf 281.3945 -283.3708 ③选取绝对值最小的元素为主元结果为： The answer is ：NaN NaN NaN -Inf -5.8976 -1.9243 -2.0291 -4.9972 23.4548 -11.1012 5.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（线性代数方程组的性态与条件数的估计）问题提出：理论上，线性代数方程组b Ax =的摄动满足矩阵的条件数确实是对矩阵病态性的刻画，但在实际应用中直接计算它显然不现实，因为计算1-A 通常要比求解方程b Ax =还困难。

实验内容：Matlab 中提供有函数“condest ”可以用来估计矩阵的条件数，它给出的是按1-范数的条件数。

首先构造非奇异矩阵A 和右端，使得方程是可以精确求解的。

再人为地引进系数矩阵和右端的摄动b A ∆∆和，使得b A ∆∆和充分小。

实验要求：（1）假设方程Ax=b 的解为x ，求解方程b b xA A ∆+=∆+ˆ)(，以1-范数，给出xx x xx -=∆ˆ的计算结果。

（2）选择一系列维数递增的矩阵（可以是随机生成的），比较函数“condest ”所需机器时间的差别.考虑若干逆是已知的矩阵，借助函数“eig ”很容易给出cond 2(A)的数值。

将它与函数“cond(A,2)”所得到的结果进行比较。

（3）利用“condest ”给出矩阵A 条件数的估计，针对（1）中的结果给出xx ∆的理论估计，并将它与（1）给出的计算结果进行比较，分析所得结果。

注意，如果给出了cond(A)和A的估计，马上就可以给出1 A的估计。

（4）估计著名的Hilbert矩阵的条件数。

5.2 实验过程如下：5.2.1.1 实验要求一的程序如下：function n=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)function x=get(A,B)[m,n]=size(A);nb=n+1;AB=[A B];for i=1:n-1pivot=AB(i,i);for k=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);for i=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为：b加扰动后的b值为：xx的值为：x0的值为：1.146958549E-06 x1的值为：6.8990e-007 5.2.1.3实验结果为：xx x xx -=∆ˆ的计算结果为：6.8990e-0075.2.2.1 实验要求二的程序如下：function cond2(A) B=A'*A;[V1,D1]=eig(B); [V2,D2]=eig(B^(-1));cond2A=sqrt(max(max(D1)))*sqrt(max(max(D2))) endfor n=10:10:100n=n A=fix(100*randn(n)); condestA=condest(A) cond2(A) condA2=cond(A,2) pause endfunction bijiao(n)a=fix(100*rand(n))+1; x=ones(n,1); b=a*x; data=rand(n)*0.00001; datb=rand(n,1)*0.00001; A=a+data; B=b+datb;xx=geshow(A,B); x1=norm(xx-x,1)/norm(x,1)x2=cond(A)/(1-norm(inv(A))*norm(xx-x))*(norm((xx-x))/(norm(A))+norm(datb)/norm(B)) datx=abs(x1-x2)5.2.3.2 实验要求三的程序运行结果如下： 给出对xx x xx -=∆ˆ的估计是：7.310559817408125e-007xx x xx -=∆ˆ的理论结果是： 3.828481757617297e-007结果相差： 3.482078059790828e-0075.2.4.1 实验要求四的程序如下： for n=4:11n=nHi=hilb(n); cond1Hi=cond(Hi,1) cond2Hi=cond(Hi,2) condinfHi=cond(Hi,inf) pause end5.2.4.2讨论：线性代数方程组的性态与条件数有着很重要的关系，既矩阵的条件数是刻画矩阵性质的一个重要的依据，条件数越大，矩阵“病态”性越严重，在解线性代数方程组的过程中较容易产生比较大的误差，则在实际问题的操作过程中，我们必须要减少对条件数来求解，把条件数较大的矩阵化成条件数较小的矩阵来进行求解。