第4章 数值积分与数值微分.ppt

1、第4章数值积分与数值微分,4.1数值积分概论4.2牛顿柯特斯公式4.3复合求积公式4.4龙贝格求积公式4.5自适应求积方法4.6高斯求积公式4.7多重积分4.8数值微分,本章基本内容,进行计算，但在工程计算和科学研究中，经常会遇到被积函数f(x)的下列一些情况：,4.1数值积分概论,实际问题当中常常要计算积分，有些数值方法，如微分方程和积分方程的求解，也都和积分计算相联系.,4.1.1数值求积的基本思想,（4）f(x)本身没有解析表达式，其函数关系由表格或图形给出，例如为实验或测量数据.,（2）f(x)的原函数不能用初等函数形式表示，例如,（3）f(x)的原函数虽然可用初等函数形式表示，但其原

2、函数表示形式相当复杂，例如,（1）f(x)复杂，求原函数困难，例如,以上的4种情况都不能用牛顿莱布尼兹公式方便地计算该函数的定积分，满足不了实际需要，因此，有必要研究定积分的数值计算问题；另外，对一些函数的求导问题，其求导、微分也相当复杂，也有必要研究求导、微分的数值计算问题。本章主要介绍数值求积分和数值求微分的方法。,由积分中值定理,对连续函数f(x),在区间a,b内至少存在一点，使,只要对平均高度f()提供一种近似算法,便可相应地获得一种数值求积方法.即所谓矩形公式.,例如,用区间a,b两端点的函数值f(a)与f(b)的算术平均值作为f()的近似值,可导出求积公式,这便是人们所熟知的梯形公

3、式.,如果改用区间a,b的中点c=(a+b)/2处的函数值f(c)近似代替f(),则又可导出所谓(中)矩形公式,一般地,在区间a,b上适当选取点xk(k=0,1,n),然后用f(xk)的加权平均值作为f()的近似值,可得到更为一般的求积公式,其中：点xk叫求积节点,系数Ak叫求积系数.Ak仅与节点xk的选取有关,而与被积函数f(x)无关.,求积公式的截断误差为,R(f)又称为求积余项.,这类数值积分方法通常称为机械求积，其特点是将积分求值问题归结为函数值的计算，这就避开了牛-莱公式寻求原函数的困难.,4.1.2代数精度的概念,定义1如果求积公式,(1)对所有次数不超过m的多项式都精确成立；(2

4、)至少对一个m+1次多项式不精确成立，则称该公式具有m次代数精度(或代数精确度).,数值求积方法的近似方法，为要保证精度，我们自然希望求积公式能对“尽可能多”的函数准确地成立，这就提出了所谓代数精度的概念.,一般来说，代数精度越高，求积公式越好。,结论一个求积公式具有m次代数精度的充要条件是该求积公式:(1)对xk(k=0,1,m)精确成立；(2)对xm+1不精确成立.,故一般地，要验证一个求积公式具有m次代数精度，只要令对于f(x)=1,x,xm求积公式精确成立等式就行.,即对于求积公式,给定n+1个互异的求积节点x0,x1,xn-1,xn,令求积公式对f(x)=1,x,xn精确成立,即得,

5、求解该方程组即可确定求积系数Ak,所得到的求积公式至少具有n次代数精度.,解当f(x)=1时,此时公式精确成立。,当f(x)=x时，,公式也精确成立.,当f(x)=x2时，,公式对x2不精确成立.,故由定理1知,梯形公式的代数精度为1次.,例2确定求积公式中的待定系数，使其代数精度尽量高，并指明求积公式所具有的代数精度.,解令f(x)=1,x,x2代入公式两端并令其相等，得,解得,得求积公式为,令f(x)=x3，得,令f(x)=x4，得,故求积公式具有3次代数精度.,如果我们事先选定求积节点xk，譬如，以区间a,b的等距分点作为节点，这时取m=n求解方程组即可确定求积系数Ak，而使求积公式至少

6、具有n次代数精度.本章第2节介绍这样一类求积公式，梯形公式是其中的一个特例.,如为了构造出上面的求积公式，原则上是一个确定参数xk和Ak的代数问题.方程组(1.4)实际上是一2n+2个参数的非线性方程组，此方程组当n1时求解非常困难，但当n=0及n=1的情形还是可以通过求解方程组得到相应的求积公式.下面对n=0讨论求积公式的建立及代数精确度.,此时求积公式为,其中，x0及A0为待定参数.根据代数精确度定义可令f(x)=1,x，由方程组知.,得,得到的求积公式就是(1.2)式的中矩形公式.再令f(x)=x2，代入(1.4)式的第三式,说明(1.2)式对f(x)=x2不精确成立，故它的代数精确度为

7、1.,方程组(1.4)是根据(1.3)式的求积公式得到的，按照代数精确度的定义，如果求积公式中除了f(xi)还有(x)在某些节点上的值，也同样可得到相应的求积公式.,例1给定形如下面的求积公式，试确定系数A0,A1,B0,，使公式具有尽可能高的代数精确度.,解根据题意可令f(x)=1,x,x2分别代入求积公式使它精确成立：,解得,当f(x)=x3时，上式右端为1/3，而左端是,于是有求积公式,故积分公式对f(x)=x3不精确成立，其代数精确度为2.,4.1.3插值型的求积公式,设给定一组节点,且已知f(x)在这些节点上的函数值f(xk),则可求得f(x)的拉格朗日插值多项式(因为Ln(x)的原

8、函数易求),其中lk(x)为插值基函数,取,由上式确定系数的公式称为插值型求积公式.,即,则f(x)Ln(x),插值型求积公式积分法几何表示,由插值余项定理,其求积余项为,其中=(x),如果求积公式(1.5)是插值型的，按照插值余项式子，对于次数不超过n的多项式f(x)，其余项R(f)等于零，因而这时求积公式至少具有n次代数精度.,反之，如果求积公式至少具有n次代数精度，则它必定是插值型的.事实上，这时求积公式对于插值基函数lk(x)应准确成立，即有,注意到lk(xj)=kj，上式右端实际上即等于Ak，因而下面式子成立.,定理1具有n+1个节点的数值求积公式(1.5),是插值型求积公式的充要条

9、件为:该公式至少具有n次代数精度.,综上所述，我们有结论为,这时令f(x)=1代入又有结论为,结论对插值型求积公式的系数必有,4.1.4求积公式的余项,若求积公式(1.3)的代数精确度为m，则由求积公式余项的表达式(1.7)可以证明余项形如,其中K为不依赖于f(x)的待定参数.这个结果表明当f(x)是次数小于等于m的多项式时，由于f(m+1)(x)=0，故此时Rf=0，即求积公式(1.3)精确成立.而当f(x)=xm+1时，f(m+1)(x)=(m+1)!，(1.8)式左端Rf0，故可求得,代入余项公式(1.8)式可以得到更细致的余项表达式.,例如梯形公式(1.1)的代数精确度为1，可以证明它

10、的余项表达式为,其中,于是得到梯形公式(1.1)的余项为,对中矩形公式(1.2)，其代数精确度为1，可以证明它的余项表达式为,其中,于是得到中矩形公式(1.2)的余项为,例2求例1中求积公式,的余项.,解由于此求积公式的代数精确度为2，故余项表达式为R=K()，令f(x)=x3，得()=3!，于是有,故得,其中h=max(xi-xi-1)，则称求积公式(1.3)是收敛的.,4.1.5求积公式的收敛性与稳定性,定义2在求积公式(1.3)中，若,在求积公式(1.3)中，由于计算f(xk)可能产生误差k，实际得到，即.记,如果对任给小正数0，只要误差|k|充分小就有,它表明求积公式(1.3)计算是稳

11、定的，由此给出,证明对任给0,若取=/(b-a),对所有k都有,故求积公式(1.3)是稳定的.证毕.,定理2若求积公式(1.3)中所有系数Ak0，则此求积公式是稳定的.,则有,定理2表明只要求积系数Ak0，就能保证计算的稳定性.,4.2牛顿柯特斯公式,为便于上机计算，通常在内插求积公式中我们通常取等距节点，即将积分区间a,b划分n等分，即令步长h=(b-a)/n，且记x0=a,xn=b，则节点记为xk=x0+kh(k=0,1,n)，然后作变换:t=(x-x0)/h,代入求积系数公式，将会简化计算.,4.2.1柯特斯系数与辛普森公式,设将积分区间a,b划分成n等分,步长h=,求积节点取为xk=a

12、+kh(k=0,1,n),由此构造插值型求积公式,则其求积系数为,引入变换x=a+th,则有,(k=0,1,n),(k=0,1,n),其中,记,于是得求积公式,称为n阶牛顿-柯特斯(Newton-Cotes)公式.,显然,柯特斯系数与被积函数f(x)和积分区间a,b无关,且为容易计算的多项式积分.,称为柯特斯系数.,常用的柯特斯系数表,当n=1时，柯特斯系数为,这时的牛顿-柯特斯公式为一阶求积公式，就是我们所熟悉的梯形公式，即,当n=2时，柯特斯系数为,相应的牛顿-柯特斯公式为二阶求积公式，就是辛普森(simpson)公式(又称为抛物形求积公式)，即,式中,(k=0,1,2,3,4),h=(b

13、-a)/4.,n=4时的牛顿-柯特斯公式就特别称为柯特斯公式.其形式是,在柯特斯系数表中(见书p104)看到n7时，柯特斯系数出现负值，于是有,特别地，假定,则有,这表明在b-a1时，初始误差将会引起计算结果误差增大，即计算不稳定，故n7的牛顿-柯特斯公式是不用的.,4.2.2偶阶求积公式的代数精度,作为插值型求积公式，n阶牛顿-柯特斯公式至少具有n次代数精度(推论1).实际的代数精度能否进一步提高呢？,先看辛普森公式，它是二阶牛顿-柯特斯公式，因此至少具有二次代数精度.进一步用f(x)=x3进行检验，按辛普森公式计算得,另一方面，直接求积得,这时有S=I，即辛普森公式对不超过三次的多项式均能

14、精确成立，又容易验证它对f(x)=x4通常是不精确的(如取a=0,b=1进行验证有，S=5/24I=1/5)，因此，辛普森公式实际上具有三次代数精度.,一般地，我们可以证明下述论断：,定理3n阶牛顿-柯特斯公式的代数精度至少为,证明由定理1已知，无论n为奇数或偶数，插值型求积公式都至少具有n次代数精度.因此我们证明n为偶数的情形，即对n+1次多项式余项为零.,令n=2k，设,为任一n+1次多项式，其最高次系数为an+1，则它的n+1阶导数为,由余项公式,有,这里变换为x=a+th，注意xj=a+jh.,下面我们证明,作变换u=t-k，则,容易验证(u)为奇函数，即(-u)=-(u)，而奇函数在

15、对称区间上的积分为零，所以,定理3说明，当n为偶数时，牛顿-柯特斯公式对不超过n+1次的多项式均能精确成立，因此，其代数精度可达到n+1.正是基于这种考虑，当n=2k与n=2k+1时具有相同的代数精度，因而在实用中常采用n为偶数的牛顿-柯特斯公式，如抛物形公式(n=2)等.,对牛顿-柯特斯求积公式通常只用n=1,2,4时的三个公式，n=1时即为梯形公式(1.1),其余项为(1.10)式.n=2时即为辛普森公式(2.3),其代数精度为3，可以证明余项可表示为,其中K由(1.9)式及(2.3)式可得,4.2.3辛普森公式的余项,从而可得辛普森公式(2.3)的余项为,也可直接积分计算得到,对n=4的

16、柯特斯公式(2.4)，其代数精度为5。故类似于求(2.3)式的余项可得到柯特斯公式的余项为,解：由梯形公式得,由辛普森公式得,例题分别用梯形公式、辛普森公式和柯特斯公式计算积分,由柯特斯公式得,积分的精确值,4.3复合求积公式,从求积公式的余项的讨论中我们看到，被积函数所用的插值多项式次数越高，对函数光滑性的要求也越高.另一方面，插值节点的增多(n的增大)，在使用牛顿-柯特斯公式时将导致求积系数出现负数(当n8时,牛顿-柯特斯求积系数会出现负数)，即牛顿-柯特斯公式是不稳定的，不可能通过提高阶的方法来提高求积精度.,为了提高精度，通常在实际应用中往往采用将积分区间划分成若干个小区间，在各小区间

17、上采用低次的求积公式(梯形公式或抛物形公式)，然后再利用积分的可加性，把各区间上的积分加起来，便得到新的求积公式，这就是复合求积公式的基本思想.为叙述方便，我们仅讨论各小区间均采用同一低次的求积公式的复合求积公式对各小区间也可分别采用不同的求积公式，也可推出新的求积公式，读者可按实际问题的具体情况讨论.,将积分区间a,bn等分,步长,xk=a+kh(k=0,1,n),则由定积分性质知,分点为,每个子区间上的积分,用低阶求积公式,然后把所有区间的计算结果求和,就得到整个区间上积分I的近似值。,所用方法:,4.3.1复合梯形公式,每个子区间xk,xk+1上的积分用梯形公式(1.1),得,将积分区间

18、a,b划分为n等分,则,得到,称为复合梯形公式，其余项可由(1.10)式得,记,由于存在f(x)C2a,b，且,所以存在(a,b)使,于是复合梯形公式的余项为,当n时，上式右端括号内的两个和式均收敛到函数的积分，所以复合梯形公式收敛.此外，Tn的求积系数均为正，由定理2知复合梯形公式是稳定的.,可以看出误差是h2阶，且由误差公式得到，当f(x)C2a,b时，则有,即复合梯形公式是收敛的.事实上只要f(x)Ca,b,则可得到收敛性，因为只要把Tn改写为,4.3.2复合辛普森求积公式,将积分区间a,b划分为n等分，在每个子区间xk,xk+1上采用辛普森公式(2.3),若记xk+1/2=xk+(1/

19、2)h,则得,记,称为复合辛普森求积公式.其余项由(2.5)式得,于是当f(x)C4a,b时,其求积余项为,由(3.6)式看出，误差阶为h4，收敛性是显然的，事实上，只要f(x)Ca,b则可得到收敛性，即,此外，由于Sn中求积系数均为正数，故知复合辛普森公式计算稳定.,例3对于函数f(x)=sinx/x，给出n=8的函数表，试用复合梯形公式和复合辛普森公式计算积分,解将积分区间0,1划分为8等分，用复合梯形公式求得,而将积分区间0,1划分为24等分，用复合辛普森公式求得,并估计误差.,比较上面两个计算结果T8与S4，它们都需要提供9个点上的函数值，然而精度却差别很大，同积分准确值I=0.946

20、0831比较，应用复合梯形公式计算的结果T8=0.9456909只有2位有效数字，而应用复合辛普森公式计算的结果S4=0.9460832却有6位有效数字.,为了利用余项公式估计误差，要求f(x)=sinx/x的高阶导数，由于,所以有,于是,复合梯形公式误差为,复合辛普森公式误差为,例4利用复合梯形公式计算使其误差限为10-4，应将区间0,1几等分?,解利用例3中的结果,取n=17可满足要求.,由复合梯形公式的余项得,对例4利用复合辛普森公式计算使其误差限为10-4，应将区间0,1几等分?,由复合辛普森公式的余项得,因此只需将区间0,1二等分，即取m=1(n=2).,解利用例3中的结果,前面用复

21、合梯形公式计算此题，满足相同的精度需要将区间0,1划分17等分，可见复合辛普森公式的精度的确比复合梯形公式精度高同样也可用|S4m-S2m|来控制计算的精度.这就是下面要介绍的龙贝格求积公式.,参见书p108的例4讲解.,4.4龙贝格求积公式,4.4.1梯形法的递推化,上节介绍的复合求积方法可提高求积精度，实际计算时若精度不够可将步长逐次分半.设将区间a,b分为n等分，共有n+1个分点，如果将求积区间再分一次，则分点增至2n+1个，我们将二分前后两个积分值联系起来加以考虑.并注意到每个子区间xk,xk+1经过二分只增加了一个分点,设h=(b-a)/n,xk=a+kh(k=0,1,n),在xk,

22、xk+1,上用梯形公式得,在xk,xk+1上用复合梯形公式得,所以,从0到n-1对k累加求和得,这就是递推的复合梯形公式.,从这一公式可以看出，将区间对分后，原复合梯形公式的值Tn作为一个整体保留.只需计算出新分点的函数值，便可得出对分后的积分值，不需重复计算原节点的函数值，从而减少了计算量.,即,例5计算积分值,它在x=0的值定义为f(0)=1，而f(1)=0.8414709，根据梯形公式计算得,解我们先对整个区间0,1使用梯形公式.对于函数,将区间二等分,再求出中点的函数值f(1/2)=0.9588510,从而利用递推公式(4.1)，有,进一步二分求积区间，并计算新分点上的函数值f(1/4

23、)=0.9896158,f(3/4)=0.9088516,再(4.1)式，有,这样不断二分下去，计算结果见表.k为二分次数，n=2k,由表可见，用复合梯形公式计算积分I要达到7位有效数字的精度需要二分区间10次，即要有分点1025个，计算量很大.,4.4.2外推技巧,上面讨论说明由梯形公式出发，将区间a,b逐次二分可提高求积公式的精度，上述加速过程还可继续下去，其理论依据是梯形公式的余项展开，设,若记Tn=T(h)，当区间a,b划分为2n等分时，则有,并且有,可以证明梯形公式余项可展开成级数形式，即,定理4设f(x)Ca,b，则有,其中I为积分值，系数l(l=1,2)与h无关.,此定理可利用f

24、(x)的泰勒展开推导得到，证略.,定理4表明T(h)I是O(h2)阶,若用h/2代替h,有,用4乘(4.3)式减去(4.2)式再除3记为S(h)，则得,这里系数l与h无关，这样构造的S(h)与积分值I近似的误差阶为O(h4).这比复合梯形公式的误差阶O(h2)提高了，容易看到S(h)=Sn，即将a,b分为n等分得到的复合辛普森公式.这种将计算I的近似值的误差阶由O(h2)提高到O(h4)的方法称为外推算法，也称为理查森(Richardson)外推算法.这是“数值分析”中一个重要的技巧，只要真值与近似值的误差能表示成h的幂级数，如(4.2)式所示，都可使用外推算法，提高精度.,与上述做法类似，从

25、(4.4)式出发，当n再增加一倍，即h减少一半时，有,从(4.6)式出发，利用外推技巧还可得到逼近阶为O(h8)的算法公式,用16乘(4.5)式减去(4.4)式再除15记为C(h)，则得,它就是把区间a,b分为n个子区间的复合柯特斯公式.C(h)=Cn，它的精度为C(h)-I=O(h6).它由辛普森法二分前后的两个积分近似值Sn与S2n=S(h/2)由(4.6)式组合得到，即,如此继续下去就可得到龙贝格(Romberg)算法.,4.4.3龙贝格算法,将上述外推技巧得到的公式(4.4),(4.6),(4.8)重新引入记号T0(h)=T(h),T1(h)=S(h),T2(h)=C(h),T3(h)

26、=R(h)等,从而可将上述公式写出统一形式,用m(h)作为I的近似值,误差量级为O(h2(m+1).,经过m(m=1,2,)次加速后，余项便取下列形式,这种处理方法通常称为理查森(Richardson)外推加速方法.,公式(4.11)也称为龙贝格求积算法.,以0(k)表示二分k次后求得的梯形值，以m(k)表示序列0(k)的m次加速值，则依递推公式(4.9)可到,龙贝格求积算法的计算过程如下：,(1)取k=0,h=b-a，求,令1k(k记区间a,b的二分次数).,(2)求值，按梯形递推公式计算0(k).,(3)求计算值，按加速公式逐个求出数表的第k行其余各元素j(k-j)(j=1,2,k).,(

27、4)若|k(0)-k-1(0)|0，计算积分,的近似值.先取步长h=b-a，应用辛普森公式有,其中,若把区间a,b对分，步长h2=h/2=(b-a)/2，在每个小区间上用辛普森公式，则得,其中,实际上(5.2)式即为,与(5.1)式比较，若f(4)(x)在区间(a,b)上变化不大，可假定f(4)()f(4)()，从而可得,与(5.2)式比较，则得,这里S1=S(a,b)，S2=S2(a,b).如果有,则可期望得到,此时可取S2(a,b)作为I(f)的近似，则可达到给定的误差精度，若不等式(5.3)不成立，则应分别对子区间a,(a+b)/2及(a+b)/2,b再用辛普森公式，此时步长h3=(1/

28、2)h2，得到S3(a,(a+b)/2)及S3(a+b)/2,b).只要分别考察下面两个不等式,是否成立.对满足要求的区间不再细分，对不满足要求的还要继续上述过程，直到满足要求为止，最后还要应用龙贝格法则求出相应区间的积分近似值.为了更直观地说明自适应积分法的计算过程及方法为何能节省计算量，看p115的例7讲解.,4.6高斯求积公式,由前面的讨论已经知道，以a=x0x11时求解是困难的.只有在节点xk(k=0,1,n)确定以后，方可利用(6.5)式求解Ak(k=0,1,n)，此时(6.5)式为关于Ak的线性方程组.下面先讨论如何选取节点xk(k=0,1,n)才能使求积公式(6.4)具有2n+1

29、次代数精度.,设a,b的n+1个节点ax0x10.由此得出h=hopt时E(h)最小.当f(x)=x1/2时,f(x)=(3/8)x-5/2,假定=(1/2)10-4，由,这与表4-10基本相符.,得,4.8.2插值型的求导公式,对于列表合适y=f(x)：,运用插值原理，可以建立插值多项式y=Pn(x)作为f(x)的近似.由于多项式的导数比较容易，我们取Pn(x)的值作为f(x)的近似值，这样建立的数值公式,统称为插值型的求导公式.,必须指出，即使f(x)与Pn(x)的值相差不多，导数的近似值Pn(x)与导数的真值f(x)仍然可能差别很大,因而在使用求导公式(8.3)时应特别注意误差的分析.,

30、利用插值多项式求导数的方法,f(x)在这些点上的函数值为f(xk),插值多项式为Pn(x),则,设,这是在节点xk处导数的近似值，误差为,因为在节点xk处有,所以在节点xk处求导有,是可以进行估计的，但是对任一非节点x处的导数误差f(x)-Pn(x)是无法估计的，因为我们无法对余项中的一项,做出进一步的说明.,下面我们仅仅考察节点处的导数.为简化讨论，假定所给的节点是等距的.,1.两点公式(n=1),设已给出两个节点x0,x1上的函数值f(x0),f(x1)，做线性插值公式,对上式两端求导，记x1-x0=h，有,于是有下列求导公式,而利用余项公式(8.4)知，带余项的两点公式是,2.三点公式(n=2),设已给出三个节点x0,x1=x0+h,x2=x0+2h上的函数值f(x0),f(x1),f(x2)，做二次插值公式,令x=x0+th，上式可表示为,上式两端对t求导，有,上式左端的导数表示对x求导数,上式分别取t=0,1,2,得到三种三点公式为