数值分析--第4章数值积分与数值微分[1]详解

1、第4章数值积分与数值微分1数值积分的根本概念实际问题当中常常需要计算定积分。在微积分中，我们熟知，牛顿一莱布尼兹公式是计算定积分的一种有效工具，在理论和实际计算上有很大作用。对定积分bf (x)dx，假设f (x)在区间aa,b上连续，且f (x)的原函数为F(x)，那么可计算定积分ba f(x)dx F(b) F(a)似乎问题已经解决，其实不然。如1) f (x)是由测量或数值计算给出数据表时，Newton-Leibnitz公式无法应用。2)许多形式上很简单的函数，例如、.3 sin x .221f (x). 1 x ,sin x ,cos x , exIn xX2等等，它们的原函数不能用初

2、等函数的有限形式表示。3)即使有些被积函数的原函数能通过初等函数的有限形式表示，但应用牛顿一莱布尼兹公式计算，仍涉及大量的数值计算，还不如应用数值积分的方法来得方便，既节省工作量，又满足精度的要求。例如以下积分-rdx-ln密1arc tg(V2x1) arc tg(V2x1) C1 x4.2 x2 2x12、2对于上述这些情况，都要求建立定积分的近似计算方法一一数值积分法。1.1数值求积分的根本思想根据以上所述，数值求积公式应该防止用原函数表示，而由被积函数的值决定。由积分中值定理：对f (x) Ca,b，存在a,b，有bf (x)dx (b a) f ()说明，定积分所表示的曲边梯形的面积

3、等于底为b a而高为f()的矩形面积(图4-1)。问题在于点 的具体位置一般是不知道的，因而难以准确算出 f ()。我们将f()称为区间a,b上的平均高度。这样，只要对平均高度f ()提供一种算法，相应地便获得一种数值求积分方法。如果我们用两端的算术平均作为平均高度f()的近似值，这样导出的求积公式b aTf(a) f(b)2(4-1)便是我们所熟悉的 梯形公式(图4-2)。而如果改用区间中点c节的“高度近似地取代平均高度f ()，那么可导出所谓 中矩形公式(简称矩形公式)(4-2)R (b a)f更一般地，我们可以在区间a,b上适中选取某些节点 Xk，然后用f(xj加权平均得到平均高度f()

4、的近似值，这样构造出的求积公式具有以下形式:图4-1图4-2式中xk称为求积节点；f (x)dxAk f (Xk)k 0(4-3)Ak成为求积系数，亦称伴随节点xk的权。权 九仅仅与节点xk的选取有关，而不依赖于被积函数f (x)的具体形式。这类由积分区间上的某些点上处的函数值的线性组合作为定积分的近似值的求积公式通常称为机械求积公式，它防止了 Newton-Leibnitz公式寻求原函数的困难。对于求积公式(4-3)，关键在于确定节点xk和相应的系数 A 。1.2代数精度的概念由Weierstrass定理可知，对闭区间上任意的连续函数，都可用多项式一致逼近。一般说来，多项式的次数越高，逼近程

5、度越好。这样，如果求积公式对m阶多项式精确成立，那么求积公式的误差仅来源于m阶多项式对连续函数的逼近误差。因此自然有如下的定义定义4.1如果某个求积公式对于次数不超过m的多项式均准确地成立，但对于 m 1次多项式就不准确成立，那么称该求积公式具有m次代数精度。例1判断求积公式11f (x)dx25f ( 一06) 8f (0)95f(丽的代数精度。解记I(f)11 f (x)dx,%f)丄5 f ( . 06)98f(0)5f( 丽因为I(1)1dx12, l%1)- (598 5)2I(x)1xdx=0，1l%x)丄5 0.6 890 5(、0.6) 02 1I(x2)1313I(x )1x

6、I(x4)12 2 2 1 2 x2dx=, %x2)(5 0.6 8 0 5 0.6)393dx=0, I%x3)丄5 (.丽)39x4dx=?, l%x4)i(5 0.36159(、06)30.36)-5l(x5)l(x6)11 x5dx=0，%x5)-5 ( 06)519：x6dx=7，%x6) 95 (0.6)3 0(、06)53(0.6) 0.24所以求积公式具有 5次代数精度。1.3插值型的求积公式最直接自然的一种想法是用原函数是容易求出的，我们以n(x)代替f (x)，由于代数多项式的 n(x)在a,b上的积分值作为所求积分I(f)的近似值，即bI (f) c n(x)dxf (

7、x)在a,b上的插值多项式这样得到的求积分公式称为插值型求积公式。通常采用Lagrange插值。设a,b上有n 1个互异节点x0,X!丄,xn， f (x)的n次Lagrange插值多项式为Ln(X)nlk(x)f (Xk)k 0其中ik(x) 二，插值型求积公式为j 0 XkXjj kbI(f) aLn(x)dxnAk f (xjk 0(4-4)其中Ak被积函数blk(x)dx, k 0,1 丄，n。可看出，Akaf (x)的形式无关。仅由积分区间a,b与插值节点求积公式(4-4)的截断误差为定义Rn(f )bf(x)dxbLn(x)dxa(n1)()a (n 1)!n l(x)dx4.2求

8、积公式xk确定，与(4-5)bf(x)dx anAkf(xQk 0如其系数bAklk(x)dx，a那么称此求积公式为插值型求积公式。定理证明4.1形如(4-3)的求积公式至少有n次代数精度的充分必要条件是插值型的。如果求积公式(4-3)是插值型的，由公式(4-5)可知，对于次数不超过n的多项式f (x)，其余项R f等于零，因而这时求积公式至少具有n次代数精度。反之，如果求积公式(4-3)至少具有n次代数精度，那么对于插值基函数 |k(x)应准确成立，并注 意到lk(xj )jk，即有lk(x)dxAj lk ( xj )j 0所以求积公式(4-3)是插值型的。1.4求积公式的收敛性与稳定性定

9、义4.3在求积公式(4-3)中，假设nf(x)dxlimAJ(Xk)n 0 k 0其中h max(xi x 1)，那么称求积公式(4-3)是收敛的。1 i n实际使用任何求积公式时，除截断误差外，还有舍入误差，因此我们必须研究其数值稳定性。 在求积公式(4-3)中，由于计算f(xQ可能产生误差k，实际得到 %，即f(Xk) % k，记nnIn(f)Akf(Xjln(%Ak%k 0k 0如果对任给正数0，只要误差 k充分小就有nIn(f) ln(% A f(Xk)%(4-6)k 0定义4.4对任给0，假设存在0，只要那么称求积公式(4-3)是稳定的。定理4.2假设求积公式(4-3)中系数Ak0

10、(k负，计算可能不稳定。证明对任给0,假设取-，对kb anIn(f)In(%Akf (Xk)它说明求积公式(4-3)计算是稳定的，由此给出:%k 0注意对任何代数精度0的求积公式均有nAkk 0In(1)f(Xk)%(k 0,1,L ,n)就有(4-6)成立,0,1,L , n)，那么此求积公式是稳定的；假设Ak有正有0,1,L ,n 都有f(xQ %，那么有Akk 0b1dx baf (Xk)%可见a 0时，有In(f) In(%nAkk 0(ba)由定义4.4可知求积公式(4-3)是稳定的。假设Ak有正有负时，假设 Ak ( f (Xk)%)0,且f (Xk)，有In(f) ln(%Ak

11、k 0nAkf (Xk)nAkAJ| f (Xk)(b a)作变换xa th，那么有n n (t j)h hdt0 j 0(k j)hj k(1)n kh nk!(n k)! 0n(tj 0 j kj)dt(1厂2 a) n n (t j)dtAk!( nk)!n 0 j 0j k令C(n)( 1)Cki kn n0 (tj)dtk!( nk)! nj 0 j k那么A(back，求积公式4-4可简化为1(f) (bna)ckn)f(Xk)(4-7)它说明初始数据的误差可能会引起计算结果误差的增大，即计算可能不稳定。2 Newton-Cotes 公式2.1 Cotes 系数被积函数在积分区间内

12、变化平缓,b a步长h，等距节点xk an可用等距节点插值公式近似。将积分区间a,b划分为n等分,kh,k 0,1丄，n。此时求积公式4-4中的积分系数可得到简化bAalk(x)dxaddxxj xxk ok b aokndxk 0称为n阶Newton-Cotes公式，简记为 N-C公式，C：n称为Cotes系数。由Ckn的表达式可看出，它不但与被积函数无关，而且与积分区间也无关。因此可将Cotes系数事先列成表格供查用见表4-1。N-C公式的截断误差为R(f)b f (n 1)()n(x x0、亠 hn 2nf(n (0n)(t j)dtj 0(4-8)a (n1)! jj(n1)!n 1时

13、11baI(f) (b a)2 f (a)-f(b)2 f(a)f(b)(4-9)为梯形公式n 2时14a b1baa bI(f)(b a)-f(a)66f(26f(b)6f (a)4f( 2 ) f(b)(4-10)为辛普生公式。n 4时b aaba3(b a).I(f)7f(a)32 f (a)12f()32 f(a)7f(b)(4-11)90424为Cotes公式。表4-1n111一22o14126661331388887162164904515451925252552889614414441993468403528010577513577132329891728017280172801

14、7280o989588892810496828350283502835028350(n) k2519962889941280358402989132335777511728017280172801728045401049692858889892835028350283502835028350790从表4-1可看出，当n 8时出现了负系数，由定理4.2可知，实际计算中将使舍入误差增大,并且往往难以估计。从而Newton-Cotes公式的收敛性和稳定性得不到保证，因此实际计算中不用高阶 Newton-Cotes 公式。2.2偶阶求积公式的代数精度作为插值型的求积公式,n阶的牛顿-柯特斯公式至少具有

15、 n次的代数精度。求积公式的代数精度能否进步提高呢?定理证明由于这里4.3当阶为偶数时，我们只要验证，当f(n1)(x) (n 1)!，从而NC公式(4-7)至少具有n为偶数时，N-C公式对n 1次代数精度。n 1f (x) X 的余项为零。按余项公式(4-8)，引进变换x a th，并注意到XjbR(f)a jh，有R(f)hn2(X(tXj )dxj)dt当n为偶数，那么n为整数，再令2进一步有尺hn2n2n2n(uj 0nj)duhnn n 22n(u j)du2 j n 2因为被积函数为奇函数。2.3几种低阶求积公式的余项梯形求积公式的余项为RtI T2!(Xa)(x b)dx由于(x

16、 a)(x b)在a,b上不变号，利用积分中值定理有Rt f Q (x a)(x b)dx)(b a)3,(a,b)2 a 12Simpson公式的余项为bb aRS I S f (x)dxf (a)4f(c) f (b)a 6a b这里c。构造次数不超过 3的多项式H (x)，使满足2H(a) f (a),H(c) f (c), H (c) f (c), H (b) f (b)由于Simps on公式具有三次代数精度，它对于这样构造的三次式H (x)是准确的，即bb aH (x)dxH(a) 4H(c) H (b)a 6所以bRSf (x) H (x) dxa由第二章的例6,可知Af(x)

17、H(x) -f(4)( )(x a)(x c)2(x b)4!因(x a)(x c)2(x b)在a,b上保号，应用积分中值定理有(4-12)1bRsf ()(x a)( x c)2 (x b)dx4!a4baba180 2(a,b)(4-13)3复化求积公式前面导出的误差估计式说明，用N-C公式计算积分近似值时，步长越小，截断误差越小。但缩小步长等于增加节点数，亦即提高插值多项式的次数，Runge现象说明，这样并不一定能提高精度。理论上已经证明，当 n 时，N-C公式所求得的近似值不一定收敛于积分的准确值，而且随着n的增大，N-C公式是不稳定的。因此，实际中不采用高阶 积函数用分段低次多项式

18、插值，由此导出复化求积公式。3.1复化梯形公式将区间a,b划分为n等分，分点xk a kh,hN-C公式，为提高计算精度，可考虑对被b a,k 0,1,L ,n，在每个区间 风公昇(k 0,1,L , n 1)上采用梯形公式，那么得ba f(X)dXk 0Xk1f (x)dxh n12f (Xk)f(Xk 1)k 0Rn(f )(4-14)称为复化梯形公式Rn(f)h n 1Tn2k 0，其余项为Tnf(Xk)f(Xk i)f (a)n 12f(Xk)k 1f(b)(4-15)由于 f(x) c2a,b,所以存在 (a,b)使于是复化梯形公式余项为n12 komkin 1f从式(4-16)可以

19、看出，余项误差是k)=(b12a)h2n 1f ( k),n k 0(Xk ,Xk 1)k)k)0mkax1f (k)f() nkk)a)h2f ()Rn(f)吟12h2阶，所以当f(x) C2a,b，有blim Tnf (x)dx，naf(x) Ca,b，那么可得收敛性，因为由b a nf (Xk)(4-15)得即复化梯形公式是收敛的。事实上只要n 11 b a .Tnf (Xk)2 n k 0n k 1所以复化梯形公式(4-15)收敛。此外，Tn的求积系数为正，由定理 4.2知复化梯形公式是稳定的。bf (x)dx (na(4-16)3.2复化辛普森公式一 1将区间a,b划分为n等分，在每

20、个区间xk,1上采用辛普森公式，记xk 12xk - h那么得2bn 1I f (x)dxak 0Xk 1f(x)dxXkhn16 k 0f (Xk) 4f (Xk 1 2)f (Xk 1)Rn( f)(4-17)Snh f(a)1f (xk 1 2)0n 12 f(Xk)k 1f(b)(4-18)称为复化辛普森求积公式，其余项由(4-13)得4RnSnh于是当f (x) C4a, b时，与复化梯形公式相似有f(4)( k),0(Xk,Xk 1)尺需f(4)()，(a,b)(4-19)可以看出误差阶是 h4，收敛性是显然的。事实上，只要f(x) Ca,b，那么有bn 1 a4f(Xk 1 2)

21、bn 1af (Xk)n k of (Xk)f (x)dx (n此外，由于Sn中求积系数均为正数，故知复化辛普森求积公式计算稳定。例2根据函数表4-1kXksin xk f(Xk)XkkXksin xk f(Xk)Xk0015580.936155611/80.997397863/40.908851621/40.989615877/80.877192533/80.9767267810.841470941/20.9588510表4-1dx的近似值，并估计误差。用复化梯形公式和复化辛普森公式计算I1 sin x0x解由复化梯形公式I 丄f(0)16f(1)kf(8)0-945691由复化辛普森公式i

22、 16f(0)f(1)f(k)44 2k 14k1f(T) 0-946084与准确值I 0.9460831L比拟，显然用复化为了利用余项公式估计误差，要求f(x)Simps on公式计算精度较高。sin x的高阶导数，由于Xf(x)sin x1o cos(xt)dt所以有f(k)(x)1 dk0 dxkcos(xt)dt1tkcos(xt2k )dtf(k)(x)1t cos(xt)02dtmax0 x 10kdt0由复化梯形误差公式(4-16)得FUf) |I T8h212 maxf (x)丄120.000434由复化辛普森误差公式(4-19)得R4(f)S4118010.271 1065例

23、3假设用复化求积分公式计算积分的近似值，要求计算结果有四位有效数字,xdxn应取多大?解因为当0 x 1时，有0.30.3xdx要求计算结果有四位有效数字，即要求误差不超过410 。又因为0,1f(k)(x)由式(4-16)得h212101即n104，开方得n 40.8。因此假设用复化梯形公式求积分,6假设用复化Simpson公式，由式4n应等于41才能到达精度。1 h180 2(4-19)(4)()h4180 16 180 16104即得n 1.62。故应取n 2。4龙贝格求积公式4.1梯形公式的逐次分半算法如前所述，复化求积公式的截断误差随着步长的缩小而减少，而且如果被积函数的高阶导数容

24、易计算和估计时，由给定的精度可以预先确定步长，不过这样做常常是很困难的，一般不值得推崇。 实际计算时，我们总是从某个步长出发计算近似值，假设精度不够可将步长逐次分半以提高近似值， 直到求得满足精度要求的近似值。设将区间a,b分为n等分，共有n 1个分点，如果将求积区间再二分一次，那么分点增至2n 1个，我们将二分前后两个积分值联系起来加以考虑。注意到每个子区间xk,xk 1经过二分只增加了1一个分点x 1 (Xk Xk1)，用复化梯形公式求得该子区间上的积分值为 k 2 2h f (xk) 2f (xk 12) f (xk 1)4b a注意，这里h代表二分前的步长，将每个子区间上的积分值相加得

25、h n 1 几 。 f(Xk)f (Xkh n 11)f (Xk 1 2)2 k 0T2n2Tnhn2k(4-20)这说明，将步长由h缩小为一时，T2n等于Tn的一半再加新增加节点处的函数值乘以当前步长。2算法4.11 .输入 a, b, f (x),bm 1,h2 对k2 .置3 .置差为5假设1,2丄h f (a) f (b)m 1,2F F f(a (2 k1)h)2T0T TohF输出IT，停机；否那么m 1hm,2h,TT0，转 3。4.2李查逊(Richards on )外推法假设用某种数值方法求量I的近似值，一般地，近似值是步长h的函数，记为4(h)，相应的误I h(h)其中卫1

26、,2丄),0 piP2LI Ii( h)1(A1hPkLh)P1计2hP2 LkhPk是与h无关的常数。假设用L2( h)P2P2 P22 hk( h)PkPPk式(4-22)减去式(4-21)乘以h)P2取满足h(2(以1P，得P1IP1)hP2P1除上式两端，得P1l1(h)11(h)(P33(R)hP3k(Pkh( h)(4-21)h代替(4-21)中的h，那么得L(4-22)R)hpk Lb2hP2 b3hP3bkhPk(4-23)1Pl)(i2,3 丄P1)：(1Pi)仍与h无关。令“1附)1 p由式(4-23)，以l2(h)作为I的近似值，其误差至少为O(hP2)，因此其中b 2(

27、l2(h)h( h)L(h)收敛于I的速度比1心)快。不断重复以上作法，可以得到一个函数序列I (h) Im1( h)Im(h)1以lm(h)近似I，误差为I lm(h) O(hPm)。随着m的增大，收敛速度越来越快，这就是RichardsonPmPm 12,3,L(4-24)外推法。4.3龙贝格求积公式由前面知道，复化梯形公式的截断误差为(Euler-Maclaurin )公式：定理4.4设f (x) C a,b，那么有2I T(h) 1h 其中系数k(k 1,2,L )与h无关。O(h2)。进一步分析，我们有如下欧拉一麦克劳林2h4 L kh2k L把李查逊外推法与欧拉一麦克劳林公式相结合

28、，可以得到求积公式的外推算法。特别地，在外1推算法式(4-24)中，取pk 2k，并记T0(h) T(h)，那么有h4jTm1(2 Tm1(h)(4-25)(4-26)Tm(h)育 ,m 1,2 丄41经过m(m 1,2,L )次加速后，余项便取以下形式：Tm(h) IME)2h2(m 2) L上述处理方法通常称为 李查逊(Richards on )外推加速方法。为研究Romberg求积方法的机器实现，引入记号：以T。表示二分k次后求得的梯形值，且以Trmk)表示序列 T0(k)的m次加速值，那么依以上递推公式得到Tmk)(mtt1 時补 1,2丄4141称为龙贝格求积算法。Romberg公式

29、的计算过程见下表 4-2表4-2khT(k)T0T(k)11T (k)T2T(k)I3T (k)1 4L0b aT0(0)1b aTT()2112b a4T0(2)T1t2(0)3b aT(3)t(2)TT(0)8T0T1T2T34b aT0(4)T1t2(2)T3t4(0)16MMMMMMMO算法4.2(1)输入 a,b, f (x),(3)计算k 12To(k)(h)1To(k 1) h f(a (i ；)h)2i i2对j 1丄,k(k j)1 jj (k j 1) (k j)4 lj 11 j 14j 1 Tk(0) Tk(01，输出 TV例4用Romberg算法计算积分If (x)d

30、x 停机;a13 2.x dx。0否那么-2h, k 1 k，返回。解 利用逐次分半算法(4-20)和Romberg算法(4-25)，计算结果见表 4-3。1T,0)f (0) f (1) 0.500000T0-T0(0)0.5 f(-)0.4267772 2AAQT,2)T。 0.25 f (-) f (-)0.407018244T0(3)1t0(2)吐 f(1) f (-) f (5) f(-) 0.4018122 2 8 8 8 8M表4-3kT0(k)T1(k)T2(k)T3(k)才)T5(k)00.50000010.4267770.40236920.40701830.4018120.

31、4000770.4000540.40005040.4004630.4000090.40000950.4001180.4000020.4000020.4000025高斯求积公式5.1 一般理论等距节点的插值型求积公式，虽然计算简单，使用方便，但是这种节点等距的规定却限制了求 积公式的代数精度。试想如果对节点不加限制，并适中选择求积系数，可能会提高求积公式的精度。Gauss型求积公式的思想也正如此，亦即在节点数n固定时，适当地选取节点xj与求积系数人讣，使求积分公式具有最高精度。设有n 1个互异节点x0,x1丄，Xn的机械求积分公式a (x)f(x)dxAk f(xk)ak0x1 (l 0,1丄，

32、m)，式(1)精确成立，即 bl(x)xldx (l0,1,L ,m)(2-27)具有m次代数精度，那么有取f(x)n Aj(xj)l j0式构成m 1阶的非线性方程组，且具有2n 2个未知数xk, Ak (k 0,1,L , n)，所以当(x)给(2-28)定后，只要 m 1 2n 2 ，即 m 2n 数精度可到达 2n 1 。另一方面，对式意的互异节点 xk k1时，方程组有解。这说明式 n 1 个节点的求积公式的代(1)，不管如何选择xk 与 Ak ，最高精度不可能超过 2n 1 。事实上，对任0 ，令p2n 2 (x)2n1(x) (x x0)2(x x1)2L (x xn)2b0 ，

33、然而an(x)P2n 1(x)dx 0 。有Ak P2n 2 (xk )k0定义 4 如果求积分公式 (4-27) 具有 2n 1 次代数精度，那么称这组节点xk 为 Gauss 点 ，相应公式(4-27)称为带权 (x) 的 高斯求积公式。定理 5 插值型求积公式的节点 a x0x1 Lxn b 是高斯点的充分必要条件是以这些节点为零点的多项式n 1(x)(x x0)(x x1)L (x xn)与任何次数不超过 n的多项式P(x)带权正交，即ba (x)P(x)a证明 必要性。设P(x) Hn，那么P(x)式 (1) 对于 f (x)n 1(x)dx 0(4-29)P(x)nn 1(x) 精

34、确成立，即有ba (x)P(x) n 1(x)dx1(x) H 2n 1 ，因此，如果x0,x1,L ,xn 是高斯点，那么nAkP(xk ) n 1(xk ) k0故式 (4-29)成立。 再证充分性。对于f(x) P(x) n1(x)f (x) H 2n 1 ，用q(x)，其中 P(x),q(x)ba (x)f(x)dxa1(x) 除 f (x) ，记商为H n ，由式 (4-29)可得 b(x)q(x)dxP(x) ，余式为 q(x) ，即(4-30)由于所给求积公式 (4-27) 是插值型的，它对于 q(x)b(x) f(x)dxa再注意到n 1(xQ0 (k 0,1,L , n)，知 q(xQbba (x) f(x)dx a (x)q(x)dxaaH n 是精确成立的，即 nAkq(xk ) k0f(xQ (k