数值计算方法实验报告Word格式文档下载.doc
- 文档编号:3618020
- 上传时间:2023-05-02
- 格式:DOC
- 页数:36
- 大小:345.50KB
数值计算方法实验报告Word格式文档下载.doc
《数值计算方法实验报告Word格式文档下载.doc》由会员分享,可在线阅读,更多相关《数值计算方法实验报告Word格式文档下载.doc(36页珍藏版)》请在冰点文库上搜索。
每步消去过程总选取按模最小或按模尽可能小的元素作为主元,观察并记录计算结果。
若每步消去过程总选取按模最大的元素作为主元,结果又如何?
分析实验的结果。
(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程序如下所示
A=hilb(n)
gauss
请输入矩阵A的阶数:
n=7
n=
7
A=
6100000
8610000
0861000
0086100
0008610
0000861
0000086
b=
15
14
pp=
3.174999999999999e+002
Ab=
61000007
861000015
086100015
008610015
000861015
000086115
000008614
r=
1
p=
i=
0.999999999998691.000000000043370.99999999964299
1.000000001211430.999999998030381.00000000152825
0.99999999954491
显然的是,该问题在主元选取与算出结果有着很大的关系,取绝对值大的元素作为主元比取绝对值小的元素作为主元时产生的结果比较准确,即选取绝对值小的主元时结果产生了较大的误差,条件数越大产生的误差就越大
实验体会:
运用高斯消去法求解线性方程组问题的时候,主元的选取和相应的消去法的选取决定了该算法的稳定性,选取绝对值大的元素比选取绝对值比较小的元素作为主元时对结果产生的误差影响比较小。
并且增加条件数反而对结果的误差产生更大的影响。
并且在运算中要尽量避免出现运用小数作为除数,使数量级加大,令大数吃掉小数的情况发生。
程序代码:
⑴
保存M文件名为:
fanshu.m
pleaseinputn:
)
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)
x0=norm(xx-x,1)/norm(x,1)
geshow.m
functionx=geshow(A,B)
AB=[AB];
pivot=AB(i,i);
AB(k,i:
nb)=AB(k,i:
nb)-(AB(k,i)/pivot)*AB(i,i:
x(n)=AB(n,nb)/AB(n,n);
x(i)=(AB(i,nb)-AB(i,i+1:
n))/AB(i,i);
⑵
cond2.m
functioncond2(A)
B=A'
*A;
[V1,D1]=eig(B);
[V2,D2]=eig(B^(-1));
cond2A=sqrt(max(max(D1)))*sqrt(max(max(D2)))
保存M文件为:
shiyan52.m
formatlong
forn=10:
10:
100
n=n
A=fix(100*randn(n));
condestA=condest(A)
cond2(A)
condA2=cond(A,2)
pause
⑶
sy5_2.m
)%输入矩阵的阶数
a=fix(100*rand(n))+1;
%随机生成一个矩阵a
x=ones(n,1);
%假设知道方程组的解全为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);
x0=norm(xx-x,1)/norm(x,1)
x00=cond(A)/(1-norm(inv(A))*norm(xx-x))*(norm((xx-x))/(norm(A))+norm(datb)/norm(B))
datx=abs(x0-x00)
(4)
forn=4:
11
n=n
Hi=hilb(n);
cond1Hi=cond(Hi,1)
cond2Hi=cond(Hi,2)
condinfHi=cond(Hi,inf)
pause
实验结果及其分析:
(1)
fanshu
n=6
n=
6
a=
822896806871
91554996764
139681667528
9297154405
641643856610
109892941883
x=
425
371
359
253
284
395
data=
1.0e-005*
Columns1through4
0.6948286229758170.7655167881490020.7093648308580730.118997681558377
0.3170994800608610.7951999011370630.7546866819823610.498364051982143
0.9502220488383550.1868726045543790.2760250769985780.959743958516081
0.0344460805029090.4897643957882310.6797026
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 计算方法 实验 报告
![提示](https://static.bingdoc.com/images/bang_tan.gif)