数值计算方法实验报告.docx
- 文档编号:10274794
- 上传时间:2023-05-24
- 格式:DOCX
- 页数:45
- 大小:207.28KB
数值计算方法实验报告.docx
《数值计算方法实验报告.docx》由会员分享,可在线阅读,更多相关《数值计算方法实验报告.docx(45页珍藏版)》请在冰点文库上搜索。
数值计算方法实验报告
重庆交通大学
学生实验报告
实验课程名称数值计算方法I
开课实验室数学实验室
学院理学院年级11专业班信息与计算科学
学生姓名李伟凯学号************
开课时间2013至2014学年第1学期
评分细则
评分
报告表述的清晰程度和完整性(20分)
程序设计的正确性(40分)
实验结果的分析(30分)
实验方法的创新性(10分)
总成绩
教师签名
邹昌文
实验五解线性方程组的直接方法
实验5.1(主元的选取与算法的稳定性)
问题提出:
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],:
)=Ab([ipi],:
);disp(Ab);pause
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);pause
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);
end
(1)⑴取矩阵A的阶数:
n=10,自动选取主元:
>>formatlong
>>gauss
请输入矩阵A的阶数:
n=10
n=10
条件数对应的范数是p-范数:
p=1
p=1
pp=2.557500000000000e+003
请输入是否为手动,手动输入1,自动输入0:
r=0
r=0
⑵取矩阵A的阶数:
n=10,手动选取主元:
①选取绝对值最大的元素为主元:
>>gauss
请输入矩阵A的阶数:
n=10
n=10
条件数对应的范数是p-范数:
p=2
p=2
pp=1.727556024913903e+003
请输入是否为手动,手动输入1,自动输入0:
r=1
r=1
ans=1111111111
②选取绝对值最小的元素为主元:
>>gauss
请输入矩阵A的阶数:
n=10
n=10
条件数对应的范数是p-范数:
p=2
p=2
pp=1.727556024913903e+003
请输入是否为手动,手动输入1,自动输入0:
r=1
r=1
ans=
1.000000000000001.000000000000001.00000000000000
1.000000000000001.000000000000001.00000000000000
0.999999999999991.000000000000010.99999999999998
1.00000000000003
(2)取矩阵A的阶数:
n=10,手动选取主元:
①选取绝对值最大的元素为主元:
>>gauss
请输入矩阵A的阶数:
n=10
n=10
条件数对应的范数是p-范数:
p=2
p=2
pp=1.727556024913903e+003
请输入是否为手动,手动输入1,自动输入0:
r=1
r=1
ans=1111111111
②选取绝对值最小的元素为主元:
>>gauss
请输入矩阵A的阶数:
n=10
n=10
条件数对应的范数是p-范数:
p=2
p=2
pp=1.727556024913903e+003
请输入是否为手动,手动输入1,自动输入0:
r=1
r=1
ans=
1.000000000000001.000000000000001.00000000000000
1.000000000000001.000000000000001.00000000000000
0.999999999999991.000000000000010.99999999999998
1.00000000000003
(3)取矩阵A的阶数:
n=20,手动选取主元:
1选取绝对值最大的元素为主元:
>>gauss
请输入矩阵A的阶数:
n=20
条件数对应的范数是p-范数:
p=1
p=1
pp=2.621437500000000e+006
ans=11111111111111111111
2选取绝对值最小的元素为主元:
>>gauss
请输入矩阵A的阶数:
n=20.
n=20
条件数对应的范数是p-范数:
p=2
p=2
pp=1.789670565881683e+006
请输入是否为手动,手动输入1,自动输入0:
r=1
r=1
ans=
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程序如下所示
functionx=gauss(n,r)
n=input('请输入矩阵A的阶数:
n=')
A=hilb(n)
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],:
)=Ab([ipi],:
);disp(Ab);pause
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);pause
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);
end
>>gauss
请输入矩阵A的阶数:
n=7
n=
7
A=
6100000
8610000
0861000
0086100
0008610
0000861
0000086
b=
7
15
15
15
15
15
14
pp=
3.174999999999999e+002
Ab=
61000007
861000015
086100015
008610015
000861015
000086115
000008614
请输入是否为手动,手动输入1,自动输入0:
r=1
r=
1
条件数对应的范数是p-范数:
p=1
p=
1
i=
1
ans=
0.999999999998691.000000000043370.99999999964299
1.000000001211430.999999998030381.00000000152825
0.99999999954491
显然的是,该问题在主元选取与算出结果有着很大的关系,取绝对值大的元素作为主元比取绝对值小的元素作为主元时产生的结果比较准确,即选取绝对值小的主元时结果产生了较大的误差,条件数越大产生的误差就越大
实验体会:
运用高斯消去法求解线性方程组问题的时候,主元的选取和相应的消去法的选取决定了该算法的稳定性,选取绝对值大的元素比选取绝对值比较小的元素作为主元时对结果产生的误差影响比较小。
并且增加条件数反而对结果的误差产生更大的影响。
并且在运算中要尽量避免出现运用小数作为除数,使数量级加大,令大数吃掉小数的情况发生。
实验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矩阵的条件数。
程序代码:
1
保存M文件名为:
fanshu.m
n=input('pleaseinputn:
n=')
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)
保存M文件名为:
geshow.m
functionx=geshow(A,B)
[m,n]=size(A);
nb=n+1;AB=[AB];
fori=1:
n-1
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
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);
end
2
保存M文件名为:
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)))
end
保存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
end
3
保存M文件为:
sy5_2.m
n=input('pleaseinputn:
n=')%输入矩阵的阶数
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)
formatlong
forn=4:
11
n=n
Hi=hilb(n);
cond1Hi=cond(Hi,1)
cond2Hi=cond(Hi,2)
condinfHi=cond(Hi,inf)
pause
end
实验结果及其分析:
(1)
>>fanshu
pleaseinputn:
n=6
n=
6
a=
822896806871
91554996764
139681667528
9297154405
641643856610
109892941883
x=
1
1
1
1
1
1
b=
425
371
359
253
284
395
data=
1.0e-005*
Columns1through4
0.6948286229758170.7655167881490020.7093648308580730.118997681558377
0.3170994800608610.7951999011370630.7546866819823610.498364051982143
0.9502220488383550.186********43790.2760250769985780.959743958516081
0.0344460805029090.4897643957882310.6797026768536750.340385726666133
0.4387443596563980.445586*********0.6550980039738410.585267750979777
0.3815584570930090.646313*********0.1626117351946310.223811939491137
Columns5through6
0.7512670593056530.547215529963803
0.2550951154592690.138********8679
0.5059570516651420.149294005559058
0.6990767226566860.257508254123737
0.8909032525357990.840717255983663
0.9592914252054450.254282178971531
datb=
1.0e-005*
0.814284826068817
0.243524968724989
0.929263623187228
0.349983765984809
0.196595250431208
0.251083857976031
A=
Columns1through4
82.00000694828622728.00000765516788196.00000709364830780.000001189976814
91.00000317099480255.00000795199900949.00000754686681996.000004983640522
13.00000950222048896.00000186872604581.00000276025076666.000009597439586
92.00000034446080197.00000489764396115.0000067970267694.000003403857266
64.00000438744359616.00000445586200843.00000655098003985.000005852677504
10.00000381558457098.00000646313010692.00000162611735394.000002238119393
Columns5through6
68.00000751267059271.000005472155294
76.0000025509511524.000001386244429
75.00000505957051228.000001492940054
40.0000069907672265.000002575082541
66.00000890903253010.000008407172560
18.00000959291425383.000002542821790
B=
1.0e+002*
4.250000081428483
3.710000024352497
3.590000092926362
2.530000034998376
2.840000019659525
3.950000025108386
xx=
1.000000428381561
0.999999831084732
1.000002591353244
0.999999600566705
0.999998226539937
0.999997826103164
x0=
1.255906711054392e-006
的计算结果为:
1.255906711054392e-0
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 计算方法 实验 报告
![提示](https://static.bingdoc.com/images/bang_tan.gif)