第四章 数值积分与数值微分 4.5 数值微分§ 4.5 数值微分§ 先看一个实例 : 已知 20世纪美国人口的统计数据为 (单位 :百万 ) 年份 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 人口 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 试计算美国 20世纪的 (相对 )年增长率 则人口的增长率为时刻的人口为若记 ),(txt )()( tx dtdxtr = 如何求dtdx 一 、 插值型求导公式 在节点处的函数值但知道不一定给出设函数 )(,)( xfxf bxxxa n ≤<<<≤ L10 nkfxf kk ,,1,0,)( L== 插值有则由阶导数存在的如果 Lagrangenxf ,1)( + )()!1( )()()( 1 )1( xnfxLxf n n n + + ++= w x 有关并与 xba ],[∈x ∏ = + ?= n j jn xxx 0 1 )()(w 插值多项式次的为 LagrangenxfxLn )()( --------(1) 对 (1)式两边求导 ,有 )()!1( )()()!1( ])([)()( 1 )1( 1 )1( xnfxnfxLxf n n n n n + + + + ′+++ ′+′=′ wxwx 将很难确定有关与由于 ])([, )1( ′+ xx nfx 可以求出时但是当 )(, kk xfxx ′= )()!1( )()()!1( ])([)()( 1 )1( 1 )1( kn n kn n knk xn fx n fxLxf + + + + ′+++ ′+′=′ wxwx )()!1( )()( 1 )1( kn n kn xn fxL + + ′++′= wx ∏ ≠ = + ?++′= n kj j jk n kn xxn fxL 0 )1( )()!1( )()( x --------(2) nk ,,1,0 L= ∏ ≠ = + ?+= n kj j jk n kn xxn fxE 0 )1( )()!1( )()( x --------(2) --------(3) (2)式称为 插值型求导公式 , (3)式为相应产生的 误差 由于公式 (2)采取的是 n次 Lagrange插值多项式 ,而高次 插值会产生 Runge现象 ,因此实际应用中多采用低次插 值型求导公式 nk ,,1,0 L= )()()( knknk xExLxf +′=′ 二 、 低阶插值型求导公式 时1=n )()()( 11 kkk xExLxf +′=′ 1.两点公式 1,0=k )(1 xL )(1 xL′ )(2 )()( )2( 1 jkk xx fxE ?= x kjk ≠= ,1,0 则若令 ,01 xxh ?= 10 1 0 xx xxf ? ?= 01 0 1 xx xxf ? ?+ 01 1 10 0 11 xxfxxf ?+?= )()()( 01010 xExLxf +′=′ )(1 01 ffh ?= )(2 )2( xfh? )()()( 11111 xExLxf +′=′ )(1 01 ffh ?= )(2 )2( xfh+ --------(4) --------(5) )(1)()( 0110 ffhxfxf ?≈′≈′ (4)(5)式称为带余项的 两点求导公式 即 精度 1阶 )(hoE =由于 2.三点公式 时2=n )()()( 22 kkk xExLxf +′=′ 2,1,0=k ))(( ))(( ))(( ))(( ))(( ))(()( 1202 10 2 2101 20 1 2010 21 02 xxxx xxxxf xxxx xxxxf xxxx xxxxfxL ?? ??+ ?? ??+ ?? ??= ))(( )()( ))(( )()( ))(( )()()( 1202 10 2 2101 20 1 2010 21 02 xxxx xxxxf xxxx xxxxf xxxx xxxxfxL ?? ?+?+ ?? ?+?+ ?? ?+?=′ ∏ ≠ = ?= 2 0 )3( 2 )(!3 )()( kj j jkk xx fxE x ， 则为等距节点 ， 即若 1201210 ,, xxxxhxxx ?=?= 22212002 2 2 2 3)( h hf h hf h hfxL ?+ ? ?+?=′ ))((!3 )()( 2010 )3( 02 xxxx fxE ??= x )43(21 210 fffh ?+?= 22212012 22)( h hf h hhf h hfxL + ? ?+?=′ )( 2 1 20 ffh +?= 22212022 2 32 2)( h hf h hf h hfxL + ?+=′ )34(2 1 210 fffh +?= )(3 )3( 2 xfh= ))((!3 )()( 2101 )3( 12 xxxx fxE ??= x )( 6 )3( 2 xfh?= ))((!3 )()( 1202 )3( 22 xxxx fxE ??= x )( 3 )3( 2 xfh= )43(21 210 fffh ?+?= )(21 20 ffh +?= )34(21 210 fffh +?= )(3 )3( 2 xfh+ )(6 )3( 2 xfh? )(3 )3( 2 xfh+ )( 0xf ′ )( 1xf ′ )( 2xf ′ --------(6) --------(7) --------(8) (6)(7)(8)式称为带余项的 三点求导公式 其中 (7)式又称为中点公式 ,其精度稍高 在分段求导公式中有着重要的地位 精度 2阶 )( 2hoE =由于 )(21)( 201 ffhxf +?=′ 3.五点公式 )(5)316364825(121)( )5( 4 432100 xf hfffff hxf +?+?+?=′ 时4=n )()()( 44 kkk xExLxf +′=′ 4,3,2,1,0=k )(20)618103(121)( )5( 4 432101 xf hfffff hxf ?+?+??=′ )(30)88(121)( )5( 4 43102 xf hffff hxf ??+?=′ )(20)310186(121)( )5( 4 432103 xf hfffff hxf ?++?+?=′ )(5)254836163(121)( )5( 4 432104 xf hfffff hxf ++?+?=′ ? ? ? ? ? ? ? --------(9) 组 (9)称为带余项的 五点求导公式 精度 4阶 综合考虑上述三种公式 ,可知五点公式的精度最高 )()()( 1 kknk xExLxf +′=′由求导公式 nk ,,1,0 L= 4,2,1=n 可发现 ∑ = =′ n j j k jkkn fAhaxL 0 )( )( 1)( nk ,,1,0 L= 4,2,1=n 且 ∑ = n j k jA 0 )( 0= 并且当步长 h越小时 ,误差会越小 但是不是 h越 小公式越好呢 ? )( 4hoE =由于 充分小时当 h 会很接近njf j ,,1,0, L= 有效数字会严重损失情形中会出现接近数相减的 , 0 )(∑ = n j j k j fA 做除数的现象中也会出现小数同时 hfAhaxL n j j k jkkn ∑ = =′ 0 )( )( 1)( 的舍入误差可能会很大因此 ∑ = =′ n j j k jkkn fAhaxL 0 )( )( 1)( 因此实际应用中步长 h不要取得太小 具体情形可参见教材 P167例 1 -3 -2 -1 0 1 2 3 -5 0 5 10 15 20 -3 -2 -1 0 1 2 3 -5 两 点 公 式 两 点 公 式 0 两 点 公 式 两 点 公 式 -3 -2 -1 -5 0 5 两 点 公 式 两 点 公 式 两 点 公 式 两 点 公 式 三 点 中 点 公 式 两点公式和三点公式的比较图 两 点 公 式 实 际 切 线 三 、 低阶插值型求导公式的分段构造 由于高次插值的 Runge现象 ,数值微分一般采用分段低 次插值公式 ,常见的就是分段两点 、 三点和五点公式 1.分段两点求导公式 bxxxa n ≤<<<≤ L10 nkfxf kk ,,1,0,)( L== 对于任取的相邻两点 1,,1,0,, 1 ?=+ nkxx kk L 由两点公式有 )( kxf ′ )(1 1 kk ffh ?≈ + 1,,1,0 ?= nk L )( nxf ′ )(1 1??≈ nn ffh --------(10) ? ? ? ??? ? 称 (10)式为 分段两点公式 2.分段三点求导公式 bxxxa n ≤<<<≤ L10 nkfxf kk ,,1,0,)( L== 对于任取的相邻三点 1,,2,1,,, 11 ?=+? nkxxx kkk L 由 三点公式 ,有 )(21)( 11 +? +?≈′ kkk ffhxf 1,,2,1 ?= nk L )43(21 210 fffh ?+?≈)( 0xf ′ )34(21 12 nnn fffh +?≈ ??)( nxf ′ --------(11) ? ? ? ??? ? 称 (11)式为 分段三点公式 )(1) 2( 1 1 kk kk ff h xxf ?≈+′ + + 用实际中下面的公式很有 )316364825(121)( 432100 fffffhxf ?+?+?≈′ )618103(121)( 432101 fffffhxf +?+??≈′ )88(121)( 2112 ++?? ?+?≈′ kkkkk ffffhxf )310186(121)( 12341 nnnnnn fffffhxf ++?+?≈′ ????? )254836163(121)( 1234 nnnnnn fffffhxf +?+?≈′ ???? ? ? ? ? ? ? ? --------(12) 3.分段五点求导公式 由 五点公式 2,,3,2,,,,, 2112 ?=++?? nkxxxxx kkkkk L 2,,3,2 ?= nk L 例 : 回到实例 (美国人口 ) 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 )(),( trtxt 人口的增长率为时刻的人口为设 )()( tx dtdxtr = 解 : 则 先用精度较高的分段五点公式求出节点处的导数值 求增长率必须先求导数 x'=21.508 13.458 16.158 11.908 12.267 25.000 27.633 23.075 22.692 28.358 r=dx/dt/x =0.0283 0.0146 0.0152 0.0097 0.0093 0.0166 0.0154 0.0113 0.0100 0.0113 1900 1920 1940 1960 1980 2000 0.005 0.01 0.015 0.02 0.025 0.03 1900 1920 1940 1960 1980 2000 0.005 0.01 0.015 0.02 0.025 0.03 将节点处的增长率作 三次样条插值 节点处增长率的值 每年增长率曲线图 年份 增长率 1900 0.0283 1901 0.0255 1902 0.0230 1935 0.0082 1936 0.0081 1937 0.0083 1953 0.0172 1954 0.0172 1979 0.0100 1980 0.0100 1981 0.0109 1989 0.0111 1990 0.0113 Renkou.m daoshu.m Renkou1.m 四 、 样条求导公式 Lagrange插值型求导公式构造比较简单 但由于误差的原因 ， 只能求出节点处的导数 作为节点必须将处的导数在某一点即若求函数 uuxf ,)( 其缺点显而易见 有根据三次样条插值若 ,],,[)( 4 baCxf ∈ )(|)()(|max 4)()( kkk bxa hoxSxf ? ≤≤ =? 2,1,0=k --------(13) )()( )()( xSxf kk ≈ 2,1,0=k 也有其实不必 ],,[)( 4 baCxf ∈ )()( xSxf 的三次样条插值多项式即先作 )(),(,)( xSxSxS ′′′二阶导数在任意一点的一阶甚至再求 )(),()( xfxfxf ′′′数值导数在相应点的一阶或二阶作为 1,,1,0),()(],[ 1 ?=+ nkxSxfxx kkk L的三次样条函数作在区间 kk k kk k yxxh xxhxS 2 13 )( )(2)( +? ?+= 1 2 3 1 )()(2 + + ???+ kk k kk yxx h xxh kk k k mxx h xx 2 12 )( )( +? ?+ 1 2 2 1 )()( + + ??+ kk k k mxx h xx Hermite 插值 kk k k yxxh xxhxS 2 13 )( )(2)( +? ?+= 1 2 3 1 )()(2 + + ???+ kk k yxx h xxh kk k mxx h xx 2 12 )( )( +? ?+ 1 2 2 1 )()( + + ??+ kk k mxx h xx 假设节点是等距的 ， 则 其中 njmxf jj ,,2,1,0,)( L==′ 1?+ = kk k k hh hl 1 1 ? ? += kk k k hh hm )(3 1 1 1 k kk k k kk kk h yy h yyg ?+?= + ? ? ml 2 1= 2 1= )(23 11 ?+ ?= kk yyh 由于 -----(14) 11 4 +? ++ kkk mmm 因此基本方程组 (三对角方程组 )为 )(3 11 ?+ ?= kk yyh 1,,2,1 ?= nk L 000 )( mfxS =′=′ nnn mfxS =′=′ )( 加上第一类边界条件 1,,2,1, ?= njmj L就可以解出 ----(15) 式代入将 )14()1,,2,1( ?= njmj L 得求导然后对 ,x ])23)(()23)([(2)( 11113 ++++ ?????+???=′ kkkkkkkkk yhxxxxxyhxxxxxhxS ])23)(()23)([(1 11112 ++++ ???+???+ kkkkkkkk mxxxxxmxxxxxh 1,,1,0],,[ 1 ?=∈ + nkxxx kk L ----(16) 上式再对 x求导 ,得 ])246()426[(1)( 1112 +++ ??+??=′′ kkkkkkk mxxxmxxxhxS ))(2(6 113 kkkk yyxxxh ??++ ++ ----(17) 从而 )()( xSxf ′≈′ )()( xSxf ′′≈′′ ---------(18) ---------(19) (18)、 (19)式组合 (16)、 (17)式统称为 样条求导公式 样条求导公式的优点 : 可以求非节点处的 1~2阶导数 精度较高 样条求导公式的缺点 : 要求预知边界条件 要解三对角方程组 比较复杂 五 、 样条求导公式的简化 介绍两种常用的简化方法 : 方程组的求解 ， 避免解三对角简化 ),,1,0()1( njmj L= ),,1,0( njmj L=以替代 首先用三点或五点公式求出节点处的导数 这种方法不但避免解三对角方程组 也不必预知边界条件 (2) 采用先求节点处导数 ,然后作三次样条插值 直接代入 (16)(17)式 ),,1,0)(( njxf j L=′求出节点处的导数即先用三点或五点公式 作三次样条插值与然后再利用 ),,1,0()( njxxf jj L=′ )(xf ′值以便求出任一点的导数 其缺点是不好求二阶导数 实验 : 找一个已知函数 ， 分别用样条求导方法和两种 简化样条求导方法求其数值导数 ， 和精确值比 较 ， 画出导数的图象 ， 并指出它们的异同