第5章-数值积分-计算方法-《代码优化》课件

第5章数值积分§1牛顿―柯特斯(Newton―Cotes)公式§2复合求积公式§3龙贝格(Romberg)积分方法第5章数值积分§1牛顿―柯特斯(Newton

1.1牛顿―柯特斯(Newton―Cotes)公式

在一元函数的积分学中,我们已经熟知,若函数f(x)在区间［a,b］上连续且其原?函数为F(x),则可用牛顿―莱布尼兹公式

(5―1)

1.1牛顿―柯特斯(Newton―Cotes)公式来求定积分。公式(5―1)虽然在理论上或在解决实际问题中都起了很大的作用,?但它并不能完全解决定积分的计算问题。因为定积分的计算常常会碰到以下三种?情况:(1)被积函数f(x)的原函数F(x)不易找到。许多很简单的函?数,例如等,其原函数都不能用初等函数表示成有限形式。来求定积分。公式(5―1)虽然在理论(2)被积函数f(x)没有具体的解析表达式。其函数关系由表格或图形表?示,无法求出原函数。(3)尽管f(x)的原函数能表示成有限形式但其表达式相当复杂。例如定?积分的被积函数的原函数就比较复杂,从数值计算角度来?看,计算量太大。

(2)被积函数f(x)没有具体的解析表如图5.1,若用?左矩形近似地代替曲边梯形,则得到左矩形公式

同样可得到右矩形公式(5―2)(5―3)如图5.1,若用?左矩形近似地代替曲边图5.1图5.1如图5.2,若用梯形的面积近似地代替曲边梯形的面积,则得到计算定积分?的梯形公式(5―4)如图5.3,若用抛物线代替曲线f(x),则可得到抛物线公式(?或辛普生公式)(5―5)如图5.2,若用梯形的面积近似地代图5.2图5.3图5.2图5.31.1牛顿―柯特斯(Newton―Cotes)公式建立数值积分公式最基本的思想是选取一个既简单又有足够精度的函数φ(x),?用φ(x)代替被积函数f(x),于是有现用第四章介绍的插值多项式Pn(x)来代替被积函数f(x),即有

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

1.1牛顿―柯特斯(Newto利用拉格朗日插值多项式(5―6)其中(5―7)利用拉格朗日插值多项式(5―6)其中(5―7这里yi=f(xi),对式(5―6)两边积分得这里yi=f(xi),对式(5―6)两边积分得为牛顿―柯特斯求积公式,Rn(f)为牛顿―柯特斯求积公式的余项。令x=x0+sh,0≤s≤n(5―8)(5―9)(5―10)我们称为牛顿―柯特斯求积公式,Rn(f(5―11)(5―11)称C(n)i为柯特斯求积系数。很显然,当n=1时,可算得此时式(5―10)为(5―12)称C(n)i为柯特斯求积系数。此时式(5―10这是梯形公式。当n=2时,可得于是(5―13)这是梯形公式。于是(5―13)这是抛物线公式。当n=3时,这是抛物线公式。代入(5―10)式得到求积公式(5―14)类似地可分别求出n=4,5,…时的柯特斯系数,从而建立相应的求积公式。具体结果见表5―1。从表中可以看出,当n≤7时,柯特斯系数为正;从n≥8开始,柯特斯系数有正有负。因此,当n≥8时,误差有可能传播扩大,牛顿―柯特斯求积公式不宜采用。代入(5―10)式得到求积公式(5―14)柯特斯系数C(n)i仅与n和i有关,与被积函数f(x)无关,且满足(5―15)事实上,式(5―10)对f(x)=1是准确成立的。柯特斯系数C(n)i仅与n和i有关,例1试分别用梯形公式和抛物线公式计算积分解利用梯形公式利用抛物线公式例1试分别用梯形公式和抛物线公式计算积原积分的准确值原积分的准确值表5―1表5―11.2误差估计现对牛顿―柯特斯求积公式所产生的误差作一个分析。由式(5―9),牛顿―柯特斯求积公式的余项为易知,牛顿―柯特斯求积公式(5―10)对任何不高于n次的多项式是准确成立的。这是因为f(n+1)(ξ)≡0故Rn(f)≡01.2误差估计一般说来,若某个求积公式对于次数不高于m的多项式都准确成立(即Rn(f)≡0),而对于某一次数为m+1的多项式并不准确成立(即Rn(f)0),则称这一求积公式的代数精确度为m。

牛顿―柯特斯求积公式的代数精确度至少为n。通常在基点个数相等的情况下,代数精确度愈高,求积公式就愈精确。定理1(梯形公式的误差)设f(x)在区间［a,b］上具有连续的二阶导数,则梯形求积公式的误差为

一般说来,若某个求积公式对于次数不高由于ω1(x)=(x-a)(x-b)

证由式(5―9)知,梯形公式的余项为(5―16)由于证由式(5―9)知,梯形公式的余项为(在区间(a,b)内不变号,f″(ξ)是x的函数且在［a,b］上连续,故根据积分第二中值定理参见有关《数学分析》教材中“一元函数积分学第二中值定理”。知,存在某一η∈(a,b)使在区间(a,b)内不变号,f″(定理2(抛物线公式的误差)设f(x)在［a,b］上有连续的四阶导数,则抛物线公式的误差为(5―17)证由式(5―9)知定理2(抛物线公式的误差)设f(§2复合求积公式2.1复合梯形公式对于定积分(5―1),将积分区间［a,b］分成n个相等的子区间［xi,xi+1］,这里步长在每一个子区间［xi,xi+1］上使用梯形公式,则§2复合求积公式2.1复合梯相加后得(5―18)(5―19)若f″(x)在［a,b］上连续,由连续函数的介值定理,存在某一ξ∈(a,b)使得相加后得(5―18)(5―19)因而于是得到复合梯形公式(5―21)其余项为因而于是得到复合梯形公式(5―21)其余例2若用复合梯形公式计算积分问积分区间要等分多少才能保证有五位有效数字解由余项(5―21)式则当0＜x＜1时,有因为又故例2若用复合梯形公式计算积分则当由于原积分的准确值具有一位整数,因此要使近似积分值有五位有效数字,只需取n满足两边取对数得整理后得到由于原积分的准确值具有一位整数,因此2.2复合抛物线公式类似复合梯形公式的做法,把区间［a,b］分成n个相等的子区间［x2i,x2i+2］(i=0,1,…,n-1),设每个子区间上的中点为x2i+1(i=0,1,…,n-1),且

在每一个子区间［x2i,x2i+2］上利用抛物线公式得(5―22)2.2复合抛物线公式(5―22)相加后得(5―23)相加后得(5―23)图5.4图5.4图5.4图5.4若f(4)(x)在［a,b］上连续,则从而得到复合抛物线公式(5―24)若f(4)(x)在［a,b］上连续,则其余项为(5―25)复合抛物线公式的计算框图见5.4。例3根据给出的函数的数据表5―2,分别用复合梯形公式和复合抛物线公式计算其余项为(5―25)复合抛物线公式的计算框图见表5―2表5―2解用复合梯形公式,这里解用复合梯形公式,这里用复合抛物线公式可得而I的准确值为0.9460831…,可见用复合抛物线公式比用复合梯形公式精确。用复合抛物线公式可得而I的准确2.3变步长公式前面介绍的复合梯形公式和复合抛物线公式的步长都是预先确定的。它的主要缺点是事先很难估计出n的大小(或步长h的大小),使结果达到预先给定的精度。在实际计算中,我们常常借助于计算机来完成积分步长h的自动选择,即采用变步长求积公式。具体地讲,就是将步长逐次折半,反复利用复合求积公式,直到满足精度要求为止。2.3变步长公式下面介绍变步长复合抛物线公式(变步长复合梯形公式留给读者作为练习)。逐次将区间［a,b］分成21,22,…,2m等分,并按复合抛物线公式逐次计算积分得到S1,S2,…,Sm,而(5―26)其中再把每个子区间分成两半,用下面介绍变步长复合抛物线公式(变步长图5.5图5.5图5.5图5.5作步长,按复合抛物线公式计算出积分的近似值S2m。对于相邻两次的积分近似值Sm、S2m,考察当｜S2m｜＜1当｜S2m｜≥1(5―27)设预先给定的精度为ε,若｜d｜＜ε则以S2m作为所要求的积分近似值,否则继续将区间分半,利用复合抛物线公式求积分,直到满足预给的精度为止。作步长,按复合抛物线公式计算出积分的§3龙贝格(Romberg)积分方法我们已经知道,当被积函数f(x)在区间［a,b］上连续时,要使得复合梯形公式或复合抛物线公式比较精确地代替定积分可将分点(即基点)加密,也就是将区间［a,b］细分,然后利用复合梯形公式或复合抛物线公式求积。§3龙贝格(Romberg)积分方法若用Tm表示把［a,b］作m等分并按复合梯形公式求积的结果,将每一小段再对分,令新的小段的长h′=h/2,则T2m与Tm之间有如下关系:(5―28)其中若用Tm表示把［a,b］作m等分并按另外,若用Sm表示把［a,b］分成m(偶数)个小段按复合抛物线公式计算的结果,那么只要把Sm中的m改为2m,h改为h′就有从Tm的定义可得到关系式(5―29)另外,若用Sm表示把［a,b］分成m我们再举一个计算上半单位圆面积的例子(它的准确面积为π/2)。现用内接正多边形的逼近方法来计算。如图5.6,图(a)、图(b)是用同样的内接正多边形计算上半单位圆的面积。图(a)是用梯形方法计算其面积,图(b)是用三角形方法计算其面积。我们再举一个计算上半单位圆面积的例子图5.6图5.6设正多边形边数为n=2k,则由图(b)利用三角形公式算得面积为同理设正多边形边数为n=2k,则由图(b如果组合一下,就会得到更精确的结果,即同理如果组合一下,就会得到更精确的结果,再以类似方法组合得这样继续下去,其值越来越接近上半单位圆面积π/2。这种方法可以用到计算定积分再以类似方法组合得这样继续下为了推广公式(5―29)和上述计算上半单位圆面积的组合方法,我们引进龙贝格求积算法。龙贝格求积算法本来是利用所谓外推法构造出的一种计算积分的方法。为了避免从外推引入而带来理论上的麻烦,我们将直接从构造一个T数表开始。首先将［a,b］依次作20,21,22,…等分,记为了推广公式(5―29)和上述计算上按复合梯形公式(5―20)算得的值相应地记为T(k)0(k=0,1,2,…);把按式(5―29)算得的S2m依次记为T(k)1(k=0,1,2,崐…),而这每一个S2m又理解为由T2m与Tm的线性组合得到的改进值,即我们可按照类似的方法继续进行改进,也即由S2m与Sm的线性组合得到改进值,依次记为T(k)2(k=0,1,2,…),即

按复合梯形公式(5―20)算得的值这样就可构造出一个数表(5-30)这样就可构造出一个数表(5-30)其中除第0列(即最左一列)的T(k)0是按复合梯形公式计算外,其余各列都按下述规则(对m)(5―31)递推地计算出来。箭头表示计算流程。其计算步骤为:(1)将区间［a,b］等分为20,用梯形公式计算T(0)0,即其中除第0列(即最左一列)的T(k)(2)将区间［a,b］等分为21,用梯形公式算出T(1)0,即再由T(0)0,T(1)0根据公式(5―31)算出T(0)1,即若｜T(0)1-T(0)0｜＜ε,(ε为预给的精度)则停止计算;否则继续往下计算;(2)将区间［a,b］等分为21,用梯形公(3)依次分别算出T(2)0,T(1)1,T(0)2,…,这一行地往下推算,每一行算完,就得验证T(0)m(m=1,2,…)是否满足预给的精度,即若则停止计算;否则继续进行下一行。为了便于在计算机上实现,可运用下列公式编制程序:(3)依次分别算出T(2)0,T(1第5章--数值积分-计算方法-《代码优化》课件图5.7图5.7图5.7图5.7例4计算积分精确到10-4。解例4计算积分精确到10-4。第5章--数值积分-计算方法-《代码优化》课件第5章--数值积分-计算方法-《代码优化》课件第5章--数值积分-计算方法-《代码优化》课件于是由于实际上于是由于实际上第5章数值积分§1牛顿―柯特斯(Newton―Cotes)公式§2复合求积公式§3龙贝格(Romberg)积分方法第5章数值积分§1牛顿―柯特斯(Newton

1.1牛顿―柯特斯(Newton―Cotes)公式

在一元函数的积分学中,我们已经熟知,若函数f(x)在区间［a,b］上连续且其原?函数为F(x),则可用牛顿―莱布尼兹公式

(5―1)

1.1牛顿―柯特斯(Newton―Cotes)公式来求定积分。公式(5―1)虽然在理论上或在解决实际问题中都起了很大的作用,?但它并不能完全解决定积分的计算问题。因为定积分的计算常常会碰到以下三种?情况:(1)被积函数f(x)的原函数F(x)不易找到。许多很简单的函?数,例如等,其原函数都不能用初等函数表示成有限形式。来求定积分。公式(5―1)虽然在理论(2)被积函数f(x)没有具体的解析表达式。其函数关系由表格或图形表?示,无法求出原函数。(3)尽管f(x)的原函数能表示成有限形式但其表达式相当复杂。例如定?积分的被积函数的原函数就比较复杂,从数值计算角度来?看,计算量太大。

(2)被积函数f(x)没有具体的解析表如图5.1,若用?左矩形近似地代替曲边梯形,则得到左矩形公式

同样可得到右矩形公式(5―2)(5―3)如图5.1,若用?左矩形近似地代替曲边图5.1图5.1如图5.2,若用梯形的面积近似地代替曲边梯形的面积,则得到计算定积分?的梯形公式(5―4)如图5.3,若用抛物线代替曲线f(x),则可得到抛物线公式(?或辛普生公式)(5―5)如图5.2,若用梯形的面积近似地代图5.2图5.3图5.2图5.31.1牛顿―柯特斯(Newton―Cotes)公式建立数值积分公式最基本的思想是选取一个既简单又有足够精度的函数φ(x),?用φ(x)代替被积函数f(x),于是有现用第四章介绍的插值多项式Pn(x)来代替被积函数f(x),即有

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

1.1牛顿―柯特斯(Newto利用拉格朗日插值多项式(5―6)其中(5―7)利用拉格朗日插值多项式(5―6)其中(5―7这里yi=f(xi),对式(5―6)两边积分得这里yi=f(xi),对式(5―6)两边积分得为牛顿―柯特斯求积公式,Rn(f)为牛顿―柯特斯求积公式的余项。令x=x0+sh,0≤s≤n(5―8)(5―9)(5―10)我们称为牛顿―柯特斯求积公式,Rn(f(5―11)(5―11)称C(n)i为柯特斯求积系数。很显然,当n=1时,可算得此时式(5―10)为(5―12)称C(n)i为柯特斯求积系数。此时式(5―10这是梯形公式。当n=2时,可得于是(5―13)这是梯形公式。于是(5―13)这是抛物线公式。当n=3时,这是抛物线公式。代入(5―10)式得到求积公式(5―14)类似地可分别求出n=4,5,…时的柯特斯系数,从而建立相应的求积公式。具体结果见表5―1。从表中可以看出,当n≤7时,柯特斯系数为正;从n≥8开始,柯特斯系数有正有负。因此,当n≥8时,误差有可能传播扩大,牛顿―柯特斯求积公式不宜采用。代入(5―10)式得到求积公式(5―14)柯特斯系数C(n)i仅与n和i有关,与被积函数f(x)无关,且满足(5―15)事实上,式(5―10)对f(x)=1是准确成立的。柯特斯系数C(n)i仅与n和i有关,例1试分别用梯形公式和抛物线公式计算积分解利用梯形公式利用抛物线公式例1试分别用梯形公式和抛物线公式计算积原积分的准确值原积分的准确值表5―1表5―11.2误差估计现对牛顿―柯特斯求积公式所产生的误差作一个分析。由式(5―9),牛顿―柯特斯求积公式的余项为易知,牛顿―柯特斯求积公式(5―10)对任何不高于n次的多项式是准确成立的。这是因为f(n+1)(ξ)≡0故Rn(f)≡01.2误差估计一般说来,若某个求积公式对于次数不高于m的多项式都准确成立(即Rn(f)≡0),而对于某一次数为m+1的多项式并不准确成立(即Rn(f)0),则称这一求积公式的代数精确度为m。

牛顿―柯特斯求积公式的代数精确度至少为n。通常在基点个数相等的情况下,代数精确度愈高,求积公式就愈精确。定理1(梯形公式的误差)设f(x)在区间［a,b］上具有连续的二阶导数,则梯形求积公式的误差为

一般说来,若某个求积公式对于次数不高由于ω1(x)=(x-a)(x-b)

证由式(5―9)知,梯形公式的余项为(5―16)由于证由式(5―9)知,梯形公式的余项为(在区间(a,b)内不变号,f″(ξ)是x的函数且在［a,b］上连续,故根据积分第二中值定理参见有关《数学分析》教材中“一元函数积分学第二中值定理”。知,存在某一η∈(a,b)使在区间(a,b)内不变号,f″(定理2(抛物线公式的误差)设f(x)在［a,b］上有连续的四阶导数,则抛物线公式的误差为(5―17)证由式(5―9)知定理2(抛物线公式的误差)设f(§2复合求积公式2.1复合梯形公式对于定积分(5―1),将积分区间［a,b］分成n个相等的子区间［xi,xi+1］,这里步长在每一个子区间［xi,xi+1］上使用梯形公式,则§2复合求积公式2.1复合梯相加后得(5―18)(5―19)若f″(x)在［a,b］上连续,由连续函数的介值定理,存在某一ξ∈(a,b)使得相加后得(5―18)(5―19)因而于是得到复合梯形公式(5―21)其余项为因而于是得到复合梯形公式(5―21)其余例2若用复合梯形公式计算积分问积分区间要等分多少才能保证有五位有效数字解由余项(5―21)式则当0＜x＜1时,有因为又故例2若用复合梯形公式计算积分则当由于原积分的准确值具有一位整数,因此要使近似积分值有五位有效数字,只需取n满足两边取对数得整理后得到由于原积分的准确值具有一位整数,因此2.2复合抛物线公式类似复合梯形公式的做法,把区间［a,b］分成n个相等的子区间［x2i,x2i+2］(i=0,1,…,n-1),设每个子区间上的中点为x2i+1(i=0,1,…,n-1),且

在每一个子区间［x2i,x2i+2］上利用抛物线公式得(5―22)2.2复合抛物线公式(5―22)相加后得(5―23)相加后得(5―23)图5.4图5.4图5.4图5.4若f(4)(x)在［a,b］上连续,则从而得到复合抛物线公式(5―24)若f(4)(x)在［a,b］上连续,则其余项为(5―25)复合抛物线公式的计算框图见5.4。例3根据给出的函数的数据表5―2,分别用复合梯形公式和复合抛物线公式计算其余项为(5―25)复合抛物线公式的计算框图见表5―2表5―2解用复合梯形公式,这里解用复合梯形公式,这里用复合抛物线公式可得而I的准确值为0.9460831…,可见用复合抛物线公式比用复合梯形公式精确。用复合抛物线公式可得而I的准确2.3变步长公式前面介绍的复合梯形公式和复合抛物线公式的步长都是预先确定的。它的主要缺点是事先很难估计出n的大小(或步长h的大小),使结果达到预先给定的精度。在实际计算中,我们常常借助于计算机来完成积分步长h的自动选择,即采用变步长求积公式。具体地讲,就是将步长逐次折半,反复利用复合求积公式,直到满足精度要求为止。2.3变步长公式下面介绍变步长复合抛物线公式(变步长复合梯形公式留给读者作为练习)。逐次将区间［a,b］分成21,22,…,2m等分,并按复合抛物线公式逐次计算积分得到S1,S2,…,Sm,而(5―26)其中再把每个子区间分成两半,用下面介绍变步长复合抛物线公式(变步长图5.5图5.5图5.5图5.5作步长,按复合抛物线公式计算出积分的近似值S2m。对于相邻两次的积分近似值Sm、S2m,考察当｜S2m｜＜1当｜S2m｜≥1(5―27)设预先给定的精度为ε,若｜d｜＜ε则以S2m作为所要求的积分近似值,否则继续将区间分半,利用复合抛物线公式求积分,直到满足预给的精度为止。作步长,按复合抛物线公式计算出积分的§3龙贝格(Romberg)积分方法我们已经知道,当被积函数f(x)在区间［a,b］上连续时,要使得复合梯形公式或复合抛物线公式比较精确地代替定积分可将分点(即基点)加密,也就是将区间［a,b］细分,然后利用复合梯形公式或复合抛物线公式求积。§3龙贝格(Romberg)积分方法若用Tm表示把［a,b］作m等分并按复合梯形公式求积的结果,将每一小段再对分,令新的小段的长h′=h/2,则T2m与Tm之间有如下关系:(5―28)其中若用Tm表示把［a,b］作m等分并按另外,若用Sm表示把［a,b］分成m(偶数)个小段按复合抛物线公式计算的结果,那么只要把Sm中的m改为2m,h改为h′就有从Tm的定义可得到关系式(5―29)另外,若用Sm表示把［a,b］分成m我们再举一个计算上半单位圆面积的例子(它的准确面积为π/2)。现用内接正多边形的逼近方法来计算。如图5.6,图(a)、图(b)是用同样的内接正多边形计算上半单位圆的面积。图(a)是用梯形方法计算其面积,图(b)是用三角形方法计算其面积。我们再举一个计算上半单位圆面积的例子图5.6图5.6设正多边形边数为n=2k,则由图(b)利用三角形公式算得面积为同理