误差分析、误差传播及算法稳定性
毕节学院实验报告
实验名称: 误差分析、误差传播及算法稳定性 实验报告序号: 1 组 别
姓 名 付制令 同组实验者 316寝室 实验项目
计算1
10
e e n x n I x dx -=?(0,1,)
n =并估计误差
实验日期
2013年4月3日
实验类别
□√ 1、验证性实验或基础性实验; □ 2、综合性实验
□ 3、设计性实验; □ 4、创新性实验和研究性实验;
教师评语
实验成绩
指导教师(签名)
赖志柱
年 月 日
实验目的:
通过本实验对求解问题的算法进行好坏判断有一个初步了解,并加强对设计一个好算法的理解,体验数值计算稳定性,从而了解数值计算方法的必要性,体会数值计算的收敛性与收敛速度。
实验任务与要求: 计算1
1
e
e n x
n
I
x dx -=?(0,1,)n =并估计误差
(1)建立若干个(不少于两个)计算公式;
(2)分析计算公式的理论误差;
(3)编写程序(推荐MATLAB )实现(1)中的计算公式、输出结果并比较实际误差;
(4)任选正整数m
n ,要求既从m I 计算n I ,又
从n
I 计算m
I ,并分析您的结果。这里0m ≠且9n ≠。
小组分工合作说明:
实验过程及内容: 解:
由分部积分可得计算n I 的递推公式
1111
01,1,2,e 1.n
n x I nI n I e dx e ---=-=???==-???……. (1) 若计算出0I ,代入(1)式,可逐次求出 1
2,,I I …
的值。
要算出0I 就要先算出1e -,若用泰勒多项式展开部分和
21
(1)(1)1(1),2!!
k
e k ---≈+-+++
…
并取k=19,用4位小数计算,则得10.3679e -≈,截断误差
14711
|0.3679|108!4
R e --=-≤
.计算过程中小数点后第5位的数字按四舍五入原则舍入,由此产生的舍入误差这里先不讨论。当初值取为
000.6321I I ≈=时,用(1)式递推的计算公式为
01
0.6321
A 1n n I I nI -?=?
=-?(),n=1,2,…。 计算结果见表1的n I 列。用0I 近似0I 产生的误差000E I I =-就是初值误差,它对后面计算结果是有影响的.
从表1中看到18I 出现负值,这与一切0n I >相矛盾。实际上,由积分估值得
111110001011
(im )(max)11
x n n n x x e e m e x dx I e x dx n n ---≤≤≤≤=<<=++?? (2) 因此,当n 较大时,用n I 近似n I 显然是不正确的。这里计算公式与
每步计算都是正确的,那么是什么原因合计算结果出现错误呢?主要就是初值0I 有误差000E I I =-,由此引起以后各步计算的误差
n n n E I I =-满足关系
1,1,2,n n E nE n -=-=….
由此容易推得
0(1)!n n E n E =-,
这说明0I 有误差0E ,则n I 就是0E 的n!倍误差。例如,n=19,若
401
||102
E -=
?,则190||19!||2E E =?>。这就说明8I 完全不能近似19I 了。它表明计算公式(A )是数值不稳定的。 我们现在换一种计算方案。由(2)式取n=19,得
11912020
e I -<<, 我们粗略取1
*9911()0.068421010
e I I -≈+==,然后将公式(1)倒过来算,
即由*9I 算出*8I ,*7I ,…,*
0I ,公式为
*
19**
10.0342()1(1),198n n I B I I n n -?=?
=?=-=??
,1, (1)
计算结果见表1的*n I 列。我们发现*
0I 与0I 的误差不超过410-。记**
n n n
E I I =-,则**01||||!
n E E n =,*0E 比*
n E 缩小了n!倍,因此,尽管*9E 较大,但由于误差逐步缩小,故可用*
n I 近似n I 。反之,当用方案(A )
计算时,尽管初值0I 相当准确,由于误差传播是逐步扩大的,因而计算结果不可靠。此例说明,数值不稳定的算法是不能使用的。
程序如下:
function x11 = facto(n)
这个函数的功能是求n 的阶乘;
x11 = 1;
if n == 0
x11 = 1;
else
end
end
function e_1 = telor(k)
这个函数的功能是求e^(-1);用泰勒多项式展开式进行计算的,k是代表展开到第k+1 项
e_1 = 1;
if k == 1
e_1 = 1;
else
for i=1:k
e_1 = e_1 +(-1)^i/facto(i);
end
end
function jifen(m)
I0=1-telor(19);
*第一种算法
I(1)=I0;
for i = 1:m
I1 = 1 - i*I0;
I(i+1)=I1;
I0=I1;
end
*第二种算法
Im=(1/2)*(1/(m+1)+telor(19)/(m+1));
B(1)=Im;
for i=1:m
In=(1/(m+1-i))*(1-Im);
B(i+1)=In;
Im=In;
end
disp(' n 第 1 种算法第2 种算法');
for i = 0:(length(B)-1)
fprintf('%4d %33.4f %12.4f |\n',i,I(i+1),B(m+1-i));
end
在Matlab命令窗口输入如下命令即可得到如表1的结果。
jifen(19)
表1 计算结果
n 第1种算法第2种算法
0 0.6321 0.6321
1 0.3679 0.3679
2 0.2642 0.2642
3 0.2073 0.2073
4 0.1709 0.1709
5 0.1455 0.1455
6 0.1268 0.1268
7 0.1124 0.1124
8 0.1009 0.1009
9 0.0916 0.0916
10 0.0839 0.0839
11 0.0774 0.0774
12 0.0718 0.0718
13 0.0669 0.0669
14 0.0627 0.0627
15 0.0590 0.0590
16 0.0555 0.0557
17 0.0572 0.0527
18 -0.0295 0.0508
19 1.5596 0.0342