数值分析实验一：误差分析、误差传播及算法稳定性

毕节学院实验报告

实验名称： 误差分析、误差传播及算法稳定性 实验报告序号： 1 组 别 实验项目 1姓 名 王朝春 同组实验者 实验日期 李亚 杨凯 2012年10月5日 计算In?e?1?xnexdx(n?0,1,?)0并估计误差 实验类别 □√ 1、验证性实验或基础性实验； □ 2、综合性实验 □ 3、设计性实验； □ 4、创新性实验和研究性实验； 指导教师（签名） 赖志柱 年 月 日 教师评语 实验成绩 实验目的： 通过本实验对求解问题的算法进行好坏判断有一个初步了解，并加强对设计一个好算法的理解，体验数值计算稳定性，从而了解数值计算方法的必要性，体会数值计算的收敛性与收敛速度。 实验任务与要求： 计算In?e?1?10xnexdx(n?0,1,?)并估计误差 （1）建立若干个（不少于两个）计算公式； （2）分析计算公式的理论误差； （3）编写程序（推荐MATLAB）实现（1）中的计算公式、输出结果并比较实际误差； （4）任选正整数m?n，要求既从Im计算In，又从In计算Im，并分析您的结果。这里m?0且n?9。 小组分工合作说明：王朝春查找资料；李亚负责编辑，杨凯协助完成。 实验过程及内容： 解： 由分部积分可得计算In的递推公式 I?1?nIn?1,n?1,2,…….??n （1） 1??1x?1I?e?edx?1?e.?0?0若计算出I0，代入（1）式，可逐次求出 I1,I2,…的值。要算出I0就要先算出e?1，若用泰勒多项式展开部分和 (?1)2(?1)ke?1?(?1)??…?2!k!?1, 并取k=19,用4位小数计算，则得e?1?0.3679，截断误差R7?|e?1?0.3679|?11??10?4.计算过程中小数点后第5位的数字按四8!4舍五入原则舍入，由此产生的舍入误差这里先不讨论。当初值取为?时，用（1）式递推的计算公式为 I0?0.6321?I0??0.6321?I0（A），n=1,2，…。 ????In?1?nIn?1?就是初值计算结果见表1的I?n列。用I?0近似I0产生的误差E0?I0?I0误差，它对后面计算结果是有影响的. 从表1中看到I?18出现负值，这与一切In?0相矛盾。实际上，由积分估值得 11e?11?e?1(mimex)?xndx?In?e?1(max)?xndx? （2） 000?x?10?x?1n?1n?1因此，当n较大时，用I?n近似In显然是不正确的。这里计算公式与每步计算都是正确的，那么是什么原因合计算结果出现错误呢？主?，由此引起以后各步计算的误差要就是初值I?0有误差E0?I0?I0?满足关系 En?In?InEn??nEn?1,n?1,2,…. 由此容易推得 En?(?1)nn!E0， 这说明I?0有误差E0，则I?n就是E0的n!倍误差。例如，n=19，若|E0|?1?完全不能近似I了。?10?4，则|E19|?19!?|E0|?2。这就说明I1982它表明计算公式（A）是数值不稳定的。 我们现在换一种计算方案。由（2）式取n=19,得 e?11?I19?， 202011e?1我们粗略取I9?(?)?0.0684?I9*，然后将公式（1）倒过来算，21010***即由I9算出I8*，I7，…，I0，公式为 *?I19?0.0342? (B)??*1*?In?1?(1?In),n?19，18，…，1；n?**计算结果见表1的In列。我们发现I0与I0的误差不超过10?4。记**，则|E0*|?En?In?In1***|En|，E0比En缩小了n!倍，因此，尽管E9*较n!*大，但由于误差逐步缩小，故可用In近似In。反之，当用方案（A）计算时，尽管初值I?0相当准确，由于误差传播是逐步扩大的，因而计算结果不可靠。此例说明，数值不稳定的算法是不能使用的。 程序如下： function x11 = facto(n) %这个函数的功能是求 n 的阶乘； x11 = 1; if n == 0 x11 = 1; else for i = 1:n x11 = x11*i; 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('M 3.4f .4f |\\n',i,I(i+1),B(m+1-i)); end 在Matlab命令窗口输入如下命令即可得到如表1的结果。 jifen(19) 表1 计算结果 n 0 1 2 第1种算法 0.6321 0.3679 0.2642 第2种算法 0.6321 0.3679 0.2642 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 0.2073 0.1709 0.1455 0.1268 0.1124 0.1009 0.0916 0.0839 0.0774 0.0718 0.0669 0.0627 0.0590 0.0555 0.0572 -0.0295 1.5596 0.2073 0.1709 0.1455 0.1268 0.1124 0.1009 0.0916 0.0839 0.0774 0.0718 0.0669 0.0627 0.0590 0.0557 0.0527 0.0508 0.0342