数值积分与数值微分课件

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

2、1km，求该卫星的轨道长度。本问题可用椭圆参数方程来描述人造地球卫星的轨道，式中a, b分别为椭圆的长短轴，该轨道的长度L就是如下参数方程弧长积分但这个积分是椭圆积分，不能用解析方法计算。5.2问题的描述与基本概念要想用计算机来计算，应对其做离散化处理。注意到定积分是如下和式的极限要离散化，做1） 去掉极限号2） 将取为具体的值3） 为减少离散化带来的误差，将用待定系数代替于是就得到定义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，使它具有尽可能高的代数精度，并指出相应的代数精度。解 本题要先求出具体的求积公式，然后再判断所求公式的代数精度。公式有3

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

5、项，有这里。两边取积分,有记 （5.3）则有 （5.4）若舍去，得求积公式求积系数的求积公式就是插值型求积公式。插值型求积公式的求积余项当为次数小于n次的多项式时，有，对应的。因此插值型求积公式的代数精度至少为n-1。若取，代入式(5.4)，可得插值型求积公式的求积系数之和为下面具体介绍常用的几个插值型求积公式。2. Newton-Cotes求积公式1) n点的Newton-Cotes公式的构造将求积节点取为a,b上的等距节点做积分变量变换: 则当时,有,于是有插值型求积公式的求积系数为记，则有常称为Cotes系数，易验证通常称 （5.6）为n点的Newton-Cotes公式。由于求积节点是等

6、距的，因此也称式（5.6）为等距节点求积公式。利用可以得出下面常用的Newton-Cotes公式A) 2 点的Newton-Cotes公式 （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 418407280-350例5.3 试分别用梯形公式和Simpson公

7、式计算解 用梯形公式计算，有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。证明 由于是插值型求积公式，故有对有记 ，易知，故s是奇函数，得，得证。3） 梯形公式与Simpson公式的余项引理 5.1 (积分中值定理)设，且在上不变号，则有梯形公式余项为在a,b不变号，有2阶连续导数，由引理 5.1，有梯形公式余项 （5.9）

8、抛物线公式的余项 （5.10）4）Newton-Cotes公式的数值稳定性设计算函数时产生舍入误差为实际在计算机中参加计算的是的近似值故Newton-Cotes公式在计算机中产生的误差为若记，则有由Cotes表5.1,当n8时, ，从而有说明此时计算舍入误差可以控制，从而Newton-Cotes公式是数值稳定的。但当n8时，有正有负， 随n增大而增大，从而导致舍入误差增加。故n8时，Newton-Cotes公式是数值不稳定的。因而一般不用n8的Newton-Cotes公式来做定积分计算。3. Gauss求积公式1）Gauss求积公式的构造与概念n点的插值型求积公式的代数精度至少是n-1 ,那么

9、是否还能提高其代数精度呢？若能，其代数精度最大能是多少？为回答这个问题，观察一下插值求积公式的构造方法，发现其至少具有n-1次代数精度的结论是在限定求积节点的前提下得出的，若让求积节点也可以自由取值，则就给提高代数精度创造了条件。为使问题讨论更具一般性，这里考虑带权的定积分求积公式 （5.11）式中是已知的非负函数，为区间a,b上的权函数，a,b可以取为。显然在式（5.11）中，若，就是前面讨论的求积公式。定理5.2 求积公式(5.11)的代数精度最大为。 证明 设是求积公式(5.11)的任意一组求积节点，用此节点构造插值型求积公式(5.11)，并令取2n次多项式代入公式(5.11), 考虑它

10、的求积余项，有因为是次多项式，由代数精度定义得的代数精度不大于。为证求积公式的代数精度可以为，设是任意一个次多项式，用去除，由多项式除法有式中都是次数不大于的多项式，于是有(5.12)因为上面的求积公式(5.11)是具有n个节点的插值型求积公式，故其代数精度不小于，从而有要让式（5.11）成为等式，由式（5.12）应有 （5.13）式（5.13）要求对固定的n次多项式和任意次数不超过的多项式都成立，这样可以用的这种任意性，选择求积节点。由正交多项式理论可知：使式（5.13）成立的点 是存在唯一的，且都在内，它就是在关于权的次正交多项式的零点。于是选取这些点作为求积公式的求积节点，并构造对应的插

11、值型求积公式，就得到具有次迭代精度的求积公式了。关于插值型求积公式的结论定理5.3 点插值型求积公式的代数精度至少是，至多为。定义5.3 若点的求积公式具有次代数精度，则称该求积公式为Gauss型求积公式，对应的求积节点和求积系数分别称为Gauss点和Gauss系数。例5.4确定参数，使求积公式成为Gauss求积公式。解 注意到被积函数中有一因式与求积公式右端无明显的关系，可将其视为权函数。为确定四个参数，依次取代入公式并将近似号取为等号，得联立方程组解出由于本题求积节点个数为，其最高代数精度，故所求的求积公式代数精度至少是3，故它是Gauss公式。2) Gauss型求积公式的求积余项设,取在

12、Gauss点上的次Hermite插值多项式,由Hermite插值余项公式，有两边积分，有由Gauss型求积公式的代数精度为及积分中值定理有，Gauss求积余项为3) Gauss型求积公式的数值稳定性对任意i, 取是关于Gauss点的次Lagrange插值函数，有,由Gauss公式的代数精度为，而是次多项式，有由i的任意性可知，Gauss求积系数。 若在Gauss公式中取，可得类似Newton-Cotes稳定性处理方法，有在计算机计算时的舍入误差这说明Gauss型求积公式是稳定的。在Gauss公式中，不同的权函数和不同积分区间，对应不同形式的Gauss公式，但最基本和常用的Gauss公式有4种。

13、4） 常用的Gauss型求积公式（1） Gauss-Legendre求积公式权函数，积分区间为, Gauss点为次Legendre正交多项式的零点，Gauss-Legendre求积公式为Gauss-Legendre求积余项为Gauss-Legendre求积公式的Gauss点与系数表 5.2nx kA k nx kA k 10.000 000 02.000 000 040.861 136 30.347 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.

14、774 596 70.555 555 60.906 179 80.236 926 90.538 469 30.478 628 7Gauss-Legendre求积公式可以计算任何有限积分区间的定积分，计算之前先作变量代换将积分区间变到，有然后再对用Gauss-Legendre求积公式。（2） Gauss-Chebyshev求积公式权函数,求积区间为，Gauss点为次Chebyshev正交多项式的零点，Gauss-Chebyshev求积公式为式中Gauss点与系数Gauss-Chebyshev求积余项（3） Gauss-Laguerre求积公式权函数,积分区间为，Gauss点为次Laguerre正

15、交多项式的零点，Gauss-Laguerre求积公式为Gauss-Laguerre求积余项为Gauss-Laguerre求积公式的Gauss点和系数也有事先计算好的表。（4）Gauss-Hermite求积公式权函数,积分区间(-,+), Gauss点为n次Hermite正交多项式的零点，Gauss-Hermite求积公式为Gauss-Hermite求积余项为Gauss-Hermite求积公式的Gauss点及系数有表。例5.5用两点Gauss公式求定积分的近似值。解 本题为有限区间的定积分，可用两点Gauss-Legendre 求积公式计算。做积分换元，将其化为-1,1上的定积分，即令有本题准确

16、值为I=3.140 52208，可见精度很高。5.4复化求积公式Newton-Cotes公式在n8时数值不稳定，因此不能用增加求积节点的方法来提高计算精度。实用中常用复化求积公式来求积区间a,b上的定积分，以获得满足给定计算精度要的定积分值。常用的复化求积公式有复化梯形公式和复化Simpson公式。 基本思想将求积区间a,b分成若干个小区间，然后在每个小区间上采用数值稳定的Newton-Cotes公式求小区间上的定积分，最后把所有小区间上的计算结果相加起来作为原定积分的近似值。1. 复化梯形公式1）复化梯形公式的构造原理取等距节点将积分区间a,b n等分,在每个小区间上用梯形公式做近似计算，就

17、有得求积公式-复化梯形公式： (5.14)2）复化梯形公式的余项记(5.15)故复化梯形公式的求积余项 由此可知,复化梯形公式的代数精度是1。若，对给定计算精度，令得出说明利用复化求积公式能得到计算误差小于的定积分近似值。2. 复化Simpson公式1） 复化Simpson公式的构造原理取等距节点将积分区间a,b n等分,在每个小区间上用Simpson公式做近似计算，再累加起来就有式中，得复化Simpson公式 (5.16)2）复化Simpson公式的余项记有复化Simpson公式的求积余项从复化simpson公式的求积余项可以看出复化simpson公式的代数精度是3，它在代数精度和计算精度上

18、都比复化梯形公式好。复化Simpson公式也称为复化抛物线公式。例5.6分别用复化梯形公式和复化Simpson公式计算，要求误差不超过。解 为较快得结果，积分区间分割数按进行。数值计算结果列表，其中代表求积余项。N复化梯形公式复化Simpson公式TnRnSnRn2-17.389 2595.32-11.592 840-0.47822-13.336 0231.27-11.984 944-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.

19、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应分成多少个子区间。以此计算积分近似值。解 复化Simpson公式的求积余项为 式中。为估计误差，要计算。注意到，故由此得从而有解出，故要求出满足计算精度要求的定积分值，只要将0,1分成4个子区间即可，此时可算出5.5 Romberg 求积方法Rom

20、berg 求积方法是对复化梯形公式用加速技术得到的一种求积方法，它也称为逐次分半加速收敛法。基本思想将Richardson 外推算法应用于复化梯形公式中，用产生的加速数列来求定积分值。1.构造原理定理 5.4 设函数逼近数的余项为式中都是与h无关的常数，且时，,则由 ，q为任意常数定义的函数也逼近，且有式中都是与h无关的非零常数。证明 用qh 代替余项式（5.18）的变元h，有(5.19)用乘式（5.18）并与式(5.19)相加，整理后，有(5.20)因为0q 1，故，用同除式(5.20)，有令 即得所证。由微积分收敛阶概念可知的收敛阶比的收敛阶大，故比逼近的速度快。通常称如上的方法为Rich

21、ardson 外推法。显然这种外推可以不断做下去以获得逼近更快的函数，一般有这样，每用一次Richardson外推，就使逼近阶提高一个等级，从而达到加速的目的。分析记 称为复化梯形公式的T形值，可以证明：此式说明逼近的阶为。利用Richardson 外推法对做加速。因为,有逼近的阶变为。若取，则有同理对再做一次Richardson 加速，有逼近I的阶为。一般，有经Richardson加速的求定积分的序列为 (5.21)逼近I的阶变为。 注意到，直接计算可知就是复化Simpson公式的，即有 (5.22)类似地可得 (5.23) (5.24)这样可以用复化梯形公式计算出序列T2k。再逐次用公式（

22、5.22），（5.23），（5.24）对其进行加工得到收敛更快的序列S2k，C2k，R2k等等。通常将序列T2k，S2k，C2k和R2k依次称为梯形序列、Simpson序列、Cotes序列和Romberg序列。若继续外推，舍入误差增大将使新序列体现不出更明显的加速效果，因此一般只外推到Romberg序列。用Romberg序列求定积分的算法称为Romberg求积方法。3Romberg求积方法的计算过程kT2kS2kC2kR2k0T1S1C1R11T2S2C2R22T4S4C43T8S84T16因为与有相重的点，计算时还可以利用下面关系式来减小计算量。在给定计算误差界后，可用R2k+1-R2k来作

23、为终止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.946 08340.946 083030.945 69090.946 083040.9845 9850由 ，故有。5.6 数值微分根据函数在若干个点处的函数值去求该函数的导数近似值称为数值微分，所求导数的近似值常称为数值导数。设是已知的个点处

24、的函数值，下面分别讨论两种求数值微分的方法。1.利用次多项式插值函数求数值导数用插值函数代替被插值函数可以用来计算被插函数的近似值；用插值函数代替定积分的被积函数可以用来计算定积分的近似值；用插值函数的导数能否用来求被插函数导数值？定义5.4 设是的次插值多项式，称用插值函数的导数来计算被插函数导数的公式 为的插值型求导公式。由插值理论，有特别有一阶数值导数的余项关系在非节点处的余项由于涉及对未知中值函数的求导，使余项不易描述和控制，但在节点处有此余项易得出误差的界限，因此插值型求导公式多用于求在节点处的数值导数。下面给出几个常用的数值微分公式。1) 两点公式给定点，有线性插值多项式记，有不带余项的数值微分公式在节点处的余项 带余项的两点数值微分公式前差公式 后差公式 2) 三点公式只给出等距节点的公式形式，对非等距节点有类似的公式。给定点集，有二次插值多项式，类似前面的做法，可得及带余项的三点数值微分公式前差公式 中心公式 后差公式 中间的节点使用中心公式，左边界节点用前插公式，右边界点用后插公式。在求数值微分公式余项表达