数值积分上机实验报告

1、数值积分上机实验报告141110038桂贤进题一：数学上已经证明了f1 4ydx = nJ。1+、2成立，所以可以通过数值积分来求71的近似值。1 .分别使用复合梯形、复合Simpson求积公式计算11的近似值。选择不同的h,对每种求 枳公式，试将误差刻画为h的函数，并比较两方法的精度。是否存在某个值，当低于这 个值之后，再继续减小h的值，计算精度不再有所改进，为什么？2 .实现Romberg求积方法，并重复上面的计算；3 .实现自适应积分方法，并重复上面的计算。解：1.1公式分析：(a)对于复合梯形公式hb - a器(D =5 Lf 3) + f 3) + 2 273 + 访)，h=f-i=

2、l离散误差为：,、 nh3 、 h2(b - a) ,、E(f) = -:Ga gb (2) JL4工乙(b)对于复合Simpson公式八mm-1%)= g / + f。) + 4W f(a + - l)h) + 2 fa + 2ih)(3)i=lf=l. b-a b-ah =2m n离散误差为：Em (T)= 一尊八旬 G) = 一 /i4(i8o n)/(4) ). a b. (4)1.2 实现算法结果：分别利用梯形公式和Simpson公式计算结果如下：(下表中Ei(f) = n-T(f)fE2(f)=忻- S(f).此处Tl为MATLAB中的数，可以认为具有 足够大的精度)1打T(f)E

3、i(f)S(f)E2(f)113.0000000000000000.141621/23.1000000000000000.04163.1333333333333330.008341/43.1311764705882350.010423.1415686274509802.4026e-0561/63.1369630664712640.004633.1415917809360438.7265e-0781/83.1389884944910890.002603.1415925024587071.5113e-07101/103.1399259889071590.001673.141592613939215

4、3.9651e-08121/123.1404352468468510.001163.1415926403053801.3284e-08201/203.1411759869541294.1667e-043.1415926529697856.2001e-10301/303.1414074684073301.8519e-043.1415926535353595.4434e-ll401/403.1414884869236121.0417e-043.1415926535801059.6878e-12501/503.1415259869232546.6667e-053.1415926535872532.5

5、402e-121001/1003.1415759869231291.6667e-053.1415926535897533.9968e-142001/2003.1415884869231304.1667e-063.1415926535897930从上表中可以看出：复合Simpson公式比复合梯形公式精度高，误差收敛的速度快不少。1.3 误差下降速度对比：误差下降速度对比0 2 4C 4- 00 1 | 11 - 7 o o o o o o O 1 1 1 1 1 1 1-e复合梯形公式0复合Simpson公式20406080100120140160180200区间数n从上图可以看出，复合Simp

6、son公式误差的收敛速度比复合梯形公式的误差的收敛速度 快不少，下面验证收敛阶。1.4 验证收敛阶：本实验的实际误差主要由离散误差和计算过程中的舍入误差组成，这里离散误差起 主导作用，故理论上实际误差的收敛阶应该与离散误差的收敛阶相同。下面利用如下公式来il算实际的收敛阶，并与理论分析所得出的离散误差的收敛阶作比较。errl /ii err 2 h2)对上表格中所列的区间长度hi值，逐次利用相邻两个小区间长l】i通过上述公式来计算收 敛阶，并绘制成图形。得到图形如下：验证收敛阶i=12345E78HD1,12(a)对复合梯形公式：由上面公式(2)可知，离散误差关于h为二阶收敛，同时由上图可知实

7、验结果的收敛阶将 近为2,故与理论分析相符。(b)对复合Simpson公式：这里却有些奇怪，由上面公式(4)可知，离散误差理论上为4阶收敛，可实验结果却是将 近6阶收敛。下面将进一步深入探究。探究如下：考虑这是由于被积函数f(x)的特殊性导致，而不是由于Simpson公式离散误 差真的能达到6阶收敛。由误差的余项公式Em(/)= 一 */(9 =工 b.考虑到f(x)=1.5 计算过程中舍入误差的影响从上表格可以看出：当区间数n从1到200时，随着h的减小，实际误差在减小。考虑如下问题：若不断减小h的值，即不断增加区间数n,是否实际误差会一直减小？ 1.5.1理论分析：我们知道影响实验结果的精

8、确度的因素主要有离散误差和舍入误差，而离散误差的 大小可以通过离散误差的余项来体现。由公式（2）和（4）,可以看出当不断增加区间 数，即区间长度h不断减小时，离散误差会越来越小。但相应地由于计算精度的限制, 当h不断减小时，舍入误差却会变大。故理论上会存在一个阈值H。肖h大于H时，由于离散误差起主导作用，随着区间长度h的减小，实际误差会变小: 而当h小于H时，此时舍入误差将在计算结果的精确度起主导作用，再减小h,会导致计 算结果的精确度基本保持不变甚至可能会有减小。1.5.2实验检验：（a）对于复合梯形公式：注意到误差收敛的速度较小，故我们首先选取区间数n= 10,注1,1O，1O3 .108

9、进行分析，得到下面图形。10复合梯形公式的误差10分段区间数nB从图中可以看出，复合梯形公式的阈值H在区间数n为10A6到10A8之间取到。下面对 n处于区间10A6到10A8进行分析。一 - 1010复合梯形公式的误差。一复合梯形107分段区间数n10由上图可以看出在n取10A6 （对应h=10A_6）左右h达到阈值，故可取阈值H=10a-6.（b）对于复合Simpson公式：注意到复合Simpson公式得收敛速度较快，取i从0到100 （对应区间数m=2i）,逐次计 算出对应于m的实际误差,并作图如下。从图中可以看出,在i=42（对应m=84,h=l/l68） 左右，h达到阈值，故可取阈值

10、H=l/l68。复a Simpson公式的1天差随区间薮m的变化情况102030405060708090100i .其中分段区间数m=2i通过上述理论分析和实验验证，我们得到如下结果：使用复合梯形公式和复合Simpson公式计算积分值时，所分成的小区间长h都存在阈值 H,当hH时，再减小h的值，计算精度不再有所改进。原因如下：当h小于H时，此时舍入误差将在计算结果的精确度起主导作用，再减小h, 虽然离散误差依旧会减小，但舍入误差会增大，这就导致计算结果的精确度基本保持不 变甚至可能会有减小。2.1公式分析对Romberg求积方法：Ti,i=y(/(a)+/W), hb-a 211b a 1Gi

11、 = 51一1，1 + %_i f(a + (k -= 2,3,.，% = -r =万儿一】k=lTmj =上与= 2,3,.分析离散误差为：对序列小力m= 1,2,3.,考虑1 = |，7(%)八一7/离散误差 的数量级为0(h2J), h为序列的小区间长。(下表中E(/)=恒一 (D|，这里冗为MATLAB中的数，可以认为是足够精确的)IruEg1130. 141621/23. 1333333333333330. 008331/43. 1421176470588235. 2499e-0441/83.1415857837618746. 8698e-0651/163. 141592665277

12、7171. 1688e-0861/323.1415926536382444. 8451e-1171/643. 1415926535897237.0166e-1481/1283. 141592653589793091/2563. 1415926535897918. 8818e-16101/5123.1415926535897928. 8818e-16111/10243. 1415926535897921. 3323e-15121/20483. 1415926535897944. 4409e-16131/40963.1415926535897921. 3323e-15141/81923. 1415

13、926535897921. 3323e-15151/163843.1415926535897914. 4409e-162.2考虑舍入误差的影响：2.2.1 理论分析：如同1.5小节中所作分析，此处Romberg积分方法同样存在阈值H。当h小于H时, 此时舍入误差将在计算结果的精确度起主导作用，再减小h,虽然离散误差依旧会减小, 但余人误差会增大，这就导致计算结果的精确度基本保持不变甚至可能会有减小。2.2.2 实验12:取 i=l,2,3,4,5,6,7,8,9,10,11,12,13,14,15,逐次计算EQ。,并作图如下：从图中可以看出，当i=9(对应h=l/256 )时，再减小h的值(即

14、增大i值)，误差基本上 不变。故可取阈值H=l/256.3.1自适应算法实现：算法思想：通过递归调用实现。(此处镇觉课本提供的算法虽然巧妙，但比较冗杂，且 会导致占用太多内存。如用递归调用算法直观上更容易理解。)实现算法如下：(自己感觉还算比较巧妙)L脚本文件zishiying.mclear;a=O;b=l;disp(,误差限为：),e=O.OOOOO(Hh=(b-a)/2;fl=f(a);f2=f(a+b)/2);f3=f(b);ticS=h/3*(fl+f2+4*f3);t=Simpson_auto(a/b/e/Szh/2/fl/f2/f3);dispC近似值为:),tdisp(误差为：)

15、,abs(pi-t)disp(,耗时为：力toe2.函数文件 Simpson_auto.mfunction y =Simpson_auto(A,B,e,S,h,Cl,C2,C3) %将已计算好的函数值传递下去，避免重复 计算fl=f(A+h);f2=flA+3*h);%利用Simpson公式分别计算左半区间和右半区间的近似值Sl=h/3*(Cl+C2+4*fl);S2=h/3*(C2+C3+4*f2);%判断是否达到所需精度，若未达到将区间分半，进行递归调用if abs(S-Sl-S2)0其次由上面推导可知：无穷积分察为dx收敛。考虑用区间截去法，将无穷积分转化为在有限区间上的积分。注意到：0

16、03I k00 x3 wdx5;00 120120=dx = .6.2899450830.103892565E- 130.3225476900.6031541041.74576110103574186924.5366202970.388879085E-19.3950709120.539294706E-340.2635603200.5217556111.4134030590.3986668113.5964257710.759424497E-17.0858100060361175868E-212.6408008440.233699724E-450.2228466040.4589646741.188

17、9321020.4170008312.9927363260.1133733825.7751435690.103991975E 19.8374674180.261017203E-315.9828739810.898547906E-6这里若取误差容限为： = 10-7 = b = 1.2 * 109但应用mathematica求解：C0043dx = 1.138566453220937 x 10-39Jioo-Tr00 x3-dx = 1.505582131320634X 10 -18Jso ” T从这里可以看出上面的放缩是很不合理的，选取b值比实际所需值大了太多，造成可后 续计算量过大和计算结果

18、失真。下面取b=50,用求枳公式计算二言心的枳分值，并与准确值总进行比较。.复合梯形求积公式 公式如下:711hb al；(f)=2 / + / Q) + 2 2+ ih)f h = ki=l离散误差为:En(/)= 一含/G)=l (b )/(2)(f),Q b这里 a = 0,b = 50 ,f(x) =且 f(0)=0.(b).复合Simpson求积公式hmm-1%(/)=： IfW + f(b) + 4+ (2i - 1)九)+ 2 W + 2访)i=li=lb a b ah =2 mn离散误差为：Em(f)= - 察-旬= 叫 Q。h-这里 a = 0,b = 50 ,f(x) =.

19、 /一1复 Simpson公式的1天差随区间数m的变化情况1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 L其中分段反间金m=2i2O1.使用Gauss-Legendre求积公式产x3I=L对上面积分,先将区间由0,50变换为-1,1,令x=25*(t+l),则 f1 (f1)3已知而前6个高斯一勒让德求积公式结点和系数为：nXkAkI0.577350269I20.774596669I00.55555555630.86113631203478548540.3399810440.65214515540.9061798460.236926885

20、0.5384693100.47862867050.9324695140.1713244920.66120938603607615730.2386191860.46791393560.9491079120.1294849660.7415311860.2797053910.405845151038183005100.417959184(d).Romberg求积公式】使用Romberg求枳序列进行计算枳分产 x3I(0 = l可得结果如N 1hEQ。116.02734E-166.493939421/27.2333E-066.4939321731/40.1293996096.3645397941/84

21、.2871052812.2068341251/167.409694370.9157549761/326.5384733880.0445339971/646.4928298410.0011095681/1286.4939423372.9344E-0691/2566.4939394074.6484E-09101/5126.4939394028.5727E-12111/10246.4939394021.3323E-14121/20486.4939394024.4409E-15131/40966.4939394023.1086E-14141/81926.4939394022.8422E-14151/1

22、63846.4939394024.1744E44其中九与准确值之间的误差见下图。(e).自适应Simpson公式在区间0,50上进行积分的结果:预先设置误差限近似值实际误差0.10.121312585637262680.016.494986150.00104670.0016.4940241488.475E-050.00016.4939174212.198E-050.000016.4939406141.212E-061E-066.4939396612.584E-071E-076.4939396242.219E-071E-086.4939394031.145E-091E-096.4939394021337E-101E