数值分析实验报告

实验五 解线性方程组的直接方法

实验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 n

i i n i i n n n n n

n n

x x x x b x x x x x x x x x x x x A 0020

10022222121102001111 ，，

其中，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:4

p=input('条件数对应的范数是p-范数：p=')

pp=cond(A,p)

end

pause

[m,n]=size(A);

nb=n+1;Ab=[A b]

r=input('请输入是否为手动，手动输入1，自动输入0：r=')

for i=1:n-1

if r==0

[pivot,p]=max(abs(Ab(i:n,i)));

ip=p+i-1;

if ip~=i

Ab([i ip],:)=Ab([ip i],:);disp(Ab); pause

end

end

if r==1

i=i

ip=input('输入i列所选元素所处的行数：ip=');

Ab([i ip],:)=Ab([ip i],:);disp(Ab); pause

end

pivot=Ab(i,i);

for k=i+1:n

Ab(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb);

end

disp(Ab); pause

end

x=zeros(n,1);x(n)=Ab(n,nb)/Ab(n,n);

for i=n-1:-1:1

x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);

end

5.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.0022

0.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:4

p=input('条件数对应的范数是p-范数：p=')

pp=cond(A,p)

end

pause

[m,n]=size(A);

nb=n+1;Ab=[A b]

r=input('请输入是否为手动，手动输入1，自动输入0：r=')

for i=1:n-1

if r==0

[pivot,p]=max(abs(Ab(i:n,i)));

ip=p+i-1;

if ip~=i

Ab([i ip],:)=Ab([ip i],:);disp(Ab); pause

end

end

if r==1

i=i

ip=input('输入i列所选元素所处的行数：ip=');

Ab([i ip],:)=Ab([ip i],:);disp(Ab); pause

end

pivot=Ab(i,i);

for k=i+1:n

Ab(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb);

end

disp(Ab); pause

end

x=zeros(n,1);x(n)=Ab(n,nb)/Ab(n,n);

for i=n-1:-1:1

x(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

②选取绝对值最大的元素为主元结果为：