计算方法第6章 数值积分

1第6章数值积分基础知识牛顿—柯特斯公式复化求积公式龙贝格求积26.1

基础知识计算连续函数在区间[a,b]上的定积分，可以用牛顿—莱布尼兹（Newton-Leibnitz）公式：

一、问题的提出=F(b)-F(a)其中F(x)是f(x)在区间[a,b]上的一个原函数。

在某些情况下，不适合用牛顿—莱布尼兹公式求定积分：

通过观测等方法得到f(x)上某些离散的点，不知道或不存在f(x)的具体解析表达式。f(x)的原函数F(x)不能用有限形式表示。

求原函数F(x)较麻烦，但利用计算机等工具，用数值积分的方法求解较方便。数值积分方法求得的结果虽然有误差，但可以控制在要求的误差范围之内。实际工程问题一般不需要理论上无误差的精确解，能得到误差满足要求的近似解就可以了。数值积分的主要思想，是避开求f(x)的原函数F(x)，只由f(x)上若干离散的点，用对应的数值积分算法，来求f(x)在[a,b]上定积分的近似值。

36.1

基础知识二、数值积分公式由定积分中值定理，在积分区间[a,b]上至少存在1个点ξ，使得=(b-a)f(ξ)f(ξ)为f(x)在区间[a,b]上的平均高度。如果知道了平均高度f(ξ)，就能够求出积分。现在的问题是，如何由f(x)上若干个点f(xi)，(i=0,1,2,……,n)，构造1个表达式，替代定积分中值定理中的平均高度f(ξ)

。

数值积分方法就是在积分区间[a,b]上选取若干节点x0,x1,……,xn，用f(x0),f(x1),……,f(xn)的线性组合作为积分的近似值。数值积分公式的一般形式为：

≈f(xi)R[f]=－

f(xi)其中xi称为求积节点，Ai称为求积系数，R[f]为数值积分公式的余项。Ai是与求积节点对应的常数，它不依赖被积函数f(x)的具体形式。这一类求积方法避开了求f(x)的原函数F(x)，易于用计算机实现，因此又称为机械求积。46.1

基础知识二、数值积分公式（续）

中矩形公式

I=(b-a)f((a+b)/2)取积分区间[a,b]的中点x=(a+b)/2近似地代替ξ，用f((a+b)/2)取代平均高度f(ξ)。中矩形公式是1点公式，又是开公式（求积节点不包含积分区间的端点）。1点公式还有左矩形公式和右矩形公式。

左矩形公式I=(b-a)f(a)取积分区间[a,b]的端点x=a近似地代替ξ，用f(a)取代平均高度f(ξ)。

右矩形公式

I=(b-a)f(b)取积分区间[a,b]的端点x=b近似地代替ξ，用f(b)取代平均高度f(ξ)。

左矩形公式和右矩形公式是半开半闭公式（求积节点仅包含积分区间1侧的端点）。

梯形公式

I=(b-a)(f(a)+f(b))/2)取f(a)、f(b)的平均值f(a)+f(b))/2近似地代替平均高度f(ξ)

。梯形公式是2点公式，又是闭公式（求积节点包含积分区间的端点）。

56.1

基础知识三、代数精度代数精度是间接地反应某一数值积分公式对被积函数逼近能力的一个参数。

定义：如果求积公式对所有不超过m次的代数多项式函数都能精确成立，而对m＋1次的代数多项式函数至少有1个不能精确成立，则称此求积公式具有m次代数精度。若被积函数为n次代数多项式函数，求积公式有m次代数精度，且m≥n，则求得的积分值误差为0。

若被积函数不是代数多项式函数，那么不能保证有限次代数精度的数值积分公式的余项为0。过高阶次代数精度的数值积分有时会出现很大误差，因此建议不使用代数精度超过6次的数值积分公式。当积分区间[a,b]较大时，可以采用复化求积公式等方法来减小误差。

定理：求积公式具有m次代数精度的充要条件是当被积函数f(x)=xr，r=0,1,2,……,m时精确成立，当被积函数f(x)=xm+1时不能精确成立。

中矩形公式和梯形公式具有1次代数精度，左矩形公式和右矩形公式具有0次代数精度。

66.1

基础知识四、插值型求积公式插值与数值积分有密切联系。我们可以先对求积节点做插值，再对插值函数求积分，如果插值函数与被积函数足够接近，就可以把插值函数的积分近似为被积函数的积分。代数插值函数具有存在唯一性，易构造，易求积分，易用计算机实现，因此本章用第5章的方法构造代数插值函数。用拉格郎日插值函数代替被积函数求积分：

≈==f(xi)式中是与求积节点xi对应的常数，不依赖被积函数的具体形式，象这样构造出的求积公式称为插值型求积公式，它的余项函数为：

R[f]=-=当x∈[a,b]时，ξ∈[a,b]，且ξ与x有关。

如果数值积分公式≈f(xi)，（n=1,2,……）是插值型求积公式，那么=b-a。76.2

牛顿—柯特斯公式一、牛顿—柯特斯公式的一般形式牛顿—柯特斯（Newton—Cotes）公式，简记为N—C公式，是一种等距节点的代数插值型求积公式。设积分区间为[a,b]，已知n＋1个求积节点x0,x1,……,xn（x0＜x1＜……＜xn），这些求积节点把积分区间n等分，步长h=(b-a)/n，即x0=a，xn=b，xi=a+i×h（i=0,1,2,……,n），则牛顿—柯特斯公式的一般形式为：≈(b-a)=其中称为柯特斯系数。与数值积分公式的一般形式对应，求积系数Ai=(b-a)，（i=0,1,2,…,n）86.2

牛顿—柯特斯公式二、柯特斯系数由柯特斯系数的定义可以知道，柯特斯系数仅由n、i决定，它是与被积函数f(x)具体形式无关的常数。下面为部分柯特斯系数。

ni=0i=1i=2i=3i=4i=5i=611/21/221/62/31/631/83/83/81/847/9016/452/1516/457/90519/28825/9625/14425/14425/9619/288641/8409/359/28034/1059/2809/3541/8407751/172803577/172801323/172802989/172802989/172801323/17280……8989/283505888/28350-928/2835010496/28350-4540/2835010496/28350……96.2

牛顿—柯特斯公式由柯特斯系数表可以写出对应的各阶牛顿—柯特斯公式。n=1时，共有2个求积节点，对应的代数插值曲线为经过这2个求积节点的直线。从几何意义上能够看出，n=1时的牛顿—柯特斯公式就是梯形公式（又称为2点公式）。n=2时，共有3个求积节点，对应的代数插值曲线为经过这3个求积节点的抛物线。此公式称为辛普生（Simpson）公式，又称为抛物线求积公式或3点公式。

≈

n=3时，共有4个求积节点。此公式称为辛普生法则。

≈(f(x0)+3f(x1)+3f(x2)+f(x3)柯特斯系数有如下性质：

=1

=

106.2

牛顿—柯特斯公式三、牛顿—柯特斯公式的代数精度牛顿—柯特斯公式是一种插值型求积公式。n阶牛顿—柯特斯公式至少具有n次代数精度。定理：牛顿—柯特斯公式的阶次n为奇数时，此求积公式至少具有n次代数精度；阶次n为偶数时，此求积公式至少具有n＋1次代数精度。四、牛顿—柯特斯公式的余项下面给出几种常用低阶牛顿—柯特斯公式余项函数的一般形式。

如果被积函数f(x)在积分区间[a,b]上有连续的2阶导数，那么梯形公式的余项函数R1[f]=－(b-a)3，其中η∈[a,b]。

如果被积函数f(x)在积分区间[a,b]上有连续的4阶导数，那么辛普生公式的余项函数R2[f]=－(b-a)5，其中η∈[a,b]。

如果被积函数在积分区间[a,b]上有连续的6阶导数，那么柯特斯公式的余项函数R4[f]=－f(6)(η)，其中η∈[a,b]。

116.2

牛顿—柯特斯公式五、牛顿—柯特斯公式的稳定性在实际应用中，求积节点不可避免地带有误差。一个稳定的算法应该能够有效地控制误差的传播。假设在求积节点x=xi处，精确值记为f(xi)，近似值记为f*(xi)，i=0,1,2,……,n，下面给出数值积分公式稳定性的定义。

定义：对于任意ε＞0，若存在δ＞0，只要|f(xi)-f*(xi)|≤δ，（i=0,1,2,……,n），就有|In(f)-In(f*)|≤ε成立，则称求积公式In(f)=f(xi)是稳定的。1～7阶牛顿—柯特斯公式的求积系数Ai＞0，（i=0,1,2,……,n），因此牛顿—柯特斯公式在阶次不大于7时是稳定的。当阶次大于7时，某些求积系数为负，这时不能保证公式的稳定性。对于较大积分区间、复杂被积函数、较高精度要求的数值积分问题，不应选择高阶牛顿—柯特斯公式，应该使用复化求积公式。

定理6.9若求积公式In(f)=f(xi)至少具有0次代数精度，而且求积系数Ai＞0，（i=0,1,2,……,n），则此求积公式是稳定的。126.3

复化求积公式对于较大积分区间、复杂被积函数、较高精度要求的数值积分问题，需要较多的求积节点。如果采用高阶插值型求积公式，那么当被积函数f(x)不是多项式函数时，求积过程可能不稳定，因此这时只能采用复化求积。

一、问题的提出复化求积的主要思想为：把积分区间[a,b]分成若干个小子区间，使被积函数f(x)在每个小子区间内足够平滑，然后用低阶数值积分公式求f(x)在每个小子区间内的积分，再对所有小子区间的积分结果求和，就得到f(x)在[a,b]上的积分。本节重点讨论等距节点复化牛顿—柯特斯公式求积：把积分区间[a,b]等分为m个小子区间，每个小子区间都被n＋1个求积节点（含小子区间端点）n等分，用n阶牛顿—柯特斯公式求每个小子区间内的积分，再对m个小子区间的积分结果求和，就得到在[a,b]上的积分。136.3

复化求积公式把积分区间[a,b]等分为m个小子区间，令步长h=(b-a)/m，求积节点xi=a＋i×h，i=0,1,2,……,m，等距节点复化梯形公式为：二、等距节点复化求积Tm=余项函数（即截断误差）为：Rm(f)=－(b-a)f"(η)，其中η

∈[a,b]。等距节点复化辛普生公式求积是指把积分区间[a,b]等分为m个小子区间，每个小子区间取子区间的中点和2个端点为求积节点，用辛普生公式求积，再对m个小子区间的积分结果求和，就得到在[a,b]上的积分。等距节点复化柯特斯公式求积是指把积分区间[a,b]等分为m个小子区间，每个小子区间被2个端点和3个内部求积节点4等分，用柯特斯公式对每个小子区间求积，再对m个小子区间的积分结果求和，就得到在[a,b]上的积分。146.3

复化求积公式等距节点复化梯形公式求积的算法输入积分区间边界a,b和求积节点的个数n。步长h=(b-a)/(n-1)输入各求积节点纵坐标f[n]。integral=(f[0]+f[n-1])/2;for(i=1;i<n-1;i++)integral+=f[i];integral*=h;输出积分结果integral。156.3

复化求积公式等距节点复化梯形公式求积的程序（1/2）#include<stdio.h>#defineMAXSIZE50voidinput(doublef[MAXSIZE],doublea,doubleh,longn);voidmain(void){ doublea,b,h,f[MAXSIZE],integral; longn,i; printf("\n请输入积分区间边界a,b："); scanf("%lf,%lf",&a,&b); printf("\n请输入求积节点的个数："); scanf("%ld",&n); h=(b-a)/(n-1); input(f,a,h,n); integral=(f[0]+f[n-1])/2; for(i=1;i<n-1;i++) integral+=f[i]; integral*=h; printf("\n积分值=%lf",integral);}166.3

复化求积公式等距节点复化梯形公式求积的程序（2/2）voidinput(doublef[MAXSIZE],doublea,doubleh,longn){ longi; printf("\n请输入求积节点纵坐标:"); for(i=0;i<=n-1;i++) { printf("\nx[%ld]=%lf,f[%ld]=",i,a+i*h,i); scanf("%lf",&f[i]); }}176.3

复化求积公式前面介绍的各种复化求积算法在求积过程中，步长固定不变，这类求积公式称为定步长求积公式。定步长求积在计算之前需要根据精度要求来确定步长，但是用余项公式估算误差时，需要对被积函数求高阶导数，实际问题经常不能满足这一要求。除此之外，还需要找出积分区间[a,b]内的η，这也难以实现。如果用[a,b]内导数的最大值代替η处的导数，可能会使估算出的误差过大。因此，实际应用中常常采用变步长求积。三、变步长求积公式变步长求积的主要思想是：先用较大的步长计算积分，在此基础上用较小步长重新积分，再用更小步长重新积分，……，反复地计算积分值，每次积分的步长都比上次积分的步长小。当步长h→0时，等距节点复化求积公式的误差→0，积分结果→精确值，相邻2次积分结果之差→0。因此，相邻2次积分结果之差能够间接地反映积分结果的误差大小。一般情况下，当相邻2次积分结果之差的绝对值足够小时，可以近似地认为精度满足要求，积分过程结束。象这样在求积过程中，步长逐步改变的求积公式称为变步长求积公式。变步长求积在计算之前不需要知道多小的步长能满足精度要求，精度要求由相邻2次积分结果之差的大小来间接地控制。

186.3

复化求积公式常用的变步长求积公式是逐次分半梯形公式求积：先用梯形公式求f(x)在积分区间[a,b]上的积分，记为T1；取积分区间的中点，把原积分区间对分为2个小积分区间，用等距节点复化梯形公式求积，记为T2；取已有小积分区间的中点，把已有小积分区间对分为2个更小的小积分区间，用等距节点复化梯形公式求积，记为T4；……，每次对分都使积分区间的个数增加1倍，把原积分区间的中点增加为新的求积节点，原求积节点在后面各次求积时仍为求积节点，这样可以避免一些重复的运算。

三、变步长求积公式（续）由Tn递推出T2n的逐次分半梯形公式为：T2n=＋h2n

196.3

复化求积公式逐次分半梯形公式求积的算法输入积分区间边界a,b和精度要求ε。步长h=b-a，子区间个数n=1。首次积分t2n=h*(f(a)+f(b))/2

暂存tn=t2nh/=2;t2n=0;for(i=0;i<=n-1;i++)x=a+h*(2*i+1);t2n+=f(x);t2n=tn/2+h*t2n;n*=2;|t2n－tn|＞ε输出积分结果t2n。206.3

复化求积公式逐次分半梯形公式求积的程序（1/2）#include<stdio.h>#include<math.h>doublef(doublex);voidmain(void){ doublex,a,b,h,tn,t2n,epsilon; longn,i; printf("\n请输入积分区间边界a,b："); scanf("%lf,%lf",&a,&b); printf("\n请输入精度要求："); scanf("%lf",&epsilon); h=b-a; n=1; t2n=h*(f(a)+f(b))/2; do{ tn=t2n; h/=2; t2n=0; for(i=0;i<=n-1;i++) {x=a+h*(2*i+1); t2n+=f(x);} t2n=tn/2+h*t2n; n*=2; }while(fabs(t2n-tn)>epsilon); printf("\n积分值=%lf",t2n);}doublef(doublex){ return(……);/*计算并返回函数值f(x)*/}216.4

龙贝格求积某些数值计算方法，是构造一个序列，去逼近精确解。可以在原序列的基础上构造一个新序列，使它能够更快地收敛于精确值。对已有近似值线性组合以求更精确的近似值的加速收敛方法为外推方法。

一、外推算法用外推方法可以对逐次分半复化求积进行加速。

①设逐次分半梯形公式求积在逐次分半的过程中，得到的积分值序列：T1,T2,T4,T8,……，梯形加速公式为I≈梯形加速公式是龙贝格积分的第1次外推。若记Rn,0=，n=0,1,2,3,……，用梯形加速公式对Rn,0,Rn+1,0加速后的结果记为Rn,1，那么梯形加速公式可以记为Rn,1=，n=0,1,2,3,……。

用梯形加速公式对逐次分半复化求积加速的效果与复化辛普生公式相同。等距节点复化辛普生公式与梯形加速公式存在如下对应关系：Rn,1=，n=0,1,2,3,……。

226.4

龙贝格求积二、龙贝格求积的一般公式②对T1,T2,T4,T8,……进行第1次外推得到的结果，即等距节点复化辛普生公式的求积结果S1,S2,S4S8,……，辛普生加速公式为I≈辛普生加速公式是龙贝格积分的第2次外推。把记为Rn,1，用辛普生加速公式对Rn,1,Rn+1,1加速后的结果记为Rn,2，那么辛普生加速公式可以记为：Rn,2≈，n=0,1,2,3,……。

用辛普生加速公式的积分结果与复化柯特斯公式求积结果相同，对应关系为：

Rn,2=，n=0,1,2,3,……。

龙