实验名称数值积分1

1、探索实验7 数值积分法一、 实验目的 了解求积公式及代数精度概念，理解并掌握求定积分的求积公式的算法构造和计算，学习用计算机求定积分的一些科学计算方法和简单的编程技术和能用程序实现这些算法。二、概念与结论1. 求积公式：计算定积分的如下形式的近似公式：称为求积公式。2代数精度： 若求积公式 对一切不高于m次的 多项都准确成立，而对于m+1次多项式等号不成立，则称此求积公式 的代数精度为m。 代数精度越高，求积公式越好。3.求积余项：4.Newton-Cotes求积公式的代数精度 n点Newton-Cotes求积公式的代数精度至少可以达到n-1,且当n为奇数时，可以达到n。5.Richardso

2、n外推定理： 设函数F1(h)逼近量F*的余项为：F*-F1(h)=a1h p1+a2h p2+····+a k p k+···式中p k>p k-1>···>p2>p1>0, F*和a i (i=1,2, ···)都是与h无关的常数,且k³1时，a k ¹0，则由：定义的函数F2(h)也逼近F*，且有F*-F2(h)= b2h p2+····+b k p k+·&#

3、183;·6. 关于复合梯形公式的展开定理设f(x)在a,b区间上无穷次可微，则有如下展开式：T(h)=I+a1h2+a2h4+a3h6+amh2m+式中T(h)是函数f(x)在a,b区间上的复化梯形值Tn,三、程序中Mathematica语句解释： 1. 随机函数 Random 随机给出闭区间0,1内的一个实数 RandomReal, xmax 随机给出闭区间0,xmax内的一个实数 RandomReal, xmin, xmax 随机给出闭区间xmin,xmax内的一个实数 RandomInteger 随机给出整数0或1 RandomInteger, xmin, xmax 随机给出

4、xmin到xmax之间的一个整数 RandomComplex 随机给出单位正方形内的一个复数2a1,a2,an 表示由元素a1,a2,an组成的一个表，元素可以是任何内容。如：1，3，4，5，1,x,2,3,x+y，1，3，1，2，3，3，2，4等3listk 表list中的第k个元素4listi,j 表list中第i个元素中的第j个元素,此时list中的第i个元素应该也是一个表。 四、方法与程序在实际问题中，往往会遇到被积函数f(x)的原函数无法用初等函数来表示，或函数只能用表格表示，或有的虽然能用初等函数表示，但过分复杂，所以这些情形都需要去建立定积分的近似计算公式来做积分计算。数值积分是

5、进行定积分计算的一种方法，它可以解决不能用定积分基本公式计算的所有定积分问题。数值积分涉及很多计算公式，这里主要介绍Newton-Cotes求积公式、复合求积公式、Romberg求积方法和Monte-Carlo方法的构造过程和算法程序。1. n点 NewtonCotes求积公式 n点 NewtonCotes求积公式又称为等距节点求积公式，它是利用被积函数f(x)在积分区间a,b的n个等分节点上的函数值构造的插值函数j(x)代替f(x)做定积分计算所构造求积公式。这个求积公式是通常做定积分近似计算的梯形公式和抛物线公式的推广，主要在理论上用的多些。1.1 n点 NewtonCotes求积公式的构

6、造过程： 将积分区间a,b 分为n-1等分，其中n个节点 xi=a+(i-1h, i=1,2,n，h=(b-a)/(n-1)，然后用f(x)在这n个节点上建立插值于f(x)的n-1次代数多项式Pn-1(x)，引入变换x=a+th, 0£t£n-1则有 带入定积分，有：Ck(n)称为Cotes(柯特斯)系数, 则得到n点 NewtonCotes求积公式:n点 NewtonCotes求积公式的求积余项为当n=2时，2点的 NewtonCotes求积公式就是如下梯形公式： 梯形求积公式求积余项为当n=3时，3点的NewtonCotes求积公式就是如下抛物线（Simpson）公式：

7、 Simpson求积公式求积余项为如果想得到其他的 NewtonCotes求积公式只要在有关书中查出Cotes系数表就可以马上得到相应的NewtonCotes求积公式。1.2 n点 NewtonCotes求积公式算法： 1 输入被积函数f(x)及积分上下限a,b2 选择Cotes系数构造求积公式3 用求积公式求定积分1.3 n点 NewtonCotes求积公式程序：Cleara,b,x,n,s;a=Input"a="b=Input"b="fx_=Input"被积函数f(x)="n= Input"求积节点个数n="c

8、=1/2,1/2,1/6,4/6,1/6,1/8,3/8,3/8,1/8,7/90,16/45,2/15,16/45,7/90,19/288,25/96,25/144,25/144,25/96,19/288,41/840,9/35,9/280,34/105,9/280,9/35,41/840,751/17280,3577/17280,1323/17280,2989/17280,1323/17280,3577/17280,751/17280,989/28350,5888/28350,-928/28350,10496/28350,-4540/28350,10496/28350,-928/28350

9、,5888/28350,989/28350;h=(b-a)/(n-1);x=Tablea+k*h,k,0,n-1;s=(b-a)*Sumcn-1,k*fxk,k,1,nPrint"定积分=",Ns,8;说明：本程序用n（n=2,3,4,5,6，7，8，9）点 NewtonCotes求积公式求a,b上的定积分近似值。程序执行后，按要求通过键盘输入积分下限a、积分上限b、被积函数f(x)和求积节点个数n后，计算机则给出定积分的近似值。程序中变量说明：a：存放积分下限b: 存放积分上限fx: 存放被积函数f(x)n: 存放求积节点个数c: 存放Cotes系数s: 存放定积分近似值

10、h: 存放节点步长x：存放节点xi注：语句c=1/2,1/2,1/6,4/6,1/6,1/8,3/8,3/8,1/8,7/90.16/45,2/15,16/45,7/90, 19/288,25/96,25/144,25/144,25/96,19/288的第i个分量表是具有i个节点的Cotes系数。1.4 例题与实验例1. 用n=3和n=4的Newton-Cotes求积公式求定积分的近似值。解：执行n点 NewtonCotes求积公式程后，在输入的窗口中按提示分别输入1、3、Exp-x/2、3，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出用n=3的Newton-Cotes求积公式计

11、算出的定积分结果： 1 2 12 (- + - + - ) 6E3/2 3 E 6 SqrtE定积分=0.76705953 再执行n点 NewtonCotes求积公式程后，在输入的窗口中按提示分别输入1、3、Exp-x/2、4，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出用n=4的Newton-Cotes求积公式计算出的定积分结果： 1 3 3 12 (- + - + - + -) 8E3/2 8E7/6 8E5/6 8 SqrtE 定积分=0.76691628因此用n=3和n=4的Newton-Cotes求积公式求本题定积分近似值分别为0.76705953和0.76691628

12、注意到本题的精确值为0.766800999.,可见n=4的Newton-Cotes求积公式计算结果较好。2. 复化求积公式 复化求积公式是把积分区间分成若干小区间，在每个小区间上采用次数不高的插值公式，如梯形公式或抛物线公式，构造出相应的求积公式，然后再把它们加起来得到整个区间上的求积公式。复化求积公式克服了高次Newton-Cotes公式计算不稳定的问题，其运算简单且易于在计算机上实现。常用的复化求积公式是复化梯形公式和复化抛物线公式，下面分别讨论。2.1复化梯形公式的构造过程: 把区间a,b n等分，取节点xk=a+kh,k=0,1,.n, h=(b-a)/n,对每个小区间xk,xk+1用

13、梯形求积公式，再累加起来得：公式就是复化梯形公式。它的求积余项为： 由复化梯形公式的余项，可以看到，当n不断变大时， Tn无限接近定积分值。因此复化梯形公式可以使定积分的计算达到任意精度。为了计算简单，提高效率，常用|T2n Tn|<e来得到满足精度要求的定积分值。2.2复化抛物线公式的构造过程在每个小区间xk,xk+1上，有 把区间a,b n等分，取节点xk=a+kh,k=0,1,.n, h=(b-a)/n,对每个小区间xk,xk+1 用抛物线求积公式，再累加起来得：公式就是复化抛物线公式。它的求积余项为： 由复化抛物线公式的余项，可以看到，当n不断变大时，Sn无限接近定积分值。因此复

14、化抛物线公式也可以使定积分的计算达到任意精度，且收敛的速度比复化梯形公式更快。为了计算简单，提高效率，常用| S2n Sn|<e来得到满足精度要求的定积分值。2.3 复化梯形求积公式算法：1.输入被积函数f(x),积分上下限a,b,和求积精度e2. nÜ1, 计算Tn3. 计算T2n4 判断|T2n Tn|<e是否成立，如果成立，输出定积分近似值，停止5 否则， Tn Ü T2n , nÜ2n,转32.4 复化抛物线求积公式算法：1.输入被积函数f(x),积分上下限a,b,和求积精度e2. nÜ1, 计算Sn3.计算S2n4.判断|S2n S

15、n|<e是否成立，如果成立，输出定积分近似值，停止5.否则， Sn Ü S2n , nÜ2n,转32.5 复化梯形求积公式程序：(*复合梯形公式*)Cleara,b,x,n,f;a=Input"a="b=Input"b="fx_=Input"被积函数f(x)="eps=Input"精度要求eps="n=1;h=b-a;t1=(fa+fb)*h/2;h=h/2;t2=t1/2+h*fa+h;er=t2-t1/N;WhileAbser>eps,Print"n=",2n

16、,"定积分值为",Nt2,10;Print"误差=",er;h=h/2;t1=t2;n=n+1;t2=t1/2+h*Sumfa+k*h,k,1,2n,2;er=t2-t1/N;Print"n=",2n,"定积分值为",Nt2,10;Print"误差=",er说明：本程序从梯形公式T1开始，用复合梯形求积公式求a,b上满足精度小于eps要求的定积分近似值。程序执行后，按要求通过键盘输入积分下限a、积分上限b、被积函数f(x)和精度要求eps 后，计算机则给出满足精度要求的定积分近似值及中间计算值和

17、误差。程序中变量说明：a：存放积分下限b: 存放积分上限fx: 存放被积函数f(x)eps: 存放求积精度eh: 存放节点步长x：为函数fx:提供变量t1: 存放复合梯形值Tnt2: 存放复合梯形值T2n和定积分近似值n: 存放复合梯形公式的节点加密次数er:存放误差注：为提高计算效率，计算T2n时使用了公式2.6 复化抛物线求积公式程序：(*复合抛物线公式*）Cleara,b,x,n,f,s;a=Input"a="b=Input"b="fx_=Input"被积函数f(x)="n=Input"分割区间数n="h=(

18、b-a)/n;a1=a+h/2;s=h/6*(fa+fb+2Sumfa+k*h,k,1,n-1+4Sumfa1+k*h,k,0,n-1);Print"n=",n;Print"定积分值为",Ns,10说明：本程序用复合抛物线求积公式求a,b上定积分近似值。程序执行后，按要求通过键盘输入积分下限a、积分上限b、被积函数f(x)和分割区间数n 后，计算机则给出定积分的近似值Sn。程序中变量说明：a：存放积分下限b: 存放积分上限fx: 存放被积函数f(x)h: 存放节点步长x：为函数fx提供变量n: 存放分割区间数s: 存放复合抛物线值Sn和定积分近似值注：M

19、athematica数学软件中也有一个求数值积分的命令，命令形式为：NIntegratefx,x,a,b 它可以求出a,b上的定积分近似值 2.7例题与实验例2. 用复合梯形公式求定积分的近似值,要求误差小于10-7。解：执行复化梯形求积公式程后，在输入的窗口中按提示分别输入：0、3、Exp-x2、10(-7)，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出的计算结果如下：n=2 定积分值为0.7313702518 误差=0.0474305n=4 定积分值为0.7429840978 误差=0.0116138n=8 定积分值为0.7458656148 误差=0.00288152n=1

20、6 定积分值为0.7465845968 误差=0.000718982n=32 定积分值为0.7467642547 误差=0.000179658n=64 定积分值为0.7468091636 误差=0.000044909n=128 定积分值为0.7468203905 误差=0.0000112269n=256 定积分值为0.7468231972 误差=2.8067´10-6n=512 定积分值为0.7468238989 误差=7.01676´10-7n=1024 定积分值为0.7468240743 误差=1.75419´10-7n=2048 定积分值为0.7468241

21、182 误差=4.38548´10-8从计算结果可以得出经过11次加密分割积分区间得到满足精度要求的定积分近似值0.7468241182。此外，通过计算过程可以明显看到复合梯形求积公式的收敛性。例3. 分别用n=1,2,4,8,16,32的复合抛物线求积公式求定积分的近似值。解：执行复化抛物线求积公式程后，在输入的窗口中按提示分别输入：0、Pi、Expx*Cosx、1，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出的计算结果如下：n=1定积分值为-11.59283955再重复如上操作，将分割区间数分别输入2，4，8，16，32，得如下输出结果：n=2定积分值为-11.98

22、494402n=4定积分值为-12.06420896n=8定积分值为-12.06995132n=16定积分值为-12.07032146n=32定积分值为-12.07034476。注意到本题的定积分值为-12.07034631639，从中可以看到复合抛物线求积公式收敛还是很快的。你能通过修改复合抛物线求积公式一次得出如上所有计算结果吗？3.Romberg求积公式 Romberg又称逐次分半加速收敛法，它是对复合梯形求积公式采用Richardson加速技术得到的一种数值求积方法。3.1 Romberg求积公式的构造过程: 由复合梯形公式的展开定理，有如下关系式：T1(h)-I=a1h2+a2h4+

23、a3h6+amh2m+这里利用Richardson外推定理对T1(h)进行加速，注意到这里p1=2,取q=1/2有于是得到收敛于I的加速公式T2(h),如果再T2(h) 进行加速，可以有继续加速下去，可以得到加速序列式中m代表对I进行的加速次数。 此外，注意到T1(h/2)=T2n ，且直接计算可以知道T2(h)就是复合抛物线公式Sn,于是有计算关系式类似地，还有如下计算格式： 和于是我们至少可以得到如下3个收敛于定积分I的数列：T型值数列：T1、T2、T4、T8 、T16 、T2n 、Simpson数列：S1、S2、S4、S8 、S16 、S2n 、Cotes数列：C1、C2、C4、C8 、

24、C16 、C2n 、Romberg数列：R1、R2、R4、R8 、R16 、R2n 、利用Romberg数列求定积分I的方法称为Romberg求积方法。 如果给定求积精度e，可以用作为终止计算的条件，并取最后算出的R值作为积分值。3.2 Romberg求积算法：1.输入被积函数f(x),积分上下限a,b,和求积精度e2.按顺序计算T1、T2、T4、T8 并做 t1Ü T1、t2Ü T2、t3Ü T4、t4Ü T8、 s1Ü (4t2-t1)/3,s2Ü (4t3-t2)/3,s3Ü (4t4-t3)/3;c1Ü (

25、16s2-s1)/15,c2Ü (16s3-s2)/15,R2Ü (64c2-c1)/63;3.nÜ1,R1Ü c2 , t1Üt4,s1Üs3 ,c1Üc24.判断|R2 R1|<e是否成立，如果成立，输出定积分近似值R2，停止5.否则， R1Ü R2 , t2Ü T16n, s2Ü (4t2-t1)/3, t1Üt2, c2Ü (16s2-s1)/15, s1Üs2 , R2Ü (64c2-c1)/63, c1Üc2 nÜ2n

26、,转43.3 Romberg求积方法程序： Cleara,b,x,n,f; a=Input"a=" b=Input"b=" fx_=Input"被积函数f(x)=" eps=Input"精度要求eps="m=1;h=b-a;t=Table0,4;t1=(fa+fb)*h/2; Dom=2m;h=h/2;tk+1=tk/2+h*Sumfa+i*h,i,1,m,2/N ,k,1,n-1s1=(4t2-t1)/3;s2=(4t3-t2)/3;s3=(4t4-t3)/3;c1=(16s2-s1)/15;c2=(16s3-s

27、2)/15;r1=c2;r2=(64c2-c1)/63;t1=t4;s1=s3;c1=c2;er=r2-r1;k=1;WhileAbser>eps,r1=r2;Print"r(",k,")=",Nr2,20," eps=",Ner,10;h=h/2;m=2m;t2=t1/2+h*Sumfa+i*h,i,1,m,2/N;s2=(4t2-t1)/3;c2=(16s2-s1)/15;r2=(64c2-c1)/63;t1=t2;s1=s2;er=r2-r1;k=k+1;c1=c2;Print"r(",k,"

28、)=",Nr2,20," eps=",Ner,10;说明：本程序从用Romberg求积方法求a,b上满足精度小于eps要求的定积分近似值。程序执行后，按要求通过键盘输入积分下限a、积分上限b、被积函数f(x)和精度要求eps 后，计算机则给出满足精度要求的定积分近似值及中间计算值和误差。程序中变量说明：a：存放积分下限b: 存放积分上限fx: 存放被积函数f(x)eps: 存放求积精度eh: 存放节点步长x：为函数fx:提供变量t,t1 ,t2: 存放复合梯形值Tns1 ,s2: 存放Simpson数列值Snc1 ,c2: 存放Cotes数列值CnR1 ,R2:

29、存放Romberg数列值Rnm: 存放复合梯形公式的节点加密次数er:存放误差注：为编程方便，程序中先算出了产生Romberg数列的值：T1、T2、T4、T8 、S1、S2、S4、C1、C2使，其中t=Table0,4提供一个4元素的表用于存放T1、T2、T4、T8。 3.4例题与实验例4.用Romberg求积方法求定积分的近似值,要求误差小于10-12。解：执行Romberg求积方法程序后，在输入的窗口中按提示分别输入：0、1、4/(1+x2）、10(-12)，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出的计算结果如下：r(1)=3.141585783761874 eps=-8

30、.310364015´10-6r(2)=3.141592638396796 eps=6.854634922´10-6r(3)=3.141592653590029 eps=1.519323334´10-8r(4)=3.141592653589794 eps=-2.358113704´10-13因此得到满足精度要求的定积分值为3.141592653589794 ，误差=-2.358113704´10-134. Monte-Carlo求积方法 Monte-Carlo求积方法以概率论大数定理为依据式，借助随机数来求定积分的近似值的一种方法，该方法特别适

31、用于求多重积分和积分区域复杂的重积分，不足之处是收敛较慢。4.1 Monte-Carlo求积方法的构造过程: 假设重积分是有限的，式中WÍRn，G(x)是定义在W上的多元函数, P(x)是定义在W上的分布函数，则有只要式中x(k)是W上关于分布P(x)的n个任意独立选取的随机点列即可，这个结论可以有概率论中的大数定律： 以概率1成立。 利用以上结果求重积分的方法就称为Monte-Carlo求积方法。对具体的积分计算有如下计算公式：1） 求定积分的计算公式式中xk是a,b均匀分布的n个任意独立选取的随机数列。2） 求二重积分的计算公式式中(xk,yk)是 平面区域D中均匀分布的n个任意

32、独立选取的随机点列，|D|表示区域D的面积。3） 求三重积分的计算公式式中(xk,yk,zk)是空间区域W中均匀分布的n个任意独立选取的随机点列，|W|表示区域W的体积。 Monte-Carlo求积方法的收敛阶为O(n1/2).4.2 Monte-Carlo求积方法算法：1.输入被积函数f(x),积分区域边界和随机点的个数n2.在积分区域中产生n个随机点列3.计算被积函数f(x)在随机点列上的函数值并相加求平均4.用平均值与积分区域的测度之积作为积分近似值，这里测度就是长度、面积、体积之类的度量。4.3 Monte-Carlo求积方法程序： （*求定积分*） Cleara,b,x,n,f; a

33、=Input"a=" b=Input"b=" fx_=Input"被积函数f(x)=" n=Input"随机点个数n="s=(b-a)*SumfRandomReal,a,b,n/nPrint"n=",n," 积分»",Ns,10;说明：本程序从用Monte-Carlo求积方法求a,b上的定积分近似值。程序执行后，按要求通过键盘输入积分下限a、积分上限b、被积函数f(x)和随机点个数n 后，计算机则给出定积分近似值。程序中变量说明：a：存放积分下限b: 存放积分上限f

34、x: 存放被积函数f(x)n: 存放随机点个数n（*求二重积分*） Cleara,b,x,y,n,f; a=Input"a="b=Input"b=" c=Input"c=" d=Input"d=" fx_,y_=Input"被积函数f(x,y)=" n=Input"随机点个数n="s=(b-a)*(d-c)*SumfRandomReal,a,b, RandomReal,c,d,n/nPrint"n=",n," 积分»",Ns,

35、10;说明：本程序从用Monte-Carlo求积方法求a,b´c,d上的定积分近似值。程序执行后，按要求通过键盘输入积分变量x的下限a、积分上限b和积分变量y的下限c、积分上限d、被积函数f(x)和随机点个数n 后，计算机则给出定积分近似值。程序中变量说明：a：存放积分变量x下限b: 存放积分变量x上限c：存放积分变量y下限d: 存放积分变量y上限fx: 存放被积函数f(x)n: 存放随机点个数n4.4例题与实验例5.用Monte-Carlo求积方法求定积分的近似值,随机点个数取为200。解：执行Monte-Carlo求积方法求定积分程序后，在输入的窗口中按提示分别输入：0、1、Sinx/x、200后，每次输入后用鼠标点击窗口的“OK”按扭，计算机在屏幕上给出的计算结果如下：n=200 积分=0.9463437539本题积分准确值为0.94608307036.例6.用Monte-Carlo求积方法求定积分的近似值,用随机点个数分别取为50，60，70，80，200的计算值进行计算实验来观察计算结果特点并提出一种从一组Monte-Carlo求积计算值推断的较好积分近似值的方法。解：修改Monte-Carlo求积方法求二重积分程序为：fx_,y_:=Exp-(x+y);Dos=SumfRandomReal,1,2,y=RandomReal,2