数值微分和数值积分

1、第第2章章 数值微分和数值积分数值微分和数值积分当函数 f (x)是以离散点给出时，或函数的表达式过于复杂时，其微分较难于求出。这时可用数值微分来近似计算f (x)的微分值 f (x)1 数值微分数值微分在微积分学中,导数 f (x)有以下三种定义:hxfhxfxfh)()()( lim0 hhxfxfh)()(lim0 hhxfhxfh2)()(lim0 差商的极限差商近似v 差商与数值微分差商与数值微分v 最简单的数值微分方法：用差商近似导数hxfhxfxf)()()( 000 向前差商hhxfxfxf)()()( 000 向后差商hhxfhxfxf2)()()( 000 中心差商误差估计

2、向前差商的截断误差由泰勒展式,)( !2)( )()(002000hxxfhxhfxfhxf 向前差商的截断误差为：hxfhxfxfxR)()()( )(000 )()( !2hOfh 向后差商的截断误差由泰勒展式hhxfxfxfxR)()()()(000 )()( !2hOfh ,)( !2)( )()(002000 xhxfhxhfxfhxf 向后差商的截断误差为：hhxfhxfxfxR2)()()()(000 )()( 6)( )( 1222212hOfhffh ,002001xhxhxx 中心差商的截断误差为：中心差商的截断误差由泰勒展式)( !3)( !2)( )()(1302000

3、 fhxfhxhfxfhxf )( !3)( !2)( )()(2302000 fhxfhxhfxfhxf 例：给出下列数据，计算例：给出下列数据，计算 f (0.02), f (0.10), f (0.08)x0.020.040.060.080.10f(x)5.065.075.0655.055.055解：解：f (0.02)(5.07-5.06) / (0.04-0.02) = 0.5 f (0.10)(5.055-5.05) / (0.10-0.08)=0.25 f (0.06)(5.05-5.07) / (0.08-0.04)=-0.5 f (0.08)(f (0.10)- f (0.0

4、6)/(0.10-0.06)=18.75用函数的差商近似导数，会产生截断误差，而且截断误差与| h |（步长）有关，当| h |越小，则截断误差也会小；但是截断误差只是误差的一部分，此外还有舍入误差，而数值微分恰好对舍入误差非常敏感，它有时会随着| h |的缩小而增大。例：对函数例：对函数 y=ex，取不同的步长计算，取不同的步长计算f (1.15), 观观察误差变化规律，察误差变化规律，(中心差商中心差商)-0.00183.16000.06-0.00253.16070.07-0.00313.16130.08-0.00403.16220.09-0.00483.16300.10误差误差f(1.1

5、5)h0.00323.15500.010.00073.15750.02-0.00013.15830.03-0.00063.15880.04-0.00083.15900.05误差误差f(1.15)h设定最佳步长方法：通常用事后估计方法选取步长h，若记D(h), D(h/2)分别为步长h, h/2 的差商计算公式，给定误差限，当| D(h) -D(h/2) | 时，h/2就是合适的步长。或以| D(hi) -D(hj) |最小确定步长。v 插值型数值微分插值型数值微分v 对于f (x),，给定n+1个点构造n次插值函数(x)，用(x) 的各阶导数近似 f (x) 的各阶导数如：给定n+1个点( x

6、i , yi )( i =0, 1, , n) 构造n次Lagrange插值多项式 niiinxfxlxLxf0)()()()( niiinxfxlxLxf0)()( )( )( )( )( )(xLxfxRn 误差估计 niinxxnfdxd0)1()(!)1()( )()(xLxfdxdn 例例2.3：给定：给定( xi , f( xi)( i = 0, 1, 2) ，并且有并且有 x2 - x1 = x1 x0 = h，计算，计算f ( x0), f ( x1), f ( x2)解：构造过点解：构造过点( xi , f( xi)( i = 0, 1, 2) 的二次多项式)()()()()

7、()()()()()(2120210121012002010212xfxxxxxxxxxfxxxxxxxxxfxxxxxxxxxL )(2)()()()(2)()(2210122002212xfhxxxxxfhxxxxxfhxxxxxL )(2)()()()(2)()( 0122022112202xxxxhxfxxxxhxfxxxxhxfxL )3(2)()( )( 20020hhxfxLxf )(2)()2()(2221hhxfhhxf )()(4)(3(21)( 2100 xfxfxfhxf 在微积分中用牛顿-莱布尼兹公式计算连续函数f (x)的定积分)()()(aFbFdxxfba 当函

8、数 f (x)是以离散点给出时，或f (x)的原函数F (x) 不易求出或过于复杂时，不宜用牛顿-莱布尼兹公式计算积分。第第2章章 数值积分数值积分在微积分中，定积分是在微积分中，定积分是imann和极限和极限 ，它是分割小区间长度趋，它是分割小区间长度趋于零时的极限，即于零时的极限，即但在数值积分中只能用有限的和近似上面的极限，通常由离散点但在数值积分中只能用有限的和近似上面的极限，通常由离散点函数值的线性组合形式给出，记函数值的线性组合形式给出，记100( )()limbniixiaf x dxf xx 0()( ),()()bnniiiaI ff x dx Ifa f x( )I f表示

9、精确积分( )nIf表示近似积分ix为积分节点ia为积分系数确定中积分系数的过程就是构造积分公式的过程( )nIfia代数精度代数精度 记a，b上以（xi，f（ xi ），i，n为积分节点的数值积分公式为：0( )( )nniiiIfa f x( )nIf若满足而则称具有m阶代数精度()()()0,0,1,.,kkknnExI xIxkm1()0mnEx( )nIf设确定以下求积公式的代数精度012,3b ahxa xah xb1293( )()()44bafx dxhfxhfx( )nIf解：取分别计算和的值是否相等2()1,.;fxx x( )I f2(1)93(1)(1)443Ibaba

10、II222221( )()2931( )()()( )4432I xbabaIxahbbaI x23322233221()()3931()()()()4433I xbabaIxahbbaI x34433332232331( )()4931( )()()(4635 )44318( )( )I xbab aI xa hba baababbI xI x所以积分公式的代数精度为v 插值型数值积分插值型数值积分 在在a, b上取上取 a x0 x1 xn b，做，做 f 的的 n 次次Lagrange插值插值多项式多项式 ，即得到，即得到 niiinxlxfxL0)()()( babaniiidxxlx

11、fdxxf0)()()(ibanijjxxxxadxjij 0)()()()(0ibainixfdxxl baniiixfadxxf0)()(积分系数是 li(x)在a, b上的积分由由节点节点 决定，决定，与与 f(x) 无关。无关。例：建立0, 2上节点为x0=0, x1=0.5, x2=2的 数值积分公式 20)(dxxf得由解 baiidxxla)(: 20200031)20)(5.00()2)(5.0()(dxxxdxxla 202011916)25.0)(05.0()2)(0()(dxxxdxxla 20202295)5.02)(02()5.0)(0()(dxxxdxxla)2(5

12、)5 .0(16)0(391)(20fffdxxf v 梯形积分和抛物线积分梯形积分和抛物线积分n=1: 把把a, b1等分，由等分，由x0=a, x1=b 构造构造1次插值多项式次插值多项式)()()()()(11001xfxlxfxlxL )()()(101iiibabaxfadxxLdxxf )(21)(00abdxbabxdxxlababa )(21)(11abdxabaxdxxlababa )()(2)(bfafabfT 记记f(x)abf(a)f(b) baiiixfadxxf10)()()(2)(2bfabafab 梯形求积公式梯形求积公式n = 2:把把a, b2等分，由等分，

13、由 x0=a, x1=(a+b)/2, x2=b构造构造2次次Lagrange插值多项式插值多项式L2(x) babadxxLdxxf)()(2dxbfxlfxlafxlbaba)()()()()()(2210 dxbabaabxbaxdxxlababa )(2/)()(2/)()(00)(61ab badxxla)(11)(64ab )()(4)(6)()(2bffafabfSdxxfbaba badxxla)(22)(61ab Simpson求积公式求积公式衡量积分公式近似程度的方法：代数精确度梯形、梯形、Simpson求积公式的误差估计求积公式的误差估计)()(2)()(bfafabfT

14、dxxfba 考察其代数精度。考察其代数精度。解：解：逐次检查公式是否精确成立逐次检查公式是否精确成立f (x) = 1： baabdx111 2 ab=f (x) = x ：=f (x) = x2 ： 222abbadxx 2baab 3233abbadxx 222baab 代数精度代数精度 = 1 考察考察Simpson求积公式代数精确度。求积公式代数精确度。解：解：逐次检查公式是否精确成立逐次检查公式是否精确成立f (x) = 1： baabdx11416 ab=f (x) = x ：=222abdxxba 246bbaaab f (x) = x2 ：3332abdxxba )2(462

15、22bbaaab =f (x) = x3 ：4443abdxxba )2(46333bbaaab =f (x) = x4 ：5554abdxxba )2(46444bbaaab 代数精确度代数精确度 = 3,)(12)(3baabf 梯形求积公式的截断误差梯形求积公式的截断误差dxbxaxfba)(!2)( dxbxaxfba)(2)( ),(ba dxabxbaxfba )(2)(2 )()(2)()(bfafabdxxffRba dxxdxxfbaba)()(1 代数精度为代数精度为1 Simpson求积公式的截断误差求积公式的截断误差(4)2( )( )()()()4!2bafabEfx

16、 xa xxb dx(4)5( )() ,2880fbaab 从梯形或Simpson求积公式的截断误差公式可知：当积分区间越小时，误差越小复化数值积分复化数值积分v复化梯形积分：复化梯形积分：),., 1, 0(,nihiaxnabhi 在每个在每个 上用梯形公式：上用梯形公式：,1 iixx 1)(iixxdxxf把把a, b分成分成n 等分等分) 1,., 0()()(211 nixfxfxxiiii在在a, b上得到上得到 11)()(2)(2niibfxfafh baniiixfxfhdxxf101)()(2)(x0 xn= T n( f ) 11)()(2)(2nibfihafafh

17、 badxxf)(复化梯形复化梯形求积公式求积公式算法过程：算法过程：1. 输入输入 a, b, n2. h= (b-a)/n3. Tn = f (a) + f (b) ; 4. 对对 i= 1, 2, . , n-1,作以下循环作以下循环 Tn=Tn+ 2*f (a+i*h)5. Tn=Tn* h/26. 输出输出 Tn的的值值用用梯梯形形法法计计算算)10(14102 ndxx 103)(12)(niinfhfE 复化梯形公式截断误差：复化梯形公式截断误差：在每个在每个 上梯形公式的截断误差：上梯形公式的截断误差：,1 iixx,12)()(131 iiiixxhffE 102)()(12

18、niifnabh 由中值定理，当由中值定理，当)( )(,102 nffbabaCfnii 有有存存在在,),(12)(2bafhab 由此可看到复化梯形公式的截断误差按照h2或者1/n2的速度下降v 复化复化Simpson积分：积分：),., 1, 0(,nihiaxnabhi 在每个在每个 上用上用Simpson公式：公式：,222 iixx把把a, b分成分成 n 等分等分)()(4)(62)(22122222 iiixxxfxfxfhdxxfii记记 m= n /2，即即n=2m为积分子区间总数为积分子区间总数 bamiiiixfxfxfhdxxf1022122)()(4)(3)(ix

19、212 ix22 ix4444x0 x1x2x3x4 )()(2)(4)(3)(1011212 mimiiibabfxfxfafhdxxf= Sn( f )注：首尾节点的系数为注：首尾节点的系数为1，奇数节点的系数为，奇数节点的系数为4，偶数节，偶数节点的系数为点的系数为24444x0 x1x2x3x4复化复化Simpson求积公式求积公式1. 输入输入 a, b;2. 输入输入n; 如果如果 n为奇数为奇数,重新输入重新输入;3. h= (b-a) /n ;4. Sn = f (a) + f (b) ; 5. 对对 i= 1, 2, . , n-1作以下循环作以下循环 如果如果 i 是偶数是

20、偶数 , 则置则置 Sn=Sn+ 2*f (a+i*h) 否则否则 置置 Sn=Sn+ 4*f (a+i*h)6. Sn= Sn*h/37. 输出输出 Sndxxfba )(用复化用复化Simpson公式求公式求 的近似值的算法的近似值的算法 )()(2)(4)(3)(1011212 mimiiibabfxfxfafhdxxf 10)4(5)(2880)2()(miinfhfE 在每个在每个 上上Simpson公式的截断误差：公式的截断误差：,222 iixx,)(2880)2()(222)4(52 iiiixxfhfE )(2880)2()4(5 mfh )()2()(901)4(5 fnn

21、ab )(180)()4(45 fnab )(180)()4(4 fhab 复化Simpson公式截断误差由此可看到复化Simpson公式的截断误差按照h4或者1/n4的速度下降)(180)()(180)()()()4(45)4(4 fnabfhabfSfIn 若给定精度 ，如何取 n ?实际计算时，总是从某个步长h出发，逐次h分半。直到满足精度v 复化积分的复化积分的自动控制误差自动控制误差算法算法复化梯形积分的复化梯形积分的自动控制误差算法自动控制误差算法nabhn 把把a, b分成分成n 等分等分 11)()(2)(2)(niinnbfxfafhfT222nnhnabh 2. 把把a,

22、b分成分成2n 等分等分 12122)()(2)(2)(niinnbfxfafhfTx0 x1x2x3x4 12122)()(2)(2)(niinnbfxfafhfT 1111222)()(2)(2)(221)(niniiinnbfxfxfafhfT niinnnxfhTfT11222)(21)( ninnnnhiafhTfT1222)12(21)(x0 x1x2x3x4Q: n =1, 2, 22, 23, 逐次分半什么时候逐次分半什么时候终止？终止？3. 复化梯形积分逐次分半复化梯形积分逐次分半终止条件终止条件)(12)()()(2nnfhabfTfI )(212)()()(222nnfh

23、abfTfI )()(41)()(2fTfIfTfInn )()(3)()(22fTfIfTfTnnn )()( |2fTfIn要要 3| )()(|2 fTfTnn只只须须)()()()(42fTfIfTfInn 复化梯形积分复化梯形积分自动控制误差算法描述自动控制误差算法描述:1. 输入输入a, b, 2. n=1; h=(b-a)/n ; Tn=h/2*(f (a)+ f (b)3. Js=0; 4. for i=1 to n Js=Js+f(a+(2i-1)*h/2);5. T2n=Tn/2+h/2*Js6. 若若 | T2n- Tn |h; 2n = n ; T2n =Tn, 转步骤

24、转步骤3 ninnnnhiafhTfT1222)12(21)(复化复化Simpson积分的积分的自动控制误差算法自动控制误差算法nabhn 把把a, b分成分成n 等分等分nabhn22 2. 把把a, b分成分成2n 等分等分 )()(2)(4)(3)(1011212 mimiiinnbfxfxfafhfSx0 xnx2nx1x2 1012)()(miinxffJ记记 112)()(miinxffO记记)(2)(4)()(3)(fOfJbfafhfSnnnn )()(2)(4)(3)(101121222 niniiinnbfxfxfafhfS)()()(2fOfJfOnnn 则则)(2)(4

25、)()(3)(2222fOfJbfafhfSnnnn )()(2)(4)()(3)(222fOfJfJbfafhfSnnnnn )(2)(4)(2)(4)()(3212fJfJfOfJbfafhnnnnn )(2)(43)(21)(222fJfJhfSfSnnnnn x0 xnx2nx1x23. 复化复化Simpson积分逐次分半终止条件积分逐次分半终止条件)(180)()()()4(4nnfhabfSfI )(2180)()()(2)4(42nnfhabfSfI )()(161)()(2fSfIfSfInn )()(15)()(22fSfIfSfSnnn )()( |2fSfIn要要 15|

26、 )()(|2 fSfSnn只只须须复化复化Simpson积分积分自动控制误差算法描述自动控制误差算法描述:1. 输入输入a, b, 2. 置置n=2; h=(b-a)/n; Jn=f(a+b)/2) Sn=h/3*( f (a)+f (b)+4*Jn)3. J2n=0; 4. for i =1 to n J2n = J2n +f (a+(2i-1)*h/2)5. S2n= Sn/2+h/6(4 *J2n- 2 *Jn)6.若若 | S2n-Sn |h; 2*n = n ; J2n =Jn, S2n=Sn , 转步骤转步骤3)(2)(43)(21)(222fJfJhfSfSnnnnn 复化复化

27、Simpson积分的积分的自动控制误差算法自动控制误差算法nabhn 把把a, b分成分成n 等分等分nabhn22 2. 把把a, b分成分成2n 等分等分 )()(2)(4)(3)(1011212 mimiiinnbfxfxfafhfSx0 xnx2nx1x2 1012)()(miinxffJ记记 112)()(miinxffO记记)(2)(4)()(3)(fOfJbfafhfSnnnn )()(2)(4)(3)(101121222 niniiinnbfxfxfafhfS)()()(2fOfJfOnnn 则则)(2)(4)()(3)(2222fOfJbfafhfSnnnn )()(2)(4

28、)()(3)(222fOfJfJbfafhfSnnnnn )(2)(4)(2)(4)()(3212fJfJfOfJbfafhnnnnn )(2)(43)(21)(222fJfJhfSfSnnnnn x0 xnx2nx1x23. 复化复化Simpson积分逐次分半终止条件积分逐次分半终止条件)(180)()()()4(4nnfhabfSfI )(2180)()()(2)4(42nnfhabfSfI )()(161)()(2fSfIfSfInn )()(15)()(22fSfIfSfSnnn )()( |2fSfIn要要 15| )()(|2 fSfSnn只只须须复化复化Simpson积分积分自动

29、控制误差算法描述自动控制误差算法描述:1. 输入输入a, b, 2. 置置n=2; h=(b-a)/n; Jn=f(a+b)/2) Sn=h/3*( f (a)+f (b)+4*Jn)3. J2n=0; 4. for i =1 to n J2n = J2n +f (a+(2i-1)*h/2)5. S2n= Sn/2+h/6(4 *J2n- 2 *Jn)6.若若 | S2n-Sn |h; 2*n = n ; J2n =Jn, S2n=Sn , 转步骤转步骤3)(2)(43)(21)(222fJfJhfSfSnnnnn 2.3.4 龙贝格（龙贝格（Romberg）求积）求积基本思想：基本思想：梯形

30、公式的逐次分半过程进行加速。梯形公式的逐次分半过程进行加速。)(322nnnTITT 3122nnnTTTI 342nnTTI )(1522nnnSISS 而而)(15122nnnSSSI 15162nnSSI Richardson外推公式外推公式),.,2.)(3 , 2(1411, 11,1,kjkRRRRjjkjkjkjk ),.,2.)(3 , 2(14411,11,1,kjkRRRjjkjkjjk R 1,1R 2,1R 2,2R 3,1R 3,2R 3,3R m,1R m,2R m,3R m,m按行进行运算按行进行运算,第一列元素先按梯形公式的逐次第一列元素先按梯形公式的逐次分半法，再向外递推分半法，再向外递推, 从而构造出精度愈来愈高从而构造出精度愈来愈高的积分，这就是的积分，这就是Romberg求积的基本思想求积的基本思想. niinnnxfhTT11222)(21 12122)()(2)(2niinnbfxfafhT或或),.,2.)(3 , 2(1411, 11,1,kjkRRRRjjkjkjkjk 例：用例：用RombergRomberg方法求下列积分：方法求下列积分：dxx 102143)24(21)1()0(21:1 ,1 ffR解解1 .3)1()21(2)0(21211 ,2 fffR13333. 3141