《数值积分与数值微分.pptx》由会员分享,可在线阅读,更多相关《数值积分与数值微分.pptx(152页珍藏版)》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、进行计算,但在工程计算和科学研究中,经常会遇进行计算,但在工程计算和科学研究中,经常会遇到被积函数到被积函数f(x)的下列一些情况:的下列一些情况:的原函数的原函数对定积分对定积分的被积函数的被积函数已知,在高等数学中可用牛顿已知,在高等数学中可用牛顿莱布尼兹公式莱布尼兹公式4.1 数值积分概论数值积分概论 实际问题当中常常要计算积分,有些数值方法,实际问题当中常常要计算积分,有些数值方法,如微分方程和积分方程的求解,也都和积分计算相如微分方程和积分方程的求解,也都和积分计算相联系联系.4.1.1 数值求积的基本思想数值求积的基本思想第1页/共152页(4)f(x)本身没有解析表达式,其函数关
2、系由表格本身没有解析表达式,其函数关系由表格或图形给出,列如为实验或测量数据或图形给出,列如为实验或测量数据.(2)f(x)的原函数不能用初等函数形式表示,例如的原函数不能用初等函数形式表示,例如(3)f(x)的原函数虽然可用初等函数形式表示,但的原函数虽然可用初等函数形式表示,但其原函数表示形式相当复杂,例如其原函数表示形式相当复杂,例如(1)f(x)复杂,求原函数困难,列如复杂,求原函数困难,列如第2页/共152页 以上的以上的 4种情况都不能用牛顿种情况都不能用牛顿莱布尼兹公莱布尼兹公式方便地计算该函数的定积分,满足不了实际需式方便地计算该函数的定积分,满足不了实际需要,因此,有必要研究
3、定积分的数值计算问题;要,因此,有必要研究定积分的数值计算问题;另外,对一些函数的求导问题,其求导、微分也另外,对一些函数的求导问题,其求导、微分也相当复杂,也有必要研究求导、微分的数值计算相当复杂,也有必要研究求导、微分的数值计算问题。本章主要介绍数值求积分和数值求微分的问题。本章主要介绍数值求积分和数值求微分的方法。方法。第3页/共152页 由积分中值定理由积分中值定理,对连续函数对连续函数f(x),在区间在区间a,b内至少存在一点内至少存在一点,使,使只要对平均高只要对平均高度度 f()提供一提供一种种近似算法近似算法,便便可相应地获得可相应地获得一种一种数值求积数值求积方法方法.即所谓
4、即所谓矩矩形公式形公式.第4页/共152页 例如例如,用区间用区间a,b两端点的函数值两端点的函数值 f(a)与与f(b)的的算术平均值作为算术平均值作为f()的近似值的近似值,可导出可导出求积公式求积公式这便是人们所熟知的这便是人们所熟知的梯形公式梯形公式.第5页/共152页 如果改用区间如果改用区间a,b的中点的中点 c=(a+b)/2 处的函数值处的函数值f(c)近似代替近似代替f(),则又可导出所谓则又可导出所谓(中中)矩形公式矩形公式第6页/共152页 一般地一般地,在区间在区间a,b上适当选取点上适当选取点xk(k=0,1,n),然后用然后用 f(xk)的的加权平均值加权平均值作为
5、作为f()的近似值的近似值,可得到可得到更为更为一般的求积公式一般的求积公式 其中:点其中:点xk叫叫求积节点求积节点,系数系数Ak叫叫求积系数求积系数.Ak仅与节仅与节点点xk的选取有关的选取有关,而与被积函数而与被积函数 f(x)无关无关.求积公式的求积公式的截断误差截断误差为为 R(f)又称为又称为求积余项求积余项.这类数值积分方法通常称为这类数值积分方法通常称为机械求积机械求积,其特点是,其特点是将积分求值问题归结为函数值的计算,这就避开了牛将积分求值问题归结为函数值的计算,这就避开了牛-莱公式寻求原函数的困难莱公式寻求原函数的困难.第7页/共152页4.1.2 代数精度的概念代数精度
6、的概念 定义定义1 如果求积公式如果求积公式(1)对所有次数不超过对所有次数不超过m的多项式都精确成立;的多项式都精确成立;(2)至少对一个至少对一个m+1次多项式不精确成立,次多项式不精确成立,则称则称该公式具有该公式具有m次代数精度次代数精度(或或代数精确度代数精确度).数值求积方法的近似方法,为要保证精度,我数值求积方法的近似方法,为要保证精度,我们自然希望求积公式能对们自然希望求积公式能对“尽可能多尽可能多”的函数准确的函数准确地成立,这就提出了所谓代数精度的概念地成立,这就提出了所谓代数精度的概念.第8页/共152页 一般来说,代数精度越高,求积公式越好。一般来说,代数精度越高,求积
7、公式越好。结论结论 一个求积公式具有一个求积公式具有m次代数精度的次代数精度的充要条充要条件件是该求积公式是该求积公式:(1)对对xk(k=0,1,m)精确成立;精确成立;(2)对对xm+1不精确成立不精确成立.故一般地,要验证一个求积公式具有故一般地,要验证一个求积公式具有m次代数次代数精度,只要令对于精度,只要令对于 f(x)=1,x,xm求积公式精确成立求积公式精确成立等式就行等式就行.第9页/共152页 即对于求积公式即对于求积公式 给定给定n+1个互异的求积节点个互异的求积节点 x0,x1,xn-1,xn,令求积公式对令求积公式对 f(x)=1,x,xn 精确成立精确成立,即得即得求
8、解该方程组即可确定求积系数求解该方程组即可确定求积系数Ak,所得到的求积公所得到的求积公式式至少具有至少具有n 次代数精度次代数精度.第10页/共152页 解解 当当 f(x)=1时时,此时公式精确成立。此时公式精确成立。例例1 验证梯形公式验证梯形公式具有一次代数精度。具有一次代数精度。当当 f(x)=x时,时,公式也精确成立公式也精确成立.当当 f(x)=x2 时,时,公式对公式对x2不精确成立不精确成立.故由定理故由定理1知知,梯形公式的代数精度为梯形公式的代数精度为1次次.第11页/共152页 例例2 确定求积公式中的待定系数,使其代数精确定求积公式中的待定系数,使其代数精度尽量高,并
9、指明求积公式所具有的代数精度度尽量高,并指明求积公式所具有的代数精度.解解 令令 f(x)=1,x,x2 代入公式两端并令其相等,得代入公式两端并令其相等,得 解得解得第12页/共152页得得求积公式求积公式为为令令 f(x)=x3,得,得令令 f(x)=x4,得,得故故求积公式求积公式具有具有3次代数精度次代数精度.第13页/共152页 如果我们事先选定求积节点如果我们事先选定求积节点xk,譬如,以区间,譬如,以区间a,b的等距分点作为节点,这时取的等距分点作为节点,这时取m=n求解方程组求解方程组即可确定求积系数即可确定求积系数Ak,而使求积公式至少具有,而使求积公式至少具有 n次次代数精
10、度代数精度.本章第本章第2节介绍这样一类求积公式,梯形节介绍这样一类求积公式,梯形公式是其中的一个特例公式是其中的一个特例.如为了构造出上面的求积公式,原则上是一个如为了构造出上面的求积公式,原则上是一个确定参数确定参数xk和和Ak的代数问题的代数问题.方程组方程组(1.4)实际上是一实际上是一2n+2个参数的非线性方程组,此方程组当个参数的非线性方程组,此方程组当n1时求时求解非常困难,但当解非常困难,但当n=0及及n=1的情形还是可以通过求的情形还是可以通过求解方程组得到相应的求积公式解方程组得到相应的求积公式.下面对下面对n=0讨论求积讨论求积公式的建立及代数精确度公式的建立及代数精确度
11、.第14页/共152页 此时求积公式为此时求积公式为其中,其中,x0及及A0为待定参数为待定参数.根据代数精确度定义可令根据代数精确度定义可令f(x)=1,x,由方程组知,由方程组知.得得得到的求积公式就是得到的求积公式就是(1.2)式的中矩形公式式的中矩形公式.再令再令f(x)=x2,代入,代入(1.4)式的第三式式的第三式第15页/共152页说明说明(1.2)式对式对f(x)=x2不精确成立,故它的代数精确不精确成立,故它的代数精确度为度为1.方程组方程组(1.4)是根据是根据(1.3)式的求积公式得到的,式的求积公式得到的,按照代数精确度的定义,如果求积公式中除了按照代数精确度的定义,如
12、果求积公式中除了f(xi)还有还有(x)在某些节点上的值,也同样可得到相应的在某些节点上的值,也同样可得到相应的求积公式求积公式.第16页/共152页 例例1 给定形如下面的求积公式,试确定系数给定形如下面的求积公式,试确定系数A0,A1,B0,,使公式具有尽可能高的代数精确度,使公式具有尽可能高的代数精确度.解解 根据题意可令根据题意可令f(x)=1,x,x2分别代入求积公式分别代入求积公式使它精确成立:使它精确成立:第17页/共152页解得解得当当f(x)=x3时,上式右端为时,上式右端为1/3,而左端是,而左端是于是有求积公式于是有求积公式故积分公式对故积分公式对f(x)=x3不精确成立
13、,其代数精确度为不精确成立,其代数精确度为2.第18页/共152页4.1.3 插值型的求积公式插值型的求积公式设给定一组节点设给定一组节点且已知且已知f(x)在这些节点上的函数值在这些节点上的函数值 f(xk),则可求得则可求得f(x)的拉格朗日插值多项式的拉格朗日插值多项式(因为因为Ln(x)的原函数易求的原函数易求)其中其中lk(x)为插值基函数为插值基函数,取取由上式确定系数的公式称为由上式确定系数的公式称为插值型求积公式插值型求积公式.即即则则 f(x)Ln(x)第19页/共152页插值型求积公式积分法几何表示第20页/共152页由插值余项定理由插值余项定理,其求积余项为其求积余项为其
14、中其中=(x)如果求积公式如果求积公式(1.5)是插值型的,按照插值余是插值型的,按照插值余项式子,对于次数不超过项式子,对于次数不超过n的多项式的多项式f(x),其余项,其余项 R(f)等于零,因而等于零,因而这时求积公式至少具有这时求积公式至少具有n次代数次代数精度精度.第21页/共152页 反之,如果求积公式至少具有反之,如果求积公式至少具有n次代数精度,次代数精度,则它必定是插值型的则它必定是插值型的.事实上,这时求积公式对于事实上,这时求积公式对于插值基函数插值基函数 lk(x)应准确成立,即有应准确成立,即有注意到注意到lk(xj)=kj,上式右端实际上即等于,上式右端实际上即等于
15、Ak,因,因而下面式子成立而下面式子成立.第22页/共152页 定理定理1 具有具有n+1个节点的数值求积公式个节点的数值求积公式(1.5)是是插值型求积公式插值型求积公式的的充要条件充要条件为为:该公式该公式至少具有至少具有n次代数精度次代数精度.综上所述,我们有结论为综上所述,我们有结论为 这时令这时令f(x)=1代入又有结论为代入又有结论为 结论结论 对插值型求积公式的系数必有对插值型求积公式的系数必有第23页/共152页4.1.4 求积公式的余项求积公式的余项若求积公式若求积公式(1.3)的代数精确度为的代数精确度为m,则由求积公,则由求积公式余项的表达式式余项的表达式(1.7)可以证
16、明余项形如可以证明余项形如其中其中K为不依赖于为不依赖于f(x)的待定参数的待定参数.这个结果表明当这个结果表明当f(x)是次数小于等于是次数小于等于m的多项式时,由于的多项式时,由于f(m+1)(x)=0,故此时故此时R f=0,即求积公式,即求积公式(1.3)精确成立精确成立.而当而当f(x)=xm+1时,时,f(m+1)(x)=(m+1)!,(1.8)式左端式左端 Rf 0,故可求得,故可求得第24页/共152页代入余项公式代入余项公式(1.8)式可以得到更细致的余项表达式式可以得到更细致的余项表达式.例如梯形公式例如梯形公式(1.1)的代数精确度为的代数精确度为1,可以证明,可以证明它
17、的余项表达式为它的余项表达式为其中其中于是得到梯形公式于是得到梯形公式(1.1)的余项为的余项为第25页/共152页对中矩形公式对中矩形公式(1.2),其代数精确度为,其代数精确度为1,可以证,可以证明它的余项表达式为明它的余项表达式为其中其中于是得到中矩形公式于是得到中矩形公式(1.2)的余项为的余项为第26页/共152页 例例2 求例求例1中求积公式中求积公式的余项的余项.解解 由于此求积公式的代数精确度为由于此求积公式的代数精确度为2,故余项,故余项表达式为表达式为R=K(),令,令f(x)=x3,得,得()=3!,于是有,于是有故得故得第27页/共152页其中其中h=max(xi-xi
18、-1),则称求积公式,则称求积公式(1.3)是是收敛的收敛的.4.1.5 求积公式的收敛性与稳定性求积公式的收敛性与稳定性 定义定义2 在求积公式在求积公式(1.3)中,若中,若 在求积公式在求积公式(1.3)中,由于计算中,由于计算f(xk)可能产生误可能产生误差差k,实际得到,实际得到 ,即,即 .记记如果对任给小正数如果对任给小正数 0,只要误差,只要误差|k|充分小就有充分小就有它表明求积公式它表明求积公式(1.3)计算是计算是稳定的稳定的,由此给出,由此给出第28页/共152页 定义定义3 对任给小正数对任给小正数 0,若存在,若存在 0,只要,只要 就有就有(1.12)式成立,则式
19、成立,则称求积公式称求积公式(1.3)是是稳定的稳定的.证明证明 对任给对任给 0,若取若取=/(b-a),对所有对所有k都有都有故求积公式故求积公式(1.3)是稳定的是稳定的.证毕证毕.定理定理2 若求积公式若求积公式(1.3)中所有系数中所有系数Ak0,则此,则此求积公式是稳定的求积公式是稳定的.则有则有定理定理2表明只要求积系数表明只要求积系数Ak0,就能保证计算的,就能保证计算的稳定性稳定性.第29页/共152页 4.2 牛顿牛顿柯特斯公式柯特斯公式 为便于上机计算,通常在内插求积公式中我们为便于上机计算,通常在内插求积公式中我们通常取等距节点,即将积分区间通常取等距节点,即将积分区间
20、a,b划分划分n等分,等分,即令步长即令步长h=(b-a)/n,且记,且记x0=a,xn=b,则节点记为,则节点记为xk=x0+kh(k=0,1,n),然后作变换,然后作变换:t=(x-x0)/h,代入代入求积系数公式,将会简化计算求积系数公式,将会简化计算.第30页/共152页4.2.1 柯特斯系数与辛普森公式柯特斯系数与辛普森公式设将积分区间设将积分区间a,b划分成划分成 n等分等分,步长步长h=求积节点取为求积节点取为xk=a+kh(k=0,1,n),由此构造插值型由此构造插值型求积公式求积公式,则其求积系数为则其求积系数为引入变换引入变换 x=a+th,则有则有(k=0,1,n)(k=
21、0,1,n)第31页/共152页其中其中记记于是得求积公式于是得求积公式称为称为n 阶牛顿阶牛顿-柯特斯柯特斯(Newton-Cotes)公式公式.显然显然,柯特斯系数与被积函数柯特斯系数与被积函数 f(x)和积分区间和积分区间a,b无关无关,且为容易计算的多项式积分且为容易计算的多项式积分.称为称为柯特斯系数柯特斯系数.第32页/共152页n11/21/221/64/61/631/83/83/81/847/9032/90 12/90 32/907/90519/28875/28850/28850/28875/28819/288641/840216/84027/840272/84027/8402
22、16/84041/840常用的柯特斯系数表常用的柯特斯系数表第33页/共152页 当当n=1时,时,柯特斯系数柯特斯系数为为这时的这时的牛顿牛顿-柯特斯公式柯特斯公式为一阶求积公式,就是我们为一阶求积公式,就是我们所熟悉的所熟悉的梯形公式梯形公式,即,即第34页/共152页 当当n=2时,时,柯特斯系数柯特斯系数为为相应的相应的牛顿牛顿-柯特斯公式柯特斯公式为二阶求积公式,就是为二阶求积公式,就是辛普辛普森森(simpson)公式公式(又称为又称为抛物形求积公式抛物形求积公式),即,即第35页/共152页式中式中(k=0,1,2,3,4),h=(b-a)/4.n=4 时的时的牛顿牛顿-柯特斯公
23、式柯特斯公式就特别称为就特别称为柯特斯公柯特斯公式式.其形式是其形式是 在在柯特斯系数表柯特斯系数表中中(见书见书p104)看到看到n 7时,时,柯特柯特斯系数斯系数出现负值,于是有出现负值,于是有第36页/共152页特别地,假定特别地,假定则有则有这表明在这表明在b-a1时,初始误差将会引起计算结果误差时,初始误差将会引起计算结果误差增大,即计算不稳定,故增大,即计算不稳定,故n 7的的牛顿牛顿-柯特斯公式柯特斯公式是是不用的不用的.第37页/共152页4.2.2 偶阶求积公式的代数精度偶阶求积公式的代数精度 作为插值型求积公式,作为插值型求积公式,n阶阶牛顿牛顿-柯特斯公式柯特斯公式至至少
24、具有少具有n次代数精度次代数精度(推论推论1).实际的代数精度能否进实际的代数精度能否进一步提高呢?一步提高呢?先看先看辛普森公式辛普森公式,它是二阶,它是二阶牛顿牛顿-柯特斯公式柯特斯公式,因此至少具有二次代数精度因此至少具有二次代数精度.进一步用进一步用f(x)=x3进行检进行检验,按验,按辛普森公式辛普森公式计算得计算得第38页/共152页另一方面,直接求积得另一方面,直接求积得这时有这时有S=I,即,即辛普森公式辛普森公式对不超过三次的多项式均对不超过三次的多项式均能精确成立,又容易验证它对能精确成立,又容易验证它对f(x)=x4通常是不精确通常是不精确的的(如取如取a=0,b=1进行
25、验证有,进行验证有,S=5/24I=1/5),因此,因此,辛普森公式辛普森公式实际上实际上具有三次代数精度具有三次代数精度.一般地,我们可以证明下述论断:一般地,我们可以证明下述论断:第39页/共152页 定理定理3 n 阶牛顿阶牛顿-柯特斯公式的代数精度至少为柯特斯公式的代数精度至少为 证明证明 由定理由定理1已知,无论已知,无论n为奇数或偶数,插值为奇数或偶数,插值型求积公式都至少具有型求积公式都至少具有n次代数精度次代数精度.因此我们证明因此我们证明n为偶数的情形,即对为偶数的情形,即对n+1次多项式余项为零次多项式余项为零.令令n=2k,设设为任一为任一n+1次多项式,其最高次系数为次
26、多项式,其最高次系数为an+1,则它的,则它的n+1阶导数为阶导数为第40页/共152页由余项公式由余项公式有有这里变换为这里变换为x=a+th,注意,注意xj=a+jh.下面我们证明下面我们证明第41页/共152页作变换作变换u=t-k,则,则容易验证容易验证(u)为奇函数,即为奇函数,即(-u)=-(u),而奇函数,而奇函数在对称区间上的积分为零,所以在对称区间上的积分为零,所以第42页/共152页 定理定理3说明,当说明,当n为偶数时,牛顿为偶数时,牛顿-柯特斯公式柯特斯公式对不超过对不超过n+1次的多项式均能精确成立,因此,其次的多项式均能精确成立,因此,其代数精度可达到代数精度可达到
27、n+1.正是基于这种考虑,当正是基于这种考虑,当n=2k与与n=2k+1时具有相同的代数精度,因而在实用中常采时具有相同的代数精度,因而在实用中常采用用n为偶数的牛顿为偶数的牛顿-柯特斯公式,如抛物形公式柯特斯公式,如抛物形公式(n=2)等等.第43页/共152页 对牛顿对牛顿-柯特斯求积公式通常只用柯特斯求积公式通常只用n=1,2,4时的三时的三个公式,个公式,n=1时即为梯形公式时即为梯形公式(1.1),其余项为其余项为(1.10)式式.n=2时即为辛普森公式时即为辛普森公式(2.3),其代数精度为其代数精度为3,可以证,可以证明余项可表示为明余项可表示为其中其中K由由(1.9)式及式及(
28、2.3)式可得式可得4.2.3 辛普森公式的余项辛普森公式的余项第44页/共152页从而可得从而可得辛普森公式辛普森公式(2.3)的余项的余项为为也可直接积分计算得到也可直接积分计算得到 对对n=4的柯特斯公式的柯特斯公式(2.4),其代数精度为,其代数精度为5。故。故类似于求类似于求(2.3)式的余项可得到式的余项可得到柯特斯公式的余项柯特斯公式的余项为为第45页/共152页 解解:由梯形公式得由梯形公式得由辛普森公式得由辛普森公式得 例题例题 分别用梯形公式、辛普森公式和柯特斯公分别用梯形公式、辛普森公式和柯特斯公式计算积分式计算积分第46页/共152页由柯特斯公式得由柯特斯公式得积分的精
29、确值积分的精确值第47页/共152页4.3 复合求积公式复合求积公式 从求积公式的余项的讨论中我们看到,被积函数从求积公式的余项的讨论中我们看到,被积函数所用的插值多项式次数越高,对函数光滑性的要求也所用的插值多项式次数越高,对函数光滑性的要求也越高越高.另一方面,插值节点的增多另一方面,插值节点的增多(n的增大的增大),在使用,在使用牛顿牛顿-柯特斯公式时将导致求积系数出现负数柯特斯公式时将导致求积系数出现负数(当当n8时时,牛顿牛顿-柯特斯求积系数会出现负数柯特斯求积系数会出现负数),即牛顿,即牛顿-柯特柯特斯公式是不稳定的,不可能通过提高阶的方法来提高斯公式是不稳定的,不可能通过提高阶的
30、方法来提高求积精度求积精度.第48页/共152页 为了提高精度,通常在实际应用中往往采用为了提高精度,通常在实际应用中往往采用将积将积分区间划分成若干个小区间,在各小区间上采用低次分区间划分成若干个小区间,在各小区间上采用低次的求积公式的求积公式(梯形公式或抛物形公式梯形公式或抛物形公式),然后再利用积,然后再利用积分的可加性,把各区间上的积分加起来,便得到新的分的可加性,把各区间上的积分加起来,便得到新的求积公式,这就是求积公式,这就是复合求积公式复合求积公式的基本思想的基本思想.为叙述为叙述方便,我们仅讨论各小区间均采用同一低次的求积公方便,我们仅讨论各小区间均采用同一低次的求积公式的复合
31、求积公式式的复合求积公式对各小区间也可分别采用不同的对各小区间也可分别采用不同的求积公式,也可推出新的求积公式,读者可按实际问求积公式,也可推出新的求积公式,读者可按实际问题的具体情况讨论题的具体情况讨论.第49页/共152页 将积分区间将积分区间a,bn等分等分,步长步长 xk=a+kh(k=0,1,n),则由定积分性质知则由定积分性质知,分点为分点为每个子区间每个子区间上的积分上的积分用用低阶求积公式低阶求积公式,然后把所有区间的然后把所有区间的计算结果求和计算结果求和,就就得到整个区间上积分得到整个区间上积分I的近似值。的近似值。所用方法所用方法:第50页/共152页4.3.1 复合梯形
32、公式复合梯形公式每个子区间每个子区间xk,xk+1上的积分用上的积分用梯形公式梯形公式(1.1),得得将积分区间将积分区间a,b划分为划分为n等分等分,则则得到得到第51页/共152页称为称为复合梯形公式复合梯形公式,其余项可由,其余项可由(1.10)式得式得 记记由于存在由于存在f(x)C C2 2a,b,且,且 所以存在所以存在 (a,b)使使 于是于是复合梯形公式的余项复合梯形公式的余项为为 第52页/共152页当当n时,上式右端括号内的两个和式均收敛到函时,上式右端括号内的两个和式均收敛到函数的积分,所以复合梯形公式收敛数的积分,所以复合梯形公式收敛.此外,此外,Tn的求积的求积系数均
33、为正,由定理系数均为正,由定理2知复合梯形公式是稳定的知复合梯形公式是稳定的.可以看出误差是可以看出误差是h2阶,且由误差公式得到,当阶,且由误差公式得到,当f(x)C2a,b 时,则有时,则有即复合梯形公式是收敛的即复合梯形公式是收敛的.事实上只要事实上只要f(x)Ca,b,则可得到收敛性,因为只要把则可得到收敛性,因为只要把Tn改写为改写为第53页/共152页4.3.2 复合辛普森求积公式复合辛普森求积公式 将积分区间将积分区间a,b 划分为划分为n等分等分,在每个子区间,在每个子区间xk,xk+1上采用上采用辛普森公式辛普森公式(2.3),若记若记xk+1/2=xk+(1/2)h,则得则
34、得记记称为称为复合辛普森求积公式复合辛普森求积公式.其余项由其余项由(2.5)式得式得第54页/共152页于是当于是当 f(x)C 4a,b时时,其其求积余项求积余项为为由由(3.6)式看出,误差阶为式看出,误差阶为h4,收敛性是显然的,收敛性是显然的,事实上,只要事实上,只要 f(x)C a,b则可得到则可得到 收敛性,即收敛性,即此外,由于此外,由于Sn中求积系数均为正数,故知复合辛中求积系数均为正数,故知复合辛普森公式计算稳定普森公式计算稳定.第55页/共152页 例例3 对于函数对于函数f(x)=sinx/x,给出,给出n=8的函数表,的函数表,试用复合梯形公式和复合辛普森公式计算积分
35、试用复合梯形公式和复合辛普森公式计算积分xf(x)01/81/43/81/25/83/47/8110.99739780.98961580.97672670.95885100.93615560.90885160.87719250.8414709 解解 将积分区间将积分区间0,1划分为划分为8等等分,用复合梯形公式求得分,用复合梯形公式求得而将积分区间而将积分区间0,1划分为划分为24等分,等分,用复合辛普森公式求得用复合辛普森公式求得并估计误差并估计误差.第56页/共152页 比较上面两个计算结果比较上面两个计算结果T8与与S4,它们都需要提供,它们都需要提供9个点上的函数值,然而精度却差别很大
36、,同积分准个点上的函数值,然而精度却差别很大,同积分准确值确值I=0.9460831比较,应用复合梯形公式计算的结比较,应用复合梯形公式计算的结果果T8=0.9456909只有只有2位有效数字,而应用复合辛普位有效数字,而应用复合辛普森公式计算的结果森公式计算的结果S4=0.9460832却有却有6位有效数字位有效数字.为了利用余项公式估计误差,要求为了利用余项公式估计误差,要求f(x)=sinx/x的的高阶导数,由于高阶导数,由于所以有所以有第57页/共152页于是于是复合梯形公式误差复合梯形公式误差为为复合辛普森公式误差复合辛普森公式误差为为第58页/共152页 例例4 利用利用复合梯形公
37、式复合梯形公式计算计算 使其误差使其误差限为限为10-4,应将区间,应将区间0,1几等分几等分?解解 利用例利用例3的结果的结果取取n=17可满足要求可满足要求.由由复合梯形公式的余项得复合梯形公式的余项得第59页/共152页 例例5 利用利用复合辛普森公式复合辛普森公式计算计算 使其误使其误差限为差限为10-4,应将区间,应将区间0,1几等分几等分?由由复化合辛普森公式的余项得复化合辛普森公式的余项得因此只需将区间因此只需将区间0,1二等分,即取二等分,即取m=1(n=2).解解 利用例利用例3的结果的结果第60页/共152页 前面用复合梯形公式计算此题,满足相同的精度前面用复合梯形公式计算
38、此题,满足相同的精度需要将区间需要将区间0,1划分划分17等分,可见复合辛普森公式的等分,可见复合辛普森公式的精度的确比复合梯形公式精度高同样也可用精度的确比复合梯形公式精度高同样也可用|S4m-S2m|来控制计算的精度来控制计算的精度.这就是下面要介绍的这就是下面要介绍的龙贝格求积龙贝格求积公式公式.参见参见书书p108的例的例4讲解讲解.第61页/共152页4.4 龙贝格求积公式龙贝格求积公式4.4.1 梯形法的递推化梯形法的递推化 上节介绍的复合求积方法可提高求积精度,实际上节介绍的复合求积方法可提高求积精度,实际计算时若精度不够可将步长逐次分半计算时若精度不够可将步长逐次分半.设将区间
39、设将区间 a,b分为分为n等分,共有等分,共有n+1个分点,如果将求积区间再个分点,如果将求积区间再分一次,则分点增至分一次,则分点增至2n+1个,我们将二分前后两个个,我们将二分前后两个积分值联系起来加以考虑积分值联系起来加以考虑.并注意到每个子区间并注意到每个子区间xk,xk+1经过二分只增加了一个分点经过二分只增加了一个分点第62页/共152页 设设h=(b-a)/n,xk=a+kh (k=0,1,n),在在xk,xk+1上用梯形公式得上用梯形公式得在在xk,xk+1上用复合梯形公式得上用复合梯形公式得所以所以第63页/共152页从从0到到n-1对对k累加求和得累加求和得 这就是这就是递
40、推的复合梯形公式递推的复合梯形公式.从这一公式可以看出,将区间对分后,原复合从这一公式可以看出,将区间对分后,原复合梯形公式的值梯形公式的值Tn作为一个整体保留作为一个整体保留.只需计算出新只需计算出新分点的函数值,便可得出对分后的积分值,不需重分点的函数值,便可得出对分后的积分值,不需重复计算原节点的函数值,从而减少了计算量复计算原节点的函数值,从而减少了计算量.即即第64页/共152页 例例5 计算积分值计算积分值它在它在x=0的值定义为的值定义为f(0)=1,而,而f(1)=0.8414709,根据,根据梯形公式计算得梯形公式计算得 解解 我们先对整个区间我们先对整个区间0,1使用梯形公
41、式使用梯形公式.对于对于函数函数将区间二等分将区间二等分,再求出中点的函数值再求出中点的函数值f(1/2)=0.9588510,从而利用递推公式从而利用递推公式(4.1),有,有第65页/共152页进一步二分求积区间,并计算新分点上的函数值进一步二分求积区间,并计算新分点上的函数值f(1/4)=0.9896158,f(3/4)=0.9088516,再再(4.1)式,有式,有这样不断二分下去,计算结果见这样不断二分下去,计算结果见表表.k为二分次数,为二分次数,n=2k kTn123456789100.93979330.94451350.94569090.94598500.94605960.94
42、607690.94608150.94608270.94608300.9460831由表可见,用复合梯形公式由表可见,用复合梯形公式计算积分计算积分I要达到要达到7位有效数字的位有效数字的精度需要二分区间精度需要二分区间10次,即要有次,即要有分点分点1025个,计算量很大个,计算量很大.第66页/共152页4.4.2 外推技巧外推技巧 上面讨论说明由上面讨论说明由梯形公式梯形公式出发,将区间出发,将区间a,b逐逐次二分可提高求积公式的精度,上述加速过程还可继次二分可提高求积公式的精度,上述加速过程还可继续下去,其理论依据是续下去,其理论依据是梯形公式的余项展开梯形公式的余项展开,设,设若记若记
43、Tn=T(h),当区间,当区间a,b划分为划分为2n等分时,则有等分时,则有并且有并且有可以证明梯形公式余项可展开成可以证明梯形公式余项可展开成级数形式级数形式,即,即第67页/共152页 定理定理4 设设f(x)Ca,b,则有,则有其中其中I为积分值,系数为积分值,系数 l(l=1,2)与与h 无关无关.此定理可利用此定理可利用f(x)的泰勒展开推导得到,的泰勒展开推导得到,证略证略.定理定理4表明表明T(h)I是是O(h2)阶阶,若用若用h/2代替代替h,有有用用4乘乘(4.3)式减去式减去(4.2)式再除式再除3记为记为S(h),则得,则得第68页/共152页这里系数这里系数 l与与h无
44、关,这样构造的无关,这样构造的S(h)与积分值与积分值I近似近似的误差阶为的误差阶为O(h4).这比这比复合梯形公式复合梯形公式的误差阶的误差阶O(h2)提高了,容易看到提高了,容易看到S(h)=Sn,即将,即将a,b分为分为n等分得等分得到的到的复合辛普森公式复合辛普森公式.这种将计算这种将计算I的近似值的误差的近似值的误差阶由阶由O(h2)提高到提高到O(h4)的方法称为的方法称为外推算法外推算法,也称为,也称为理查森理查森(Richardson)外推算法外推算法.这是这是“数值分析数值分析”中中一个重要的技巧,只要真值与近似值的误差能表示一个重要的技巧,只要真值与近似值的误差能表示成成h
45、的幂级数,如的幂级数,如(4.2)式所示,都可使用外推算法,式所示,都可使用外推算法,提高精度提高精度.与上述做法类似,从与上述做法类似,从(4.4)式出发,当式出发,当n再增加一再增加一倍,即倍,即h减少一半时,有减少一半时,有第69页/共152页从从(4.6)式出发,利用外推技巧还可得到逼近阶为式出发,利用外推技巧还可得到逼近阶为O(h8)的算法公式的算法公式用用16乘乘(4.5)式减去式减去(4.4)式再除式再除15记为记为C(h),则得,则得它就是把区间它就是把区间a,b分为分为n个子区间的复合个子区间的复合柯特斯公式柯特斯公式.C(h)=Cn,它的精度为,它的精度为C(h)-I=O(
46、h6).它由辛普森法二它由辛普森法二分前后的两个积分近似值分前后的两个积分近似值Sn与与S2n=S(h/2)由由(4.6)式组式组合得到,即合得到,即如此继续下去就可得到龙贝格如此继续下去就可得到龙贝格(Romberg)算法算法.第70页/共152页4.4.3 龙贝格算法龙贝格算法将上述外推技巧得到的公式将上述外推技巧得到的公式(4.4),(4.6),(4.8)重新引重新引入记号入记号T0(h)=T(h),T1(h)=S(h),T2(h)=C(h),T3(h)=R(h)等等,从而可将上述公式写出统一形式从而可将上述公式写出统一形式用用m(h)作为作为I 的近似值的近似值,误差量级为误差量级为O
47、(h2(m+1).经过经过m(m=1,2,)次加速后,余项便取下列形式次加速后,余项便取下列形式 这种处理方法通常称为这种处理方法通常称为理查森理查森(Richardson)外推外推加速方法加速方法.第71页/共152页公式公式(4.11)也称为也称为龙贝格求积算法龙贝格求积算法.以以0(k)表示二分表示二分k次后求得的梯形值,以次后求得的梯形值,以m(k)表示序列表示序列0(k)的的m次加速值次加速值,则依递推公式,则依递推公式(4.9)可到可到第72页/共152页 龙贝格求积算法龙贝格求积算法的的计算过程计算过程如下:如下:(1)取取k=0,h=b-a,求,求令令1k(k记区间记区间 a,
48、b的二分次数的二分次数).(2)求值求值 ,按梯形递推公式计算,按梯形递推公式计算0(k).(3)求计算值,按加速公式逐个求出求计算值,按加速公式逐个求出数表数表的第的第k行其余各元素行其余各元素j(k-j)(j=1,2,k).(4)若若|k(0)-k-1(0)|0,计算积分,计算积分的近似值的近似值.先取步长先取步长h=b-a,应用辛普森公式有,应用辛普森公式有其中其中若把区间若把区间a,b对分,步长对分,步长h2=h/2=(b-a)/2,在每个小区,在每个小区间上用辛普森公式,则得间上用辛普森公式,则得第80页/共152页其中其中实际上实际上(5.2)式即为式即为第81页/共152页与与(
49、5.1)式比较,若式比较,若f(4)(x)在区间在区间(a,b)上变化不大,可假上变化不大,可假定定f(4)()f(4)(),从而可得,从而可得与与(5.2)式比较,则得式比较,则得这里这里S1=S(a,b),S2=S2(a,b).如果有如果有则可期望得到则可期望得到第82页/共152页此时可取此时可取S2(a,b)作为作为I(f)的近似,则可达到给定的误差的近似,则可达到给定的误差精度精度 ,若不等式,若不等式(5.3)不成立,则应分别对子区间不成立,则应分别对子区间a,(a+b)/2 及及(a+b)/2,b再用辛普森公式,此时步长再用辛普森公式,此时步长h3=(1/2)h2,得到,得到S3
50、(a,(a+b)/2)及及S3(a+b)/2,b).只要分只要分别考察下面两个不等式别考察下面两个不等式是否成立是否成立.对满足要求的区间不再细分,对不满足要对满足要求的区间不再细分,对不满足要求的还要继续上述过程,直到满足要求为止,最后还求的还要继续上述过程,直到满足要求为止,最后还要应用龙贝格法则求出相应区间的积分近似值要应用龙贝格法则求出相应区间的积分近似值.为了为了更直观地说明自适应积分法的计算过程及方法为何能更直观地说明自适应积分法的计算过程及方法为何能节省计算量,节省计算量,看看p115的的例例7讲解讲解.第83页/共152页4.6 高斯求积公式高斯求积公式 由前面的讨论已经知道,