数值积分上机实验报告

数值积分上机实验报告 桂贤进题一：数学上已经证明了0141+x2dx=成立，所以可以通过数值积分来求的近似值。1.分别使用复合梯形、复合Simpson求积公式计算的近似值。选择不同的h,对每种求积公式，试将误差刻画为h的函数，并比较两方法的精度。是否存在某个值，当低于这个值之后，再继续减小h的值，计算精度不再有所改进，为什么？2.实现Romberg求积方法，并重复上面的计算；3.实现自适应积分方法，并重复上面的计算。解：1.1公式分析：(a) 对于复合梯形公式Tnf=h2fa+fb+2i=1n-1fa+ih ,h=b-an (1)离散误差为：Enf=-nh312f(2)=-h2b-a12f(2),ab (2)(b) 对于复合Simpson公式Smf=h3fa+fb+4i=1mfa+2i-1h+2i=1m-1f(a+2ih) (3)h=b-a2m=b-an离散误差为： Emf=-mh590f(4)=-h4b-a180f(4),ab. (4) 1.2实现算法结果：分别利用梯形公式和Simpson公式计算结果如下：（下表中E1f=-Tf,E2f=-Sf. 此处为MATLAB中的数,可以认为具有足够大的精度）ihiT(f)E1(f)S(f)E2(f)113.00000.141621/23.00000.04163.33330.008341/43.82350.010423.09802.4026e-0561/63.12640.004633.60438.7265e-0781/83.10890.002603.87071.5113e-07101/103.71590.001673.92153.9651e-08121/123.68510.001163.53801.3284e-08201/203.41294.1667e-043.97856.2001e-10301/303.73301.8519e-043.53595.4434e-11401/403.36121.0417e-043.01059.6878e-12501/503.32546.6667e-053.72532.5402e-121001/1003.31291.6667e-053.97533.9968e-142001/2003.31304.1667e-063.97930从上表中可以看出：复合Simpson公式比复合梯形公式精度高，误差收敛的速度快不少。1.3误差下降速度对比：从上图可以看出，复合Simpson公式误差的收敛速度比复合梯形公式的误差的收敛速度快不少，下面验证收敛阶。1.4验证收敛阶：本实验的实际误差主要由离散误差和计算过程中的舍入误差组成，这里离散误差起主导作用，故理论上实际误差的收敛阶应该与离散误差的收敛阶相同。下面利用如下公式来计算实际的收敛阶，并与理论分析所得出的离散误差的收敛阶作比较。err1err2=h1h2p对上表格中所列的区间长度 hi 值，逐次利用相邻两个小区间长 hi 通过上述公式来计算收敛阶,并绘制成图形。得到图形如下：(a) 对复合梯形公式：由上面公式(2)可知，离散误差关于h为二阶收敛，同时由上图可知实验结果的收敛阶将近为2，故与理论分析相符。(b) 对复合Simpson公式：这里却有些奇怪，由上面公式(4)可知，离散误差理论上为4阶收敛，可实验结果却是将近6阶收敛。下面将进一步深入探究。探究如下：考虑这是由于被积函数f(x) 的特殊性导致，而不是由于Simpson公式离散误差真的能达到6阶收敛。由误差的余项公式Emf=-mh590f(4)=-h4b-a180f(4),ab.考虑到f(x)=1.5计算过程中舍入误差的影响从上表格可以看出：当区间数n从1到200时，随着h的减小，实际误差在减小。考虑如下问题：若不断减小h的值，即不断增加区间数n，是否实际误差会一直减小?1.5.1理论分析：我们知道影响实验结果的精确度的因素主要有离散误差和舍入误差 ，而离散误差的大小可以通过离散误差的余项来体现。由公式（2）和（4），可以看出当不断增加区间数，即区间长度h不断减小时，离散误差会越来越小。但相应地由于计算精度的限制，当h不断减小时，舍入误差却会变大。故理论上会存在一个阈值H。当h大于H时，由于离散误差起主导作用，随着区间长度h的减小，实际误差会变小；而当h小于H 时，此时舍入误差将在计算结果的精确度起主导作用，再减小h,会导致计算结果的精确度基本保持不变甚至可能会有减小。1.5.2实验检验：（a）对于复合梯形公式：注意到误差收敛的速度较小，故我们首先选取区间数 n=100，101,102,103108 进行分析，得到下面图形。从图中可以看出，复合梯形公式的阈值H在区间数n为106到108之间取到。下面对n处于区间106到108进行分析。由上图可以看出在n取106 (对应h=10-6)左右h达到阈值，故可取阈值H=10-6.（b） 对于复合Simpson公式：注意到复合Simpson公式得收敛速度较快，取i从0到100（对应区间数m=2i）,逐次计算出对应于m的实际误差，并作图如下。从图中可以看出，在i=42(对应m=84,h=1/168)左右，h达到阈值，故可取阈值H=1/168。通过上述理论分析和实验验证，我们得到如下结果：使用复合梯形公式和复合Simpson公式计算积分值时，所分成的小区间长h都存在阈值H，当hH时，再减小h的值，计算精度不再有所改进。原因如下：当h小于H 时，此时舍入误差将在计算结果的精确度起主导作用，再减小h，虽然离散误差依旧会减小，但舍入误差会增大，这就导致计算结果的精确度基本保持不变甚至可能会有减小。2.1公式分析对Romberg求积方法:T1,1=h12fa+fb, h1=b-aTi,1=12Ti-1,1+hi-1k=12i-1f(a+k-1hi- 1),i=2,3, hi=b-a2i-1=12hi-1Tm,j=4j-1Tm,j-1-Tm-1,j-14j-1-1,j=2,3,(mj)分析离散误差为：对序列Tm,j，m=1,2,3 ， 考虑Em,j=abfxdx-Tm,j, 离散误差的数量级为O(h2j)，h为序列的小区间长。(下表中Ef=-Ti,if , 这里为MATLAB中的数，可以认为是足够精确的)IhiTi,iE(Ti,i)1130.141621/23.33330.008331/43.88235.2499e-0441/83.18746.8698e-0651/163.77171.1688e-0861/323.82444.8451e-1171/643.97237.0166e-1481/1283.9793091/2563.97948.8818e-16101/5123.97928.8818e-16111/10243.97921.3323e-15121/20483.97944.4409e-16131/40963.97921.3323e-15141/81923.97921.3323e-15151/163843.97944.4409e-16 2.2考虑舍入误差的影响：2.2.1理论分析：如同1.5小节中所作分析，此处Romberg积分方法同样存在阈值H。当h小于H 时，此时舍入误差将在计算结果的精确度起主导作用，再减小h，虽然离散误差依旧会减小，但舍入误差会增大，这就导致计算结果的精确度基本保持不变甚至可能会有减小。2.2.2实验验证：取i=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15，逐次计算E(Ti,i),并作图如下：从图中可以看出，当 i=9 (对应h=1/256 ) 时，再减小h的值（即增大i值），误差基本上不变。故可取阈值H=1/256.3.1自适应算法实现：算法思想：通过递归调用实现。（此处感觉课本提供的算法虽然巧妙，但比较冗杂，且会导致占用太多内存。如用递归调用算法直观上更容易理解。）实现算法如下：（自己感觉还算比较巧妙）1.脚本文件zishiying.mclear;a=0;b=1;disp(误差限为：),e=0.h=(b-a)/2;f1=f(a);f2=f(a+b)/2);f3=f(b);tic S=h/3*(f1+f2+4*f3);t=Simpson_auto(a,b,e,S,h/2,f1,f2,f3);disp(近似值为：),tdisp(误差为：),abs(pi-t)disp(耗时为：),toc2.函数文件Simpson_auto.mfunction y =Simpson_auto(A,B,e,S,h,C1,C2,C3) % 将已计算好的函数值传递下去，避免重复计算f1=f(A+h);f2=f(A+3*h);% 利用Simpson公式分别计算左半区间和右半区间的近似值S1=h/3*(C1+C2+4*f1); S2=h/3*(C2+C3+4*f2); % 判断是否达到所需精度，若未达到将区间分半，进行递归调用if abs(S-S1-S2)0其次由上面推导可知：无穷积分 0x3ex-1dx 收敛。考虑用区间截去法，将无穷积分转化为在有限区间上的积分。注意到：bx3ex-1dxbx3x55!dx=b120x2dx=120b120 .这里若取误差容限为：=10-7b=1.2*109但应用mathematica求解：100x3ex-1dx=1.093710-3950x3ex-1dx=1.063410-18从这里可以看出上面的放缩是很不合理的，选取b值比实际所需值大了太多，造成可后续计算量过大和计算结果失真。下面取 b=50，用求积公式计算 050x3ex-1dx 的积分值，并与准确值 415 进行比较。(a).复合梯形求积公式 公式如下：Tnf=h2fa+fb+2i=1n-1fa+ih ,h=b-an 离散误差为： Enf=-nh312f(2)=-h2b-a12f(2),ab 这里 a=0,b=50 ,f(x)=x3ex-1，且f(0)=0.求解结果为：(b).复合Simpson 求积公式Smf=h3fa+fb+4i=1mfa+2i-1h+2i=1m-1f(a+2ih) h=b-a2m=b-an离散误差为： Emf=-mh590f(4)=-h4b-a180f(4),ab. 这里 a=0,b=50 ,f(x)=x3ex-1.(c).使用Gauss-Legendre求积公式If=050x3ex-1dx对上面积分，先将区间由 0,50 变换为-1,1，令x=25*(t+1)，则If=-11254*t+13e25t+1-1dt已知而前6个高斯勒让德求积公式结点和系数为：nxkAk10.120.100.555 555 55630.0.0.0.40.0.0.0.50.0.0.0.0.0.60.0.0.0.0.0.00.(d).Romberg求积公式使用Romberg求积序列进行计算积分If=050x3ex-1dx可得结果如下：ihiTi,iE(Ti,i)116.02734E-166.21/27.2333E-066.31/40.6.41/84.2.51/167.0.61/326.0.71/646.0.81/1286.2.9344E-0691/2566.4.6484E-09101/5126.8.5727E-12111/10246.1.3323E-14121/20486.4.4409E-15131/40966.3.1086E-14141/81926.2.8422E-14151/163846.4.1744E-14其中Ti,i与准确值之间的误差见下图。(e).自适应Simpson公式在区间0,50上进行积分的结果：预先设置误差限近似值实际误差0.10.6.0.016.0.0.0016.8.475E-050.00016.2.198E-050.000016.1.212E-061E-066.2.584E-071E-076.2.219E-071E-086.1.145E-091E-096.1.337E-101E-106.1.077E-111E-116.1.585E-121E-126.6.59E-131E-136.1.776E-141E-146.2.665E-151E-156.1.776E-15结论：从上面分析可以看出，如果不用截断法的话，我们可以直接用Gauss-Leguerre积分公式，且在6个结点的情况下，精度就已经达到了10-6的数量级，考虑到假设结点和积分系数已知，我们所做的只需要计算f(x)在这六个结点的函数值而已。这样看的话，计算量很小就能达到较高的精度。注意到考虑用截断法进行求积分的话，一方面要对被积函数做适当的放大以对于给定的误差限求得一个合适的b值，是很困难的一件事情，我之前尝试的放大被证明是完全不合