数值积分课件 (《计算方法》)

第7章数值积分§1插值型求积公式§2复化求积公式§3龙贝格(Romberg)求积方法2/6/20231

§1插值型求积公式在一元函数的积分学中,我们已经熟知,若函数f(x)在区间［a,b]上连续且其原函数为F(x),则可用牛顿―莱布尼兹公式(7―1)来求定积分。2/6/20232

公式(7―1)虽然在理论上或在解决实际问题中都起了很大的作用,但它并不能完全解决定积分的计算问题。因为定积分的计算常常会碰到以下三种情况:

(1)被积函数f(x)的原函数F(x)不易找到。许多很简单的函数,例如

其原函数都不能用初等函数表示成有限形式。2/6/20233(2)被积函数f(x)没有具体的解析表达式。其函数关系由表格或图形表示,无法求出原函数。其被积函数的原函数就比较复杂,从数值计算角度来看,计算量太大。(3)尽管f(x)的原函数能表示成有限形式但其表达式相当复杂。例如定积分2/6/20234图7.1如图7.1,若用左矩形近似地代替曲边梯形,则得到左矩形公式

(7―2)2/6/20235同样可得到右矩形公式：(7―3)2/6/20236图7.2如图7.2,若用梯形的面积近似地代替曲边梯形的面积,则得到计算定积分的梯形公式(7―4)2/6/20237如图7.3,若用抛物线代替曲线f(x),则可得到抛物线公式(或辛普生公式)(7―5)图7.32/6/20238此外,众所周知的梯形公式:

I(f)≈(b-a)[f(a)+f(b)]/2和Simpson公式:

I(f)≈(b-a)[f(a)+4f((a+b)/2)+f(b)]/6则分别可以看作用a,b,c=(a+b)/2,

三点高度的加权平均值[f(a)+f(b)]/2和[f(a)+4f(c)+f(b)]/6作为平均高度f(ξ)的近似值.2/6/20239

更一般地,取区间[a,b]内n+1个点{xi},(i=0,1,2,…n)处的高度{f(xi)}(i=0,1,…,n)通过加权平均的方法近似地得出平均高度f(ξ),这类求积方法称为机械求积:2/6/202310或写成:数值积分公式求积系数

求积节点

(1)2/6/202311记称(2)为数值求积公式,(3)为求积公式余项(误差).构造或确定一个求积公式，要讨论解决的问题有(i)

确定求积系数Ak和求积节点xk；(ii)

求积公式的误差估计和收敛性为了构造形如式(2)的求积公式,需要提供一种判定求积方法精度高低准则2/6/202312求积公式的代数精度定义1

称求积公式(2)具有m次代数精度,如果它满足如下两个条件:

(i)对所有次数≤m次的多项式,有

(ii)存在m+1次多项式,使得定义1中的条件(i),(ii)等价于:2/6/202313插值型求积公式

在积分区间[a,b]上取n+1个节点xi，i=0,1,2,…,n,作f(x)的n次代数插值多项式（拉格朗日插值公式）:则有

为插值余项于是有2/6/202314取称(4)式为插值型求积公式,其中求积系数Ak由(5)式确定.(4)(5)Ak由节点决定，与f(x)无关。2/6/2023152/6/202316推论1

求积系数满足:误差定理1形如的求积公式至少有

n

次代数精度该公式为插值型（即：）2/6/202317现用第六章介绍的插值多项式Pn(x)来代替被积函数f(x),即有1.1牛顿―柯特斯公式(Newton―Cotes)

取节点为等距,即a=x0＜x1＜…＜xn=b

建立数值积分公式最基本的思想是选取一个既简单又有足够精度的函数φ(x),用φ(x)代替被积函数f(x),于是有

2/6/202318利用拉格朗日插值多项式(7―6)其中(7―7)2/6/202319这里yi=f(xi),对式(7―6)两边积分得

2/6/202320为牛顿―柯特斯(Newton-Cotes)求积公式,Rn(f)为牛顿―柯特斯求积公式的余项。我们称2/6/202321令x=x0+sh,0≤s≤ndx=hds=(b-a)/nds(7―11)2/6/202322Newton-Cotes公式的误差为:与x有关注意:由(7-11)式确定的Cotes系数只与i和n有关,与f(x)和积分区间[a,b]无关,且满足:(7-9)2/6/202323称Ci(n)为柯特斯求积系数。很显然,当n=1时,可算得此时式(7―10)为(7―12)这是梯形公式。2/6/202324

当n=2时,可得于是(7―13)这是抛物线(Simpson)公式。2/6/202325当n=3时,代入(7―10)式得到求积公式

2/6/202326类似地可分别求出n=4,5,…时的柯特斯系数,从而建立相应的求积公式。具体结果见表7―1。从表中可以看出,当n≤7时,柯特斯系数为正;从n≥8开始,柯特斯系数有正有负。因此,当n≥8时,误差有可能传播扩大,牛顿―柯特斯求积公式不宜采用。柯特斯系数Ci

(n)

仅与n和i有关,与被积函数f(x)无关,且满足(7―15)事实上,式(7―10)对f(x)=1是准确成立的。2/6/202327表7―12/6/202328定理当阶数n为偶数时,Newton-Cotes公式(8)至少具有n+1次代数精度.证明

只需验证当n为偶数时,Newton-Cotes公式对f(x)=xn+1的余项为零.由于f(x)=xn+1,所以f(n+1)(x)=(n+1)!.由式(7-9)得引进变换t=u+n/2,因为n为偶数,故n/2为整数,于是有据此可断定R(f)=0,因为上述被积函数是个奇函数.2/6/202329Newton-Cotes公式的数值稳定性

现在讨论舍入误差对计算结果产生的影响.设用公式

近似计算积分时,其中计算函数值f(xj)有误差εj

(j=0,1,2,…,n).设计算Cj(n)没有误差,中间计算过程中的舍入误差也不考虑,则在式(10)的计算中,由εj引起的误差为(10)2/6/202330如果Cj(n)都是正数,并设故en是有界的,即由εj引起的误差受到控制,不超过ε的(b-a)倍,保证了数值计算的稳定性.而当n>7时,Cj(n)将出现负数,保证数值稳定性.因此高阶公式不宜采用,有实用价值的仅仅是几种低阶的求积公式.将随n增大,因而不能则有2/6/202331解利用梯形公式利用抛物线公式例1试分别用梯形公式和抛物线公式计算积分:原积分的准确值

2/6/202332现对牛顿―柯特斯求积公式所产生的误差作一个分析。由式(7―9),牛顿―柯特斯求积公式的余项为

1.2误差估计

易知,牛顿―柯特斯求积公式(7―10)对任何不高于n次的多项式是准确成立的。这是因为

f(n+1)(ξ)≡0

故Rn(f)≡02/6/202333牛顿―柯特斯求积公式的代数精确度至少为n。通常在基点个数相等的情况下,代数精确度愈高,求积公式就愈精确。一般说来,若某个求积公式对于次数不高于m的多项式都准确成立(即Rn(f)≡0),而对于某一次数为m+1的多项式并不准确成立,则称这一求积公式的代数精确度为m。定理1(梯形公式的误差)设f(x)在区间［a,b］上具有连续的二阶导数,则梯形求积公式的误差为2/6/202334由于ω1(x)=(x-a)(x-b)证由式(7―9)知,梯形公式的余项为ω1(x)在区间(a,b)内不变号,f″(ξ)是x的函数且在［a,b］上连续,故根据积分第二中值定理知,存在某一η∈(a,b)使2/6/202335

定理2(抛物线公式的误差)设f(x)在［a,b］上有连续的四阶导数,则抛物线公式的误差为(7―17)2/6/202336n=1:TrapezoidalRule/*令x=a+th,h=ba,用中值定理*/代数精度=1n=2:Simpson’sRule代数精度=3n=4:CotesRule,代数精度=5,2/6/202337复合求积公式高次插值有Runge现象,高阶Newton-Cotes公式会出现数值不稳定,低阶Newton-Cotes公式有时又不能满足精度要求.解决这个矛盾的办法是将积分区间[a,b]分成若干小区间,在每个小区间上用低阶求积公式计算,然后将它们加起来,这就是复合求积方法.2/6/202338§2复合求积公式

2.1复合梯形公式对于定积分(7―1),将积分区间［a,b］分成n个相等的子区间［xi,x

i+1］,这里步长在每一个子区间［xi,x

i+1］上使用梯形公式,则2/6/202339

复化梯形公式积分法2/6/202340相加后得

(7―18)(7―19)若f″(x)在［a,b］上连续,由连续函数的介值定理,存在某一ξ∈(a,b)使得2/6/202341因而

于是得到复合梯形公式(7―21)其余项为2/6/202342例2若用复合梯形公式计算积分

则当0＜x＜1时,有因为又解由余项(7―21)式问积分区间要等分多少才能保证有五位有效数字?2/6/202343由于原积分的准确值具有一位整数,因此要使近似积分值有五位有效数字,只需取n满足

两边取对数得整理后得到

取n=68.2/6/202344类似复合梯形公式的做法,把区间［a,b］分成n个相等的子区间［x2i,x2i+2］(i=0,1,…,n-1),设每个子区间上的中点为x2i+1(i=0,1,…,n-1),且(7―22)2.2复合抛物线公式在每一个子区间［x2i,x2i+2

］上利用抛物线公式得2/6/202345

复化Simpson公式积分法2/6/202346相加后得(7―23)2/6/202347若f(4)(x)在［a,b］上连续,则从而得到复合抛物线公式(7―24)其余项为

(7―25)2/6/202348图7.4复合抛物线公式框图2/6/202349的数据表,分别用复合梯形公式和复合抛物线公式计算例3已知函数

x

f(x)011/80.99739782/80.98961583/80.97672674/80.95885105/80.93615566/80.90885167/80.877192510.84147092/6/202350解用复合梯形公式,这里2/6/202351用复合抛物线公式可得比较上面两个结果T8和S4,它们都需要提供9个点上的函数值工作量基本相同,然而精度却差别很大.同积分的准确值I(f)=0.9460831比较,复化梯形法的结果T8=0.9456909只有两位有效数字,而复化Simpson法的结果S4=0.9460832却有六位有效数字.2/6/202352

复化梯形公式：在每个上用梯形公式：=

Tn/*中值定理*/2/6/202353

复化Simpson公式：2/6/2023542.3变步长公式前面介绍的复合梯形公式和复合抛物线公式的步长都是预先确定的。它的主要缺点是事先很难估计出n的大小(或步长h的大小),使结果达到预先给定的精度。在实际计算中,我们常常借助于计算机来完成积分步长h的自动选择,即采用变步长求积公式。具体地讲,就是将步长逐次折半,反复利用复合求积公式,直到满足精度要求为止。

2/6/202355下面介绍变步长复合抛物线公式(变步长复合梯形公式留给读者作为练习)。

(7―26)其中再把每个子区间分成两半,用逐次将区间［a,b］分成2,4,…,2m等分,并按复合抛物线公式逐次计算积分得到S1,S2,…,Sm,而2/6/202356作步长,按复合抛物线公式计算出积分的近似值S2m。对于相邻两次的积分近似值Sm、S2m,考察当｜S2m｜＜1当｜S2m｜≥1(7―27)设预先给定的精度为ε,若｜d｜＜ε则以S2m作为所要求的积分近似值,否则继续将区间分半,利用复合抛物线公式求积分,直到满足预给的精度为止。2/6/202357图7.5变步长复合抛物线公式2/6/202358

图7.5

变步长复合抛物线公式2/6/202359§3龙贝格(Romberg)积分方法我们已经知道,当被积函数f(x)在区间［a,b］上连续时,要使得复合梯形公式或复合抛物线公式比较精确地代替定积分可将分点(即基点)加密,也就是将区间［a,b］细分,然后利用复合梯形公式或复合抛物线公式求积。2/6/202360若用Tm表示把［a,b］作m等分并按复合梯形公式求积的结果,将每一小段再对分,令新的小段的长h′=h/2,则T2m与Tm之间有如下关系:(5―28)

其中

2/6/202361另外,若用Sm表示把［a,b］分成m(偶数)个小段按复合抛物线公式计算的结果,那么只要把Sm中的m改为2m,h改为h′就有从Tm的定义可得到关系式(5―29)2/6/202362我们再举一个计算上半单位圆面积的例子(它的准确面积为π/2)。现用内接正多边形的逼近方法来计算。如图5.6,图(a)、图(b)是用同样的内接正多边形计算上半单位圆的面积。图(a)是用梯形方法计算其面积,图(b)是用三角形方法计算其面积。2/6/202363图5.62/6/202364设正多边形边数为n=2k,则由图(b)利用三角形公式算得面积为同理2/6/202365如果组合一下,就会得到更精确的结果,即同理2/6/202366再以类似方法组合得这样继续下去,其值越来越接近上半单位圆面积π/2。这种方法可以用到计算定积分2/6/202367为了推广公式(5―29)和上述计算上半单位圆面积的组合方法,我们引进龙贝格求积算法。龙贝格求积算法本来是利用所谓外推法构造出的一种计算积分的方法。为了避免从外推引入而带来理论上的麻烦,我们将直接从构造一个T数表开始。首先将［a,b］依次作20,21,22,…等分,记2/6/202368按复合梯形公式(5―20)算得的值相应地记为T(