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

1、PAGE PAGE 163第5章 数值积分与数值微分方法关于定积分计算，已经有较多方法，如公式法、分步积分法等，但实际问题中，经常出现不能用通常这些积分方法计算的定积分问题。怎样把这些通常方法失效的定积分在一定精度下快速计算出来，特别是通过计算机编程计算出来就是本章研究的内容。此外，怎样根据函数在若干个点处的函数值去求该函数的导数近似值也是本章介绍的内容。本章涉及的方法有Newton-Cotes求积公式、Gauss求积公式、复化求积公式、Romberg求积公式和数值微分。5.1 引 例人造地球卫星轨道可视为平面上的椭圆。我国的第一颗人造地球卫星近地点距离地球表面439km，远地点距地球表面23

2、84km，地球半径为6371km，求该卫星的轨道长度。本问题可用椭圆参数方程来描述人造地球卫星的轨道，式中a, b分别为椭圆的长短轴，该轨道的长度L就是如下参数方程弧长积分但这个积分是椭圆积分，不能用解析方法计算。5.2问题的描述与基本概念要想用计算机来计算，应对其做离散化处理。注意到定积分是如下和式的极限要离散化，做去掉极限号将取为具体的值为减少离散化带来的误差，将用待定系数代替于是就得到定义5.1 若存在实数且任取都有 （5.1）则称式(5.1)为一个数值求积公式。式中称为求积系数，称为求积节点，而称 （5.2）为求积余项或求积公式（5.1）的截断误差。从定义可以看到，数值求积公式依赖于求

3、积节点个数n、求积节点和求积系数，这三个量有一个发生变化，则产生不同的求积公式。定义5.2 若求积公式对所有不超过m次的多项式有求积余项，而对某一个m+1次多项式有，则称该求积公式的代数精度为m。一般，一个求积公式的代数精度越大，则该求积公式越好。确定代数精度的方法依次取代入公式并验证是否成立。若第一个使不成立的k值为m，则对应的代数精度为m-1。例 5.1确定求积公式的代数精度。解 取代入求积公式有易验证，但，故本题求积公式代数精度为3。例 5.2确定下面求积公式的参数A，B，C，使它具有尽可能高的代数精度，并指出相应的代数精度。解 本题要先求出具体的求积公式，然后再判断所求公式的代数精度。

4、公式有3个待定参数，h不是求积公式的参数，故利用3个条件得到的3个等式关系就可以解决求出具体求积公式的问题。依次取代入求积公式并取等号，有A+B+C=4hA-2C=0A+4C=163h解之得A=169h , B=43h , C=89 h故所求的求积公式为为确定其代数精度，再取代入求出的公式继续计算，有，故所求的求积公式具有二阶代数精度。5.3 插值型求积公式借助多项式插值函数来构造的求积公式称为插值型求积公式。一般选用不同的插值公式就可以得到不同的插值型求积公式。基本思想利用被积函数 QUOTE fx 的插值函数x代替 QUOTE fx 做定积分的近似计算来构造求积公式。1.构造原理考虑 QU

5、OTE fx 在n个节点 QUOTE x1,x2,xn 上的n-1次Lagrange插值多项式 QUOTE Ln-1(x) 与 QUOTE fx 的余项，有这里 QUOTE lin-1x=k=1kinx-xkxi-xk, nx=k=1n(x-xk) 。两边取积分,有记 （5.3）则有 （5.4）若舍去，得求积公式求积系数的求积公式就是插值型求积公式。插值型求积公式的求积余项当为次数小于n次的多项式时，有 QUOTE fn(x)0 ，对应的。因此插值型求积公式的代数精度至少为n-1。若取，代入式(5.4)，可得插值型求积公式的求积系数之和为下面具体介绍常用的几个插值型求积公式。2. Newton

6、-Cotes求积公式1) n点的Newton-Cotes公式的构造将求积节点 QUOTE xi 取为a,b上的等距节点做积分变量变换: 则当 QUOTE xa,b 时,有 QUOTE t0,n-1 ,于是有插值型求积公式的求积系数为记 QUOTE ci(n-1)=1n-10n-1k=1k1nt-k+1i-kdt ，则有 QUOTE cin-1 常称为Cotes系数，易验证通常称 （5.6）为n点的Newton-Cotes公式。由于求积节点 QUOTE xi 是等距的，因此也称式（5.6）为等距节点求积公式。利用可以得出下面常用的Newton-Cotes公式A) 2 点的Newton-Cotes

7、公式 （5.7）这正是我们熟悉的梯形公式。B) 3点的Newton-Cotes公式为 （5.8）称它为Simpson公式或抛物线公式。表5.1 部分Cotes系数n2345678912 1216 46 1618 38 38 18790 1645 215 1645 79019288 2596 25144 25144 2596 1928841840 935 9280 34105 9280 935 41840751172803577172801323172802989172802989172801323172803577172807511728098928350588828350-928283501

8、049628350-4540283501049628350-9282835058882835098928350例5.3 试分别用梯形公式和Simpson公式计算解 用梯形公式计算，有I1-02f0+f1=121+sin10.920 735 49用Simpson公式计算，有I1-06f0+4f12+f1=161+8sin12+sin10.946 145 882）n点Newton-Cotes公式的代数精度 定理5.1 当求积节点个数n为奇数时，对应的Newton-Cotes求积公式的代数精度至少为n。证明 由于是插值型求积公式，故有对 QUOTE xn 有记 QUOTE s=k=12l+1s+l-

9、k+1=m=-lls+m, m=l-k+1 ，易知，故s是奇函数，得 QUOTE Rxn=0 ，得证。3） 梯形公式与Simpson公式的余项引理 5.1 (积分中值定理)设，且在上不变号，则有梯形公式余项为在a,b不变号，有2阶连续导数，由引理 5.1，有梯形公式余项 （5.9）抛物线公式的余项 （5.10）4）Newton-Cotes公式的数值稳定性设计算函数 QUOTE fxk 时产生舍入误差为 QUOTE k 实际在计算机中参加计算的是 QUOTE fxk 的近似值故Newton-Cotes公式在计算机中产生的误差为若记 QUOTE =max1kn|k| ，则有由Cotes表5.1,当

10、n8时, ， QUOTE ck(n-1)0 从而有说明此时计算舍入误差可以控制，从而Newton-Cotes公式是数值稳定的。但当n8时， QUOTE ckn-1 有正有负， QUOTE k=1n|ckn-1| 随n增大而增大，从而导致舍入误差增加。故n8时，Newton-Cotes公式是数值不稳定的。因而一般不用n8的Newton-Cotes公式来做定积分计算。3. Gauss求积公式1）Gauss求积公式的构造与概念n点的插值型求积公式的代数精度至少是n-1 ,那么是否还能提高其代数精度呢？若能，其代数精度最大能是多少？为回答这个问题，观察一下插值求积公式的构造方法，发现其至少具有n-1次

11、代数精度的结论是在限定求积节点 QUOTE xi 的前提下得出的，若让求积节点 QUOTE xi 也可以自由取值，则就给提高代数精度创造了条件。为使问题讨论更具一般性，这里考虑带权的定积分求积公式 （5.11）式中 QUOTE x 是已知的非负函数，为区间a,b上的权函数，a,b QUOTE a,b 可以取为 QUOTE 。显然在式（5.11）中，若 QUOTE x1 ，就是前面讨论的求积公式。定理5.2 求积公式(5.11)的代数精度最大为 QUOTE 2n-1 。 证明 设 QUOTE x1,x2,xn 是求积公式(5.11)的任意一组求积节点，用此节点构造插值型求积公式(5.11)，并令

12、取2n次多项式代入公式(5.11) QUOTE fx=wn2x , 考虑它的求积余项，有因为 QUOTE wn2x 是 QUOTE 2n 次多项式，由代数精度定义得的代数精度不大于 QUOTE 2n-1 。为证求积公式的代数精度可以为 QUOTE 2n-1 ，设 QUOTE f(x) 是任意一个 QUOTE 2n-1 次多项式，用 QUOTE wnx 去除 QUOTE f(x) ，由多项式除法有 QUOTE fx=qxwnx+rx,fxi=rxi (i=1,2,n) 式中 QUOTE qx，rx 都是次数不大于 QUOTE n-1 的多项式，于是有(5.12)因为上面的求积公式(5.11)是具

13、有n个节点的插值型求积公式，故其代数精度不小于 QUOTE n-1 ，从而有要让式（5.11）成为等式，由式（5.12）应有 （5.13）式（5.13）要求对固定的n次多项式和任意次数不超过 QUOTE n-1 的多项式 QUOTE qx 都成立，这样可以用 QUOTE qx 的这种任意性，选择求积节点 QUOTE x1,x2,xn。 。由正交多项式理论可知：使式（5.13）成立的点 QUOTE x1,x2,xn 是存在唯一的，且都在 QUOTE a,b 内，它就是在 QUOTE a,b 关于权 QUOTE x 的 QUOTE n 次正交多项式的零点。于是选取这些点作为求积公式的求积节点，并构

14、造对应的插值型求积公式，就得到具有 QUOTE 2n-1 次迭代精度的求积公式了。关于插值型求积公式的结论定理5.3 QUOTE n 点插值型求积公式的代数精度至少是 QUOTE n-1 ，至多为 QUOTE 2n-1 。定义5.3 若点的求积公式具有 QUOTE 2n-1 次代数精度，则称该求积公式为Gauss型求积公式，对应的求积节点 QUOTE xi 和求积系数分别称为Gauss点和Gauss系数。例5.4确定参数 QUOTE x1,x2,A1,A2 ，使求积公式成为Gauss求积公式。解 注意到被积函数中有一因式与求积公式右端无明显的关系，可将其视为权函数。为确定四个参数，依次取代入公

15、式并将近似号取为等号，得联立方程组解出由于本题求积节点个数为 QUOTE n=2 ，其最高代数精度 QUOTE 2n-1=3 ，故所求 QUOTE x1,x2,A1,A2 的求积公式代数精度至少是3，故它是Gauss公式。2) Gauss型求积公式的求积余项设 QUOTE f(x)C2na,b ,取 QUOTE f(x) 在Gauss点 QUOTE x1,x2,xn 上的 QUOTE 2n-1 次Hermite插值多项式 QUOTE H2n-1(x) ,由Hermite插值余项公式，有两边积分，有由Gauss型求积公式的代数精度为 QUOTE 2n-1 及积分中值定理有，Gauss求积余项为3

16、) Gauss型求积公式的数值稳定性对任意i, 取 QUOTE fx=lin-12x, lin-1x 是关于Gauss点 QUOTE x1,x2,xn 的 QUOTE n-1 次Lagrange插值函数，有 QUOTE lin-1xj=ij ,由Gauss公式的代数精度为 QUOTE 2n-1 ，而 QUOTE fx 是 QUOTE 2n-2 次多项式，有由i的任意性可知，Gauss求积系数。 QUOTE Ai0 若在Gauss公式中取 QUOTE fx1 ，可得类似Newton-Cotes稳定性处理方法，有在计算机计算时的舍入误差这说明Gauss型求积公式是稳定的。在Gauss公式中，不同的

17、权函数 QUOTE x 和不同积分区间，对应不同形式的Gauss公式，但最基本和常用的Gauss公式有4种。4） 常用的Gauss型求积公式Gauss-Legendre求积公式权函数 QUOTE x1 ，积分区间为 QUOTE -1,1 , Gauss点为 QUOTE n 次Legendre正交多项式的零点，Gauss-Legendre求积公式为Gauss-Legendre求积余项为Gauss-Legendre求积公式的Gauss点 QUOTE xi 与系数 QUOTE Ai 表 5.2nx kA k nx kA k 10.000 000 02.000 000 040.861 136 30.3

18、47 854 820.577 350 31.000 000 00.339 981 00.652 145 230.000 000 00.888 888 950.000 000 00.568 888 90.774 596 70.555 555 60.906 179 80.236 926 90.538 469 30.478 628 7Gauss-Legendre求积公式可以计算任何有限积分区间的定积分，计算之前先作变量代换将积分区间 QUOTE a,b 变到 QUOTE -1，1 ，有然后再对用Gauss-Legendre求积公式。Gauss-Chebyshev求积公式权函数 QUOTE x=11-

19、x2 ,求积区间为 QUOTE -1，1 ，Gauss点为 QUOTE n 次Chebyshev正交多项式的零点，Gauss-Chebyshev求积公式为式中Gauss点 QUOTE xi 与系数Gauss-Chebyshev求积余项Gauss-Laguerre求积公式权函数 QUOTE x ,积分区间为 QUOTE 0，+ ，Gauss点为 QUOTE n 次Laguerre正交多项式的零点，Gauss-Laguerre求积公式为Gauss-Laguerre求积余项为Gauss-Laguerre求积公式的Gauss点和系数也有事先计算好的表。（4）Gauss-Hermite求积公式权函数 Q

20、UOTE x=e-x2 ,积分区间(-,+), Gauss点为n次Hermite正交多项式的零点，Gauss-Hermite求积公式为Gauss-Hermite求积余项为Gauss-Hermite求积公式的Gauss点及系数 QUOTE Ak 有表。例5.5用两点Gauss公式求定积分的近似值。解 本题为有限区间的定积分，可用两点Gauss-Legendre 求积公式计算。做积分换元，将其化为-1,1上的定积分，即令有本题准确值为I=3.140 52208，可见精度很高。5.4复化求积公式Newton-Cotes公式在n8时数值不稳定，因此不能用增加求积节点的方法来提高计算精度。实用中常用复化

21、求积公式来求积区间a,b上的定积分，以获得满足给定计算精度要的定积分值。常用的复化求积公式有复化梯形公式和复化Simpson公式。 基本思想将求积区间a,b分成若干个小区间，然后在每个小区间上采用数值稳定的Newton-Cotes公式求小区间上的定积分，最后把所有小区间上的计算结果相加起来作为原定积分的近似值。复化梯形公式1）复化梯形公式的构造原理取等距节点 QUOTE xk=a+kh, h=b-an(k=0,1,n) 将积分区间a,b n等分,在每个小区间 QUOTE xk,xk+1(k=0,1,n-1) 上用梯形公式做近似计算，就有得求积公式复化梯形公式： (5.14)2）复化梯形公式的余

22、项记(5.15)故复化梯形公式的求积余项 由此可知,复化梯形公式的代数精度是1。若 QUOTE |fx|M2 ，对给定计算精度 QUOTE ，令得出说明利用复化求积公式能得到计算误差小于的定积分近似值。2. 复化Simpson公式复化Simpson公式的构造原理取等距节点 QUOTE xk=a+kh, h=b-an(k=0,1,n) 将积分区间a,b n等分,在每个小区间 QUOTE xk,xk+1(k=0,1,n-1) 上用Simpson公式做近似计算，再累加起来就有式中 QUOTE xk+12=xk+h2 ，得复化Simpson公式 (5.16)2）复化Simpson公式的余项记有复化Si

23、mpson公式的求积余项从复化simpson公式的求积余项可以看出复化simpson公式的代数精度是3，它在代数精度和计算精度上都比复化梯形公式好。复化Simpson公式也称为复化抛物线公式。例5.6分别用复化梯形公式和复化Simpson公式计算，要求误差不超过 QUOTE 0.510-3 。解 为较快得结果，积分区间分割数按 QUOTE 2n(n=0,1,) 进行。数值计算结果列表，其中 QUOTE Rn 代表求积余项。N复化梯形公式复化Simpson公式TnRnSnRn2-17.389 2595.32-11.592 840-0.47822-13.336 0231.27-11.984 944

24、-0.85410-123-12.382 1620.312-12.064 209-0.61410-224-12.148 0040.77710-1-12.069 951-0.39510-325-12.089 7420.19410-126-12.075 1940.48510-227-12.071 5580.12110-228-12.070 6490.30310-3本题积分的准确值为I=-12.070 346 316 4，可见复化梯形公式和复化Simpson公式能求出精度较高的解。例5.7考虑用复化Simpson公式计算要使误差小于0.510-6，那么求积区间0,1应分成多少个子区间。以此计算积分近似

25、值。解 复化Simpson公式的求积余项为 QUOTE Rf, Sn =-1-02590h4f4() QUOTE (0,1) 式中 QUOTE h=1-0n=1n, fx=sinxx 。为估计误差，要计算 QUOTE f4(x) 。注意到，故由此得从而有解出 QUOTE n10472 ，故要求出满足计算精度要求的定积分值，只要将0,1分成4个子区间即可，此时可算出5.5 Romberg 求积方法Romberg 求积方法是对复化梯形公式用加速技术得到的一种求积方法，它也称为逐次分半加速收敛法。基本思想将Richardson 外推算法应用于复化梯形公式中，用产生的加速数列来求定积分值。1.构造原理

26、定理 5.4 设函数 QUOTE F1(h) 逼近数 QUOTE F* 的余项为 QUOTE F*-F1h=1hP1+2hP2+khPk+ (0P1P2) 式中 QUOTE Pk,k 都是与h无关的常数，且 QUOTE k1 时，,则由 ，q为任意常数定义的函数 QUOTE F2h 也逼近 QUOTE F* ，且有式中 QUOTE k1(k2) 都是与h无关的非零常数。证明 用qh 代替余项式（5.18）的变元h，有(5.19)用 QUOTE -qP1 乘式（5.18）并与式(5.19)相加，整理后，有(5.20)因为0q 1，故 QUOTE 1-qP10 ，用 QUOTE 1-qP1 同除式

27、(5.20)，有令 QUOTE F2h=F1qh-qP1F1h1-qP1 QUOTE k1=kqPk-qP11-qP1 即得所证。由微积分收敛阶概念可知 QUOTE F2h 的收敛阶 QUOTE P2 比的收敛阶 QUOTE P1 大，故 QUOTE F2h 比 QUOTE F1h 逼近的速度快。通常称如上的方法为Richardson 外推法。显然这种外推可以不断做下去以获得逼近更快的函数，一般有 QUOTE Fm+1h=Fmqh-qPmFmh1-qPm (m=1,2,3) QUOTE F*-Fm+1h=m+1(m)hPm+1+m+1(m)hPm+2+ 这样，每用一次Richardson外推，

28、就使逼近阶提高一个等级，从而达到加速的目的。分析记 称 QUOTE I=abf(x)dx 为复化梯形公式的T形值，可以证明：此式说明 QUOTE T2h 逼近 QUOTE I 的阶为 QUOTE O(h4) 。利用Richardson 外推法对做加速。因为,有 QUOTE T2h=T1qh-q2T1h1-q2 QUOTE T2h 逼近的阶变为 QUOTE O(h4) 。若取 QUOTE q=12 ，则有同理对 QUOTE T2h 再做一次Richardson 加速，有逼近I的阶为。一般，有经Richardson加速的求定积分的序列为 (5.21) QUOTE Tmh=4m-1Tm-1h2-Tm

29、-1(h)4m-1-1 QUOTE T2h 逼近 QUOTE I I的阶变为 QUOTE O(h4) 。 注意到 QUOTE T1h2=T2n ，直接计算可知 QUOTE T2h 就是复化Simpson公式的 QUOTE Sn ，即有 QUOTE Sn=4T2n-Tn4-1 (5.22)类似地可得 QUOTE Cn=42S2n-Sn42-1 (5.23) QUOTE Rn=43C2n-Cn43-1 (5.24)这样可以用复化梯形公式计算出序列T2k。再逐次用公式（5.22），（5.23），（5.24）对其进行加工得到收敛更快的序列S2k，C2k，R2k等等。通常将序列T2k，S2k，C2k和R

30、2k依次称为梯形序列、Simpson序列、Cotes序列和Romberg序列。若继续外推，舍入误差增大将使新序列体现不出更明显的加速效果，因此一般只外推到Romberg序列。用Romberg序列求定积分的算法称为Romberg求积方法。3Romberg求积方法的计算过程kT2kS2kC2kR2k0 eq oac(,1)T1 eq oac(,3)S1 eq oac(,6)C1 eq oac(,10)R11 eq oac(,2)T2 eq oac(,5)S2 eq oac(,9)C2 eq oac(,14)R22 eq oac(,4)T4 eq oac(,8)S4 eq oac(,13)C43 e

31、q oac(,7)T8 eq oac(,12)S84 eq oac(,11)T16因为 QUOTE Tn 与 QUOTE T2n 有相重的点，计算时还可以利用下面关系式来减小计算量。在给定计算误差界 QUOTE 后，可用R2k+1-R2k来作为终止Romberg计算的条件。例5.8 用Romberg算法计算计算结果精确到0.510-6。解 由Romberg求积方法，得计算结果如下表kT2kS2kC2kR2k00.920 73550.946 14590.946 08300.946 083110.939 79330.946 08690.946 08310.946 083020.944 51350.

32、946 08340.946 083030.945 69090.946 083040.9845 9850由 ，故有。5.6 数值微分根据函数在若干个点处的函数值去求该函数的导数近似值称为数值微分，所求导数的近似值常称为数值导数。设是已知的个点处的函数值，下面分别讨论两种求数值微分的方法。1.利用次多项式插值函数求数值导数用插值函数代替被插值函数可以用来计算被插函数的近似值；用插值函数代替定积分的被积函数可以用来计算定积分的近似值；用插值函数的导数能否用来求被插函数导数值？定义5.4 设是的次插值多项式，称用插值函数的导数来计算被插函数导数的公式 为的插值型求导公式。由插值理论，有特别有一阶数值导数的余项关系在非节点处的余项由于涉及对未知中值函数的求导，使余项不易描述和控制，但在节点处有此余项