5数值积分与数值微分课件

第5章数值积分与数值微分方法分步积分法等，但实际问题中，经常出现不能用通常这些积分方法计算的定积分问题。怎样把这些通常方法失效的定积分在一定精度下快速计算出来，特别是本章涉及的方法有Newton-Cotes求积公式、5.1引例人造地球卫星轨道可视为平面上的椭圆。我国的第一颗人造地球卫星近地点距离地球表面439km，远地点距地球表2384km，地球半径为6371km，求该卫星的轨道长度。本问题可用椭圆参数方程该轨道的长度L就是如下参数方程弧长积分0但这个积分是椭圆积分，不能用解析方法计算。aaa入0iii=1要离散化，做1)去掉极限号lim2)将i取为具体的xi值i3)为减少离散化带来的误差，将xi用待定系数Ai代替于是就得到12n12n取f(x)C[a,b],都有iii=1则称式(5.1)为一个数值求积公式。 (5.1)iR(f)=jbf(x)dxnAf(x)(5.2)aiiai=1为求积余项或求积公式(5.1)的截断误差。个数n、求积节点{x}和求积系数{A}，这三个量有一ii定义5.2若求积公式ii对所有不超过m对所有不超过m次的多项式p(x)有求积余m项R(p)=0，而对某一个m+1次多项式mp(x)有R(p)0，则称该求积公式的代数mm+1精度为m。iii=1a一般，一个求积公式的代数精度越大，则该求积公式越好。确定代数精度的方法依次取f(x)=xk(k=0,1,)代入公式iii=1aRxk=0不成立的k值为m，例5.1确定求积公式jfxdxff3)133解取f(x)=xk代入求积公式有133=[()k](1+(1)k)k+13易验证故本题求积公式代数精度为3。例5.2确定下面求积公式2h精度，并指出相应的代数精度。解本题要先求出具体的求积公式，然后代数精度。就可以解决求出具体求积公式的问题。一2h939为确定其代数精度，再取f(x)=x3代入求出的公式继35.3插值型求积公式借助多项式插值函数来构造的求积公一般选用不同的插值公式就可以得到基本思想利用被积函数f(x)的插值函数代替f(x)做定积分的近似计算来构造求积公式。1.构造原理2考虑f(x)在n个节点x,2n上的n-1次Lagrange插值多项式L(x)与f(x)的余项，有f(x)=f(x)l(x)+f(n)()(x)x[a,b]i=1iin1n!nai=1iain-1an!n (5.4)记则有AAbl(x)dxiin1a (5.3)ai1iian!n若舍去bf(n)()(x)dx，得求积公式an!nbf(x)dxnAf(x)iii1求积系数Abl(x)dx的求积公式就是插值型aiin1a求积公式。插值型求积公式的求积余项Rfjbfn)o(x)dxan!n当f(x)为次数小于n次的多项式时，有因此插值型求积公式的代数精度至少为n-1。若取f(x)=1，代入式(5.4)，可得插值型求积公式的求积系数之和为k下面具体介绍常用的几个插值型求积公式。2.Newton-Cotes求积公式1)n点的Newton-Cotes公式的构造i将求积节点x取为[a,b]上的等距节点iin-1求积公式的求积系数为--iiix xi通常称 (通常称 (5.6)为n点的Newton-Cotes公式。jjbf(x)dx~(b-a)xnC(n)f(x)iii=1a由于求积节点x是等距的，因此也称式i (5.6)为等距节点求积公式。Newton-Cote公s式A)2点的Newton-Cotes公式a2我们熟悉的梯形公式。 (5.7)B)3点的Newton-Cotes公式为称它为Simpson公式或抛物线公式。表5.1部分Cotes系数nnC(n)i223456789例5.3试分别用梯形公式和Simpson公式计算解用梯形公式计算，有用Simpson公式计算，有2)n点Newton-Cotes公式的代数精度定理5.1当求积节点个数n为奇数时，对应的Newton-Cotes求积公式的代数精度至少为n。证明由于是插值型求积公式，故有an!k0kk=1记(s)=(1)m=k+1nl(s+m)，易知(-s)=-(s)，-l3)梯形公式与Simpson公式的余项引理5.1(积分中值定理)有a引理5.1，有12!a12 (5.10)梯形公式余项12抛物线公式的余项 (5.9)2880fx)=ffx)=f(x)-ck=1,2,nkkk4)Newton-Cotes公式的数值稳定性设计算函数f(x)时产生舍入误差为ckk实际在计算机中参加计算的是f(x)的近似值k((故Newton-Cotes公式在计算机中产生的误差为niiiii=1i=1iinii=1由Cotes表5.1,当时,C(n)>0，从而有inini说明此时计算舍入误i=1差可以控制，从而iiii增大而增大，从而导致舍入误差i=增1加。故n>8时，Newton-Cotes公式是数值不稳Newton-Cotes公式来做定积分计算。3.Gauss求积公式1)Gauss求积公式的构造与概念n点的插值型求积公式的代数精度至少n数精度呢？为回答这个问题，观察一下插值求积公式的构造方法，发现其至少具有n-1次代数精度的结论是在限定求积节点x的前提下得i出的，若让求积节点x也可以自由取值，则i就给提高代数精度创造了条件。 (5.11)为使问题讨论更具一般性，这里考虑带权的定积分求积公式jjbp(x)f(x)dxnAf(x)kka式中p(x)是已知的非负函数，为区间[a,b]显然在式(5.11)中，若p(x)=1，就是定理5.2求积公式(5.11)的代数精度最大为证明设x,x,,x是求积公式(5.11)的任意一12n组求积节点，用此节点构造插值型求积公式(5.11)，取2n次多式f(x)=12(x)代公式(5n.11),考虑它n的求积余项，有nknknkn因为2(x)是2n次多式，由代数精度定义得nnakka为证求积公式jbp(x)f(x)dx则xnAf(x)的代数精度akka用o(x)去除f(x)，由多项式除法有nnkkaaaa(5.12)因为上面的求积公式(5.11)是具有n个节点的插值型akkkka要让式(5.11)成为等式，由式(5.12)应有jbpxqxxdx0(5.13)na式(5.13)要求对固定的n次多项式(x)和任意次数不超过n1的多项式q(x)都成立，这可以用q(x)的这种任意性，选择求积节点x,x,,x。12n由正交多项式理论可知：使式(5.13)成立的点度的求积公式了关于插值型求积公式的结论定理5.3n点插值型求积公式的代数精定义5.3若n点的求积公式具有2n1次代数精度，则称该求积公式为Gauss型求积公式，对应的求积节点x和求积系数Akk分别称为Gauss点和Gauss系数。例5.4确定参数x,x,A,A，使求积公式1212成为Gauss求积公式。解注意到被积函数中有一因式与求积公式右端无明显的关系，可将其视为权函数。为确定四个参数，依次取代入公式并将近似号取为等号，得联立方程组1919解出22于本题求积节点个数为n=2，其最高2)Gauss型求积公式的求积余项fxCnabfxGaussxxx12n上的2n-1次Hermite插值多项式H(x),由2n-1Hermite插值余项公式，有两边积分，有由Gauss型求积公式的代数精度为2n-1及积分中值定理有，Gauss求积余项为R(f)=jbp(x)f(x)dx-xnAf(x)3)Gauss型求积公式的数值稳定性对任意i,取f(x)=l2(x),l(x)是关于in-1in-1nge函数，有l1()=6n,由Gauss公式的代数精ain-1kin-1kiaikann类似Newton-Cotes稳定性处理方法，有在计算机计算时的舍入误差n=kkAkckkkka公式，但最基本和常用的Gauss公式有44)常用的Gauss型求积公式 (1)Gauss-Legendre求积公式Gauss点为n次Legendre正交多项式的零点，Gauss-Legendre求积公式为jj1f(x)dx如xnAf(x)-1kkGauss-Legendre求积余项为表5.2nnxknxkAkAkGauss-Legendre求积公式可以计算任何有限积分区间的定积分，计算之前先作变量代换22然后再对用Gauss-Legendre求积公式。 (2)Gauss-Chebyshev求积公式ussnj1f(x)dx~xnAf(x)-11-x2k=1kk式中Gauss点xk与系数Ak-"=k-"="xcos,A,k1,,nk2nknGauss-Chebyshev求积余项GaussGauss-Laguerre求积余项为 (3)Gauss-Laguerre求积公式为n次Laguerre正交多项式ndxn的零点，Gauss-Laguerre求积公式为0kkRffn0,Gauss-Laguerre求积公式的Gauss点和系数也有 (4)Gauss-Hermite求积公式ex次Hermite正交多项式,Gauss点为n=_ee_H(x)(1)nx2d=_ee_nnGaussHermite求积公式为jj+we_x2f(x)dx如xnAf(x)kk_wGauss-Hermite求积余项为Gauss-Hermite求积公式的Gauss点及系数有表例5.5用两点Gauss公式求定积分0解本题为有限区间的定积分，可用两点Gauss-Legendre求积公式计算。做积分换元，将其化为[-1,1]上的定积分，即令01常用复化求积公式来求积区间[a,b]上的定积分，以获得满足给定计算精度要的定积分值。常用的复化求积公式有复化梯形公式和复化Simpson公式。将求积区间[a,b]分成若干个小区间，然后在每个小区间上采用数值稳定的Newton-Cotes公式求小区间上的定积分，最后把所有小区间上的计算结果相1.复化梯形公式1)复化梯形公式的构造原理njbf(x)dx=jxk+1f(x)dx~h[f(x)+f(x)]ax2kk+1k=0kk=0得求积公式---复化梯形公式：2)复化梯形公式的余项jbfxdxThf,(n)n=[x,x]an12kkkk+1k=0nan12()22()2222说明利用复化求积公式能得到计算误差小于c的定积2分近似值。(5.16)(5.16)2.复化Simpson公式1)复化Simpson公式的构造原理取等距节点x=a+ih,h=b-a,i=0,1,2,,n将积分区in间[a,b]n等分,在每个小区间[x,x],k=0,1,,n-1上用Simpson公式做近似计算，再累加起来就有式中x=x+h，得复化Simpson公式kk22)复化Simpson公式的余项有复化Simpson公式的求积余项n从复化simpson公式的求积余项可以看出复化0例5.6分别用复化梯形公式和复化Simpson公式计算nN复化梯形公式复化Simpson公式2-17.3892595.3222-13.3360231.27592840-0.47823-12.3821620.3122.148004-12.08974226-12.07519427-12.071558-12.070649题积分的准确值为，可见复化梯形例5.7考虑用复化Simpson公式计算nxdx0x要使误差小于，那么求积区间[0,1]应分成多少解复化Simpson公式的求积余项为nnxx0x出45.5Romberg求积方法Romberg求积方法是对复化梯形公式用加速技术得到的一种求积方法，它也称为逐次分半加速收敛思想将Richardson外推算法应用于复化梯形公式中，F(qh)F(qh)qp1F(h)F(h)=00,0<q<11qp1k1.构造原理定理5.4设函数F(h)逼近数F*的余项为0F*F*F(h)01=+1=+++<<< p 22 p 33pp,02p,p,都是与h无关的常数，且k1时，0,kkkF(h)也逼近F*，且有1FF*F(h)=(1)hp2+(1)hp3+123F*FF*F(qh)=a(qh)p1+a(qh)p2+012(5.19)k证明用qh代替余项式(5.18)的变元h，有k用qp1乘式(用qp1乘式(5.18)并与式(5.19)相加，整理后，有有qhqpFhqhqpFhF*00(h),(h),由微积分收敛阶概念可知F(h)的收敛阶p比2F(h)的收敛阶p大，故F(h)比F(h)逼近F*的速度快。0110通常称如上的方法为Richardson外推法。显然这种外推可以不断做下去以获得逼近F*更快FF(qh)qpmF(h)F(h)=mm,m=0,1,2,mqpmF*F(h)=a(m)hpm+1+a(m)hpm+2+这样，每用一次Richardson外推，就使逼近阶提高称为复化梯形公式的T形值，可以证明：0123此式说明(h)逼近I的阶为O(h2)。利用Richardson外推法对T(h)做加速。因为p=2,有(h)逼近I的阶变为O(h4)。若取q=12，则有T(h/2)-(12)2T0(h)4T(h2((5.22)同理对T(h)再做一次Richardson加速，有1T(h)=11242124212一般，有经Richardson加速的求定积分的序列为T(h)=m1m1,m=1,2,m4m10n(5.21)n1n1SS=2nn4TTC=42S2nSn(5.23)n421R=43C2nCn(5.24)n431这样可以用复化梯形公式计算出序列。再逐次用公式(5.22)，(5.23)，(5.24)对其进行加工得到若继续外推，舍入误差增大将使新序列体现不出更明显的加速效果，因此一般只外推到Romberg序用Romberg序列求定积分的算法称为Romberg3．Romberg求积方法的计算过程k02343因为T与T有相重的点，计算T时还可以利用下面n2n2n关系式来减小计算量。例5.8用Romberg算法计算。解由Romberg求积方法，得计算结果如下表kk000.9200.9460.9460.946735514590830083110.9390.9460.9460.946793308690831083020.94451350.94608340.946083030.94569090.946083040.98459850由RR2125.6数值微分根据函数在若干个点处的函数值去求该函数的导数近似值称为数值微分，所求导数的近似值常称为设y=f(x)(k=0,1,...,n)是已知的n+1个点kknx<x<nn值导数用插值函数代替被插值函数可以用来计算被插函数的近似值；用插值函数代替定积分的被积函数可以用来计算定积分的近似值；用插值函数的导数能否用来求被插函数导数定义5.4设P(x)是f(x)的n次插值多项式，称用n插值函数的导数来计算被插函数导数的公式n为f(x)的插值型求导公式。n在非节点处的余项由于涉及对未知中值函数x述和控制，但在节点处有R(x)=f(x)P(x)=f(n+1)().(x)k=0,1,,n1kknk(n+1)!n+1k此余项易得出误差的界限，因此插值型求导公式多用在节点处的数值导数。f(f(x)=[f(x)f(x)]f()0h102f(x)=[f(x)f(x)]+f()1h1021)两点公式给定点x,f(x)(k=0,1)，有线性插值多项式kkP(x)=L(x)=xx1f(x)+xx0f(x)11xx0xx1P(x)=f(x1)f(x0),k=0,11kh在节点处的余项1k2!2kk012!01R(x)=f()(x)=(2xxx)f(),(x,x1k2!2kk012!01f,f,(x)=[一3f(x)+4f(x)一f(x)]+f,(飞)02h01231h2f,(x)=[一f(x)+f(x)]一f,(飞)12h0261h2f,(x)=[f(x)一4f(x)+3f(x)]+f,(飞)22h01231h22)三点公式只给出等距节点的公式形式，对非等距节点有kkk022f,(x)必P,(x)及带余项的三点数值微分公式k2k前差公式02中间的节点使用中心公式，左边界节点用前插公式，在求数值微分公式余项表达式时，有时从插值函数的余项表达式不易得出数值微分公式的余项，但借例5.9推导计算二阶导数的中心差分公式f,(x)=1(f(x)一2f(x)+f(x))一h2f(4)(飞)1h20121202k0解用数据集{x,f(x)}做二次插值多项式,有kk0210122021022020122两边对t求二阶导数20012故有220h2012ff(x)P(x)=1(f(x)一2f(x)+f(x))121h2012导，故不能用求导的方式得出余项表达式，下面用ffxfxhfxfxhfxhfx1)h3+f(4)(1)h401112!3!4!ff(x)=f(x+h)=f(x)+f(x)h+f(x1)h2+f(x1)h3+f(4)(2)h421112!3!4!101212有ff(x)2f(x)+f(x)=f,(x)h2+h4(f(4)()+f(4)())01214!1注意到f(4)()+f(4)()=2f(4)