重庆交通大学学生实验报告
实验课程名称数值计算方法I
开课实验室数学实验室
学院理学院年级11专业班信息与计算科学学生姓名李伟凯学号631122020203
开课时间2013 至2014 学年第 1 学期
实验五 解线性方程组的直接方法
实验5.1 （主元的选取与算法的稳定性）
问题提出：Gauss 消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss 消去法作为数值算法的稳定性呢？Gauss 消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。
实验内容：考虑线性方程组
n n n R b R A b Ax ∈∈=⨯,,
编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss 消去过程。
实验要求：
（1）取矩阵⎥⎥
⎥
⎥⎥
⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816 b A ，则方程有解T x )1,,1,1(* =。
取n=10计算矩阵的条件数。让程序自动选取主元，结果如何？
（2）现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。
（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。
（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。
实验5.2（线性代数方程组的性态与条件数的估计） 问题提出：理论上，线性代数方程组b Ax =的摄动满足
⎪⎪⎭
⎫
⎝⎛∆+∆∆-≤∆-b b A A A
A A cond x x
11)( 矩阵的条件数确实是对矩阵病态性的刻画，但在实际应用中直接计算它显然不现实，因为计算1-A 通常要比求解方程b Ax =还困难。
实验内容：MATLAB 中提供有函数“condest ”可以用来估计矩阵的条件数，
它给出的是按1-范数的条件数。首先构造非奇异矩阵A 和右端，使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动b A ∆∆和，使得b A ∆∆和充
分小。
实验要求：
（1）假设方程Ax=b 的解为x ，求解方程b b x
A A ∆+=∆+ˆ)(，以1-范数，给出
x
x x x
x -=
∆ˆ的计算结果。
（2）选择一系列维数递增的矩阵（可以是随机生成的），比较函数“condest ”
所需机器时间的差别.考虑若干逆是已知的矩阵，借助函数“eig ”很容易给出cond 2(A)的数值。将它与函数“cond(A,2)”所得到的结果进行比较。
（3）利用“condest ”给出矩阵A 条件数的估计，针对（1）中的结果给出
x
x ∆的理论估计，并将它与（1）给出的计算结果进行比较，分析所得结果。注意，如果给出了cond(A)和A 的估计，马上就可以给出1-A 的估计。 （4）估计著名的Hilbert 矩阵的条件数。
n j i j i h h H j i n n j i ,,2,1,,
1
1
,
)(,, =-+=
=⨯
思考题一：（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
100222221211020011
1
1
，，
其中，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）你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因吗？
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)
p=input('条件数对应的范数是p-范数：p=')
pp=cond(A,p)
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
