计算方法与计算 实验一误差分析

1、计算方法与分析实验报告实验一、误差分析一、上机前的准备工作1、复习和掌握与本次实验有关的教学内容。2、根据本次实验要求,在纸上编写算法及上机的程序,并经过人工模拟运 行检验,减少不必要的错误,提高上机效率。切忌不编程序、不作人工检查就进 行程序输入,这只能使上机调试的难度增加,甚至可能带来学习自信心的下降, 影响后续课程的学习。二、上机实验步骤1、启动开发环境;2、建立源程序文件,输入源程序;3、编译产生目标程序,连接生成可执行程序,运行程序,输出结果;4、对数值计算结果进行误差分析,讨论数值算法的收敛性与稳定性;5、整理实验报告。三、实验报告实验报告是记录实验工作全过程的技术文档,实验报告的

2、撰写是科学技术工 作的一个组成部分。 数值分析实验报告包括下列要求:1、实验原理;2、实验内容和要求;3、数值算法描述,包括数据输入、数据处理和数据输出;4、算法的实现(1 给出具体的计算实例,(2 经调试正确的源程序清单,(3 对具体的数值例子给出数值结果;5、计算结果的误差分析,算法的收敛性与稳定性的讨论;6、实验心得。误差分析误差问题是数值分析的基础, 又是数值分析中一个困难的课题。 在实际计算 中,如果选用了不同的算法,由于舍入误差的影响,将会得到截然不同的结果。 因此, 选取算法时注重分析舍入误差的影响, 在实际计算中是十分重要的。 同时, 由于在数值求解过程中用有限的过程代替无限的

3、过程会产生截断误差, 因此算法 的好坏会影响到数值结果的精度。 一、实验目的1、通过上机编程,复习巩固以前所学程序设计语言及上机操作指令; 2、通过上机计算,了解误差、绝对误差、误差界、相对误差界的有关概念;3、 通 过上机计算,了解舍入误差所引起的数值不稳定性。二、算法实例 例 1.1 用差商 ha f h a f a f ( ( (-+'求 x x f ln (=在 3=x 处导数的近似值。取1. 0=h , 1000. 0=h , h =0.000 000 000 000 001和 h =0.000 000 000 000 000 1分别用 MATLAB 软件计算,取十五位数字计

4、算。解: 在 MATLAB 工作窗口输入下面程序>>a=3;h=0.1;y=log(a+h-log(a;yx=y/h 运行后得yx = 0.32789822822991. 将此程序中 h 改为 0.000 1,运行后得yx = 0.33332777790385.后者比前者好。再取 h = 0.000 000 000 000 001,运行后得yx = 0.44408920985006,不如前者好。取 h = 0.000 000 000 000 000 1,运行后得yx = 0,算出的结果反而毫无价值。例 1.2 分别求方程组 b AX =在下列情况时的解,其中 A =01. 1111

5、. (1 =22b ; (2 =01. 22b . 解: (1 首先将方程组 b AX =化为同解方程 b A X 1-=,然后在 MATLAB 工作窗口输入程序>> b=2,2'A=1,1;1,1.01; X=Ab运行后输出当 =22b 时, b AX =的解为 =02X ;(2同理可得,当 =01. 22b 时, b AX =的解为 =11X .例 1.3 计算 e 的近似值。 解:泰勒级数e +=! ! 4! 3! 21432 n x x x x x nx(<<-x ,取 1=x ,得e +=! 1! 41! 31! 2111 n . (1.1 这是一个无

6、限过程,计算机无法求到精确值。只能在(1.1取有限项时计算,再估计误差。如果取有限项! ! ! ! ( n s n 1413121111+=作为 e 的值必然会有误差,根据泰勒余项定理可知其截断误差为e ! 1( 1(+=-n e s n 10(<<.如果取(1.1的前九项,输入程序>> n=8; s=1;S=1; for k=1:n s=s*k; S=S+1/s, end s,S,R=3/(s*(n+1>>n=8; s=1;S=1; for k=1:n s=s*k; S=S+1/s; end s,S,R=3/(s*(n+1运行结果:s =40320 S =

7、2.7183R =8.2672e-006或>>S1=1+1+1/2+1/(1*2*3+1/(1*2*3*4+1/(1*2*3*4*5+1/(1*2*3*4*5*6+1/(1*2*3*4*5*6*7+1/(1*2*3*4*5*6*7*8,R1=3/(1*2*3*4*5*6*7*8*9运行后结果S =8.267195767195768e-006 R =2.71827876984127 颠倒!因为截断误差为e , 10(101968.267! 93! 18( 1(6-8<<<+=- e s 所以 e 的近似值 e +=!81! 71! 61! 51! 41! 31! 2

8、111 1(8 s 2.718 28. 例 1.4 取 282.718作为 e 的四舍五入近似值时,求其绝对误差和相对误差。 解:在 MATLAB 工作窗口输入程序>>juewu=exp(1-2.71828运行后输出结果为juewu = 1.828 459 045 505 326e-006例 1.5 计算 2sin xxd x 的近似值,并确定其绝对误差和相对误差。 解 因为被积函数xxsin 的原函数不是初等函数,故用泰勒级数求之。 +-+-=! ! ! 1sin 429 75 386x x x x x x (<<-x , (1.5 这 是 一 个 无 限 过 程 ,

9、 计 算 机 无 法 求 到 精 确 值 。 可 用 (1.5 的 前 四 项! ! 14275 36x x x -+-代替被积函数 x x sin ,得=2sin x x y d 20x ! ! 14275 36x x x -+-d x =! 7 (! 5 (! 3 (275375 3-+-=y . 根据泰勒余项定理和交错级数收敛性的判别定理,得到绝对误差!992(9<-=yy R = WU, 在 MATLAB 命令窗口输入计算程序如下:syms xf=1-x2/(1*2*3+x4/(1*2*3*4*5-x6/(1*2*3*4*5*6*7 y=int(f,x,0,pi/2,y1=dou

10、ble(yy11=pi/2-(pi/23/(3*3*2+(pi/25/(5*5*4*3*2-(pi/27/(7*7*6*5*4*3*2inf=int(sin(x/x,x,0,pi/2 ,infd=double(inf WU =(pi/29/(9*9*8*7*6*5*4*3*2, R =infd-y11因为运行后输出结果为: =y 1.370 762 168 154 49, y=1.370 744 664 189 38, =R 1.750 396 510 491 47e-005, WU= 1.782 679 830 970 664e-005410-<. 所以, y 的绝对误差为 =410-

11、,故 =20sin xxy d 7 1.370x 。 y的相对误差为 =r 71.370104-=y <0.007 3%. 例 1.6 设计三种算法求方程 01522=-+x x 在 3, 2(的一个正根 *x 的近似值,并研究每种算法的误差传播情况 .解:为解已知方程,我们可以设计如下三种算法,然后将计算结果与此方程 的精确解 5. 2*=x 比较,观察误差的传播 .算法 1 将已知方程化为同解方程 =x 2215x -. 取初值 20=x ,按迭代公式21215k k x x -=+ 依次计算 , , , , 21n x x x ,结果列入表 1 3中。算法 2 将已知方程化为同解方

12、程 1215+=x x . 取初值 20=x ,按迭代公式12151+=+k k x x依次计算 , , , , 21n x x x ,结果列入表 1 1中。算法 3 将已知方程化为同解方程 141522+-+-=x x x x x . 取初值 20=x , 按迭代公式为 1415221+-+-=+k k kk k x x x x x 依次计算 , , , , 21n x x x ,结果列入表 1 1中。我们为这三种算法的计算编写两套 MATLAB 程序如下:(1 MATLAB 主程序function k,juecha,xiangcha,xk= liti112(x0,x1,limax % 输入

13、的量 -x0是初值 , limax是迭代次数和精确值 x;% 输出的量 -每次迭代次数 k 和迭代值 xk,% -每次迭代的绝对误差 juecha 和相对误差 xiangcha , x(1=x0;for i=1:limaxx(i+1=fl(x(i;%程序中调用的 fl.m juecha = abs(x(i-x1;xiangcha = juecha /( abs(x(i+eps;xk=x(i;k=i-1;(i-1,juecha,xiangcha,xk end(2 MATLAB 调用函数程序及其计算结果 算法 2的 MATLAB 调用函数程序function y1=fl(x y1=15/(2*x+

14、1; 将 MATLAB 主程序和调用函数程序分别命名 liti112.m 和 fl.m ,分别保存为 M 文件,然后在 MATLAB 工作窗口输入命令>> k,juecha,xiangcha,xk= liti112(2,2.5,100 运行后输出计算结果列入表 1 1和表 1-2中。 将算法 2的 MATLAB 调用函数程序的函数分别用 y1=15-2*x2和 y1=x-(2*x2+x-15/(4*x+1代替,得到算法 1和算法 3的调用函数程序,将其保 存,运行后将三种算法的前 8个迭代值 821, , , x x x 列在一起(见表 1-1 ,进行 比较 . 将三种算法的 82

15、1, , , x x x 对应的绝对误差和相对误差的值列在一起(见表 1-2 ,进行比较。 例 1.7 求数 18(71915-+=-x 的近似值。 解 (1直接用 MATLAB 命令>> x=(715*(sqrt(1+8(-19-1运行后输出结果x = 0.问题出现在两个相近的数 198-+与 1相减时,计算机运行程序>>sqrt(1+8(-19-1运行后输出结果ans = 0.由于计算机硬件只支持有限位机器数的运算, 因此在计算中可能引入和传播 舍入误差 . 因为有效数字的严重损失, 导致输出 1819-+-的结果为 0, 计算机不 能再与数 157继续进行真实的计

16、算,所以,最后输出的结果与 x 的精确值不符。(2如果化为1887 18(71919151915+=-+=-x ,再用 MATLAB 命令>> x=(715*( (8(-19/(sqrt(1+8(-19+1运行后输出结果x = 1.6471e-005这是因为 1819-+-化为1881919+-后,计算机运行程序>> x= (8(-19/(sqrt(1+8(-19+1运行后的结果为x =3.4694e-018由于有效数字的损失甚少, 所以运算的结果 -18103.4694再与 157继续计算, 最后 输出的结果与 x 的精确值相差无几。例 1.8 求数 13030ln(

17、2-=y 的近似值。 解 (1直接用 MATLAB 程序>> x=30;x1= sqrt(x2-1运行后输出结果x1 = 29.9833输入 MATLAB 程序>> x=30; x1=29.9833;y=log(x-x1运行后输出结果y = -4.0923(2因为 13030ln(2-中的 30=x 很大,如果采用倒数变换法111221-+=-=x x x x z ,即130301ln13030ln(22-+=- 190030ln(-+-=.输入 MATLAB 程序>> x=30;y=-log(x+sqrt(x2-1运行后输出结果y = -4.0941(3输

18、入 MATLAB 程序>> x=30; y=log(x-sqrt(x2-1运行后输出结果y = -4.0941可见, (2计算的近似值比(1的误差小。参加计算的数,有时数量级相差很大 . 如果不注意采取相应的措施,在它们的加减法运 算中,绝对值很小的那个数经常会被绝对值较大的那个数“吃掉” ,不能发挥其作用, 造成计算结果失真。例 1.9 请在 16位十进制数值精度计算机上利用软件 MA TLAB 计算下面的两个数0.30.1111111111111111*+=x 和 0.30.11111111111111111*+=y将计算结果与准确值比较,解释计算结果。 解 在 MA TLAB

19、 工作窗口输入下面程序>> x=111111111111111+0.1+0.3, y=1111111111111111+0.1+0.3运行后输出结果x = 1.111111111111114e+014, y =1.111111111111111e+015从输出的结果可以看出, x *x =, 而 y *y . 为什么 *y 仅仅比 *x 多一位 1, 而 y *y 呢? 这是因为 计算机进行运算时,首先要把参加运算的数写成绝对值小于 1而“阶码”相 同的数, 这一过程称为数的 “对阶” 。 在 16位十进制数值精度计算机上利用软件 MATLAB 计算这两个数,把运算的数 *x 写成

20、浮点规格化形式为,151515*1030000000000000.0001010000000000000.000100111111111111111. 0+=x在 16位十进制数值精度计算机上, 三项的数都表示为小数点后面 16位数字的数与 1510之积,所以,计算机没有对数进行截断,而是按原来的三个数进行计算。因此,计算的 结果 x *x =。而161616*10030000000000000.00010010000000000000.000101111111111111111. 0+=y三项的数都表示写成绝对值小于 1而“阶码”都为 1610的数以后,第一项的纯小数的小数点后面有 16位数

21、字 . 但是,后两项数的纯小数的小数点后面有 17位数字,超过了 16位十进制数值精度计算机的存储量,计算机对后两项的数都进行截断最后一位,即 后两项的数都是 16位机上的零,再进行计算,所以计算结果与实际不符。三 、 实验任务对 20, , 2, 1, 0 =n , 计算定积分+=15dx x x y nn .算法 1:利用递推公式151-=n n y ny , 20, , 2, 1 =n , 取 -=+=10182322. 05ln 6ln 51dx x y . >> a=0.182322;n=1;b=(1/n-(5*a算法 2:利用递推公式n n y n y 51511-=-

22、 1, , 19, 20 =n . 注意到=+=101202010201051515611261dx x dx x x dx x , 取 008730. 0 12611051(20120+y . 思考:从计算结果看,哪个算法是不稳定的,哪个算法是稳定的。>> a=0.008730;n=20;b=(1/5*1/n-(1/5*a解: 由递推公式:y n = 1 -5y n 1 n 可得 当 n =1,2,L ,20 时,以及y 0 0.182322=0.182322 (小数点后面取 6位 使用Matlab计算 ,得出下列各值: Ni 编号 1 2 3 4 5 6 7 8 9 10 11

23、 12 13 14 15 16 17 18 19 20 0.0884 0.0580 0.0433 0.0335 0.0325 0.0042 0.1219 -0.4845 2.5336 -12.5680 62.9309 -314.5712 1.5729e+003 -7.8644e+003 3.9322e+004 -1.9661e+005 9.8305e+005 4.9153e+006 -2.4576e+007 1.2288e+008 计算结果 由以上计算结果可以看出，当 n8 时，出现负值，完全失真，但与原题中： xn n min x + 5 x dx < x + 5dx < max

24、 0 x 1 0 x 1 0 0 1 1 6 ( n + 1 1 1 xn x + 5dx x+5 0 1 1 < y n < 1 5 ( n + 1 完全不符， | E n |= | y n y n |= | 5( y n 1 y n 1 | | E n |= 5 n | E 0 | E 20误 差 被 放 大 5 20 因此第一种算法是不稳定的。 我们来看第 2 种方法： 由递推公式 :y n -1 = 当 y n 可得 5 n = 20， ， ， 时 ,以及 y 20 0.008730=0.008730 19 L 1 5n 1 - 1 ( 小数点后面取 6 位 使用Matlab计算 ,得出下列各值 : 编号 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 0