(精品)第5章数值积分.ppt

《(精品)第5章数值积分.ppt》由会员分享，可在线阅读，更多相关《(精品)第5章数值积分.ppt（85页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第第5章章 数值积分数值积分 5.1 引言引言 利用牛顿利用牛顿莱布尼兹（莱布尼兹（NewtonLeibniz）公式公式 (5.1)解决函数解决函数 在在 上的积分问题在理论和应用上都有重大的意上的积分问题在理论和应用上都有重大的意义。然而，在实际问题中，往往会遇到一些困难。有些形式上较简单的义。然而，在实际问题中，往往会遇到一些困难。有些形式上较简单的函数，其原函数函数，其原函数 不易求出或不能用初等函数表示成有限形式；有些不易求出或不能用初等函数表示成有限形式；有些被积函数的原函数过于复杂；而有些函数的函数值是由实验、观测等方被积函数的原函数过于复杂；而有些函数的函数值是由实验、观测等方法

2、得出，并没有给出具体的解析表达式。这些情形说明公式法得出，并没有给出具体的解析表达式。这些情形说明公式(5.1)在应用在应用上是有局限性的，因此研究定积分的数值计算问题就显得十分必要。上是有局限性的，因此研究定积分的数值计算问题就显得十分必要。本章主要介绍一些常用的数值积分方法，包本章主要介绍一些常用的数值积分方法，包括梯形积分法、辛卜生括梯形积分法、辛卜生积分法、变步长积分法、牛顿积分法、变步长积分法、牛顿柯特斯积分法、高斯积分法、龙贝格积柯特斯积分法、高斯积分法、龙贝格积分法以及分法以及高振荡函数的积分法高振荡函数的积分法。5.2 5.2 梯形积分法梯形积分法 5.2.1 5.2.1 梯形

3、积分法的基本思想梯形积分法的基本思想 梯形积分法的基本思想：在积分区间梯形积分法的基本思想：在积分区间 上，根上，根据给定的插值条件据给定的插值条件 和和 ，构造，构造一个一次二项式一个一次二项式 ，并以，并以 的积分值近似地的积分值近似地代替代替 。从几何角度而言，是以梯形面积近。从几何角度而言，是以梯形面积近似地代替曲边梯形的面积。似地代替曲边梯形的面积。5.2.2 5.2.2 梯形求积公式梯形求积公式 依据梯形积分法的基本思想，将区间 分成 个 相 等 的 小 区 间,则 每 个 小 区 间 的 长 度 为 ,对每个小区间均实施如下的梯形求积：将这些小梯形的求积值加起来，可以得到如下梯形

4、求积公式：其中，5.2.3 5.2.3 实现梯形积分法的基本步骤实现梯形积分法的基本步骤 (1)(1)输入区间输入区间 的端点的端点 值以及分割数值以及分割数 ；(2)(2)将区间将区间 等分成等分成 个小区间，每一个小区间的个小区间，每一个小区间的 长度长度 ；(3)(3)计算每一个等分点的函数值计算每一个等分点的函数值 (4)(4)计算计算 (5)(5)输出输出 的值；的值；(6)(6)结束。结束。图图 5.2 5.2 梯形积分法的梯形积分法的N-SN-S图描述图描述 例例5.1 5.1 使用梯形求积公式求下列定积分的值。使用梯形求积公式求下列定积分的值。#define N 16 /*#d

5、efine N 16 /*等分数等分数*/*/float float funcfunc(float x)(float x)float y;float y;y=4.0/(1+x*x);y=4.0/(1+x*x);return(y);return(y);void void gedianzhigedianzhi(float y,float a,float h)(float y,float a,float h)intint i;i;for(i=0;i=N;i+)for(i=0;i=N;i+)yi=yi=funcfunc(a+i*h);(a+i*h);float trapeze(float y,floa

6、t h)float trapeze(float y,float h)float s;float s;intint i;i;s=(y0+yN)/2.0;s=(y0+yN)/2.0;for(i=1;iN;i+)for(i=1;iN;i+)s+=yi;s+=yi;return(s*h);return(s*h);main()main()float a,b,h,s,fN+1;float a,b,h,s,fN+1;clrscrclrscr();();printfprintf(input a,b=);(input a,b=);scanfscanf(%f,%f,&a,&b);(%f,%f,&a,&b);h=(

7、b-a)/(float)N;h=(b-a)/(float)N;gedianzhigedianzhi(f,a,h);(f,a,h);s=trapeze(f,h);s=trapeze(f,h);printfprintf(s=%fn,s);(s=%fn,s);程序运行结果：程序运行结果：input a,b=0,1input a,b=0,1s=s=3.1409423.140942 辛卜生积分法的基本思想：在积分区间辛卜生积分法的基本思想：在积分区间 上，根据上，根据给定的插值条件给定的插值条件 、和和 ，构造，构造个二次插值求积多项式个二次插值求积多项式 ，并以，并以 的的积分值近似地代替积分值近似地

8、代替 。从几何角度而言，是用过三。从几何角度而言，是用过三点的抛物线面积近似地代替积分的曲边面积。点的抛物线面积近似地代替积分的曲边面积。5.3 5.3 辛卜生（辛卜生（Simpson）积分法积分法 5.3.1 5.3.1 辛卜生积分法的基本思想辛卜生积分法的基本思想5.3.2 5.3.2 辛卜生求积公式辛卜生求积公式 依依据据辛辛卜卜生生积积分分法法的的基基本本思思想想，将将区区间间 分分成成 (必必须须是是 偶偶 数数）个个 相相 等等 的的 小小 区区 间间，则则 每每 个个 小小 区区 间间 的的 长长 度度 为为 ，在小区在小区间间 均均实实施如下的辛卜生求施如下的辛卜生求积积：将这

9、些求积值加起来，可以得到如下辛卜生求积公式：将这些求积值加起来，可以得到如下辛卜生求积公式：其中其中:为寄数项的函数为寄数项的函数 值之和。值之和。为偶数项的函数为偶数项的函数 值之和。值之和。5.3.3 5.3.3 实现辛卜生积分法的基本步骤实现辛卜生积分法的基本步骤(1)(1)输入区间输入区间 的端点的的端点的 值以及分割数值以及分割数 ；(2)(2)将区间将区间 等分成等分成 个小区间，每一个小区间的长度个小区间，每一个小区间的长度 ；(3)(3)计算每一个等分点的函数值计算每一个等分点的函数值 ；(4)(4)计算计算:（计算奇数项的函数值之和）（计算奇数项的函数值之和）（计算偶数项的函

10、数值之和）（计算偶数项的函数值之和）(5)(5)计算计算 ；(6)(6)输出输出 的值；的值；(7)(7)结束。结束。图图 5.4 5.4 辛卜生积分法的辛卜生积分法的N-SN-S图描述图描述例例5.25.2 使用辛卜生求积公式求下列定积分的值。使用辛卜生求积公式求下列定积分的值。#include include.h#define N 16 /*#define N 16 /*等分数等分数*/*/float float funcfunc(float x)(float x)float y;float y;y=4.0/(1+x*x);y=4.0/(1+x*x);return(y);return(y)

11、;void void gedianzhigedianzhi(float y,float a,float h)(float y,float a,float h)intint i;i;for(i=0;i=N;i+)for(i=0;i=N;i+)yi=yi=funcfunc(a+i*h);(a+i*h);float float simpsonsimpson(float y,float h)(float y,float h)float s,s1,s2;float s,s1,s2;intint i;i;s1=y1;s1=y1;s2=0.0;s2=0.0;for(i=2;i=N-2;i=i+2)for(i

12、=2;i=while(p=epseps)/*)/*判断是否达到精度要求判断是否达到精度要求,若没有达到若没有达到,继续循环继续循环*/*/s=0.0;s=0.0;for(i=0;i=n-1;i+)for(i=0;i=n-1;i+)x=a+(i+0.5)*h;x=a+(i+0.5)*h;s=s+s=s+funcfunc(x);(x);t2=(t1+h*s)/2.0;/*t2=(t1+h*s)/2.0;/*计算计算*/*/p=p=fabsfabs(t1-t2);/*(t1-t2);/*计算精度计算精度*/*/t1=t2;t1=t2;n=n+n;n=n+n;h=h/2.0;h=h/2.0;retur

13、n(t2);return(t2);void main()void main()float a,b,t;float a,b,t;clrscrclrscr();();printfprintf(input a,b=);(input a,b=);scanfscanf(%f,%f,&a,&b);(%f,%f,&a,&b);t=t=btrapezebtrapeze(a,b);(a,b);printfprintf(t=%fn,t);(t=%fn,t);程序运行结果：程序运行结果：input a,b=0,1input a,b=0,1t=0.746824t=0.746824 5.4.4 5.4.4 变步长辛卜生

14、求积分法变步长辛卜生求积分法变步长辛卜生求积分法的基本过程变步长辛卜生求积分法的基本过程:（1 1）利用梯形公式，将积分区间）利用梯形公式，将积分区间 一等分，一等分，（2 2）将其中的每一个求积小区间再二等分一次）将其中的每一个求积小区间再二等分一次（3 3）根据上面两式）根据上面两式 和和 的余项的余项 、，可以，可以 推导出如下的变步长辛卜生求积公式。推导出如下的变步长辛卜生求积公式。进一步得到再二次等分一次后的变步长辛卜生求积公式为进一步得到再二次等分一次后的变步长辛卜生求积公式为（4 4）若）若 ，二等分后的积分值，二等分后的积分值 就是最后的结果；否则保存当就是最后的结果；否则保存

15、当前的变步长梯形积分值、等分数、积分值与步长，转到第（前的变步长梯形积分值、等分数、积分值与步长，转到第（2 2）步继续做）步继续做二等分处理。二等分处理。5.4.5 5.4.5 实现变步长辛卜生积分法的基本步骤实现变步长辛卜生积分法的基本步骤变变步步长长梯梯形形积积分分法法的的N N-S S图图描描述述例例5.4 5.4 使用变步长辛卜生求积分法求下列定积分的值。使用变步长辛卜生求积分法求下列定积分的值。#include include.h#include#include#define#define epseps 0.000001 /*0.000001 /*容许误差容许误差*/*/float

16、 float funcfunc(float x)(float x)float y;float y;y=y=sqrtsqrt(1-x*x);(1-x*x);return(y);return(y);float float bsimpsonbsimpson(float a,float b)(float a,float b)intint i,n;i,n;float h,p,e,s;float h,p,e,s;float t1,t2,s1,s2,x;float t1,t2,s1,s2,x;n=1;n=1;h=b-a;h=b-a;t1=h*(t1=h*(funcfunc(a)+(a)+funcfunc(b

17、)/2.0;(b)/2.0;s1=t1;/*s1=t1;/*用代替用代替*/*/e=e=epseps+1.0;+1.0;while(e=while(e=epseps)s=0.0;s=0.0;for(i=0;i=n-1;i+)for(i=0;i=n-1;i+)x=a+(i+0.5)*h;x=a+(i+0.5)*h;s=s+s=s+funcfunc(x);(x);t2=(t1+h*s)/2.0;/*t2=(t1+h*s)/2.0;/*计算计算*/*/s2=(4*t2-t1)/3.0;/*s2=(4*t2-t1)/3.0;/*计算计算*/*/e=e=fabsfabs(s2-s1);/*(s2-s1)

18、;/*计算精度计算精度*/*/t1=t2;t1=t2;s1=s2;s1=s2;n=n+n;n=n+n;h=h/2.0;h=h/2.0;return(s2);return(s2);void main()void main()float a,b,s;float a,b,s;clrscrclrscr();();printfprintf(input a,b=);(input a,b=);scanfscanf(%f,%f,&a,&b);(%f,%f,&a,&b);s=s=bsimpsonbsimpson(a,b);(a,b);printfprintf(s=%fn,s);(s=%fn,s);程序运行结果：

19、程序运行结果：input a,b=0,1input a,b=0,1s=0.785398s=0.7853985.5 5.5 牛顿牛顿柯特斯柯特斯(NewtonCotes)积分法积分法 5.5.1 5.5.1 牛顿牛顿柯特斯积分法的基本思想柯特斯积分法的基本思想 牛牛顿顿柯柯特特斯斯积积分分法法的的基基本本思思想想：用用高高次次的的插插值值求求积积多多项项式式 去去逼逼近近被被积积函函数数 ，以获得高精度的积分值。，以获得高精度的积分值。事事实实上上，梯梯形形积积分分是是当当 时时的的牛牛顿顿柯柯特特斯斯积积分分，辛辛卜卜生生积积分分是是当当 时时的的牛牛顿顿柯柯特斯积分，它们都是牛顿特斯积分，它

20、们都是牛顿柯特斯积分的特例。柯特斯积分的特例。5.5.2 5.5.2 牛顿牛顿柯特斯求积公式柯特斯求积公式 下面给出三到五阶牛顿下面给出三到五阶牛顿柯特斯求积公式。柯特斯求积公式。实现三阶牛顿实现三阶牛顿柯特斯求积公式的基本步骤如下：柯特斯求积公式的基本步骤如下：三阶牛顿三阶牛顿柯特斯求积公式的柯特斯求积公式的N-SN-S图描述图描述例例5.5 5.5 使用牛顿使用牛顿柯特斯求积公式求下列定积分的值。柯特斯求积公式求下列定积分的值。#include include.h#define N 5#define N 5float float funcfunc(float x)(float x)floa

21、t y;float y;y=4.0/(1+x*x);y=4.0/(1+x*x);return(y);return(y);void void gedianzhigedianzhi(float y,float a,float b,(float y,float a,float b,intint n)n)float h,s;float h,s;intint i;i;h=(b-a)/(float)n;h=(b-a)/(float)n;for(i=0;i=n;i+)for(i=0;i=n;i+)yi=yi=funcfunc(a+i*h);(a+i*h);float nc3(float y,float a,

22、float b,float nc3(float y,float a,float b,intint n)n)float h,s,s0,s1;float h,s,s0,s1;intint i;i;h=(b-a)/(float)n;h=(b-a)/(float)n;s0=s1=0.0;s0=s1=0.0;for(i=0;i=n-3;i=i+3)for(i=0;i=n-3;i=i+3)s0+=yi+yi+3;s0+=yi+yi+3;s1+=yi+1+yi+2;s1+=yi+1+yi+2;s=s0+3.0*s1;s=s0+3.0*s1;return(s*h*3.0/8.0);return(s*h*3.0

23、/8.0);main()main()float a,b,h,s,fN+1;float a,b,h,s,fN+1;float n3,n4,n5;float n3,n4,n5;printfprintf(input a,b=);(input a,b=);scanfscanf(%f,%f,&a,&b);(%f,%f,&a,&b);gedianzhigedianzhi(f,a,b,3);(f,a,b,3);n3=nc3(f,a,b,3);n3=nc3(f,a,b,3);printfprintf(n 3-(n 3-nc nc 4-4-nc nc 5-5-ncncn);n);printfprintf(%f%

24、f%fn,n3,n4,n5);(%f%f%fn,n3,n4,n5);程序运行结果：程序运行结果：input a,b=0,1input a,b=0,1 3-3-nc nc 4-4-nc nc 5-5-ncnc 3.138461 3.142118 3.1418783.138461 3.142118 3.1418785.6 5.6 高斯积分法高斯积分法 5.6.1 5.6.1 高斯积分法的基本思想高斯积分法的基本思想 前前面面介介绍绍的的几几种种数数值值积积分分方方法法，都都是是先先寻寻找找一一个个插插值值求求积积多多项项式式 ，并以，并以 近似代替函数近似代替函数 进行积分而获得积分的近似值，即进

25、行积分而获得积分的近似值，即 (5.2)(5.2)表明，函数表明，函数 在区间在区间 上的积分可以用函数上的积分可以用函数 在该区间上在该区间上 的的 个点的函数值的线性组合来近似代替。但是，由于插值求积公式个点的函数值的线性组合来近似代替。但是，由于插值求积公式是利用插值多项式的积分得到的，因此，如果被积函数是利用插值多项式的积分得到的，因此，如果被积函数 为次数不超为次数不超过过 的多项式，则利用插值求积公式计算得到的积分值是准确的。的多项式，则利用插值求积公式计算得到的积分值是准确的。在实际应用中，为了提高数值求积公式的精度，一般要求数值求积在实际应用中，为了提高数值求积公式的精度，一般

26、要求数值求积公式对于次数尽可能高的多项式能准确成立。由此提出了高斯积分法。公式对于次数尽可能高的多项式能准确成立。由此提出了高斯积分法。即：如果插值求积公式即：如果插值求积公式(5.2)(5.2)具有具有 次代数精度，那么称该插值求次代数精度，那么称该插值求积公式为高斯求积公式。其节点积公式为高斯求积公式。其节点 称为高斯点，称为高斯点，称为高斯求积系数。称为高斯求积系数。下面具体介绍几个常用的高斯求积公式。下面具体介绍几个常用的高斯求积公式。5.6.2 5.6.2 勒让德勒让德高斯高斯(LegendreLegendre-Gauss-Gauss)求积公式求积公式 勒勒让让德德高高斯斯求求积积公

27、公式式，特特别别适适合合于于计计算算区区间间-1,1-1,1的积分，其求积公式可以表示为以下形式：的积分，其求积公式可以表示为以下形式：而而对对于于一一般般区区间间 ，通通过过变变换换可可以以得得到到以以下下形形式式的的求积公式：求积公式：高高斯斯点点 及及高高斯斯求求积积系系数数 ，参参见见表表5.15.1（表表示示阶阶数）。数）。例例5.6 5.6 使使用用勒勒让让德德高高斯斯求求积积公公式式求求下下列列定定积分的值。积分的值。#include#include double double funcfunc(double x)(double x)double y;double y;y=x*x

28、+sin(x);y=x*x+sin(x);return(y);return(y);double double legendrelegendre_gauss(double a,double b,_gauss(double a,double b,intint m,m,intint n)n)double h,double h,hxhx,y,s,y,s,dxdx,x0;,x0;intint i,k;i,k;static double x=-0.9061798459,-0.5384693101,0.0000000000,static double x=-0.9061798459,-0.538469310

29、1,0.0000000000,0.5384693101,0.9061798459;0.5384693101,0.9061798459;static double w=0.2369268851,0.4786286705,0.5688888889,static double w=0.2369268851,0.4786286705,0.5688888889,0.4786286705,0.2369268851;0.4786286705,0.2369268851;dxdx=(b-a)/(double)m;=(b-a)/(double)m;hxhx=dxdx/2;/2;s=0.0;s=0.0;for(k=

30、0;km;k+)for(k=0;km;k+)x0=a+(double)k*x0=a+(double)k*dxdx+hxhx;for(i=0;in;i+)for(i=0;in;i+)y=x0+xi*y=x0+xi*hxhx;s=s+wi*s=s+wi*funcfunc(y);(y);return(s*return(s*hxhx););main()main()double a,b,s;double a,b,s;intint m,n;m,n;clrscrclrscr();();printfprintf(input(input duandianduandian a,b=);a,b=);scanfsca

31、nf(%f,%f,&a,&b);(%f,%f,&a,&b);printfprintf(input(input jieshujieshu m=);m=);scanfscanf(%d,&m);(%d,&m);printfprintf(input(input fengeshufengeshu n=);n=);scanfscanf(%d,&n);(%d,&n);s=s=legendrelegendre_gauss(a,b,m,n);_gauss(a,b,m,n);printfprintf(s=%lfn,s);(s=%lfn,s);程序运行结果：程序运行结果：input input duandiandu

32、andian a,b=2.5,8.4 a,b=2.5,8.4input input jieshujieshu m=10 m=10input input fengeshufengeshu n=5 n=5s=192.077774s=192.0777745.6.3 5.6.3 埃尔米特埃尔米特高斯高斯(HermiteHermite-Gauss-Gauss)求积公式求积公式 埃埃尔尔米米特特高高斯斯求求积积公公式式，特特别别适适合合于于计计算算如下形式的积分：如下形式的积分：其求积公式可以表示为以下形式：其求积公式可以表示为以下形式：高高斯斯点点 及及高高斯斯求求积积系系数数 ，参参见见表表5.25.

33、2（表表示阶数）。示阶数）。埃尔米特埃尔米特高斯求积公式的高斯求积公式的N-SN-S图描述图描述 例例5.7 5.7 使用埃尔米特使用埃尔米特高斯求积公式求下列定积分的值。高斯求积公式求下列定积分的值。#include#include double double funcfunc(double x)(double x)double y;double y;y=x*x;y=x*x;return(y);return(y);double double hermitehermite_gauss(_gauss(intint n)n)double s;double s;intint i,k;i,k;stat

34、ic double x=-2.02018287046,-0.95857246461,0.00000000000,static double x=-2.02018287046,-0.95857246461,0.00000000000,0.95857246461,2.02018287046;0.95857246461,2.02018287046;static double w=0.01995324206,0.39361932315,0.94530872048,static double w=0.01995324206,0.39361932315,0.94530872048,0.3936193231

35、5,0.01995324206;0.39361932315,0.01995324206;s=0.0;s=0.0;for(i=0;in;i+)for(i=0;in;i+)s=s+wi*s=s+wi*funcfunc(xi);(xi);return(s);return(s);main()main()double s;double s;intint n;n;clrscrclrscr();();printfprintf(input n=);(input n=);scanfscanf(%d,&n);(%d,&n);s=s=hermitehermite_gauss(n);_gauss(n);printfp

36、rintf(s=%lfn,s);(s=%lfn,s);程序运行结果：程序运行结果：input n=5input n=5s=0.886227s=0.8862275.6.4 5.6.4 拉盖尔拉盖尔高斯高斯(LaguerreLaguerre-Gauss-Gauss)求积公式求积公式 拉拉盖盖尔尔高高斯斯求求积积公公式式，特特别别适适合合于于计计算算如下形式的积分：如下形式的积分：其求积公式可以表示为以下形式：其求积公式可以表示为以下形式：高高斯斯点点 及及高高斯斯求求积积系系数数 ，参参见见表表5.35.3（表示阶数）。表示阶数）。拉盖尔拉盖尔高斯求积公式的高斯求积公式的N-SN-S图描述图描述例

37、例5.8 5.8 使用拉盖尔使用拉盖尔-高斯求积公式求下列定积分的值。高斯求积公式求下列定积分的值。#include include.hdouble double funcfunc(double x)(double x)double y;double y;y=x;y=x;return(y);return(y);double double hermitehermite_gauss(_gauss(intint n)n)double s;double s;intint i,k;i,k;static double x=0.26356031972,1.41340305911,3.59642577104,

38、static double x=0.26356031972,1.41340305911,3.59642577104,7.08581000586,12.64080084428;7.08581000586,12.64080084428;static double w=0.52175561058,0.39866681108,0.07594244968,static double w=0.52175561058,0.39866681108,0.07594244968,0.00361175868,0.00002336997;0.00361175868,0.00002336997;s=0.0;s=0.0;

39、for(i=0;in;i+)for(i=0;in;i+)s=s+wi*s=s+wi*funcfunc(xi);(xi);return(s);return(s);main()main()double s;double s;intint n;n;clrscrclrscr();();printfprintf(input n=);(input n=);scanfscanf(%d,&n);(%d,&n);s=s=hermitehermite_gauss(5);_gauss(5);printfprintf(s=%lfn,s);(s=%lfn,s);程序运行结果：程序运行结果：input n=5input

40、n=5s=1.000000s=1.0000005.7 5.7 龙贝格（龙贝格（Romberg）积分法积分法 5.7.1 5.7.1 龙贝格积分法的基本思想龙贝格积分法的基本思想 前面讲述的各种求积方法是插值求积的思想，前面讲述的各种求积方法是插值求积的思想，而龙贝格积分法的基本思想是，使用一个诸如而龙贝格积分法的基本思想是，使用一个诸如梯形求积法等代数精度较低的求积公式，相继梯形求积法等代数精度较低的求积公式，相继以步长以步长 和和 求得定积分的两个近似结果，求得定积分的两个近似结果，然后再做它们适当的线性组合，就可以得到一然后再做它们适当的线性组合，就可以得到一个代数精度更高的公式。个代数精

41、度更高的公式。5.7.2 5.7.2 实现龙贝格积分法的基本步骤实现龙贝格积分法的基本步骤（1 1）输入区间）输入区间 的端点的端点 的值的值,最大迭代次数最大迭代次数 以及容许误差以及容许误差；（2 2）计算区间计算区间 的长度的长度 ；（3 3）用梯形积分法计算积分近似值）用梯形积分法计算积分近似值 ；（4 4）对）对 计算计算 对对 计算计算 ，如果如果 ，则退出循环。，则退出循环。（5 5）如果）如果 ，则继续；否则输出无解信息，转（，则继续；否则输出无解信息，转（7 7）；）；（6 6）输出）输出 的值；的值；（7 7）结束。）结束。表表5.15.1龙贝格求积算法元素进行运算的顺序龙

42、贝格求积算法元素进行运算的顺序实现龙贝格积分法的实现龙贝格积分法的NSNS图，如图图，如图5.115.11所示。所示。例例5.9 5.9 使用龙贝格求积公式求下列定积分的值。使用龙贝格求积公式求下列定积分的值。#include include.h#include#include#define DFS_N 20 /*#define DFS_N 20 /*等分数等分数*/*/#define MAX_N 10 /*define MAX_N 10 /*最大循环次数最大循环次数*/*/#define define epseps 0.00001 /*0.00001 /*容许误差容许误差*/*/double

43、 double funcfunc(double x)(double x)double y;double y;y=4.0/(1+x*x);y=4.0/(1+x*x);return(y);return(y);double sum(double double sum(double aaaa,double bb,long,double bb,long intint n)n)double h,s;double h,s;intint i;i;h=(bb-h=(bb-aaaa)/n;)/n;s=0.0;s=0.0;for(i=1;in;i+)for(i=1;in;i+)s+=s+=funcfunc(aaaa

44、+i*h);+i*h);s=s+(s=s+(funcfunc(aaaa)+)+funcfunc(bb)/2.0;(bb)/2.0;return(h*s);return(h*s);void void rombergromberg(double a,double b)(double a,double b)double s,tMAX_N+12;double s,tMAX_N+12;intint i,flag=0;i,flag=0;long long intint n=DFS_N,m;n=DFS_N,m;t01=sum(a,b,n);t01=sum(a,b,n);n*=2;n*=2;for(m=1;m

45、MAX_N;m+)for(m=1;mMAX_N;m+)for(i=0;im;i+)for(i=0;im;i+)ti0=ti1;ti0=ti1;t01=sum(a,b,n);t01=sum(a,b,n);n*=2;n*=2;for(i=1;i=m;i+)for(i=1;i=m;i+)ti1=ti-11+(ti-11-ti-10)/(ti1=ti-11+(ti-11-ti-10)/(powpow(2,2*m)-1);(2,2*m)-1);if(if(fabsfabs(tm1-tm-11)(tm1-tm-11)epseps)printfprintf(t%ld0=%lfn,m,tm1);(t%ld0=

46、%lfn,m,tm1);flag=1;flag=1;break;break;if(flag=0)if(flag=0)printfprintf(Return no(Return no solovedsolovedn);n);main()main()double a,b;double a,b;clrscrclrscr();();printfprintf(input a,b=);(input a,b=);scanfscanf(%lf,%lf,&a,&b);(%lf,%lf,&a,&b);rombergromberg(a,b);(a,b);程序运行结果：程序运行结果：input a,b=0,1inpu

47、t a,b=0,1t20=3.141570t20=3.1415705.8 5.8 高振荡函数的积分法高振荡函数的积分法 5.8.1 5.8.1 高振荡函数的积分法的基本思想高振荡函数的积分法的基本思想 在工程实际问题中，经常会遇到如下形如在工程实际问题中，经常会遇到如下形如 的积分，当的积分，当m m充分大时为高振荡函数的积分。对于高振荡函数的积分，如果采充分大时为高振荡函数的积分。对于高振荡函数的积分，如果采用插值求积法进行积分，则在建立被积函数用插值求积法进行积分，则在建立被积函数 或或 的插值的插值多项式多项式 时，为了使时，为了使 能够很好地逼近它们，就要求能够很好地逼近它们，就要求

48、也要振荡也要振荡得厉害，即要求插值多项式得厉害，即要求插值多项式 的次数足够高。但是，高次插值实际的逼近的次数足够高。但是，高次插值实际的逼近性质很不好，实用价值不大。即使采用分段低次插值，效果也不会很理想。性质很不好，实用价值不大。即使采用分段低次插值，效果也不会很理想。因此，引进计算高振荡函数的积分的重要方法分部积分法。因此，引进计算高振荡函数的积分的重要方法分部积分法。分部积分法的基本思想是，令：其中：则有：反复利用分部积分法可以得到 分离出实部和虚部后就得到以下分部积分公式。5.8.2 5.8.2 分部积分公式分部积分公式当当积积分区分区间为间为 时时，则变为则变为 实实现现高高振振荡

49、荡函函数数积积分分的的N NS S图图例例5.10 5.10 用分部积分法计算下列高振荡积分的值。用分部积分法计算下列高振荡积分的值。其中其中 ，。取取 ，则有则有#include#include void part(double a,double b,void part(double a,double b,intint m,m,intint n,double n,double fafa,double,double fbfb,double s)double s)intint mm,i,j;mm,i,j;double double smasma,smbsmb,cmacma,cmbcmb;doub

50、le double sasa4,4,sbsb4,ca4,4,ca4,cbcb4;4;smasma=sin(m*a);=sin(m*a);smbsmb=sin(m*b);=sin(m*b);cmacma=coscos(m*a);(m*a);cmbcmb=coscos(m*b);(m*b);sasa0=0=smasma;sasa1=1=cmacma;sasa2=-2=-smasma;sasa3=-3=-cmacma;sbsb0=0=smbsmb;sbsb1=1=cmbcmb;sbsb2=-2=-smbsmb;sbsb3=-3=-cmbcmb;ca0=ca0=cmacma;ca1=-;ca1=-sm