数值分析实验报告

  • 格式:doc
  • 大小:3.58 MB
  • 文档页数:41

下载文档原格式

  / 41
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

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

实验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)选取其他你感兴趣的问题或者随机生成矩阵,计算其条件数。重复上述实验,观察记录并分析实验结果。

思考题一:(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

10022

222121102001111 ,, 其中,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)你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因吗?

相关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程序:

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,手动选取主元:

选取绝对值最大的元素为主元: