第四章 数值积分与数值微分 4.2 复合求积法§ 4.2 复合求积法§ 固定时而节点个数的长度较大当积分区间 1,],[ +nba 直接使用 Newton-Cotes公式的余项将会较大 增加时即而如果增加节点个数 1, +n 公式的舍入误差又很难得到控制 为了提高公式的精度 ,又使算法简单易行 ,往往使用复合方法 分成若干个子区间即将积分区间 ],[ ba 然后在每个小区间上使用低阶 Newton-Cotes公式 最后将每个小区间上的积分的近似值相加 一 、 复合求积公式 等份分割为的积分区间将定积分 nbadxxfb a ],[)(∫ nkkhaxk ,,1,0, L=+= n abh ?=各节点为 公式上使用在子区间 CotesNewtonnkxx kk ??=+ )1,,1,0](,[ 1 L 节点为步长为等份分割为将 ,,],[ 1 lhlxx kk + 1,, 2,, +=+++ kkkkk xl lhx l hx l hxx L 记为 121 ,,,, ++++ = k l lk lklk k xxxxx L )(1 )( k l x x Idxxfk k∫ + ≈ ∑ = + + ?= l i l ik l ikk xfCxx 0 )( 1 )()( 求积公式阶的上作在 CotesNewtonlxfxx kk ?+ )(],[ 1 ∑ = + = l i l ik l i xfCh 0 )( )( ∫ba dxxf )( ∑∫ ? = += 1 0 1 )( n k x x k k dxxf ∑? = ≈ 1 0 )( n k k lI 由积分的区间可加性 ,可得 ∑∑? = = + = 1 0 0 )( )( n k l i l ik l i xfCh复合求积公式 nI= 可得复合梯形求积公式时 ,1=l n b a Tdxxf ≈∫ )( ∑∑? = = += 1 0 1 0 )1( )( n k i iki xfCh ∑? = ++= 1 0 1 )]()([2 1n k kk xfxfh ])()(2)([2 1 1 ∑? = ++?= n k kn bfxfafn abT复合梯形公式 求积公式可得复合时 Simpsonl ,2= ∑∑? = = + = 1 0 2 0 2 )2( )( n k i iki xfChn b a Sdxxf ≈∫ )( 复合 Simpson公式 复合抛物线公式 n b a Sdxxf∫ ≈)( ∑ ? = ++ ++= 1 0 1 2 1 )]()(4)([6 1n k kkk xfxfxfh )]()(2)(4)([6 1 1 1 0 2 1 bfxfxfafn ab n k k n k k +++?= ∑∑ ? = ? = + 求积公式可得复合时 Cotesl ,4= ∑∑? = = + = 1 0 4 0 4 )4( )( n k i iki xfChn b a Cdxxf ≈∫ )( )](7)(32)(12)(32)(7[90 1 1 0 4 3 4 2 4 1 + ? = +++ ++++= ∑ k n k kkk k xfxfxfxfxf h )](7)(14)](32)(12)(32[)(7[90 1 1 1 0 4 3 4 2 4 1 bfxfxfxfxfafn ab n k k n k kkk +++++?= ∑∑ ? = ? = +++ 复合 Cotes公式 )]()([2 )( bfafabT +?= ],[ 10 xx )]()([2 1)0( xfafhT += ],[ 21 xx 2 )1( hT = )]()([ 21 xfxf + ],[ 32 xx 2 )2( hT = )]()([ 32 xfxf + ],[ 43 xx 2)3( hT = )]()([ 43 xfxf + ],[ 12 ?? nn xx 2)2( hT n =? )]()([ 12 ?? + nn xfxf ],[ 1 nn xx ? 2 )1( hT n =? )]()([ 1 bfxf n +? 2 hT n = +)([ af +)(2 1xf +)(2 2xf +L +? )(2 1nxf )](bf 复合梯形公式分解 )]()2(4)([6 bfbafafabS +++?= ],[ 10 xx 6)0( hS = )]()(4)([ 1 2 10 xfxfaf ++ + ],[ 21 xx 6)1( hS = )]()(4)([ 2 2 111 xfxfxf ++ + ],[ 1 nn xx ? 6)1( hS n =? )]()(4)([ 2 111 bfxfxf nn ++ +?? 6 hS n = +)([ af +∑ ? = + 1 0 2 1 )(4 n k k xf +∑ ? = 1 1 )(2 n k kxf )](bf 复合 Simpson公式分解 例 1. ∫= 10 sin dxx xI计算定积分使用各种复合求积公式 解 : 为简单起见 ,依次使用 8阶复合梯形公式 、 4阶 复合 Simpson公式和 2阶复合 Cotes公式 可得各节点的值如右表 0 1 0.125 0.99739787 0.25 0.98961584 0.375 0.97672674 0.5 0.95885108 0.625 0.93615564 0.75 0.90885168 0.875 0.87719257 1 0.84147098 )( ii xfx 8 7 6 5 4 3 2 1 0 x x x x x x x x x Trapz 4 2 13 3 2 12 2 2 11 1 2 10 0 . x x x x x x x x x Simp + + + + 2 4 31 2 11 4 11 1 4 30 2 10 4 10 0 x x x x x x x x x Cotes + + + + + + 复合求积 公式的程序 newtoncotes.m 函数程序 func.m 8T ])1()(2)0([16 1 7 1 ∑ = ++= k k fxff 分别由复合 Trapz、 Simpson、 Cotes公式有 94569086.0= 4S )]1()(2)(4)0([24 1 3 1 3 0 2 1 fxfxff k k k k +++= ∑∑ == + 94608331.0= 2C )]1(7)(14)](32)(12)(32[)0(7[1801 1 1 1 0 4 3 4 2 4 1 fxfxfxfxff k k k kkk +++++= ∑∑ == +++ 94608307.0= 8T 94569086.0= 4S 94608331.0= 2C 94608307.0= 原积分的精确值为 ∫= 1 0 sin dx x xI L671839460830703.0= 精度最高 精度次高 精度最低 比较三个 公式的结果 那么哪个复合求积公式的收敛最快呢 ？ 二 、 复合求积公式的余项和收敛的阶 我们知道 ,三个求积公式的余项分别为 )(TR )(12 )( 3 hfab ′′??= )(SR )() 2(180 )4(4 hfabab ???= )(CR )() 4(945 )(2 )6(6 hfabab ???= 单纯的求积公式 复合求积公式的每个小区间 )(12 2 kfhh h′′???= )(2180 )4( 4 kf hh h? ? ?? ? ??= )(49452 )6( 6 kf hh h? ? ?? ? ??= ])(12[ 1 0 3∑? = ′′?= n k kf h h ],,[)(.1 2 baCxf ∈设被积函数 则复合梯形公式的余项为 nTI ? ∑ ? = ′′?= 1 0 3 )(12 n k kf h h )(max)()}({min 1 0 xfnfxf bxa n k k bxa ′′≤′′≤′′ ≤≤ ? =≤≤ ∑ h由于 使得由介值定理 ],,[, ba∈?h )()( 1 0 hh fnf n k k ′′=′′∑ ? = nTI ? ∑? = ′′?= 1 0 3 )( 12 n k k n fnh h )( 12 3 hfnh ′′?=即有 )( )(12 )( 2 nTR fhab = ′′??= h nTI ? 2h TI n? ∑? = ′′?= 1 0 3 )(12 n k kf h h又由 ∑? = ′′?= 1 0 )(121 n k k hf h ∑ ? = ?′′?= 1 0 )(121 n k kk xf h ∫ ′′?→ ba dxxf )(121 )]()([121 afbf ′?′?= ???????? ∞→→nh 0 复合梯形公式的余项为足够大时因此当 ,n nTI ? )]()([12 2 afbfh ′?′?≈ )]()([121 2 afbfh ′?′?= ∑? = ??? ? ??? ? ??= 1 0 )4( 4 5 )(2180 n k kf h h ],[)(.2 4 baCxf ∈若被积函数 nSI ? 公式的余项为复合足够大时则 Simpsonn , )(2180 )4( 4 hfhab ????????= )]()([2180 4 4 afbfh ′′′?′′′??≈ ∑? = ??? ? ??? ? ??= 1 0 )6( 6 7 )(49452 n k kf h h 公式的余项同样可得复合若 CotesbaCxf ],,[)(.3 6∈ nCI ? )(4945 )(2 )6(6 hfhab ? ? ?? ? ???= )]()([49452 )5()5(6 6 afbfh ???≈ )]()([49452 )5()5( 6 afbfh ????????= )]()([21801 4 afbfh ′′′?′′′???????= nTI ? )]()([12 1 2 afbfh ′?′?≈ nSI ? )]()([2180 1 4 afbfh ′′′?′′′? ? ?? ? ??≈ nCI ? )]()([4945 2 )5()5(6 afbfh ?? ? ?? ? ??≈ 比较三种复合公式的的余项 的速度依次更快趋于定积分即 ICST nnn ,, 阶无穷小量，，的分别是 642h 为此介绍收敛阶的概念 )( 2ho= )( 4ho= )( 6ho= 定义 1. 满足使其余项及若存在对于复合求积公式 nn IIcpI ?≠> ,00 ch II p n h =? →0 lim 阶收敛的是则称复合求积公式 pIn )( pn hII o=? 阶收敛的概念也等价于显然 p, 不难知道 ,复合梯形 、 Simpson、 Cotes公式的收敛阶分别为 2阶 、 4阶和 6阶 通常情况下 ,定积分的结果只要满足所要求的精度即可 精度越高越大分割的小区间数而积分区间 nInba ,],[ 运算量也很大太大但 ,n 但精度可能又达不到运算量虽较小太小 ,,n 取多大值合理呢 ？那么 n 三 、 复合求积公式步长的自动选取 复合梯形公式的余项为 nTI ? )]()([12 2 afbfh ′?′?≈ 个时小区间数量增加到即将步长缩小一倍 nhh 2,21, 1 = nTI 2? )]()([12 2 1 afbfh ′?′?≈ )]()([) 2(12 1 2 afbfh ′?′?= n n TI TI ? ? 2 4 1≈因此有 nn TITI ?≈? 244 )(31 22 nnn TTTI ?≈?即 为的近似值的截断误差约作为因此 IT n2 )(31 22 nnn TTTI ?=?≈? 为的近似值的截断误差约作为 IS n2 )(151 22 nnn SSSI ?=?≈? 为的近似值的截断误差约作为 IC n2 )(631 22 nnn CCCI ?=?≈? 依此类推 e若预先给定的误差限为 有因此对一般的复合积分 nI )(1 22 nnn IIpII ?=?≈? epII nn <?2只要 e<? nII 2就有 的近似值即为满足要求的 II n2 步长自动选取的步骤 : 11 ,1.1 Iabhn 计算，取 ?== 取不同的值 不同的方法 p ||1,212, 12212 IIpIhhn ?=?== 和计算，取步长折半 否则停止计算若 ,,, 2II =<? e ||1,214,.2 24424 IIpIhhn ?=?== 和计算，取步长折半 否则停止计算若 ,,, 4II =<? e 依此类推 kII 2,, =<? 停止计算直到 e 以上这种方法称为 自适应求积法 有时也去掉 精度会更高 nS )]()(2)(4)([6 1 1 1 0 2 1 bfxfxfafn ab n k k n k k +++?= ∑∑ ? = ? = + 以复合 Simpson求积公式的特点为例 具有以下特点 : 1)()( 的系数总是bfaf + 2)( 的和的系数总是kxf 4)( 2 1 的和的系数总是+kxf 旧节点 新节点 步长折半 =0s =2s =1s 01),()(0,/)(,1,1 =+=?=== sbfafsnabhkn 2 , 2 10 10 hax bxax += == + )()(0,/)(,2,2 bfafsnabhkn +=?=== ,23,2 , 2 11 2 10 2 101 hathat xhat +=+= =+= ++ + 0s→ )1(2s→ )1(211 sss +?→ )2(2s→ ))1(24120(6 ssshI ++= )]2(24120[6 ssshI ++= ,1, +?+? kknnn )()(0,/)(,3,4 bfafsnabhkn +=?=== 2 13 2 12 2 11 2 10 2 11312 2 101 ,,, ,, ++++ ++ === yyyy tytyty )2(211 sss +?→ )3(2s→ )]3(24120[6 ssshI ++= ,1, +?+? kknnn 四 、 复合自适应求积法的算法设计 规定步长的求积法的算法比较简单 ,这里只以自动选取步 长的复合 Simpson求积公式为例介绍自适应法的算法设计 (一 ) 算法名称 ∫= ba dxxfI )(求定积分 ),,,( EPSbafunnautosimpso (二 ) 存储方式 aa 存积分下限 bb 存积分上限 e存预先给定的误差限EPS )(1 afy 存函数值 )(2 bfy 存函数值 存前一次积分值1I 存后一次积分值2I 一维向量 ， 存分段节点x 数值一维向量 ， 存节点处函f 存旧节点处函数值之和1S 和存新增节点处函数值之2S (三 ) 自然语言 e误差限输入积分上下限 ,,.1 ba 01,01,1,1.2 ==== Iskn )(2),(1.3 bfyafy == 210.4 yys += nabh /)(.5 ?= 0)(2.6 =ks nj ,,2,1.7 L= 2)12().1(7 hjax ?+= )()(2)(2).2(7 xfksks += )1(211,1.8 ?+=> kssskif 6)](24120[2.9 hksssI ++= 15/)12(.10 IIabsd ?= EPSdif >.11 否则,5,,21,1 gotonnnIIkk +==+= hI ,2.12 输出 停机.13 (四 ) 程序实例 程序 名 ： Autosimpson.m 命令 格式 ： autosimpson('fun',a,b,eps) 例 2. 用自适应 Simpson公式计算下列定积分 ,并比较 ∫ ?= 10 2 dxeI x Li42.m 解 : m = n=8 trapz n=20 trapz n=50 trapz n=100 trapz n=4 simpson n=20 simpson n=4 cotes autosimp 1e-4 autosimp 1e-6 autosimp 1e-10 z = 0.745865614845695 0.746670836939873 0.746799607189351 0.74681800146797 0.746826120527467 0.746824136005348 0.746824133229615 0.746826120527467 0.746824140606985 0.74682413281433 t = 0 0 0 0.06 0 0 0 0.11 0.16 0.22 d = -9.6e-04 -1.5e-04 -2.5e-05 -6.1e-06 2.0e-06 3.2e-09 4.2e-10 2.0e-06 7.8e-09 1.9e-12 h =0.25 4 h =0.0625 16 h =0.0078125 128 z1= 0.746824132812427 区间数及方法 积分结果 运行时间 误差