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

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

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

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

4、节点 Xk,然后用f (Xk)加权平均得到平均高度f()的近似值，这样构造出的求积公式具有下列形式:图4-1图4-2f (x)dxAkf(Xk)k 0(4-3)式中Xk称为求积节点;Ak成为求积系数，亦称伴随节点Xk的权。权 人仅仅与节点Xk的选取有关,而不依赖于被积函数f (X)的具体形式。这类由积分区间上的某些点上处的函数值的线性组合作为定积分的近似值的求积公式通常称为机械求积公式，它避免了 Newton-Leibnitz公式寻求原函数的困难。对于求积公式(4-3),关键在于确定节点xk和相应的系数Ak 。1.2代数精度的概念由Weierstrass定理可知，对闭区间上任意的连续函数，都可

5、用多项式一致逼近。一般说来，多 项式的次数越高，逼近程度越好。这样，如果求积公式对 m阶多项式精确成立，那么求积公式的误 差仅来源于m阶多项式对连续函数的逼近误差。因此自然有如下的定义定义4.1如果某个求积公式于次数不超过m的多项式均准确地成立，但对于 m 1次多项式就不准确成立，则称该求积公式具有m次代数精度。例1判断求积公式1 f(x)dx -5f ( .0.6) 8f (0) 5f( .0.6)19的代数精度。解记1I1f(x)dx，I%f) 95f( 0.6) 8f(0) 5f(。6)因为I(1)1 dx 2, I%1) 1(5 8 5) 21 9I(x) 1 xdx=0, I%x)

6、15 0.6 8 0 5 (、0.6) 01 921I(x2)15 ( 1迸)35 0.36)- 55 (、,06)53_(0.6) 0.24最直接自然的一种想法是用f (x)在a,b上的插值多项式n(x)代替f (x),由于代数多项式的原函数是容易求出的，我们以n(x)在a,b上的积分值作为所求积分I (f)的近似值，即bI (f) a n(x)dxa这样得到的求积分公式称为插值型求积公式。通常采用设a,b上有n 1个互异节点x0,x1,L ,xLagrange 插值。f (x)的n次Lagrange插值多项式为Ln(x)lk(x)f(xk)k 0其中lk(x)x x ,插值型求积公式为0

7、xkxikbI(f) Ln(x)dxanAkf (xk)k 0(4-4)其中Akblk(x)dx, k 0,1,L ,n。可看出， Aa仅由积分区间a,b与插值节点xk确定，与被积函数f (x)的形式无关。求积公式(4-4)的截断误差为定义Rn(f )bf(x)dxbLn(x)dx ab f(n 1)()n 1(x)dx a (n 1)!(4-5)4.2求积公式ba f(x)dxnAkf (xk)k 0如其系数bAklk(x)dx,a则称此求积公式为插值型求积公式。22 n/ 212x2dx=3，%x2)§ (5 0.6 8 0 5 0.6)-I (x3)1 x3dx=0, I%x3

8、) 15 ( .06)3 019I (x4)1 x4dx=-, I%x4) 1(5 0.36 015' '9I(x5)1 x5dx=0, %x5) 15 ( 06)5 01 9I (x6)11 x6dx= 7, I%x6) 15 (0.6)3 0 5所以求积公式具有 5次代数精度。1.3插值型的求积公式定理 证明4.1形如(4-3)的求积公式至少有n次代数精度的充分必要条件是插值型的。如果求积公式(4-3)是插值型的，由公式(4-5)可知，对于次数不超过 n的多项式f(x),其余项Rf等于零，因而这时求积公式至少具有n次代数精度。反之，如果求积公式(4-3)至少具有n次代数精度

9、，那么对于插值基函数1k (x)应准确成立，并注意到lk(xj) jk ,即有所以求积公式(4-3)是插值型的。lk(X)dX aAjlk(Xj) Ak1.4求积公式的收敛性与稳定性定义4.3在求积公式(4-3)中，若limn h 0nAkf(Xk)k 0bf(X)dXa(4-3)是收敛的。其中h maX(x X 1),则称求积公式1 i n实际使用任何求积公式时，除截断误差外，还有舍入误差，因此我们必须研究其数值稳定性。在求积公式(4-3)中，由于计算f(Xk)可能产生误差k ,实际得到 的，即f(Xk) f% k,记In(f)Akf (Xk),In(%)Ak%如果对任给正数0,只要误差k

10、0充分小就有In(f) In(的nAk f(Xk)的k 0(4-6)它表明求积公式(4-3)计算是稳定的，由此给出:定义4.4对任给 0 ,若存在 0,只要f(Xk) 的(k 0,1,L ,n)就有(4-6)成立,则称求积公式定理4.2(4-3)是稳定的。若求积公式(4-3)中系数Ak0 (k0,1,L ,n),则此求积公式是稳定的；若Ak有正有负，计算可能不稳定。证明对任给0,若取In(f)InM,对kb anAk f(Xk)k 00,1,L ,n 都有Akk 0注意对任何代数精度0的求积公式均有nAkk 0Inb1dX baf(Xk)的f (Xk)的可见Ak0时，有In(f) In(的nA

11、kk 0(ba)0,且f (Xk)由定义4.4可知求积公式(4-3)是稳定的。若Ak有正有负时，假设 Ak(f (Xk)In(f) In(的Ak k 0nAkf (Xk)Ak|f(Xk)nAk(b a)它表明初始数据的误差可能会引起计算结果误差的增大，即计算可能不稳定。2 Newton-Cotes 公式2.1 Cotes 系数被积函数在积分区间内变化平缓,可用等距节点插值公式近似。将积分区间a,b划分为n等分,步长hba ,等距节点 nXk akh, k0,1,L ,n。此时求积公式（4-4）中的积分系数可得到简化作变换th ,则有则Ak(bblk(x)dx an x xj -dxxk为O-J

12、hdxk (k j)h0(k j)hk!(n k)! 0(t0 kj)dt(1)nk(b a)k!(n k)! n(t0 kj)dt(n)( 1)nkCkk!(nk)! n(t0kj)dta）Ckn）,求积公式（4-4）可简化为I(f) (bna)k 0Ckn)f(Xk)(4-7)称为n阶Newton-Cotes公式，简记为 N-C公式， Ckn）称为Cotes系数。由Ckn）的表达式可看出，它不但与被积函数无关，而且与积分区间也无关。因此可将Cotes 系数事先列成表格供查用（见表4-1)。N-C公式的截断误差为R(f)b f (n 1) ( ) n(n 1)! j(X0Xj)dxhn 2(

13、n1)!n / f(n1)(n)(t j)dtj 0(4-8)I(f)1(b a) - f (a)12 f(b)f(a)f(b)(4-9)为梯形公式I(f)(ba) -f(a) 6f(16 f(b)f(a)a b4f() f(b)(4-10)为辛普生公式。I(f)b 、 90ab a7f(a) 32 f (a 4b)12f( 2a-)32 f(a3(b a)4) 7f(b)(4-12)为Cotes公式。表4-1nC(n)Ck11122214166631331888847162167904515459051925252525192889614414496288641993499418403528

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

15、积公式,n阶的牛顿-柯特斯公式至少具有 n次的代数精度。求积公式的代数精度能否进步提高呢?定理证明4.3当阶为偶数时,我们只要验证，当NC公式(4-7)至少具有n为偶数时，N-C公式对n 1次代数精度。n 1, “ 一 一f (x) x 的余项为零。按余项公式(4-8),由于这里f(n1)(x) (n 1)!,从而bR(f)a(xxj )dx引进变换x a th ,并注意到xjR(f)hn 2(tj)dt当n为偶数，则n为整数，再令2进一步有R(f) hn2n2 n2(u0n2 j)duhnn n 22n(u j)du2 j n 2因为被积函数为奇函数。2.3 几种低阶求积公式的余项梯形求积公

16、式的余项为b f ()RTI T (x a)(x b)dxa 2!由于(x a)(x b)在a,b上不变号，利用积分中值定理有(4-12)Rtf ( ) (x a)(x b)dx f ( ) (b a)3,(a,b)2 a12Simpson公式的余项为_b .ba.,RSI S f(x)dx f(a)4f (c)f(b)a6a b 这里c 。构造次数不超过 3的多项式H (x),使满足2H(a) f(a),H(c) f (c),H (c) f (c), H (b) f (b)由于Simpson公式具有三次代数精度，它对于这样构造的三次式H(x)是准确的，即bb aH (x)dx 一 H(a)

17、4H(c) H (b)a6所以bRsf(x) H (x) dxa由第二章的例6,可知f(x) H(x) ；f(4)( )(x a)(x c)2(x b)4!2因(x a)(x c) (x b)在a,b上保号，应用积分中值定理有1小 boRsf ( ) (x a)( x c) (x b)dx4!b a180f(4)(),(a,b)(4-13)3复化求积公式前面导出的误差估计式表明，用N-C公式计算积分近似值时，步长越小，截断误差越小。但缩小步长等于增加节点数，亦即提高插值多项式的次数，Runge现象表明，这样并不一定能提高精度。理论上已经证明，当 n 时，N-C公式所求得的近似值不一定收敛于积分

18、的准确值，而且随着n的增大，N-C公式是不稳定的。因此，实际中不采用高阶N-C公式，为提高计算精度，可考虑对被积函数用分段低次多项式插值，由此导出复化求积公式。3.1复化梯形公式将区间a,b划分为n等分，分点xka kh,hb a 一,k 0,1,L ,n，在每个区间xk ,xk 1(k 0,1,L ,n 1)上采用梯形公式，则得称为复化梯形公式Rn(f)ba f (x)dxTn2TnXk k 0 kf(Xk)h3 n12 k1f (x)dxf(Xki)k)=(bf(Xk) f(Xki)f(a)a)h212Rn(f)(4-14)由于 f(X) C2a,b,0mkin 1fk)k)所以存在 （a

19、,b）使于是复化梯形公式余项为n 12f(Xk)k 1n 1f ( k)，max0 k n 1f (k)k)f(b)(4-15)(Xk , Xk 1)(4-16)从式（4-16）可以看出，余项误差是即复化梯形公式是收敛的。事实上只要22h2阶，所以当f (x) C2a,b,有blim Tnf (x)dx ,naf(x) Ca,b,则可得收敛性，因为由(4-15)得a)h2f ()n 11 b aTn - f(Xk)2 n k 0b a n rb rf(Xk)f (x)dx (nn k 1a所以复化梯形公式（4-15）收敛。此外，Tn的求积系数为正，由定理 4.2知复化梯形公式是稳定的。3.2复

20、化辛普森公式1将区间a,b划分为n等分，在每个区间xk,xk 1上米用辛普森公式，记 xk 12 xk ；1 h则得bI f (x)dxan 1 x<1f(x)dxXkh n 1hf(Xk) 4f (Xk12) f(Xk 1)Rn(f)6 k 0(4-17)Sn 6 f(a)n 14 f(Xk12)k 0n 12 f(Xk)k 1f(b)(4-18)称为复化辛普森求积公式，其余项由（4-13）得4f(4)( k),(Xk,Xk 1)Rn1Sn180 2于是当f （x） C4a, b时，与复化梯形公式相似有4f(4)()，(a,b)(4-19)可以看出误差阶是 h4,收敛性是显然的。事实上

21、，只要 f (x) Ca,b,则有n 1f(xk12)k 0b a n 1 b af (xk) nf(xk)k 1R(f) |i T83由复化辛普森误差公式(4-19)得f(x)dx (n ) a此外，由于Sn中求积系数均为正数，故知复化辛普森求积公式计算稳定。例2根据函数表4-1表4-1kxksin xk f(xk) xkk00151”80.9973978621/40.9896158733/80.9767267841/20.9588510Xk5 83 47 81sin xk f(xk)0.93615560.90885160.87719250.8414709用复化梯形公式和复化辛普森公式计算I

22、1 sin xsni dx的近似值，并估计误差。解由复化梯形公式1 .I f(0)16kf ( ) 0.9456918由复化辛普森公式1 ,I -f(0)f(1)16f(k)44 2k 14 f () 0.946084k 18与准确值I 0.9460831L比较，显然用复化为了利用余项公式估计误差，要求f(x)Simpson公式计算精度较高。sin x的高阶导数，由于xf(x)sin x1o cos(xt)dt所以有f(k)(x)1 dk0 dxkcos(xt)dt1tk0k - cos(xt )dt于是max f (k)(x)tk cos(xt7)dt1tkdt0由复化梯形误差公式(4-16

23、)得h 一、行 maxf (x)1120.000434R(f)| |IS41180160.271 1065例3若用复化求积分公式计算积分xdx的近似值，要求计算结果有四位有效数字，n应取多大?解因为当0 x 1时，有0.3 e 1 ex 1于是1 x0.3 e xdx 10要求计算结果有四位有效数字，即要求误差不超过-10 4。又因为2f 的(x) ex 1 x 0,1由式(4-16)得工h212h2121040.8。因此若用复化梯形公式求积分，n应等于41才能达到精度。若用复化Simpson公式，由式(4-19)180 2f(4)()h4180 16411180 16 n104即得n 1.6

24、2。故应取n 2。4龙贝格求积公式4.1梯形公式的逐次分半算法如前所述，复化求积公式的截断误差随着步长的缩小而减少，而且如果被积函数的高阶导数容 易计算和估计时，由给定的精度可以预先确定步长，不过这样做常常是很困难的，一般不值得推崇。 实际计算时，我们总是从某个步长出发计算近似值，若精度不够可将步长逐次分半以提高近似值， 直到求得满足精度要求的近似值。设将区间a,b分为n等分，共有n 1个分点，如果将求积区间再二分一次，则分点增至2n 1个，我们将二分前后两个积分值联系起来加以考虑。注意到每个子区间xk,xk 1经过二分只增加了1一个分点x 1 (xk xk 1),用复化梯形公式求得该子区间上

25、的积分值为k 22hf(xk) 2f(xk 12) f(xk1)4 ba汪意，这里h 代表二分前的步长，将每个子区间上的积分值相加得T2nf(Xk)f (Xk i)f (xk 1 2 )T2n2Tnf ( Xk 12)(4-20)h这表明，将步长由h缩小为一时，T2n等于Tn的一半再加新增加节点处的函数值乘以当前步长。 2算法4.11 .输入 a, b, f (x),ba2 .置 m 1,h ,T0 h f (a) f(b)2m 13 .置 F 0,对 k 1,2,L ,2F F f (a (2k 1)h)14 . T -T0 hF2_. h. _1.1 TTO3,输出I T ,停机；否则 m

26、 1m,-h,TT0,转3。1.2 李查逊(Richardson )外推法h的函数，记为I1(h),相应的误假设用某种数彳1方法求量I的近似值，一般地，近似值是步长差为I L(h)其中 i(i 1,2,L ),0 P1 P2 L1hPkI Ii( h)h) nhLP1Pi式(4-22)减去式(4-21)乘以取满足其中b由式(4-23)1,I1(2(以1h)P2P，,得P1IP1)hP2I1(h)P33(P1除上式两端，得11( h)2(Pi1P1).(1RPl(h)Pl)(i 2,3,LI2(h)2hP2L khPk2k是与h无关的常数。若用2( h)P2P2 P22 hR)hP3b2hP2

27、b3hP3)仍与h无关。令I1( h)"(h)k( h)PkP/Pkk(PkbkhPk以I2S)作为I的近似值，其误差至少为O(hP2),因此不断重复以上作法，可以得到一个函数序列I m 1 ( h)Im(h) 71PmPm 11Im1(h)，m(4-21)h代替(4-21)中的h ,则得 LR)hPk L(4-22)(4-23)I2(h)收敛于I的速度比I(h)快。2,3,L(4-24)以Im(h)近似I ,误差为I Im(h) O(hPm)。随着m的增大，收敛速度越来越快，这就是Richardson外推法。1.3 龙贝格求积公式由前面知道，复化梯形公式的截断误差为O(h2)。进一

28、步分析，我们有如下欧拉一麦克劳林(Euler-Maclaurin )公式：定理4.4设f(x) C a,b,则有2h4Lkh2kL_2I T(h) 1h其中系数k(k 1,2,L)与h无关。把李查逊外推法与欧拉麦克劳林公式相结合，可以得到求积公式的外推算法。特别地，在外1 一推算法式(4-24)中，取 -,pk 2k ,并记T0(h) T(h),则有j h4jTm1(-) Tm1(h)(4-25)(4-26)Tm(h)象二,m 1,2,L41经过m(m 1,2,L )次加速后，余项便取下列形式： 2(m 1)2( m 2)Tm(h) I 1h2h L上述处理方法通常称为 李查逊(Richard

29、son )外推加速方法。为研究Romberg求积方法的机器实现，引入记号：以T0(k)表示二分k次后求得的梯形值，且以Tm?)表示序列To(k)的m次加速值，则依以上递推公式得到Tmk)(h)-4Tmk11)-m1-Tmk1,k 1,2,L4141称为龙贝格求积算法。Romberg公式的计算过程见下表4-2表4-2kh丁(k)TO丁(k)T1丁 (k)T2丁(k)T3丁 (k)T 4L0b aT(0)1b aTT(0)21。 12b a4To(2)工T2(0)3b a丁(3)T(2)T丁(0)8TOT1T2T34b aT(4)t13)t2(2)T3(1)T4(0)16MMMMMMMO算法4.2

30、输入a,b, f (x),(0) h(2)置 h b a,To-f(a) f(b) ,k 1(3)计算ok 12To(k)(h)1To(k1) h f(a (i 1)h)2i i2对 j 1,L ,kj (k j 1) (k j)(k j) 4 Tj 1Tj 1j4j1bh, k 1 k,返回(3)。(4) TkTk(°1),输出 Tk(0)a f(x)dx 停机；否则-例4用Romberg算法计算积分I0x32dx。解 利用逐次分半算法(4-20)和Romberg算法(4-25),计算结果见表 4-3。(0)1.X -f (0) f (1) 0.500000琛)1To(0) 0.5

31、 f(-) 0.42677722(2)1 (1)13T()T00.25 f (-) f (-) 0.407018244T1T。 015 f(1) f(3) f(5) f(7) 0.401812228888M表4-3kT0(k)T1(k)T2(k)T3(k)T(k)T5(k)00.50000010.4267770.40236920.40701830.4018120.4000770.4000540.40005040.4004630.4000090.40000950.4001180.4000020.4000020.4000025高斯求积公式5.1 一般理论等距节点的插值型求积公式，虽然计算简单，使用

32、方便，但是这种节点等距的规定却限制了求 积公式的代数精度。试想如果对节点不加限制，并适当选择求积系数，可能会提高求积公式的精度。Gauss型求积公式的思想也正如此，亦即在节点数n固定时，适当地选取节点xj与求积系数Ak,使求积分公式具有最高精度。设有n 1个互异节点 ,x1,L ,4的机械求积分公式(x) f(x)dxAk f(xk)ak0(2-27)具有 m 次代数精度，那么有取f (x)nAj(xj)lj0xl (l0,1,L ,m) ，式 (1)精确成立，即bl(x)x dx (l 0,1,L , m)a式(2)构成 m 1 阶的非线性方程组， 且具有 2n 2个未知数xk, Ak (k

33、 0,1,L , n) ， 所以当(2-28)(x)给定后，只要m 1 2n 2 ，即 m 2n 1 时，方程组有解。这表明式n 1 个节点的求积公式的代数精度可达到 2n 1。另一方面，对式(1)，不管如何选择xk 与 Ak ，最高精度不可能超过2n 1 。事实上，对任意的互异节点 xk kn 0 ，令p2n 2(x)n2 1(x) (x x0)2(x x1)2L (x xn)2nb有AkP2n 2(xk) 0，然而(x)P2n 1(x)dx 0。ak0定义 4 如果求积分公式 (4-27)具有 2n 1次代数精度，则称这组节点xk 为 Gauss 点 ，相应公式 (4-27) 称为带权 (

34、x) 的 高斯求积公式。定理 5 插值型求积公式的节点 a x0x1Lxnb 是高斯点的充分必要条件是以这些节点为零点的多项式n 1(x) (x x0)(x x1)L (x xn )与任何次数不超过 n的多项式P(x)带权正交，即(x)P(x) n 1(x)dx(4-29)证明 必要性。设P(x)H n ，则 P(x) n 1 (x) H 2n 1 ，因此，如果x0 ,x1 ,L , xn 是高斯点，则式(1)对于f (x)P(x)n 1 (x) 精确成立，即有ba (x)P(x) n 1(x)dxnAkP(xk) n 1(xk) k0故式 (4-29)成立。再证充分性。对于对于f(x) H2n 1 ，用1 (x) 除 f (x) ，记商为P(x) ，余式为 q(x) ，即f(x) P(x) n 1(x)q(x),其中 P(x),q(x)b(x) f (x)dxaH n ，由式 (4-29)可得(x)q(x)dx(4-30)(4-27) 是插值型的，它对于 q(x)b(x) f (x)dxa再注意到 n 1(xk)0 (k 0,1,L ,n) ，知q(xk)H n 是精确成立的，即nAkq(xk)k0f (xk ) ( k 0,1,L , n) ，从而由 (4-30)有b(x)f(x)dxan(x)q(x)dxAk f (xk)k0可见求积公式(4-27)对一切次