数值计算方法实验报告Word格式.doc
- 文档编号:13036572
- 上传时间:2022-10-03
- 格式:DOC
- 页数:36
- 大小:345.50KB
数值计算方法实验报告Word格式.doc
《数值计算方法实验报告Word格式.doc》由会员分享,可在线阅读,更多相关《数值计算方法实验报告Word格式.doc(36页珍藏版)》请在冰豆网上搜索。
问题提出:
Gauss消去法是我们在线性代数中已经熟悉的。
但由于计算机的数值运算是在一个有限的浮点数集合上进行的,如何才能确保Gauss消去法作为数值算法的稳定性呢?
Gauss消去法从理论算法到数值算法,其关键是主元的选择。
主元的选择从数学理论上看起来平凡,它却是数值分析中十分典型的问题。
实验内容:
考虑线性方程组
编制一个能自动选取主元,又能手动选取主元的求解线性方程组的Gauss消去过程。
实验要求:
(1)取矩阵,则方程有解。
取n=10计算矩阵的条件数。
让程序自动选取主元,结果如何?
(2)现选择程序中手动选取主元的功能。
每步消去过程总选取按模最小或按模尽可能小的元素作为主元,观察并记录计算结果。
若每步消去过程总选取按模最大的元素作为主元,结果又如何?
分析实验的结果。
(3)取矩阵阶数n=20或者更大,重复上述实验过程,观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异,说明主元素的选取在消去过程中的作用。
(4)选取其他你感兴趣的问题或者随机生成矩阵,计算其条件数。
重复上述实验,观察记录并分析实验结果。
实验5.2(线性代数方程组的性态与条件数的估计)
理论上,线性代数方程组的摄动满足
矩阵的条件数确实是对矩阵病态性的刻画,但在实际应用中直接计算它显然不现实,因为计算通常要比求解方程还困难。
MATLAB中提供有函数“condest”可以用来估计矩阵的条件数,它给出的是按1-范数的条件数。
首先构造非奇异矩阵A和右端,使得方程是可以精确求解的。
再人为地引进系数矩阵和右端的摄动,使得充分小。
(1)假设方程Ax=b的解为x,求解方程,以1-范数,给出的计算结果。
(2)选择一系列维数递增的矩阵(可以是随机生成的),比较函数“condest”所需机器时间的差别.考虑若干逆是已知的矩阵,借助函数“eig”很容易给出cond2(A)的数值。
将它与函数“cond(A,2)”所得到的结果进行比较。
(3)利用“condest”给出矩阵A条件数的估计,针对
(1)中的结果给出的理论估计,并将它与
(1)给出的计算结果进行比较,分析所得结果。
注意,如果给出了cond(A)和的估计,马上就可以给出的估计。
(4)估计著名的Hilbert矩阵的条件数。
思考题一:
(Vadermonde矩阵)设
,
其中,,
(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矩阵
实验程序:
M文件程序为:
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)
p=input('
条件数对应的范数是p-范数:
p='
pp=cond(A,p)
pause
[m,n]=size(A);
nb=n+1;
Ab=[Ab]
r=input('
请输入是否为手动,手动输入1,自动输入0:
r='
fori=1:
n-1
ifr==0
[pivot,p]=max(abs(Ab(i:
n,i)));
ip=p+i-1;
ifip~=i
Ab([iip],:
)=Ab([ipi],:
);
disp(Ab);
pause
end
end
ifr==1
i=i
ip=input('
输入i列所选元素所处的行数:
ip='
Ab([iip],:
end
pivot=Ab(i,i);
fork=i+1:
n
Ab(k,i:
nb)=Ab(k,i:
nb)-(Ab(k,i)/pivot)*Ab(i,i:
nb);
end
disp(Ab);
end
x=zeros(n,1);
x(n)=Ab(n,nb)/Ab(n,n);
fori=n-1:
-1:
1
x(i)=(Ab(i,nb)-Ab(i,i+1:
n)*x(i+1:
n))/Ab(i,i);
(1)⑴取矩阵A的阶数:
n=10,自动选取主元:
>
formatlong
gauss
n=10
n=10
p=1
p=1
pp=2.557500000000000e+003
r=0
r=0
⑵取矩阵A的阶数:
n=10,手动选取主元:
①选取绝对值最大的元素为主元:
p=2
p=2
pp=1.727556024913903e+003
r=1
r=1
ans=1111111111
②选取绝对值最小的元素为主元:
pp=1.727556024913903e+003
ans=
1.000000000000001.000000000000001.00000000000000
1.000000000000001.000000000000001.00000000000000
0.999999999999991.000000000000010.99999999999998
1.00000000000003
(2)取矩阵A的阶数:
(3)取矩阵A的阶数:
n=20,手动选取主元:
①选取绝对值最大的元素为主元:
n=20
pp=2.621437500000000e+006
ans=11111111111111111111
②选取绝对值最小的元素为主元:
n=20.
n=20
pp=1.789670565881683e+006
1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000010.999999999999971.00000000000006
0.999999999999891.000000000000230.99999999999955
1.000000000000900.999999999998211.00000000000352
0.999999999993181.000000000012730.99999999997817
1.00000000002910
(4)该题目的M程序如下所示
n=input(
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 计算方法 实验 报告