第四章N-C公式

1、第四章 数值积分与数值微分Newton-CotesNewton-Cotes公式公式nkknknxfCabI0)()()(badxxfA)(120001)2)12(2)1(21)(kjkkabjafabkTkT微积分学微积分学- “- “人类精神的卓越胜利人类精神的卓越胜利”微积分就是微分运算和积分运算这两种互逆运算方法微积分就是微分运算和积分运算这两种互逆运算方法的合称，就像加法与减法，乘法与除法是互逆运算一的合称，就像加法与减法，乘法与除法是互逆运算一样，但微积分的运算法则要比加减乘除，乘方，开方样，但微积分的运算法则要比加减乘除，乘方，开方等运算复杂得多，现在已成为高等数学的核心内容。等运

2、算复杂得多，现在已成为高等数学的核心内容。 为什么要数值积分？为什么要数值积分？在微积分里，按在微积分里，按Newton-LeibnizNewton-Leibniz公式求定积分公式求定积分( )( )( )( )baI ff x dxF bF aWhy do we do numerical integral? f(x)没有解析表达式，只有数表形式没有解析表达式，只有数表形式x12345f(x)44.5688.5 f(x)有表达式，但原函数不是初等函数有表达式，但原函数不是初等函数210 xedx10(arctan)x x dx它们的原函它们的原函数都不是初数都不是初等函数等函数. . f(x)

3、原函数表达式很复杂，计算量很大原函数表达式很复杂，计算量很大讨论数值积讨论数值积分的必要性分的必要性更一般地，我们可以在区间更一般地，我们可以在区间a,ba,b上选取某些节点上选取某些节点4.-(1.3)P98-99)()()(2221bfafdxxfAAabbaab 考察其代数精度。考察其代数精度。f(x)abf(a)f(b)梯形公式梯形公式/* trapezoidal rule*/解：解：逐次检查公式是否精确成立逐次检查公式是否精确成立代入代入 P0 = 1： baabdx111 2 ab=代入代入 P1 = x ：=代入代入 P2 = x2 ： 222abbadxx 2baab 3233

4、abbadxx 222baab 代数精度代数精度 = 1( )1( )1baf xI fdxba时，3( )(14 1)()6baIfba22( )( )2baf xxI f时，223( )(22)62babaIfaabb332( )( )3baf xxI f时，332223( )()63babaIfaabb443( )( )4baf xxI f时，344333()( )()624baabbaIfab554( )( )5baf xxI f时，4443()( )()( )64baabIfabI f试确定下面积分公式中的参数使其代数精确度尽量高.220)()0()()0(2)(Ihffahhffh

5、dxxfIhhdxxI00222hI 202232hahhI0)(xxf对于hI 2hhdxxI011)(xxf对于22hhdxxI022)(xxf对于33h3)221(ha2II 令121a3022242hahhIhdxxI033)(xxf对于44h44h4023252hahhIhdxxI044)(xxf对于55h65h3 ,2 , 1 , 0)()(2jxIxIjj)()(424xIxI,)(baCxf设函数等份分割为将积分区间nba,nkkhaxk, 1 ,0,为步长其中nabh各节点为Cotes系数系数)(niC 注：注：Cotes 系数仅取决于系数仅取决于 n 和和 i，可查表得到。

6、与可查表得到。与 f (x) 及区及区间间a, b均无关。均无关。)()(2bfafab)(1fI10)1()()(kkkxfCab)()(210 xfxfab)(1fI即132321883838( )()( )()()( )baa babf x dxb af afff b 3n 时，时，3/8公式公式例例 用用n=6的牛顿柯特斯公式计算定积分值的牛顿柯特斯公式计算定积分值101.dxIx 解：解：将积分区间将积分区间0,1划分为划分为n份，得到节点列为份，得到节点列为100 1 266, , ,ixii在这些节点处的函数值为在这些节点处的函数值为160 1 2616(), , ,iif xi

7、xi 则则n=6的牛顿柯特斯公式为的牛顿柯特斯公式为01230124199341 084035280105419984035280() ()( )()()()( )()If xf xf xf xf xf xf x 41 69 69634 64169696840 6 35 7 280 8 105 9 840 10 35 11 280 12 0 6931. (1) ( ) 0 bnnaf xn(x) I(f) (x)dx , nfP当为次多项式时从而至少有次代数精确度。2,2,2210abhbxabxaxn则取Cotes系数为系数为dtttC20)2(0)2)(1(4161dtttC20)2(1)

8、2(2164dtttC20)2(2)1(4161求积公式为求积公式为2I20)2()()(kkkxfCab)(61)(64)(61)(210 xfxfxfab)()2(4)(6bfbafafab)(2fI上式称为Simpson求积公式，也称三点公式或抛物线公式记为)(2fIS Simpson公式的余项为公式的余项为)()(2IRSRbadxxR)(2)()2(180)4(4fababSimpson公式具有公式具有3次代数精度次代数精度4,4 , 1 , 0, 4abhkkhaxnk则取dtttttC)4)(3( )2)(1(! 44140)4(0907dtttttC)4)(3( )2(! 34

9、140)4(19032dtttttC)4)(3( )1(! 2! 24140)4(29012dtttttC)4)(2( )1(! 34140)4(39032dtttttC)3)(2( )1(! 44140)4(4907)(4fI40)4()()(kkkxfCab)(907)(9032)(9012)(9032)(907)(43210 xfxfxfxfxfab)(7)(32)(12)(32)(79043210 xfxfxfxfxfab上式称为上式称为Cotes求积公式，也称也称五点公式记为)(4fIC )()(4IRCRbadxxR)(4)()4(945)(2)6(6fabab使用使用n次次Lag

10、range插值多项式的插值多项式的Newton-Cotes公式至少具有公式至少具有n次代数精度次代数精度,并且并且n为偶数时至为偶数时至少具有少具有n+1次代数精度次代数精度.n=n=偶数时偶数时Newton-Cotes Newton-Cotes 求积公式的代数精确度求积公式的代数精确度bqx bn)(nknkk)!(n(x)q(x)n11101 1则次多项式对dt)t ()t ( t)thax(dx)x(wdx)x(wdx)x(dx)x(nnnbanbabanbahbb)!n()(qPq)n(n1 102111令积分误差：次代数精确度。有为偶数时至少求积公式当即即：，为偶数时，上式积分为当

11、1n n CotesNewton )()( 0n dxxdxxbanbaPqk)dt)(tk(t)kk)(t(t)(ttn)dt(t)(tt kn212111200k)du)(uk(u)u(u)k(uk)(ukt ukk111 令分为零！点上的对称区间上的积奇函数在以零点为对称令H(u)H(u)(k)u)(ku()uu)()kuk)(u(u) H(k)(uk(u)u(u)kk)(u(uH(u)k1111111 12dtjtknknCnkjnjknnk 00)()()!( !)1(无关与函数的划分有关的节点只与积分区间)(,xfxbaj因此用因此用Newton-CotesNewton-Cotes

12、公式计算积分的舍入误差主要由公式计算积分的舍入误差主要由的计算引起函数值)(kxf其值可以精确给定其值可以精确给定响的舍入误差对公式的影只需讨论)(kxf)()()(,)(计算值的近似值作为而以为精确值假设kkkxfxfxf为误差)()(kkkxfxfnInkknkxfCab0)()()(记)(计算值的近似值为nI而理论值为nInkknkxfCab0)()()(的误差为与nnIInnIInkkknkxfxfCab0)()()()(nnIInkkknkxfxfCab0)()()()(nkknkCab0)()(nkkCab0)()()(kkkxfxfkmax设10nkkC, 1, 0, 70nkk

13、kCCn)(ab的增大而增大且随有正有负nCCnnkkk, 1, 80定义定义2 2 在机械求积公式中，若在机械求积公式中，若00lim()( )nbkkankhA f xf x dx其中其中11max()iii nhxx 则称机械求积公式是则称机械求积公式是的。的。使用机械求积公式计算使用机械求积公式计算()kf x得到的近似值记为得到的近似值记为kf记记kkkfxf)(为误差为误差0( )()nnkkkIfA f x0( )nnkkkIfA f舍入误差舍入误差nkkknkkkknnAfxfAfIfI00)()()(这表明求积公式计算是稳定的。这表明求积公式计算是稳定的。定义定义3 3 对任

14、给对任给只要只要成立，就称机械求积公式是成立，就称机械求积公式是的。的。0若若0kkfxf)(就有就有nkkknkkkknnAfxfAfIfI00)()()(),1 , 0(0nkAk知，由稳定性定义的条件取ab, 0kkfxf)(nkkknkkkknnAfxfAfIfI00)()()()(0abAnkkHavent we had enough formulae? Whats up now?Oh come on, you dont seriously consider h=(b a)/2 acceptable, do you?Why cant you simply refine the par

15、tition if you have to be so picky?Dont you forget the oscillatory nature of high-degree polynomials!Uh-oh固定时而节点个数的长度较大当积分区间1,nba直接使用直接使用Newton-CotesNewton-Cotes公式的余项将会较大公式的余项将会较大增加时即而如果增加节点个数1,n公式的舍入误差又很难得到控制公式的舍入误差又很难得到控制为了提高公式的精度为了提高公式的精度, ,又又使算法简单易行使算法简单易行, ,往往使用往往使用复化方法复化方法分成若干个子区间即将积分区间,ba然后然后在

16、每个小区间上使用低阶在每个小区间上使用低阶Newton-CotesNewton-Cotes公式公式最后将每个小区间上的积分的近似值最后将每个小区间上的积分的近似值相加相加 复化复化 Simpson 公式：公式：),., 0(,nkhkaxnabhk )()(4)(6)(1211 kkkxxxfxfxfhdxxfkkkx21 kx1 kx44444 )()(2)(4)(6)(1010121 nknkkkbabfxfxfafhdxxf= Sn注：注：为方便编程，可采用另一记法：令为方便编程，可采用另一记法：令 n = 2n 为偶数，为偶数， 这时这时 ，有，有hkaxhnabhk ,2 )()(2

17、)(4)(3 koddkevenkknbfxfxfafhS复化求积公式的余项和收敛的阶复化求积公式的余项和收敛的阶我们知道,三个求积公式的余项分别为)(TR)(12)(3fab )(SR)()2(180)4(4fabab)(CR)()4(945)(2)6(6fabab单纯的求积公式单纯的求积公式复化求积公式的每个小区间)(122kfhh )(2180)4(4kfhh)(49452)6(6kfhh )(12103 nkkfh,)(. 12baCxf设被积函数则复合梯形公式的余项为nTI 103)(12nkkfh)(max)()(min10 xfnfxfbxankkbxa 由于使得由介值定理,ba

18、)()(10fnfnkk 即有即有)(12)(2fhab nTI 103)(12nkknfnh)(123fnh 10)4(45)(2180nkkfh,)(. 24baCxf若被积函数nSI 公式的余项为复合足够大时则Simpsonn,)(2180)4(4fhabnTI nSI 阶无穷小量，的分别是42h)(2ho)(4ho)(12)(2fhab )(2180)4(4fhab 例例1 对于函数对于函数f(x)=sinx/x，给出，给出n=8的函数表，的函数表，试用复化梯形公式和复化辛普森公式试用复化梯形公式和复化辛普森公式 计算积分计算积分 10.sindxxxIxf(x)01/81/43/81

19、/25/83/47/8110.99739780.98961580.97672670.95885100.93615560.90885160.87719250.8414709 解解 将积分区间将积分区间0,1划分为划分为8等分，用复化梯形公式求得等分，用复化梯形公式求得.9456909. 08 T而将积分区间而将积分区间0, 1划分为划分为24等等分，用复化辛普森公式求得分，用复化辛普森公式求得.9460832. 04 S 比较上面两个计算结果比较上面两个计算结果T8与与S4，它们都需要提供，它们都需要提供9个点上的函数值，然而精度却差别很大，同积分准个点上的函数值，然而精度却差别很大，同积分准确

20、值确值I=0.9460831比较，应用复化梯形公式计算的结比较，应用复化梯形公式计算的结果果T8=0.9456909只有只有2位有效数字，而应用复化辛普位有效数字，而应用复化辛普森公式计算的结果森公式计算的结果S4= 0.9460832却有却有6位有效数字位有效数字. 为了利用余项公式估计误差，要求为了利用余项公式估计误差，要求f(x)=sinx/x的的高阶导数，由于高阶导数，由于.)cos(sin)(10 dtxtxxxf所以有所以有.2cos)cos()(1010)( dtkxttdtxtdxdxfkkkk 于是于是.112cos)(max10)(10 kdtkxttxfkkx 复化梯形公

21、式误差复化梯形公式误差为为.000434. 03181121)(max12)(210288 xfhTIfRx复化辛普森公式误差复化辛普森公式误差为为.10271. 0514128801)(6444 SIfR 例例2 利用复利用复化梯形公式化梯形公式计算计算 使其误使其误差限为差限为0.5*10-6，应将区间，应将区间0, 1几等分几等分?,sin10 dxxxI 解解 利用例利用例1的结果的结果取取n=236可满足要求可满足要求.由由复化梯形公式的余项得复化梯形公式的余项得.11)(max)(10 kxfkx222410111()10.12122136nR Thfhh 6因此只需将区间因此只需

22、将区间0, 1四等分，即取四等分，即取9个点个点42 h444611( )101802900 22nbahhISf n3.4取取n=4 使用复合使用复合Cotes公式，只需要将区间公式，只需要将区间1等分，即四小等等分，即四小等分，五个点。分，五个点。 前面用复化梯形公式计算此题，满足相同的精度前面用复化梯形公式计算此题，满足相同的精度需要将区间需要将区间0, 1划分划分17等分，可见复化辛普森公式的等分，可见复化辛普森公式的精度的确比复化梯形公式精度高同样也可用精度的确比复化梯形公式精度高同样也可用| S4m- -S2m |来控制计算的精度来控制计算的精度. 这就是下面要介绍的这就是下面要介

23、绍的龙贝格求积龙贝格求积公式公式.例例 用复化用复化Simpson公式计算积分公式计算积分 的近似值，的近似值，并估计误差。（取并估计误差。（取n=5）101dxIx 解：解：n=5，h=(1-0)/n=0.2，节点列为，节点列为0 10 110., ,ixii 则复化则复化Simpson公式为公式为11111121 0 1 11 02 1 04 1 06 1 080264444441 01 1 03 1 05 1 07 1 09.I 截断误差估计：截断误差估计：44450 145454242412100 1241 3 10180180( )( ) , ( ),max( )()()( / )(

24、).xfxmfxxba hmR 1 052 083333 071429 062500 055560033334 090909 076923 066667 058824 052632(. )( .).( .) 0 03333 20 79450 69315. 理查德森理查德森外推法外推法 /* Richardsons extrapolation */利用利用低低阶公式产生阶公式产生高高精度的结果。精度的结果。设对于某一设对于某一 h 0，有公式，有公式 T0(h) 近似计算某一未知值近似计算某一未知值 I。由。由Taylor展开得到：展开得到： T0(h) I = 1 h + 2 h2 + 3 h

25、3 + i 与与 h 无关无关现将现将 h 对分，得：对分，得：( () )( () )( () ).)(3232222120 hhhhIT Q：如何将公式精度由如何将公式精度由 O(h) 提高到提高到 O(h2) ？.432112)()(23322020 hhIhTTh 即：即：.12)()(2)(32210201 hhIhTThTh .)(42312 hhIhT 12)()(221212 hTTh.)(2211 mmmhhIhT 12)()(2121 mmhmmhTT求积公式 (1)当求积系数 、求积节点 都可以自由选取时,其代数精确度最高可以达到多少次? 1( )( ) ( )()( ,

26、)nbkkakI fx f x dxA f xRf1nkkA1nkkx4.5 高斯型积分高斯型积分 /* Gaussian Quadrature */bankkkxfAdxxfx0)()()(构造具有构造具有2n+1次代数精度的求积公式次代数精度的求积公式将节点将节点 x0 xn 以及系数以及系数 A0 An 都作为待定系数。都作为待定系数。令令 f (x) = 1, x, x2, , x2n+1 代入可求解，得到的公式代入可求解，得到的公式具有具有2n+1 次代数精度。这样的节点称为次代数精度。这样的节点称为Gauss 点点，公式称为公式称为Gauss 型求积公式型求积公式。例：例：求求 的

27、的 2 点点 Gauss 公式。公式。dxxfx)(10解：解：设设 ，应有，应有 3 次代数精度。次代数精度。 101100)()()(xfAxfAdxxfx代入代入 f (x) = 1, x, x2, x3 31130092211200721100521032xAxAxAxAxAxAAA2776. 03891. 02899. 08212. 01010 AAxx不是线性方程组，不是线性方程组，不易求解。不易求解。正交多项式正交多项式证明：证明： “” 对任意次数对任意次数不大于不大于n 的多项式的多项式 Pm(x)， Pm(x) w(x)的次数的次数不大于不大于2n+1，则代入公式应则代入公

28、式应精确成立精确成立： nkkkmkbamxwxPAdxxwxPx0)()()()()( 0= 0 x0 xn 为为 Gauss 点点, 则公式则公式 bankkkxfAdxxfx0)()()( 至少有至少有2 2n n+1+1次代数精度。次代数精度。与任意次数不大于与任意次数不大于n 的多项式的多项式 P(x) （带权）正交（带权）正交。 nkkxxxw0)()(P119;Th5 x0 xn 为为 Gauss 点点 0)()()(1dxxxwxPnba“” 要证明要证明 x0 xn 为为 Gauss 点，即要证公式对点，即要证公式对任意次数任意次数不大于不大于2n+1 的多项式的多项式 Pm

29、(x) 精确成立，精确成立，即证明：即证明： nkkmkbamxPAdxxPx0)()()( 设设)()()()(xrxqxwxPm bababamdxxrxdxxqxwxdxxPx)()()()()()()( 0 nkkkxrA0)( nkkmkxPA0)( n n次多项式次多项式)( 0 )()()()(kkkkkmxrxrxqxwxPTH5TH5表明，在表明，在a,ba,b上带权的上带权的n+1n+1次正交多项式次正交多项式的零点就是求积公式的高斯点。的零点就是求积公式的高斯点。 如何通过正交多项式求高斯点？如何通过正交多项式求高斯点？& 正交多项式族正交多项式族 0, 1, , n,

30、有性质：有性质：任意次数不大于任意次数不大于n 的多项式的多项式 P(x) 必与必与 n+1 正交。正交。若取若取 w(x) 为其中的为其中的 n+1则则 n+1的根的根就是就是 Gauss 点。点。再解上例：再解上例： 101100)()()(xfAxfAdxxfxStep 1：构造正交多项式构造正交多项式 2设设cbxxxaxxx 2210)(,)(, 1)( 53 a0)(10 dxaxx0),(10 1021102100)(53(0),(0)(0),(dxcbxxxxdxcbxxx 215910 cb即：即：215910)(22 xxx Step 2：求求 2 = 0 的的 2 个根，即为个根，即为 Gauss 点点 x0 ，x1221/20)9/10(9/1021;0 xStep 3：代入代入 f (x) = 1, x