数值微分和数值积分

1、第第2章章 数值微分和数值积分数值微分和数值积分2.1 数值微分2.2 数值积分2.3 复化数值积分2.4 重积分计算2.5 高斯型积分公式2.6 程序示例2.1 数值微分 导数与差商的关系 数值微分取成导数的近似值，即差商。000()( ) ( )()( ) (2.1)()() limlimlimhhhf xhf xhf xf xhfxhf xhf xhh2.1.1 差商型数值微分 数值微分用向前差商表示 截断误差：用Taylor展式， 存在x0,x0+h， 数值微分用向后差商表示000()()()f xhf xfxh0 0 +xxh200000( )()()()(),2!ff xhf xf

2、xxhxh000()()( )( )()( ).2f xhf xfR xfxhO hh000()()()f xf xhfxh00 xhx 同理得截断误差 数值微分用中心差商表示 截断误差：由Taylor展式，000()()( )( )()( )2f xf xhfR xfxhO hh000()()()2f xhf xhfxh000 xhxxh23010002302000()()()()() (1)2!3!()()()()() (2)2!3!fxff xhf xfx hhhfxff xhf xfx hhh2001202200()()()()( )()262 ( )(), .(mean theore

3、m)6f xhf xhffhR xfxhhfO hxhxh 数值微分的几何意义000()()()f xhf xfxh000()()limhf xhf xh 设定最佳步长 计算数值微分时产生的误差有截断误差和舍入误差两部分组成。由数值微分公式产生截断误差，有原始数据产生舍入误差。 一般情况，步长 h 越小，误差也越小，但步长太小，会引起误差的增长。因此，实际计算时，需要选择一个最佳步长。 用中心差商分析数值微分的误差。 截断误差： 223max( ).66hhfM f(x0-h)f(x0+h)有效数字严重损失 舍入误差： 设 e 为计算 f(x0-h) 和 f(x0+h) 时的最大舍入误差，可以

4、证明，舍入误差量不超过 e/h。于是得到计算数值微分的总误差为实际计算中，采用事后估计方法来选去步长h/2：若23( )6heg hMhh 越小，舍入误差越大，并产生病态323opt3From ( )0, we know that min( ).33So, the optimal is =. Meg hhg hhehhM( )( /2).D hD hD(h)表示数值微分计算公式例2.2 对函数 y = ex,选取不同的步长进行计算 f (1.15),观察 误差的变化规律，确定最佳步长。解 用中心差商表示的数值微分计算公式得到： h f (1.15) error h f (1.15) error

5、 0.1 3.1630 -0.0048 0.05 3.1590 -0.0008 0.09 3.1622 -0.0040 0.04 3.1588 -0.0006 0.08 3.1613 -0.0031 0.03 3.1583 -0.0001 0.07 3.1607 -0.0025 0.02 3.1575 0.0007 0.06 3.1600 -0.0018 0.01 3.1550 0.0032实际计算时，准确误差难以得到，|D(0.05)-D(0.04)|0, 证明用梯形积分公式计算积分 所得到的数值计算结果比准确值大，并说明其几何意义。 证明 由公式的截断误差知，几何意义：因为 f (x)0，

6、所以 f(x)为凹函数，因此梯形面积 大于曲边梯形的面积，如图。( )dbaf xx3() ( )( )d( )( )0( )0)12 ( )d( ).bababaE ff xxT fffxf xxT f 例2.6 试确定求积公式中的系数 A, B, C 时期代数精度尽量高，并确定代数精度。解 1111( )d(0),22f xxAfBfCf201211, 0, . Taking ( )1, ,we have22xxxf xx x 1111114 d2 31120d0 22311240d 4433AABCxABCx xBABCxC 1141241 ( )d(0). (1)32332f xxff

7、f3421 Take ( ), leftright0; ( ), left,right,56 (1) does not hold. 3.f xxf xxm例 2.7 计算中矩形求积公式的误差：解 ( )( )d().2baabI ff xxfba2( )( ),2222!2fababababf xffxx22( )d()d2221 ( )d221 ()( )d ,222bbaababaabababf xxfbafxxabfxxababfbafxx2231 ( )d()( )d222( )( ) d() .2224bbaabaababf xxfbafxxffabxxba Newton-Cotes

8、积分系数 a, b n 等分, 节点 xi = a +ih (i = 0,1,2,n), h = (b - a)/n, 作变换，令 x = a +th，则 0tn 。于是，li(x)的分子分母：于是我们得到( )d .biial xx02101112()()()(1)() (1)(),()()()()()() (1)2() ( 2 )( () ) nniiiiiiiiinxxxxxxth thtn ht ttn hxxxxxxxxxxxxih ihh hhhni h ( 1)! ()!,n ininih 02011100( )d()()() d()()()()()(1)() d()( 1)!

9、()!( 1)(1)() d! ()!()( 1) ()! (biiabnaiiiiiiinnnn inn inn il xxxxxxxxxxxxxxxxxxxt ttn hh tti hinihht ttntinitiban i0( )(1)()d)!() ()nnit ttntnitiba c ci(n)称为Newton-Cotes系数，它与积分节点和积分区间无直接关系，只与插值的节点数有关。表 2.1 Newton-Cotes系数(n =1,2,8)( )( )( )( )( )( )( )( )( )012345678 111 221412 666133 88nnnnnnnnnnccc

10、cccccc31 8871621674 90451545901925252525195 288961441449628841993496 84035280105280941 3584098958889281049645401049692858889898 283502835028350283502835028350283502835028350( )0( )00 1.().nniinnniiiicbacba Newton-Cotes积分公式的代数精度 由截断误差知， Newton-Cotes积分公式至少有 n 阶代数精度。 可以证明：当 n 为偶数时，Newton-Cotes公式的代数精度可达

11、 n +1，即对 n+1次多项式，求积公式也精确成立。 2.3 复化数值积分 求积公式的收敛性和稳定性 考虑求积公式： 定义1 若 则称求积公式(*)是收敛的。 0()( )d ,0limnbiianibaf xf xxhn0( )d() (*)nbiiaif xxf x 计算 f(xi)时可能产生误差i而实际得到的是 定义2 如果对任意的0，存在0,只要 就有 则称求积公式(*)是稳定的。, i.e. ().iiiiff xf (),iif xf000()( (),nnniiiiiiiiiif xff xf定理 如果N.-C.求积系数i0(i=0,1,2,n), 则求积公式是 稳定的。证明

12、如果对任意的0，若取(存在)= /(b - a), 都有 (), 0,1,2,iif xfin00000then , by the definition of stability, we have ()() , (nnniiiiiiiiiinniiiif xff xfb). a因为i0由表2.1知，当n 7时，Newton-Cotes积分公式是不稳定的。因此不可能用高阶N.-C.积分公式来提高求积精度。 复化求积法 通常把积分区间等分成若干个子区间，在每个子区间上用低阶的求积公式（如梯形积分公式、Simpson积分公式），对所有的子区间求和即得整个区间a, b上的积分公式，这种方法称为复化求积

13、法。2.3.1 复化梯形积分 将a, b分成若干小区间，在每个区间xi, xi+1上用梯形积分公式，再将这些小区间上的数值积分累加起来，就得到区间a, b上的数值积分。这种方法称为复化梯形积分。 计算公式 a, b n等分, h = xi+1- xi= (b -a)/n, xi = a + ih, i = 0,1,2,n,11031111012131100( )( )d( )d() ()()()2121111 ( )()()( )()222212iinbxaxiniiiiiiiinnniiiiiiI ff xxf xxxxxxf xf xfhhf af xf xf bf复化梯形公式，记为 T(

14、h) 或 Tn( f ):截断误差：1111( )( )( )()( ) (2.11)22nniiT hTfhf af xf b133032( )( )( )()( )1212() ( ), , (2.12)12nnniihhEfI fTffnfbafa bn 10()()niifnfFrom (2.12), it is easy to see that lim( )( )( )d .bnanTfI ff xx2332322Let max |( )|, then 0, when ()() 1, ( ).1212a x bnMfxbaMbanEfMn Romberg2.3.2 复化Simpson

15、积分 计算公式 a, b 2m 等分, m 为积分子区间数，记 n = 2m，n+1为节点总数 ，h = xi+1- xi= (b -a)/n, xi = a + ih, i = 0,1,2,n,2221051(4)221220112120151(4)0( )( )d( )d(2 )2 ()4 ()()()62880 ( )4()2()( )3(2 ) ()2880iimbxaximiiiiimmiiiimiiI ff xxf xxhhf xf xf xfhf af xf xf bhf复化Simpson积分为1121201( )( )4()2()( ) (2.13)3mmniiiihSff a

16、f xf xf b系数首尾为1，奇数点为4，偶数点为2 截断误差显然，51(4)05(4)5(4)45(4)44(4)(2 )( )( )( )()2880(2 ) ( ) ( , )2880() ( ) (,2 )2880() ( )180 ( ), 180mnniihEfI fSffhmfa bbabafhnmnmbafnbah f , (2.14)a blim( )( )( )d .bnanEfI ff xx记 5(4)444()max( ), then ( ) .2880na x bbaMfxEfMm 544544 In order to let ( ), only need to l

17、et () , 2880()i.e. when 1, we have ( ).2880nnEfbaMmbaMmEf例2.8 计算 使计算结果有5位有效数字。 分别用复化梯形和复化辛普森求积公式，求分点数。 解由复化梯形误差公式得由复化辛普森误差公式得10( )e d ,xI fx4(4)(4)110 , ( )e , ( )( )e ,2 |( )| |( )|ee, 01.xxxf xfxfxfxfxx34222()1( )( )10 ,2121267.3, take 68.nbaeI fTfMnnnn441 ( )( )10 22880 2.084 13, 26.neI fSfmmnm 2

18、.3.3 自动控制误差算法 在实际计算中，复化积分公式的误差通常用事后估计的方法，如用 |T2n( f )-Tn( f )| 来估计误差 |I( f )-T2n( f )| ，即 先计算 T2n( f ) 将a, b n等分, 记 hn = (b - a)/n, 由 (2.11) 式,对a, b 2n等分，即对xi, xi+1 在等分一次，记其中点为22( )( )( )( ) (Later we know 1/3)nnnI fTfC TfTfC1111( )( )()( ) .22nnniiTfhf af xf b1/2(21),222nnniinhhhxxaihai,inxaih11/ 2

19、11/ 21100 ( )( )d( )d( )d( )d iiiiiinnbxxxaxxxiiI ff xxf xxf xxf xx11/221/2011/21/21111/2101 ( )()()2 ()()211 ( )()()( ) 2221 ( )()22niiniiiiiiinnniiiinniixxTff xf xxxf xf xhf af xf xf bhf af x111/2011/2011( )() 221 ( ) ( ) (2.15)2where ( )()(21) .2nniinnnnnnniniihf bf xTfHfhHfhf xhf aiRomberg1/22ni

20、ihxx12111For example, from ( ), we know when 1, ( )( ( )( );211 when 2, ( )( )()( )2221 ();22nTfbanTff af bbanTff af xf bbaTf x1( )Tf2( )Tf 误差估计 由复化梯形公式截断误差公式(2.12)222 ( )( )( ), 12( )( )( ),122nnbaI fTfh fbahI fTff 1210011 ( )(), ( )(), i.e. ( ),( ) 2are mean values on nodes and 2 nodes respectivel

21、y, thus we can regard them as ( )( ) .Therefore we have nniiiiffffffnnnnff222 ( )( )4( )( ) ,1 ( )( )( )( ) (2.17)3nnnnnI fTfI fTfI fTfTfTf22Thus , 0,provided ( )( )3 , we have ( )( ).nnnTfTfI fTf 注注 1 1：由(2.17), 复化梯形求积公式 T2n( f ) 的误差大致为 我们将这个误差值看成公式 T2n( f ) 的一种补偿，可以得到注注 2 2： 类似地，辛普森求积公式的误差可用 S2n(

22、f ) 和 Sn( f ) 来表示，即21( )( )3nnTfTf2221( )( )( )( ) 31 = 4( )( ) ( ) ).3nnnnnnI fTfTfTfTfTfSf2211( )( )( )( ) (2.17)15nnnI fSfSfSf2.3.4 龙贝格(Romberg)积分 理查森（Richardson）外推法 我们已知: (2.11) , (2.15) 对a, b的 n 等分区间 xi, xi+1 应用辛普森公式，有11211/20 ( )( )2()( ) (2.11)21 ( )( ) ( ) ) (2.15)2where ( )() .nniinnnnnniih

23、Tff af xf bTfTfHfHfhf x 由此可见，复化梯形公式经过线性组合就可以得到复化辛普森公式，其结果是截断误差由 O(h2) 提高到 O(h4) 。这个过程称为理查森(Richardson)外推法。11111/2100111/21022( )( )d()4 ()()612 ( )2()( )()3 2312 ( )2( )( )331 4( )( )3iinnxiiniiixiinniiiinnnnnxxSff xxf xf xf xhf af xf bhf xTfTfTfTfTf21( ) 4( )( ) ( ) ).3nnnI fTfTfSf由注1知，From (2.15)例

24、2.9 将Richardson外推法应用于复化Simpson 公式，可 以得到复化柯特斯公式： (截断误差O(h6) 解 由复化辛普森公式的截断误差得 S2n(f ) 和 Sn(f ) 的误差： 对充分大的n, f (4)()和 f (4)()几乎是a, b上 f (4)(x)的均值，因此可视为f (4)() f (4)() 。于是由上述二式得2161( )( )( )( ) (2.19)1515nnnI fSfSfCf4(4)24(4)( )( )( ), 1802 ( )( )( ),180nnbahI fSffbaI fSfh f 22222 ( )( )16( )( )4( )( )1

25、61 ( )( )( )( ). 151541nnnnnnnI fSfI fSfSfSfI fSfSfCf 同理，再对复化柯特斯公式进行组合就可以得到龙贝格(Romberg)公式: (截断误差O(h8) 为便于在计算机上实现龙贝格算法，将Tn, Sn, Cn, Rn,用统一的符号 Rk, j 表示, k, j = 0,1,2, 。 列标 j 表示外推次数： j = 0, 无外推（梯形公式）； j = 1，外推一次（Simpson公式）； j = 2，外推2次（Cotes公式）； j = 3，外推3次（Romberg公式）；2323641( )( )( )63634( )( ) . (1.20)

26、41nnnnnRfCfCfCfCf行标 k 表示区间a, b等分的次数： k a, b等分 小区间个数 0 不等分 20 = 1 1 1次 21 = 2 2 2次 22 = 4 3 3次 23 = 8 于是，龙贝格计算公式：停止计算：0,011,022,041,112,123,14, , , , ,RTRTRTRSRSRS,11,1,1, 1,2,; 1,2,41k jkjk jk jjRRRRkjk,1,10, .k kkkRR Rk,k的计算顺序 (同差商)如下表： 外推方向 区 间 细 分 0,01,0 1,12,0 2,12,23,0 3,13,23,34,0 4,14,24,34,4

27、 nnRRRRRRRRRRRRRRRTS nnCR不同公式的值2005.11.4 2.6 程序示例 分别给出用复化梯形公式和龙贝格积分公式的自动控制误差算法。 复化梯形公式算法描述 11211/2011 ( )( )()( ) (2.11)221 ( )( ) ( ) ) (2.15)2where ( )() ,nniinnnnniiTfhf af xf bTfTfHfbHfhf xh1/2,.2iiahxaih xaihn算法描述： 1. 输入 m, a, b, 定义 f(x); n = m, h = (b -a)/n; 1111/20112. 2( )( )()( ) ;22 12 100

28、; while | 12| 12; (); 2( 1)/2,/2,2 ; end while3. output 2nniiniiTTfhf af xf bTTTTTTHhf xTTHhhnnT void main () int n=m; int i;double T_n,H_n,T1,T2;double h=(b-a)/n;T_n=(f(a)+f(b)/2;for (i=1;in;i+)T_n+=f(a+h*i);T_n*=h;T2=T_n;printf(At the time,T2=%16.10lfn,T2);T1=T2+100; printf(At the time,T1=%16.10lf

29、n,T1);用复化梯形公式 (2.11)计算TndoT1=T2; printf(At the time,T1=%16.10lfn,T1); for (i=0,H_n=0;iepsilon);printf(T=%16.10lfnn,T2);计算T2n龙贝格积分算法 输入数据： a, b, e(误差), m (区间数), 定义 f(x), n =1, h0 = b-a 计算梯形面积： R0,0 = ( f(a) + f(b)*h/2 for k =1 to m,11,1,141k jkjk jk jjRRRR12,01,011(21) /2;/2;kkkkkikkRRhf aihhh Romber

30、g公式： 输出 Rk, k,1,11,1,1,1for 1 to /(41) if exit cyclejk jk jk jkjk kkkjkRRRRRR主要程序 double computeT(double aa,double bb,long int n) int i;double sum,h=(bb-aa)/n;sum=0;for (i=1;in;i+) sum+=f(aa+i*h);sum+=(f(aa)+f(bb)/2; printf(TT=%16.10lfn,h*sum);return (h*sum); /实现复化梯形公式void main() int i;long int n=N_H, m=0;double T2MAXREPT+1;T10= computeT (a, b, n);printf(T=%16.10lfn,T10);n*=2;for (m=1;mMAXREPT; m+) for (i=0;im; i+