各种数值积分

1、各种数值积分第1页，共85页，2022年，5月20日，7点25分，星期二5.1 引言 利用牛顿莱布尼兹（NewtonLeibniz）公式 (5.1) 解决函数 在 上的积分问题在理论和应用上都有重大的意义。然而，在实际问题中，往往会遇到一些困难。有些形式上较简单的函数，其原函数 不易求出或不能用初等函数表示成有限形式；有些被积函数的原函数过于复杂；而有些函数的函数值是由实验、观测等方法得出，并没有给出具体的解析表达式。这些情形说明公式(5.1)在应用上是有局限性的，因此研究定积分的数值计算问题就显得十分必要。 本章主要介绍一些常用的数值积分方法，包括梯形积分法、辛卜生积分法、变步长积分法、牛顿

2、柯特斯积分法、高斯积分法、龙贝格积分法以及高振荡函数的积分法。第2页，共85页，2022年，5月20日，7点25分，星期二5.2 梯形积分法 5.2.1 梯形积分法的基本思想 梯形积分法的基本思想：在积分区间 上，根据给定的插值条件 和 ，构造一个一次二项式 ，并以 的积分值近似地代替 。从几何角度而言，是以梯形面积近似地代替曲边梯形的面积。第3页，共85页，2022年，5月20日，7点25分，星期二5.2.2 梯形求积公式 依据梯形积分法的基本思想，将区间 分成 个 相等的小区间,则每个小区间的长度为 ,对每个小区间均实施如下的梯形求积： 将这些小梯形的求积值加起来，可以得到如下梯形求积公式

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

4、*/float func(float x) float y; y=4.0/(1+x*x); return(y);void gedianzhi(float y,float a,float h) int i; for(i=0;i=N;i+) yi=func(a+i*h);float trapeze(float y,float h) float s; int i; s=(y0+yN)/2.0; for(i=1;iN;i+) s+=yi; return(s*h);第8页，共85页，2022年，5月20日，7点25分，星期二main() float a,b,h,s,fN+1; clrscr(); pri

5、ntf(input a,b=); scanf(%f,%f,&a,&b); h=(b-a)/(float)N; gedianzhi(f,a,h); s=trapeze(f,h); printf(s=%fn,s);程序运行结果：input a,b=0,1s=3.140942第9页，共85页，2022年，5月20日，7点25分，星期二 辛卜生积分法的基本思想：在积分区间 上，根据给定的插值条件 、 和 ，构造个二次插值求积多项式 ，并以 的积分值近似地代替 。从几何角度而言，是用过三点的抛物线面积近似地代替积分的曲边面积。 5.3 辛卜生（Simpson）积分法 5.3.1 辛卜生积分法的基本思想第

6、10页，共85页，2022年，5月20日，7点25分，星期二5.3.2 辛卜生求积公式 依据辛卜生积分法的基本思想，将区间 分成 ( 必须是偶数）个相等的小区间，则每个小区间的长度为 ，在小区间 均实施如下的辛卜生求积： 将这些求积值加起来，可以得到如下辛卜生求积公式： 其中: 为寄数项的函数 值之和。为偶数项的函数 值之和。第11页，共85页，2022年，5月20日，7点25分，星期二5.3.3 实现辛卜生积分法的基本步骤(1) 输入区间 的端点的 值以及分割数 ；(2)将区间 等分成 个小区间，每一个小区间的长度 ；(3) 计算每一个等分点的函数值 ；(4) 计算: （计算奇数项的函数值之

7、和） （计算偶数项的函数值之和）(5) 计算 ；(6) 输出 的值；(7) 结束。第12页，共85页，2022年，5月20日，7点25分，星期二图 5.4 辛卜生积分法的N-S图描述第13页，共85页，2022年，5月20日，7点25分，星期二例5.2 使用辛卜生求积公式求下列定积分的值。第14页，共85页，2022年，5月20日，7点25分，星期二#include #define N 16 /* 等分数 */float func(float x) float y; y=4.0/(1+x*x); return(y);void gedianzhi(float y,float a,float h)

8、 int i; for(i=0;i=N;i+) yi=func(a+i*h);第15页，共85页，2022年，5月20日，7点25分，星期二float simpson(float y,float h) float s,s1,s2; int i; s1=y1; s2=0.0; for(i=2;i=eps) /* 判断是否达到精度要求,若没有达到,继续循环 */ s=0.0; for(i=0;i=n-1;i+) x=a+(i+0.5)*h; s=s+func(x); t2=(t1+h*s)/2.0; /* 计算 */ p=fabs(t1-t2); /* 计算精度 */ t1=t2; n=n+n;

9、h=h/2.0; return(t2);第22页，共85页，2022年，5月20日，7点25分，星期二void main() float a,b,t; clrscr(); printf(input a,b=); scanf(%f,%f,&a,&b); t=btrapeze(a,b); printf(t=%fn,t);程序运行结果：input a,b=0,1t=0.746824第23页，共85页，2022年，5月20日，7点25分，星期二 5.4.4 变步长辛卜生求积分法变步长辛卜生求积分法的基本过程:（1）利用梯形公式，将积分区间 一等分，（2）将其中的每一个求积小区间再二等分一次（3）根据上

10、面两式 和 的余项 、 ，可以 推导出如下的变步长辛卜生求积公式。 进一步得到再二次等分一次后的变步长辛卜生求积公式为（4）若 ，二等分后的积分值 就是最后的结果；否则保存当前的变步长梯形积分值、等分数、积分值与步长，转到第（2）步继续做二等分处理。第24页，共85页，2022年，5月20日，7点25分，星期二5.4.5 实现变步长辛卜生积分法的基本步骤第25页，共85页，2022年，5月20日，7点25分，星期二第26页，共85页，2022年，5月20日，7点25分，星期二变步长梯形积分法的N-S图描述第27页，共85页，2022年，5月20日，7点25分，星期二例5.4 使用变步长辛卜生求

11、积分法求下列定积分的值。#include #include #define eps 0.000001 /* 容许误差 */float func(float x) float y; y=sqrt(1-x*x); return(y);第28页，共85页，2022年，5月20日，7点25分，星期二 float bsimpson(float a,float b) int i,n; float h,p,e,s; float t1,t2,s1,s2,x; n=1; h=b-a; t1=h*(func(a)+func(b)/2.0; s1=t1; /* 用代替 */ e=eps+1.0; while(e=e

12、ps) s=0.0; for(i=0;i=n-1;i+) x=a+(i+0.5)*h; s=s+func(x); t2=(t1+h*s)/2.0; /* 计算 */ s2=(4*t2-t1)/3.0; /* 计算 */ e=fabs(s2-s1); /* 计算精度 */ t1=t2; s1=s2; n=n+n; h=h/2.0; return(s2); 第29页，共85页，2022年，5月20日，7点25分，星期二void main() float a,b,s; clrscr(); printf(input a,b=); scanf(%f,%f,&a,&b); s=bsimpson(a,b);

13、 printf(s=%fn,s);程序运行结果：input a,b=0,1s=0.785398第30页，共85页，2022年，5月20日，7点25分，星期二5.5 牛顿柯特斯(NewtonCotes)积分法 5.5.1 牛顿柯特斯积分法的基本思想 牛顿柯特斯积分法的基本思想：用高次的插值求积多项式 去逼近被积函数 ，以获得高精度的积分值。 事实上，梯形积分是当 时的牛顿柯特斯积分，辛卜生积分是当 时的牛顿柯特斯积分，它们都是牛顿柯特斯积分的特例。第31页，共85页，2022年，5月20日，7点25分，星期二5.5.2 牛顿柯特斯求积公式 下面给出三到五阶牛顿柯特斯求积公式。第32页，共85页，

14、2022年，5月20日，7点25分，星期二实现三阶牛顿柯特斯求积公式的基本步骤如下：第33页，共85页，2022年，5月20日，7点25分，星期二三阶牛顿柯特斯求积公式的N-S图描述第34页，共85页，2022年，5月20日，7点25分，星期二例5.5 使用牛顿柯特斯求积公式求下列定积分的值。第35页，共85页，2022年，5月20日，7点25分，星期二#include #define N 5float func(float x) float y; y=4.0/(1+x*x); return(y);void gedianzhi(float y,float a,float b,int n) fl

15、oat h,s; int i; h=(b-a)/(float)n; for(i=0;i=n;i+) yi=func(a+i*h);第36页，共85页，2022年，5月20日，7点25分，星期二float nc3(float y,float a,float b,int n) float h,s,s0,s1; int i; h=(b-a)/(float)n; s0=s1=0.0; for(i=0;i=n-3;i=i+3) s0+=yi+yi+3; s1+=yi+1+yi+2; s=s0+3.0*s1; return(s*h*3.0/8.0);main() float a,b,h,s,fN+1; f

16、loat n3,n4,n5; printf(input a,b=); scanf(%f,%f,&a,&b); gedianzhi(f,a,b,3); n3=nc3(f,a,b,3); printf(n 3-nc 4-nc 5-ncn); printf(%f %f %fn,n3,n4,n5);第37页，共85页，2022年，5月20日，7点25分，星期二程序运行结果：input a,b=0,1 3-nc 4-nc 5-nc 3.138461 3.142118 3.141878第38页，共85页，2022年，5月20日，7点25分，星期二5.6 高斯积分法 5.6.1 高斯积分法的基本思想 前面介

17、绍的几种数值积分方法，都是先寻找一个插值求积多项式 ，并以 近似代替函数 进行积分而获得积分的近似值，即 (5.2) 表明，函数 在区间 上的积分可以用函数 在该区间上 的 个点的函数值的线性组合来近似代替。但是，由于插值求积公式是利用插值多项式的积分得到的，因此，如果被积函数 为次数不超过 的多项式，则利用插值求积公式计算得到的积分值是准确的。 在实际应用中，为了提高数值求积公式的精度，一般要求数值求积公式对于次数尽可能高的多项式能准确成立。由此提出了高斯积分法。即：如果插值求积公式(5.2)具有 次代数精度，那么称该插值求积公式为高斯求积公式。其节点 称为高斯点， 称为高斯求积系数。 下面

18、具体介绍几个常用的高斯求积公式。 第39页，共85页，2022年，5月20日，7点25分，星期二5.6.2 勒让德高斯(Legendre-Gauss)求积公式 勒让德高斯求积公式，特别适合于计算区间-1,1的积分，其求积公式可以表示为以下形式： 而对于一般区间 ，通过变换可以得到以下形式的求积公式： 高斯点 及高斯求积系数 ，参见表5.1（ 表示阶数）。第40页，共85页，2022年，5月20日，7点25分，星期二第41页，共85页，2022年，5月20日，7点25分，星期二第42页，共85页，2022年，5月20日，7点25分，星期二第43页，共85页，2022年，5月20日，7点25分，星

19、期二第44页，共85页，2022年，5月20日，7点25分，星期二第45页，共85页，2022年，5月20日，7点25分，星期二第46页，共85页，2022年，5月20日，7点25分，星期二第47页，共85页，2022年，5月20日，7点25分，星期二第48页，共85页，2022年，5月20日，7点25分，星期二例5.6 使用勒让德高斯求积公式求下列定积分的值。第49页，共85页，2022年，5月20日，7点25分，星期二#include double func(double x) double y; y=x*x+sin(x); return(y);double legendre_gauss(

20、double a,double b,int m,int n) double h,hx,y,s,dx,x0; int i,k; static double x=-0.9061798459,-0.5384693101,0.0000000000, 0.5384693101,0.9061798459; static double w=0.2369268851,0.4786286705,0.5688888889, 0.4786286705,0.2369268851; dx=(b-a)/(double)m; hx=dx/2; s=0.0; for(k=0;km;k+) x0=a+(double)k*dx+

21、hx; for(i=0;in;i+) y=x0+xi*hx; s=s+wi*func(y); return(s*hx);第50页，共85页，2022年，5月20日，7点25分，星期二main() double a,b,s; int m,n; clrscr(); printf(input duandian a,b=); scanf(%f,%f,&a,&b); printf(input jieshu m=); scanf(%d,&m); printf(input fengeshu n=); scanf(%d,&n); s=legendre_gauss(a,b,m,n); printf(s=%lfn

22、,s);程序运行结果：input duandian a,b=2.5,8.4input jieshu m=10input fengeshu n=5s=192.077774第51页，共85页，2022年，5月20日，7点25分，星期二5.6.3 埃尔米特高斯(Hermite-Gauss)求积公式 埃尔米特高斯求积公式，特别适合于计算如下形式的积分： 其求积公式可以表示为以下形式： 高斯点 及高斯求积系数 ，参见表5.2（ 表示阶数）。第52页，共85页，2022年，5月20日，7点25分，星期二第53页，共85页，2022年，5月20日，7点25分，星期二第54页，共85页，2022年，5月20日

23、，7点25分，星期二埃尔米特高斯求积公式的N-S图描述 第55页，共85页，2022年，5月20日，7点25分，星期二例5.7 使用埃尔米特高斯求积公式求下列定积分的值。第56页，共85页，2022年，5月20日，7点25分，星期二#include double func(double x) double y; y=x*x; return(y);double hermite_gauss(int n) double s; int i,k; static double x=-2.02018287046,-0.95857246461,0.00000000000, 0.95857246461,2.02

24、018287046; static double w=0.01995324206,0.39361932315,0.94530872048, 0.39361932315,0.01995324206; s=0.0; for(i=0;in;i+) s=s+wi*func(xi); return(s);第57页，共85页，2022年，5月20日，7点25分，星期二main() double s; int n; clrscr(); printf(input n=); scanf(%d,&n); s=hermite_gauss(n); printf(s=%lfn,s);程序运行结果：input n=5s=

25、 0.886227第58页，共85页，2022年，5月20日，7点25分，星期二5.6.4 拉盖尔高斯(Laguerre-Gauss)求积公式 拉盖尔高斯求积公式，特别适合于计算如下形式的积分： 其求积公式可以表示为以下形式： 高斯点 及高斯求积系数 ，参见表5.3（ 表示阶数）。第59页，共85页，2022年，5月20日，7点25分，星期二第60页，共85页，2022年，5月20日，7点25分，星期二第61页，共85页，2022年，5月20日，7点25分，星期二拉盖尔高斯求积公式的N-S图描述第62页，共85页，2022年，5月20日，7点25分，星期二例5.8 使用拉盖尔-高斯求积公式求下

26、列定积分的值。第63页，共85页，2022年，5月20日，7点25分，星期二#include double func(double x) double y; y=x; return(y);double hermite_gauss(int n) double s; int i,k; static double x=0.26356031972,1.41340305911,3.59642577104, 7.08581000586,12.64080084428; static double w=0.52175561058,0.39866681108,0.07594244968, 0.003611758

27、68,0.00002336997; s=0.0; for(i=0;in;i+) s=s+wi*func(xi); return(s); 第64页，共85页，2022年，5月20日，7点25分，星期二main() double s; int n; clrscr(); printf(input n=); scanf(%d,&n); s=hermite_gauss(5); printf(s=%lfn,s);程序运行结果：input n=5s= 1.000000第65页，共85页，2022年，5月20日，7点25分，星期二5.7 龙贝格（Romberg）积分法 5.7.1 龙贝格积分法的基本思想 前面

28、讲述的各种求积方法是插值求积的思想，而龙贝格积分法的基本思想是，使用一个诸如梯形求积法等代数精度较低的求积公式，相继以步长 和 求得定积分的两个近似结果，然后再做它们适当的线性组合，就可以得到一个代数精度更高的公式。 第66页，共85页，2022年，5月20日，7点25分，星期二5.7.2 实现龙贝格积分法的基本步骤（1）输入区间 的端点 的值,最大迭代次数 以及容许误差；（2）计算区间 的长度 ；（3）用梯形积分法计算积分近似值 ；（4）对 计算 对 计算 ， 如果 ，则退出循环。（5）如果 ，则继续；否则输出无解信息，转（7）；（6）输出 的值；（7）结束。第67页，共85页，2022年，

29、5月20日，7点25分，星期二 表5.1龙贝格求积算法元素进行运算的顺序实现龙贝格积分法的NS图，如图5.11所示。 第68页，共85页，2022年，5月20日，7点25分，星期二例5.9 使用龙贝格求积公式求下列定积分的值。第69页，共85页，2022年，5月20日，7点25分，星期二#include #include #define DFS_N 20 /* 等分数 */#define MAX_N 10 /* 最大循环次数 */#define eps 0.00001 /* 容许误差 */double func(double x) double y; y=4.0/(1+x*x); return

30、(y);double sum(double aa,double bb,long int n) double h,s; int i; h=(bb-aa)/n; s=0.0; for(i=1;in;i+) s+=func(aa+i*h); s=s+(func(aa)+func(bb)/2.0; return(h*s);第70页，共85页，2022年，5月20日，7点25分，星期二void romberg(double a,double b) double s,tMAX_N+12; int i,flag=0; long int n=DFS_N,m; t01=sum(a,b,n); n*=2; for

31、(m=1;mMAX_N;m+) for(i=0;im;i+) ti0=ti1; t01=sum(a,b,n); n*=2; for(i=1;i=m;i+) ti1=ti-11+(ti-11-ti-10)/(pow(2,2*m)-1); if(fabs(tm1-tm-11)eps) printf(t%ld0=%lfn,m,tm1); flag=1; break; if(flag=0) printf(Return no solovedn);第71页，共85页，2022年，5月20日，7点25分，星期二main() double a,b; clrscr(); printf(input a,b=);

32、scanf(%lf,%lf,&a,&b); romberg(a,b);程序运行结果：input a,b=0,1t20=3.141570第72页，共85页，2022年，5月20日，7点25分，星期二5.8 高振荡函数的积分法 5.8.1 高振荡函数的积分法的基本思想 在工程实际问题中，经常会遇到如下形如 的积分，当m充分大时为高振荡函数的积分。对于高振荡函数的积分，如果采用插值求积法进行积分，则在建立被积函数 或 的插值多项式 时，为了使 能够很好地逼近它们，就要求 也要振荡得厉害，即要求插值多项式 的次数足够高。但是，高次插值实际的逼近性质很不好，实用价值不大。即使采用分段低次插值，效果也不会

33、很理想。 因此，引进计算高振荡函数的积分的重要方法分部积分法。 第73页，共85页，2022年，5月20日，7点25分，星期二分部积分法的基本思想是， 令： 其中： 则有： 反复利用分部积分法可以得到 分离出实部和虚部后就得到以下分部积分公式。第74页，共85页，2022年，5月20日，7点25分，星期二5.8.2 分部积分公式当积分区间为 时，则变为 第75页，共85页，2022年，5月20日，7点25分，星期二实现高振荡函数积分的NS图第76页，共85页，2022年，5月20日，7点25分，星期二例5.10 用分部积分法计算下列高振荡积分的值。 其中 ， ， 。取 ， ， ， 则有 第77

34、页，共85页，2022年，5月20日，7点25分，星期二#include void part(double a,double b,int m,int n,double fa,double fb,double s) int mm,i,j; double sma,smb,cma,cmb; double sa4,sb4,ca4,cb4; sma=sin(m*a); smb=sin(m*b); cma=cos(m*a); cmb=cos(m*b); sa0=sma; sa1=cma; sa2=-sma; sa3=-cma; sb0=smb; sb1=cmb; sb2=-smb; sb3=-cmb; ca0=cma; ca1=-sma; ca2=-cma; ca3=sma; cb0=cmb; cb1=-smb; cb2=-cmb; cb3=smb; s0=0.0; s1=0.0;