各种数值积分

5.1引言利用牛顿—莱布尼兹（Newton—Leibniz）公式(5.1)解决函数在上的积分问题在理论和应用上都有重大的意义。然而，在实际问题中，往往会遇到一些困难。有些形式上较简单的函数，其原函数不易求出或不能用初等函数表示成有限形式；有些被积函数的原函数过于复杂；而有些函数的函数值是由实验、观测等方法得出，并没有给出具体的解析表达式。这些情形说明公式(5.1)在应用上是有局限性的，因此研究定积分的数值计算问题就显得十分必要。本章主要介绍一些常用的数值积分方法，包括梯形积分法、辛卜生积分法、变步长积分法、牛顿—柯特斯积分法、高斯积分法、龙贝格积分法以及高振荡函数的积分法。)()()d(aFbFxxfba)(xf],[ba)(xF5.2梯形积分法5.2.1梯形积分法的基本思想梯形积分法的基本思想：在积分区间上，根据给定的插值条件和，构造一个一次二项式，并以的积分值近似地代替。从几何角度而言，是以梯形面积近似地代替曲边梯形的面积。],[1kkxx])(,[kkxfx])(,[11kkxfx)(1xP1d)(1kkxxxxP1d)(kkxxxxf5.2.2梯形求积公式依据梯形积分法的基本思想，将区间分成个相等的小区间,则每个小区间的长度为,对每个小区间均实施如下的梯形求积：将这些小梯形的求积值加起来，可以得到如下梯形求积公式：其中，],[bannabh/)())()((2)()d(111kkkkxxxfxfxxxxfkk11110)(2)()()(2)()(d)(niniinbahiafbfafhxfxfxfhxxfnabh/)(5.2.3实现梯形积分法的基本步骤(1)输入区间的端点值以及分割数；(2)将区间等分成个小区间，每一个小区间的长度；(3)计算每一个等分点的函数值(4)计算(5)输出的值；(6)结束。][ba,ba,N][ba,NN/abh)(),2,,1,0()(Nihiafyi1102NiiNyyyhss图5.2梯形积分法的N-S图描述例5.1使用梯形求积公式求下列定积分的值。102d14xx#defineN16/*等分数*/floatfunc(floatx){floaty;y=4.0/(1+x*x);return(y);}voidgedianzhi(floaty[],floata,floath){inti;for(i=0;i=N;i++)y[i]=func(a+i*h);}floattrapeze(floaty[],floath){floats;inti;s=(y[0]+y[N])/2.0;for(i=1;iN;i++)s+=y[i];return(s*h);}main(){floata,b,h,s,f[N+1];clrscr();printf(inputa,b=);scanf(%f,%f,&a,&b);h=(b-a)/(float)N;gedianzhi(f,a,h);s=trapeze(f,h);printf(s=%f\n,s);}程序运行结果：inputa,b=0,1↙s=3.140942辛卜生积分法的基本思想：在积分区间上，根据给定的插值条件、和，构造—个二次插值求积多项式，并以的积分值近似地代替。从几何角度而言，是用过三点的抛物线面积近似地代替积分的曲边面积。][2kkx,x])([kkxf,x)]([11kkxf,x)]([22kkxf,x)(2xP2d)(2kkxxxxP2d)(kkxxxxf5.3辛卜生（Simpson）积分法5.3.1辛卜生积分法的基本思想5.3.2辛卜生求积公式依据辛卜生积分法的基本思想，将区间分成(必须是偶数）个相等的小区间，则每个小区间的长度为，在小区间均实施如下的辛卜生求积：将这些求积值加起来，可以得到如下辛卜生求积公式：其中:],[bannnabh/)(],[2kkxx))(2)(4)((3)d(212kkkkkxxxfxfxfhxxfshxxfkkxx3d)(22s*21s*4xfxfsba)()(ks:1ks:2为寄数项的函数值之和。为偶数项的函数值之和。)(kxf)(kxf5.3.3实现辛卜生积分法的基本步骤(1)输入区间的端点的值以及分割数；(2)将区间等分成个小区间，每一个小区间的长度；(3)计算每一个等分点的函数值；(4)计算:（计算奇数项的函数值之和）（计算偶数项的函数值之和）(5)计算；(6)输出的值；(7)结束。][ba,ba,N][ba,NN/abh)(),2,,1,0()(Nihiafyi/2)2N(1i1i21yy1s2/)2N(1ii2y2s2s2140*s*yysn3/s*h图5.4辛卜生积分法的N-S图描述例5.2使用辛卜生求积公式求下列定积分的值。102d14xx#includestdio.h#defineN16/*等分数*/floatfunc(floatx){floaty;y=4.0/(1+x*x);return(y);}voidgedianzhi(floaty[],floata,floath){inti;for(i=0;i=N;i++)y[i]=func(a+i*h);}floatsimpson(floaty[],floath){floats,s1,s2;inti;s1=y[1];s2=0.0;for(i=2;i=N-2;i=i+2){s1+=y[i+1];/*计算奇数项的函数值之和*/s2+=y[i];/*计算偶数项的函数值之和*/}s=y[0]+y[N]+4.0*s1+2.0*s2;return(s*h/3.0);}main(){floata,b,h,s,f[N+1];scanf(%f,%f,&a,&b);h=(b-a)/(float)N;gedianzhi(f,a,h);s=simpson(f,h);printf(s=%f\n,s);}程序运行结果：inputa,b=0,1↙s=3.1415935.4变步长求积分法5.4.1变步长求积分法的基本思想变步长求积分法是以梯形公式为基础，逐步改变步长，以达到预先所要求的精度。变步长求积分法主要有变步长梯形求积分法和变步长辛卜生求积分法，本节我们将介绍这两种方法。5.4.2变步长梯形求积分法变步长梯形求积分法的基本过程:（1）利用梯形公式，将积分区间一等分，即其中（2）将每一个求积小区间再二等分一次（即由原来的等分变成等分），则有][ba,101)()(2niiinxfxfhTabhnn2101050110150502)(2)]()([42)()(2)()(2nini.iiinii.i.iinxfhxfxfhxfxfxfxfhT其中，为再二等分一次后新增加的结点，它们都是原来各小区间的中点。由上式可以看出，在对每一个小区间在二等分后，在积分值的第一项中只包含再二等分之前的各结点上的函数值，并且第一项的值正好是再二等分之前积分值的一半，显然，这一项中所包含的函数值就不必计算了。再二等分后需要计算的函数都包含在第二项中，它们都是二等分后出现的新的结点。因此有（3）若，二等分后的积分值就是最后的结果；否则保存当前的等分数、积分值与步长，即转到第（2）步继续做二等分处理。50.ixnT2nT10502)(221ni.innxfhTTnnTT2nT2n2nnnTT22hh5.4.3实现变步长梯形积分法的基本步骤（1）输入区间的端点的值以及容许误差ε；（2）计算区间的长度；（3）首先将区间进行一等分，并设；（4）再将每一个求积小区间，由原来的等分变成等分，即，；（5）计算（6）成立，则继续；否则，转（4）；（7）输出的值；（8）结束。][b,ab,a][b,aabh-=][b,a1=n2/))()((bfafhtnnn2nn22/h=h2/))0.5)(((102hiafhttninn|tt|nn-2nntt2nt2图5.5变步长梯形积分法的N-S图描述floatbtrapeze(floata,floatb){inti,n;floath,p,t1,t2,s,x;n=1;h=b-a;t1=h*(func(a)+func(b))/2.0;p=eps+1.0;while(p=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;h=h/2.0;}return(t2);}voidmain(){floata,b,t;clrscr();printf(inputa,b=);scanf(%f,%f,&a,&b);t=btrapeze(a,b);printf(t=%f\n,t);}程序运行结果：inputa,b=0,1↙t=0.7468245.4.4变步长辛卜生求积分法变步长辛卜生求积分法的基本过程:（1）利用梯形公式，将积分区间一等分，（2）将其中的每一个求积小区间再二等分一次（3）根据上面两式和的余项、，可以推导出如下的变步长辛卜生求积公式。进一步得到再二次等分一次后的变步长辛卜生求积公式为（4）若，二等分后的积分值就是最后的结果；否则保存当前的变步长梯形积分值、等分数、积分值与步长，转到第（2）步继续做二等分处理。][b,a101)()(2niiinxfxfhT1n0i5.0inn2xf2hT21T)(nTnT2)(1212fhab)(21222fhab342nnnTTS34242nnnTTSnnSS2nS25.4.5实现变步长辛卜生积分法的基本步骤变步长梯形积分法的N-S图描述例5.4使用变步长辛卜生求积分法求下列定积分的值。#includestdio.h#includemath.h#defineeps0.000001/*容许误差*/floatfunc(floatx){floaty;y=sqrt(1-x*x);return(y);}102d1xxfloatbsimpson(floata,floatb){inti,n;floath,p,e,s;floatt1,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=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;/*计算*/s2=(4*t2-t1)/3.0;/*计算*/e=fabs(s2-s1);/*计算精度*/t1=t2;s1=s2;n=n+n;h=h/2.0;}return(s2);}voidmain(){floata,b,s;clrscr();printf(inputa,b=);scanf(%f,%f,&a,&b);s=bsimpson(a,b);printf(s=%f\n,s);}程序运行结果：inputa,b=0,1↙s=0.7853985.5牛顿—柯特斯(Newton—Cotes)积分法5.5.1牛顿—柯特斯积分法的基本思想牛顿—柯特斯积分法的基本思想：用高次的插值求积多项式去逼近被积函数，以获得高精度的积分值。事实上，梯形积分是当时的牛顿—柯特斯积分，辛卜生积分是当时的牛顿—柯特斯积分，它们都是牛顿—柯特斯积分的特例。)(xPn)(xf1n2n5.5.2牛顿—柯特斯求积公式下面给出三到五阶牛顿—柯特斯求积公式。实现三阶牛顿—柯特斯求积公式的基本步骤如下