第七章数值微分与数值积分初次修改稿1

1、1第七章第七章 数值微分与数值积分数值微分与数值积分1 数值微分数值微分2 NewtonCotes求积公式求积公式3 复化求积公式复化求积公式4 Romberg求积公式求积公式5 Gauss型求积公式型求积公式2 利用离散点上函数的信息求函数导数近似值利用离散点上函数的信息求函数导数近似值的方法的方法, 称为称为数值微分数值微分. 差商型求积公式差商型求积公式 插值型求积公式插值型求积公式1 数值微分数值微分3由导数定义由导数定义hxfhxfxfh)()(lim)(0 当当h很小时很小时, 可用可用差商差商近似导数近似导数.4 差商型求导公式差商型求导公式 (3)中心差商公式中心差商公式0,)

2、()()( hhxfhxfxf,)()()(hhxfxfxf .2)()()(hhxfhxfxf (1) 向前差商公式向前差商公式 (2) 向后差商公式向后差商公式5 几何意义几何意义hx xhx ABChxfhxfkBC)()( hhxfxfkAB)()( hhxfhxfkAC2)()( B点切线斜率点切线斜率)(xf 从几何直观看从几何直观看: 中心差商效果最好中心差商效果最好6 截断误差截断误差)(2)( )()()( 1hOhhxfhxfhxfxf )(2)( )()()( 2hOhhxfhhxfxfxf )(12)()(2)()()( 223)3(3)3(hOhhxfhxfhhxfh

3、xfxf 其中其中1,0321 由由Taylor公式可得公式可得7 二阶导数的中心差商公式二阶导数的中心差商公式2)()(2)()( hhxfxfhxfxf 截断误差截断误差)(12)()(2)()( )4(22 fhhhxfxfhxfxf 8近似计算近似计算 badxxfI)(数值积分数值积分9 依据微积分基本定理依据微积分基本定理, 只要找到被积函数只要找到被积函数 f (x)的的原函数原函数 F (x), F (x)=f (x), 便有便有)()()(aFbFdxxfba 为什么还要对积分进行近似计算为什么还要对积分进行近似计算 大量的被积函数找不到用初等函数表示的原函数大量的被积函数找

4、不到用初等函数表示的原函数 实验测量或数值计算给出的通常是一张函数表实验测量或数值计算给出的通常是一张函数表, 即被积函数的表达式未知即被积函数的表达式未知.数值积分数值积分10 依据积分中值定理依据积分中值定理 就是说，底为就是说，底为b a 而高为而高为 f ( ) 的矩形面积恰恰等的矩形面积恰恰等于所求曲边梯形于所求曲边梯形 f (x)的面积的面积.)()(abfdxxfba 取取a, b内若干个节点内若干个节点xk 处的高度处的高度 f (xk ), 通过加通过加权平均的方法生成平均高度权平均的方法生成平均高度 f ( ), 这类求积公式称这类求积公式称机械求积公式机械求积公式式中式中

5、 xk 称为称为求积节点求积节点, Ak 称为称为求积系数求积系数, 亦称伴随亦称伴随节点的权节点的权. nkkkbaxfAdxxf0)()(数值积分基本思想数值积分基本思想112 Newton-Cotes 公式公式基本思想基本思想: 利用利用插值多项式插值多项式).()(xfxLn banbadxxLdxxffI.)()()(其中其中Ln(x)是是n阶阶Lagrange插值多项式，用插值多项式，用Ln (x)的的积分近似积分近似 f (x)的积分，即的积分，即插值型求积公式插值型求积公式12 bankjjjkjkdxxxxxA,0)()(由由 决定决定,与与 无关无关.节点节点 f (x)

6、在在a, b上取上取 a x0 x1 0, 使得使得|,|max0knk 则称该求积公式是则称该求积公式是稳定稳定的的. 求积公式的稳定性求积公式的稳定性39 若求积公式是稳定的若求积公式是稳定的, 则则 f (x)的观察值的较小的的观察值的较小的误差引起的求积结果的误差也是较小的误差引起的求积结果的误差也是较小的. 求积公式求积公式没有把没有把 f (x)的误差的误差“放大放大”很多很多.40), 1 , 0()(nkfxfkk nkkkknnfxfAfIfI0)()()(证明证明因此复化梯形公式是数值稳定的因此复化梯形公式是数值稳定的.当当).(2)1(2abhhnh 定理定理 复化梯形公

7、式是复化梯形公式是数值稳定数值稳定的的. nkkkkfxfA0)( nkkA0 41x0 x2xf (x)x4hhxn 2hxnmnnabh2, .hx3x1xn 1 复化复化Simpson公式公式分片二次多项式近似分片二次多项式近似42 将积分区间将积分区间a, b划分为划分为n=2m等分等分, 步长步长 h=( b a )/n, 分点分点 xk= a+kh ( k=0, 1, , n). 在每个在每个小区间小区间 x2k 2 , x 2k ( k=1, , m)上用上用Simpson公式：公式： kkxxkkkkkxfxfxfxxdxxf222)()(4)(6)(21222222 )()(

8、4)(321222kkkxfxfxfh 复化复化Simpson公式公式k=1, , m43 111122)(4)(2)()(3mkmkkkxfxfbfafh= Sn( f ) mkxxbadxxfdxxffIkk1222)()()( )()(4)(3212221kkkmkxfxfxfh 复化复化Simpson公式公式44)()()(fSfIfRnn 当当 f (x)在在a, b上具有四阶连续导数时上具有四阶连续导数时, ),(),()(1)4()4(bamffmkk ),(222kkkxx 故得故得),(),(180)()(90)()4(4)4(5bafhabfmhfRn 复化复化Simpso

9、n公式的截断误差公式的截断误差 mkkfh1)4(5),(90 45 由复化由复化Simpson公式的截断误差知公式的截断误差知, 误差阶为误差阶为 h4, 收敛性是显然的收敛性是显然的, 事实上事实上,只要只要 f (x) Ca, b则则可得到可得到收敛性收敛性, 即即.)()(limdxxffSbann 由于求积系数均为正由于求积系数均为正, 与复化梯形公式一样的与复化梯形公式一样的证法可得复化证法可得复化 Simpson公式是公式是数值稳定数值稳定的的.46例：例：计算计算dxx 10214 解：解： )1()(2)0(161718fxffTkk8kxk 其中其中= 3.13898849

10、4 )1()(2)(4)0(241oddeven8fxfxffSkk8kxk 其中其中= 3.141592502运算量运算量基本相同基本相同 显然用复化显然用复化Simpson公式计算精度较高公式计算精度较高, 这与它们这与它们的误差阶的结论是相符的的误差阶的结论是相符的.47例例 对于函数对于函数,sin)(xxxf 给出给出n=8的函数表的函数表, 试用试用复化梯形公式及复化复化梯形公式及复化Simpson公式计算积分公式计算积分 10.sindxxxI解解ix)(ixf0.00.1250.250.3750.50.6250.750.8751.01.00.99739780.98961580.

11、97672670.95885100.93615560.90885160.84147090.8771925应用复化梯形公式求得应用复化梯形公式求得T8=0.9456909应用复化应用复化Simpson公式求得公式求得S8=0.9460832准确值准确值 I=0.9460831两者运算量基本相同两者运算量基本相同48trapz: 复化梯形公式求积分复化梯形公式求积分.用法用法: trapz(X, Y), 其中其中X, Y为相同维数的向量为相同维数的向量.例例: X=0.125:0.125:1.0;Y=sin(X)./X;X=0,X;Y=1,Y;trapz(X,Y)ans = 0.945690863

12、58270Matlab函数函数49例例 若用复化求积公式计算积分若用复化求积公式计算积分dxeIx 10的近似值的近似值, 若要求计算结果有若要求计算结果有4位有效数字位有效数字, n应取多大应取多大?解解, 1110 dxeIex.105 . 04 1 , 0, 1| )(|)( xexfxk复化梯形公式的误差复化梯形公式的误差)( )(12|2 fabhRT .83.40 n若用复化梯形公式求积分若用复化梯形公式求积分, n取取41能达到精度要求能达到精度要求.2121n 41021 50故应取故应取n=4. 该例表明该例表明, 为达到相同的精度为达到相同的精度, 用复化用复化Simpso

13、n公式所需的计算量比复化梯形公式要少公式所需的计算量比复化梯形公式要少, 这也说明这也说明了复化了复化Simpson公式的精度高公式的精度高.复化复化Simpson公式公式的误差的误差)()(180|)4(4 fabhRS .25. 3 n41801n 41021 51 复化复化梯形公式的逐次分半算法梯形公式的逐次分半算法将区间将区间a, b分成分成n=2m等分等分, 记记,2mmabh .)(2)()(21212 mkmmkhafbfafhTm), 2 , 1 , 0( m称称 为为梯形值序列梯形值序列.2mT52 mhTTmm2122所有新增加节点的函数值之和所有新增加节点的函数值之和其中

14、其中.2mmabh 复化复化梯形公式的逐次分半算法梯形公式的逐次分半算法53以以n=8, m=3为例为例. 记记 fk= f (xk)x0 x2x4x6x3x1x5x7x8 )(2167654321808fffffffffabT )(21664280fffffab 75318ffffab 24T 所有新增加节点的函数值之和所有新增加节点的函数值之和. 3h54 复化梯形公式余项的后验估计复化梯形公式余项的后验估计);,(),(12112bafhabTIn );,(),(2122222bafhabTIn f ( 1 ), f ( 2 ) 分别是分别是 f (x) 在在a, b上的上的n个点与个点

15、与 2n 个点处的算术平均值个点处的算术平均值 (每个小区间上取一个点每个小区间上取一个点). 当当n较大时较大时, 有有.)( 1)()(21dxxfabffba 55因此因此, 若事先给定误差限若事先给定误差限 , 则当则当.3|2 nnTT时时, 就可停止计算就可停止计算, 并认为并认为 T2n是满足精度要求的近是满足精度要求的近似值似值.;412 nnTITI);()(21 ff ).(3122nnnTTTI 56 复化复化Simpson公式的逐次分半算法公式的逐次分半算法将区间将区间a, b分成分成 n=2m 等分等分, 记记,2mmabh evenodd2)(2)(4)()(3km

16、kmmkhafkhafbfafhSm, 2 , 1 m称称 为为Simpson序列序列.2mS57;1612 nnSISI);,(),(18011)4(4bafhabSIn );,(),(218022)4(42bafhabSIn );()(2)4(1)4( ff ).(15122nnnSSSI 因此因此, 若事先给定误差限若事先给定误差限 , 则当则当.15|2 nnSS时可停止计算时可停止计算, 取取 S2n为满足精度要求的近似值为满足精度要求的近似值. 复化复化Simpson公式余项的后验估计公式余项的后验估计58nnnnnTTTTTI3134)(31222 4 Romberg求积公式求积

17、公式启示启示: 是否用是否用 复化梯形公式余项的后验估计表明复化梯形公式余项的后验估计表明nnTT31342 逼近逼近 I ( f ) 比用比用 T2n要好要好. 事实上有事实上有.313422nnnSTT 即梯形值序列的巧妙线性组合得到即梯形值序列的巧妙线性组合得到Simpson序列序列!59以以n=4为例加以说明为例加以说明. 记记 fk= f (xk), )(2167654321808fffffffffabT )(28642804fffffabT 8abkaxk )(422847654321808fffffffffabT .3)4(488TTS )( 2)( 48642753180fff

18、ffffffab 484TT 60).(15122nnnSSSI .1516)(151222nnnnnSSSSSI 逼近逼近 I ( f ) 比用比用 S2n要好要好.回答回答: 是的是的, 记记,15)16(22nnnSSC 15)16(2nnSS 则它恰为复化则它恰为复化Cotes公式公式; 且有如下误差估计式且有如下误差估计式).(| )(|62hOfICn 复化复化Simpson公式余项的后验估计表明公式余项的后验估计表明 问题问题: 是否用是否用61.631636422nnnRCC 类似地可以得到类似地可以得到2nR其中其中被称为被称为Romberg序列序列.).(| )(|82hO

19、fIRn 截断误差截断误差:621T2T4T8T16T32T2S4S8S16S32S4C8C16C32C8R16R32R3422nnnTTS 151622nnnSSC 636422nnnCCR 停机准则停机准则:梯形值序列梯形值序列Simpson序列序列Cotes序列序列Romberg序列序列 Romberg求积公式求积公式63例例 计算计算.sin10dxxxI 2)1()0(1ffT =0.9207355)5 . 0(212112fTT =0.9397933 )43()41(412124ffTT=0.9445135=0.9456909 )87()85()83()81(812148ffffT

20、T解解先求梯形值序列先求梯形值序列64nnTnSnCnR24180.92073550.93979330.94451350.94569090.94614590.95608690.94608330.94608300.94608310.94608313422nnnTTS 151622nnnSSC .636422nnnCCR 利用只有两三位有效数字利用只有两三位有效数字的的T1, ,T8 经过三次外推得经过三次外推得到到7位有效数字位有效数字. 可见加速的可见加速的效果十分显著效果十分显著.用用Romberg算法计算如下算法计算如下65 理论依据理论依据: 复化梯形公式的余项展开复化梯形公式的余项展开

21、.记记),(hTTn 定理定理 设设,)(baCxf 则则 kkhhhIhT24221)( 其中系数其中系数 k ( k=0, 1, )是与是与 h 无关的常数无关的常数. T (h) 逼近逼近 I 的速度是的速度是 O ( h2 )阶阶.66,3)()2(4)(262411 kkhhhIhThThT 当区间当区间a, b 2n等分时等分时, 则有则有),2(2hTTn 在定理中以在定理中以 h/2 代替代替 h 得得,2164224221 kkhhhIhT 上式乘以上式乘以4减去减去 T(h) 再除以再除以3, 记之为记之为 T1(h), 得得T1 (h) 逼近逼近 I 的速度是的速度是 O

22、 ( h4 )阶阶, 效果比效果比 T (h)好好, 它它不是别的不是别的, 就是就是Simpson序列序列.67 类似地类似地 kkhhhIhThThT2826111215)()2(16)( ,162262411 kkhhhIhT 上式乘以上式乘以16减去减去 T1(h) 再除以再除以15, 记之为记之为 T2(h), 得得T2 (h) 逼近逼近 I 的速度是的速度是 O ( h6 )阶阶, 效果比效果比 T1 (h)好好, 它不是别的它不是别的, 就是就是Cotes公式序列公式序列.68 对对Cotes公式序列进行同样处理得到公式序列进行同样处理得到Romberg公式序列公式序列.Rich

23、ardson外推加速方法外推加速方法也称为也称为Romberg求积算法求积算法 收敛性说明收敛性说明: 如果如果 f (x) 充分光滑充分光滑, 那么梯形公那么梯形公式序列式序列, Simpson公式序列公式序列, Cotes公式序列公式序列, Romberg公式序列均收敛到所求的积分值公式序列均收敛到所求的积分值. 对于对于 f (x)不充分光滑的函数也可用不充分光滑的函数也可用Romberg算算法计算法计算, 只是收敛慢一些只是收敛慢一些. 也可以直接使用复化也可以直接使用复化Simpson公式计算公式计算.69例例 用用Bomberg算法计算积分算法计算积分.1023dxxI 解解23)

24、(xxf 在在0, 1上仅是一次连续可微上仅是一次连续可微用用Romberg算法计算结果见下表算法计算结果见下表nnTnSnCnR241816320.50.4267770.4070180.4018120.4004630.4001180.4023690.4004320.4000770.4000140.4000020.4003020.4000540.4000090.4000020.4000500.4000090.40000270. )()(1 nkkkbaxfAdxxf5 Gauss型求积公式型求积公式 基本思想基本思想 设计求积公式设计求积公式:在节点数在节点数 n 固定时固定时, 适当地选取求

25、积节点适当地选取求积节点 xk 与求与求积系数积系数 Ak , 使求积公式具有使求积公式具有最高的代数精确度最高的代数精确度.71例例 确定确定x1, x2, A1, A2, 使求积公式使求积公式)()()(221111xfAxfAdxxf 具有最高次的代数精确度具有最高次的代数精确度.x2x11 选取选取 (A1 , A2 , x1 , x2)使该求积公式对使该求积公式对 f (x) = 1, x, x2, x3 时等号成立时等号成立. 172)()()(221111xfAxfAdxxf 313111 0 32 0 21 12121322311113322221111222211112111

26、xxAAxAxAdxxxfxAxAdxxxfxAxAxdxxfAAdxf)31()31()(11ffdxxfI 对对 f = 1, x, x2, x3 积分精确成立积分精确成立 四个方程四个未知数四个方程四个未知数73梯形公式与梯形公式与Gauss求积公式的比较求积公式的比较 对对1, x求积公式精确求积公式精确成立成立(1(1次代数精确度次代数精确度) ) 对对1, x, x2, x3求积公式精求积公式精确成立确成立(3 次代数精确度次代数精确度)74)()()()( :333221111xfAxfAxfAdxxfn x3x11x2 1 选取选取(A1, A2 , A3 , x1, x2 ,

27、 x3) 使该求积公式对使该求积公式对 f (x) = x0, x1, x2, x3, x4, x5 时等号成立时等号成立.区间区间 1, 1上的上的Gauss求积公式求积公式75 3/503/55/98/95/9 0 52 0 32 0 21 132132153352112511554334211241144333321123113323322112211223321121131121xxxAAAxAxAxAdxxxfxAxAxAdxxxfxAxAxAdxxxfxAxAxAdxxxfxAxAxAxdxxfAAAdxf)53(95)0(98)53(95)(11fffdxxfI 对对 f = x

28、0, x1, x2, x3, x4, x5 求积公式等号成立求积公式等号成立76syms A1 A2 A3 x1 x2 x3eq1=A1+A2+A3-2eq2=A1*x1+A2*x2+A3*x3eq3=A1*x12+A2*x22+A3*x32-2/3eq4=A1*x13+A2*x23+A3*x33eq5=A1*x14+A2*x24+A3*x34-2/5eq6=A1*x15+A2*x25+A3*x35A1,A2,A3,x1,x2,x3=solve(eq1,eq2,eq3,eq4,eq5,eq6)Matlab求解求解77)()()()(221111nnxfAxfAxfAdxxf 选取选取(A1,

29、A2, , An , x1, x2 , , xn)使该求积公式对使该求积公式对f (x) = x0, x1, x2, , x2n 1时等号成立时等号成立. .区间区间 1, 1上的上的Gauss求积公式求积公式78 0 121 32 0 21 112122112121112122222211222112222222112211222112111121 nnnnnnnnnnnnnnnnnnnxAxAxAdxxxfxAxAxAndxxxfxAxAxAdxxxfxAxAxAxdxxfAAAdxf 2n个方程个方程2n个未知数个未知数, 非线性方程非线性方程, 解是否存在唯解是否存在唯一？若存在一？若

30、存在, 如何求解？如何求解？ 对对 f = x0, x1, x2, , x2n 1求积公式等号成立求积公式等号成立79. )()()(1 nkkkbaxfAdxxfx 研究最一般情形的带权积分研究最一般情形的带权积分:具有具有2n1次代数精确度次代数精确度, 则称这组节点则称这组节点 xk为为 Gauss 点点, 上述公式称为带权函数上述公式称为带权函数 (x)的的Gauss型型求求积公式积公式.定义定义 如果一组节点如果一组节点x1, x2, , xn a, b能使求积公式能使求积公式 nkkkbaxfAdxxfx1)()()( Gauss型型求积公式求积公式80. 12, 1 , 0,)(

31、1 nmxAdxxxmknkkbam 由上述由上述 2n个方程确定全部个方程确定全部 的的2n个待定参数个待定参数 xk , Ak ( k=1, , n), 使求积公式至少具有使求积公式至少具有 2n 1次代数精次代数精确度确度. 但上述方程组是非线性方程组但上述方程组是非线性方程组, 求解十分困求解十分困难难. 一般利用一般利用 正交多项式正交多项式来求出来求出Gauss点与求积系点与求积系数数. Gauss型型求积公式求积公式 对对 f = x0, x1, x2, , x2n 1求积公式等号成立求积公式等号成立81定理定理 一求积公式的节点一求积公式的节点 x1, x2, , xn是高斯点

32、的充是高斯点的充分必要条件是以这些节点为根分必要条件是以这些节点为根(零点零点)的多项式的多项式与任何次数不超过与任何次数不超过n1的多项式的多项式P(x)带权正交带权正交)()()(21nnxxxxxxx bandxxxPx. 0)()()( 高斯求积基本定理高斯求积基本定理82 bajdxxlx)()(02 高斯求积公式的稳定性定理高斯求积公式的稳定性定理定理定理 高斯求积公式的系数高斯求积公式的系数Ak (k=1, 2, , n)全是正全是正的，因而高斯求积公式是稳定的的，因而高斯求积公式是稳定的.证明证明 对于拉格朗日插值基函数对于拉格朗日插值基函数 lj(x) (j=1, , n),

33、 有有 nkkjkxlA12)(.jA 83.)()()(1lim bankkkndxxfxxfA 高斯求积公式的收敛性定理高斯求积公式的收敛性定理定理定理 设设 f (x) C a, b, 则高斯求积公式是收敛的则高斯求积公式是收敛的.84高斯求积公式中高斯点和求积系数的确定高斯求积公式中高斯点和求积系数的确定确定高斯求积公式中的两组不同性质的待定系数确定高斯求积公式中的两组不同性质的待定系数 xk , Ak (k=1, 2, , n)采用不同的方法：采用不同的方法： 由由a, b上带权上带权(x)正交的正交的n次多项式的零点确定次多项式的零点确定高斯点高斯点 xk (k=1, 2, n);

34、 由由n 1次代数精确度条件确定求积系数次代数精确度条件确定求积系数Ak (k=1, , n).85在在 1,1上权函数为上权函数为1的积分的积分 的高斯求积的高斯求积公式的求积节点公式的求积节点(高斯点高斯点)即为即为n次勒让德多次式次勒让德多次式Pn (x)的零点的零点. 11)(dxxf).0(2)(11fdxxf 高斯高斯-勒让德求积公式勒让德求积公式 n=1时，时，P1(x)=x, 高斯点高斯点x1=0, 求积系数求积系数A1=2, 一点高斯一点高斯-勒让德求积公式为中矩形公式勒让德求积公式为中矩形公式86).31)(31(21)13(21)(22 xxxxP 11.3131)(ffdxxf n=2两点高斯两点高斯-勒让德求积公式为：勒让德求积公式为：高斯高斯-勒让德求积公式勒让德求积公式高斯点高斯点31,3121 xx求积系数求