南开大学计算物理讲义5

1、6 Numerical integration数值积分数值积分Integration)()()(aFbFdxxfba我们知道，如果函数我们知道，如果函数f (x)在区间在区间a, b上连续，上连续，且原函数为且原函数为F(x)，则可用牛顿，则可用牛顿莱布尼兹公莱布尼兹公式式来求得定积分。例如：来求得定积分。例如：0)sin(？dxx02)0cos(cos)sin(dxxdxeax02然而，对有些函数来说，找到原函数往往很困然而，对有些函数来说，找到原函数往往很困难，有些原函数不能用初等函数表示出来，例难，有些原函数不能用初等函数表示出来，例如：如：（a为有限数）为有限数） 还有些原函数尽管能用

2、初等函数的有限形式还有些原函数尽管能用初等函数的有限形式表示出来，但表达式过于复杂，也不便使用。表示出来，但表达式过于复杂，也不便使用。特别是在实际问题中，更多的函数是用表格或特别是在实际问题中，更多的函数是用表格或图形表示的，对这种函数，更无法用牛顿图形表示的，对这种函数，更无法用牛顿莱布尼兹公式求积分。因此，有必要研究用数莱布尼兹公式求积分。因此，有必要研究用数值方法求定积分的问题。这种数值积分方法也值方法求定积分的问题。这种数值积分方法也是微分方程数值解法的基础。是微分方程数值解法的基础。6.1 Manual method; 6.2 Constant rule6.3 Trapezium

3、rule; 6.4 Mid-point rule6.5 Simpsons rule; 6.6 Romberg integration6.7 Gauss quadrature6.8 Example of numerical integration6.9 Numerical Integration of Multi- dimensions6.1 Manual methodSuppose we need to integrate from x0 to x1. We shall subdivide this interval into n steps of size x=(x1x0)/n as sho

4、wn in figure.6.2 Constant rulex)f(xf(x)dxI(f)xxx000Integration of a Taylor Series expansion of f(x) shows the error in this approximation to bexxxf(x)dxI(f)00 xxxdxxx)(xfxx)(xf)f(x00)( 21)(20000030200 6121x)(xfx)(xfx)f(x)x(x)f(x20O6.3 Trapezium rule （梯形公式梯形公式 ）xxf(x)f(xf(x)dxI(f)xxx)210000 xxxf(x)dx

5、I(f)00300020000 121)21) 312121x)(xfxxf(x)f(xxx)(xfx)(xf)f(x)f(x（30200200000 6121)( 21)(00 x)(xfx)(xfx)f(xdxxx)(xfxx)(xf)f(xxxxxxf(x)f(xI(f)21006.4 Mid-point rulex)f(xf(x)dxxxxx0212100 x0 xxxxxxxxdxxx)(xfxx)(xf)f(xf(x)dx212120000021210000)( 21)(300 241x)(xfx)f(x xxxf(x)dx00 Trapezium rule3000 121)21x

6、)(xfxxf(x)f(xbafabdxxf)()()()(f积分中值定理积分中值定理: :如果函数如果函数 f(x) 在区间在区间a,b上连续上连续, ,则在积分区间则在积分区间a,b内存在一点使内存在一点使下式下式成立。成立。由于 的具体位置一般是未知的,因而难以准确地计算出 .Constant ruleTrapezium ruleMid-point ruleExampleComputing the approximate value for the integralTrapezium rule Mid-point rule102dxxTrapezium rule Mid-point ru

7、lefrom which we conclude that the error in M is about 1/12, and the error in T is about 1/6.31|31103102xdxxSuppose we need to integrate from x0 to x1. We shall subdivide this interval into n steps of size x=(x1x0)/n as shown in figure复化复化 Compound 复化梯形公式复化梯形公式 Compound Trapezium rulexxf(x)f(xf(x)dxI

8、(f)xxx)210000The Compound Trapezium Rule approximation to the integral is thereforeWhile the error for each step is O( x3), the cumulative error is n times this or O( x2) O(n-2) x=h ( x=(x1x0)/n )1. Choose the number of mesh points and fix the step. 2. Calculate f(a) and f(b) and multiply with h/2.

9、3. Perform a loop over i=1 to n-1 (f(a) and f(b) are known) and sum up the terms f(a+x)+ f(a+2x)+ f(a+3x)+ + f(b-x)4. Multiply the final result by h and add (f(a)+ f(b)* h/2. Mid-point rulex)f(xf(x)dxxxxx0212100There are two further advantages of the Mid-point Rule over the Trapezium Rule. The first

10、 is that it is requires one fewer function evaluations for a given number of subintervals.The second is that it can be used more effectively for determining the integral near an integrable singularity.bandxxP)(badxxf)(数值积分的基本思想是构造一个简单函数（例数值积分的基本思想是构造一个简单函数（例如多项式）如多项式）Pn(x) 来近似代替被积分函数来近似代替被积分函数 f(x)，

11、然后通过求然后通过求 求得求得的近似值。的近似值。插值型求积公式插值型求积公式 由拉格朗日插值公式由拉格朗日插值公式)()()(0knkknxfxlxP)()()()(00knkbakknkbakxfdxxldxxfxlI6.5 Simpsons rule1两点公式两点公式构造以构造以a，b为结点的线性插值多项式为结点的线性插值多项式)()()(1bfabaxafbabxxPdxbfabaxafbabxdxxPTbaba)()()(1 y(x)yxxxxyxxxxy(x)21211212)()()(1bfabaxafbabxxP)()()(21)(21)()(21)()()()()()()()

12、(221bfafabababbfbabaafdxaxabbfdxbxbaafdxbfabaxafbabxdxxPTbabababa f(b)f(a)a)(b21dxf(x)TbaTrapezium rule2一点公式一点公式构造以构造以 (x0 x/2)和和(x0+ x/2)的中点为结点的线性的中点为结点的线性插值多项式插值多项式)()(00 xfxP x)f(x)dxf(xdxxPTxxxxxxxx021210212100000)( Mid-point rule2bac)()()( )()()()()()()(2bfcbabcxaxcfbcacbxaxafbacabxcxxP)()()()(

13、)()( )()()()()()()(2102bfcfafdxcbabcxaxbfdxbcacbxaxcfdxbacabxcxafdxxPSbabababa3三点公式三点公式，b为节点，构造二次插值多项式为节点，构造二次插值多项式 以以a, )(61)()(2)(31)(1)()(1)()(1)()(223320ababbccbababbacadxbcxcbxbacadxbxcxbacadxbacabxcxbababa)(64)()(1abdxbcacbxaxba)(61)()(2abdxcbabcxaxba)()(4)(6bfcfafabS称为辛卜生（称为辛卜生（Simpson）求积公式。）

14、求积公式。我们也可以类似地构我们也可以类似地构造五点、七点，造五点、七点，n n点求积公式，在实际点求积公式，在实际运算中，以梯形和运算中，以梯形和SimpsonSimpson公式为常用。公式为常用。x0 x0+ xx0+2 xacb2bac，b为节点，构造二次插值多项式为节点，构造二次插值多项式以以a, Simpsons rulex0 x0+ xx0+2 xacb)()(4)(6bfcfafabS)()(020 xxxxf(x),)2()()(000 xxfxxfxfx，20002)2()(2)(xxxfxxfxfxxxfxxfxf2)2()(4)(3000)(0 xf)()(4)(6bfc

15、fafabSx0 x0+ xx0+2 xacbinterpolation at one, two, and three equally spaced points on the interval a, b gives the first three Newton-Cotes quadrature rules:)(30212100 xOx)f(xf(x)dxxxxx1. Interpolating the function value at the midpoint of the interval by a constant gives a one-point quadrature rule k

16、nown as the midpoint rule2. Interpolating the function values at the two endpoints of the interval by a straight line (i.e., a polynomial of degree one) gives a two-point quadrature rule known as the trapezoid rule:)()2130000 xOxxf(x)f(xf(x)dxI(f)xxx3. Interpolating the function values at three poin

17、ts (the two endpoints and the midpoint) by a quadratic polynomial gives a three-point quadrature rule known as Simpsons rule:)()()(4)(65xObfcfafabSTrapezium rule Mid-point rule Simpsons rulex)f(xf(x)dxxxxx0212100 xxf(x)f(xf(x)dxT(f)xxx)210000)()(4)(6bfcfafabSTrapezium rule Mid-point rule Simpsons ru

18、lexx)f(xf(x)dxfMxxx2)(0200 xxf(x)f(xf(x)dxT(f)xxx2)22100200)(32)(31)(fMfTfS243)(000200 x)f(xx)f(x)f(xxf(x)dxfSxxxExampleError Estimation. We illustrate error estimation by computing the approximate value for the integralTrapezium rule Mid-point rule Simpsons rule102dxxTrapezium rule Mid-point rulefr

19、om which we conclude that the error in M is about 1/12, and the error in T is about 1/6.31|31103102xdxx)()(4)(6bfcfafabS Simpsons rule3/14/14161S31|31103102xdxx15.0dxxAs an example, we approximate the integralTrapezium rule Mid-point rule Simpsons rule15 .0dxx精确值为精确值为 4267767. 0)15 . 0(25 . 0115 . 0

20、dxx梯形公式梯形公式 43093403. 0)175. 045 . 0(65 . 0115 . 0dxx Simpsons rule Mid-point rule43301.075.0*)5.01(15.0dxx4309644. 0|3210.53/2xTrapezium rule Mid-point rule Simpsons ruleThe correctly rounded result for this problem is 0.746824Compound Simpsons rule复化辛卜生公式复化辛卜生公式 Compound Simpsons rule复化辛卜生公式复化辛卜生公

21、式 n=2mand the corresponding error O(nx5) or O(x4).xk01/81/43/81/25/83/47/81f (xk)43.938463.764703.506853.200002.876402.460002.265492利用数据表计算积分利用数据表计算积分dxxI102*14用复化梯形公式 用复化辛卜生公式 dxxI102*141415926. 3|arctg410* xI 13899. 31872432852212832412812)0(21818fffffffffT这个问题有明显的答案这个问题有明显的答案取取n = 8用复化梯形公式用复化梯形公式

22、 14159. 31874432854212834412814)0(81314fffffffffS取取n=4, 用辛卜生公式用辛卜生公式 With the Trapezium Rule6.6 Romberg integration龙贝格求积公式龙贝格求积公式 )()21300 xOxxf(x)f(xxxxf(x)dxI(f)00With the Compound Trapezium RuleWhile the error for each step is O( x3), the cumulative error is n times this or O( x2) O(n-2) x=hWith

23、the Compound Trapezium Rule we know from above section the error in some estimate T( x) of the integral I using a step size x goes like c x2 as x0, for some constant c. Likewise the error in T( x/2) will be c x2/4. From this we may construct a revised estimate T(1)( x/2) for I as a weighted mean of

24、T( x) and T( x/2): Romberg integration龙贝格求积公式龙贝格求积公式 T( x) c x2T( x/2) c x2/42/2T( x) c x2T( x/2) c x2/4By choosing the weighting factor = 4/3 we eliminate the leading order (O(x2) error terms, relegating the error to O(x4). Thus we have /4+1- =04/3 I+c(/4+1- x2+O (x4 )()(4)(62bffafabba 3134 112STT)

25、()(231)()(2)(4342bfafabbffafabba123134TT 144 222nnnCSS 3134 112STT 3134 224STT 3134 448STT 144 2nnnSTT 龙贝格公式同理可得同理可得 144 222nnnCSS 一般地，有一般地，有 144 2nnnSTT 144 323nnnRCC 龙贝格公式龙贝格公式注：注：(1)(1)上述加速技巧称为上述加速技巧称为龙贝格求积算法龙贝格求积算法； (2)(2)每加速一次，计算精度提高每加速一次，计算精度提高二阶二阶； (3)(3)该技巧可以不断继续下去，该技巧可以不断继续下去， 但通常但通常最多用到龙贝格

26、公式最多用到龙贝格公式。Romberg公式计算方法 T1 T2 S1 T4 S2 C1 T8 S4 C2 R1 T16 S8 C4 R2 T32 S16 C8 R4. 上面是Romberg的计算表若 则计算停止；否则用i+1代i，转入下一步。RRnn1221442TTSnnn144323CCRnnn144222SSCnnnabhn0, 1)(21)(2100bfafhT22, 201abhhn221210001hafhTT2222,2abhn02312012212221hiafhTTi按照区间逐次分半的方法，计算梯形和序列按照区间逐次分半的方法，计算梯形和序列 )()(2)(4T21bffaf

27、abbakkkabhn2,2021012122211hiafhTTkikkkk ,ba)43(),4(abafabaf),(),(bfaf1T(1)计算算出1S,ba),2(baf2T(2)把2等分计算算出与(3)把4等分,计算4T2S1C算出, 与;龙贝格求积的计算步骤如下：算出 T8 , S4 , C2 ,算出 R1用龙贝格求积法计算积分用龙贝格求积法计算积分10214dxxI的近似值的近似值, , xk01/81/43/81/25/83/47/81f (xk)43.938463.764703.506853.200002.876402.460002.265492214)(xxf; 2) 1

28、 (, 4)0(ff. 3)1 ()0(211ffT; 2 . 3)21(f; 1 . 3)21(212112fTT.13333. 33134121TTS 10214dxxI214)(xxf; 2) 1 (, 4)0(ff; 2 . 3)21(f)( fS13333. 3f(1)4f(0.5)(f(0)61s;56. 2)43(,76471. 3)41(ff;13118. 3)43()41(412124ffTT;14157. 33134242TTS.14212. 31511516121SSC 214)(xxf;50685. 3)83(,93846. 3)81(ff;26549. 2)87(,8

29、7640. 2)85(ff;13899. 3)87()85()83()81(812148ffffTT;14159. 33134484TTS;14159. 31511516242SSC.14158. 36316364121CCR ;14159. 3,14094. 3816ST.14159. 3,14159. 324RC;14159. 3,14143. 31632ST.14159. 3,14159. 348RC14159. 314102dxx 所以I 的准确值为14159265. 3arctan41410102xdxxI6.7 Gauss quadrature在插值型求积公式中，插值节点是事先固定

30、的，在插值型求积公式中，插值节点是事先固定的，有时还进一步限定是等距的。现在进一步考虑，有时还进一步限定是等距的。现在进一步考虑，在节点个数一定的情况下，是否可以在在节点个数一定的情况下，是否可以在a, ba, b上自由选择节点的位置，使求积公式的精度提上自由选择节点的位置，使求积公式的精度提得更高得更高? ? 高斯型求积公式高斯型求积公式代数精确度代数精确度为一般求积公式。这里为一般求积公式。这里Ak为不依赖为不依赖f (x)的常数的常数. .)()(0inikbaxfAdxxf（1）若若（1）式）式对任意不高于对任意不高于m次的多项式次的多项式精确精确成立，成立，而对于而对于xm+1不能精

31、确成立，就说不能精确成立，就说（1）式具有式具有m次代数精确度。次代数精确度。mkdxxxAbaknikii, 2 , 1 , 0,0称：称： 例例 ：形如形如 111100)()()(xfAxfAdxxf的两点求积公式。的两点求积公式。 用梯形公式（即以用梯形公式（即以x0 = -1，x1 = 1为节点的插为节点的插值型求积公式）立即可得值型求积公式）立即可得 ) 1 () 1()(11ffdxxfmkdxxxAbaknikii, 2 , 1 , 0,0) 1 () 1()(11ffdxxf一次代数精确度。一次代数精确度。 mkdxxxAbaknikii, 2 , 1 , 0,0cxf)(c

32、dxxf2)(11左左= = 右右= =c+c=2cc+c=2cxxf)(0|21)(11211xdxxf左左= = 右右= =-1+1=0-1+1=02)(xxf32|31)(11311xdxxf左左= = 右右= =1+1=21+1=2)1 ()0(4) 1(31)(11fffdxxfmkdxxxAbaknikii, 2 , 1 , 0,0cxf)(cdxxf2)(11左左= = 右右=(=(c+4c+c)/3c+4c+c)/3xxf)(0|21)(11211xdxxf左左= = 右右= =-1+1=0-1+1=02)(xxf32|31)(11311xdxxf左左= = 右右=(=(1+1

33、)/31+1)/33)(xxf0|41)(11411xdxxf左左= = 右右= =-1+1=0-1+1=04)(xxf52|51)(11511xdxxf左左= = 右右=(=(1+1)/31+1)/3作为插值型求积公式它至少具有作为插值型求积公式它至少具有n次代数精确度次代数精确度 mkdxxxAbaknikii, 2 , 1 , 0,0)()(0inikbaxfAdxxf可以证明可以证明其最大可能的代数精确度是其最大可能的代数精确度是?可以证明可以证明: :对于对于n+1个节点的求积公式，其最大个节点的求积公式，其最大可能的代数精确度是可能的代数精确度是2n + 1次。次。)()(0ini

34、kbaxfAdxxf可以证明可以证明, ,含有含有n+1个节点而代数精确度为个节点而代数精确度为2n + 1的插值型求积公式是存在的的插值型求积公式是存在的, ,它所用的节点是它所用的节点是a,ba,b上的第上的第n+1次正交多项式的零点。次正交多项式的零点。称这一称这一 类求积公式为高斯类求积公式为高斯(Gauss)(Gauss)公式公式, ,其节点其节点为高斯点。为高斯点。 梯形公式梯形公式Simpson公式公式可以证明可以证明: :对于对于n+1个节点的求积公式，其最大个节点的求积公式，其最大可能的代数精确度是可能的代数精确度是2n + 1次。次。对于对于2 个节点的求积公式，其最大可能

35、的代数精个节点的求积公式，其最大可能的代数精确度是确度是3 3次。次。 111100)()()(xfAxfAdxxf若对求积公式中的四个待定系数若对求积公式中的四个待定系数A0, A1, x0, x1适适当选取，使求积公式的代数精确度是当选取，使求积公式的代数精确度是3 3次次若对求积公式中的四个待定系数若对求积公式中的四个待定系数A0, A1, x0, x1适当选取，适当选取，使求积公式使求积公式对对f (x) = 1，x，x2，x3都都准确成立准确成立03202311300211200110010 xAxAxAxAxAxAAA 111100)()()(xfAxfAdxxf032023113

36、00211200110010 xAxAxAxAxAxAAA求解得110 AA3333)(11ffdxxf)()()(110011xfAxfAdxxf333310 xx此例告诉我们，对两点求积公式，只要适当选取此例告诉我们，对两点求积公式，只要适当选取求积节点求积节点x0, x1和求积系数和求积系数A0, A1，可使其代数精，可使其代数精确度达到三次。由此可见，求积公式的代数精确确度达到三次。由此可见，求积公式的代数精确度不仅与积分节点数目有关，而且与这些这点的度不仅与积分节点数目有关，而且与这些这点的所在位置有关。适当调整这些点的分布和求积系所在位置有关。适当调整这些点的分布和求积系数，能使求

37、积公式达到最高的代数精确度。本节数，能使求积公式达到最高的代数精确度。本节就是要讨论这种具有最高代数精确度的求积公式，就是要讨论这种具有最高代数精确度的求积公式，它们也叫做高斯型求积公式。它们也叫做高斯型求积公式。 构造高斯型求积公式的方法就是去找构造高斯型求积公式的方法就是去找a, b上上的的n + 1次正交多项式，再把它的次正交多项式，再把它的n + 1个零点个零点求出来，由于正交多项式具有性质求出来，由于正交多项式具有性质: :在在a, ba, b上的上的n + 1次正交多项式一定有次正交多项式一定有n + 1个不同零个不同零点，且全部位于点，且全部位于a, ba, b内，所以只要将此内

38、，所以只要将此n + 1个零点作为个零点作为n + 1次插值多项式的节点，构造出次插值多项式的节点，构造出的插值多项式即为高斯型求积公式。的插值多项式即为高斯型求积公式。3333)(11ffdxxf不失一般性，假定积分区间为不失一般性，假定积分区间为（-1, 1），），因为因为总可以利用变换总可以利用变换将区间将区间（a, b）变成变成（-1, 1）而积分变为：而积分变为：tabbx22adttgabdxxfba11)(2)(tababftg22)()()()(110011xfAxfAdxxf3333)(11ffdxxf3333)(11ffdxxf高斯高斯勒让德计算积分公式勒让德计算积分公式d

39、ttgabdxxfba11)(2)(tababftg22)(例例 399529. 25 . 111dxx运用高斯运用高斯勒让德公式计算积分勒让德公式计算积分两点梯形公式两点梯形公式三点辛卜生公式：三点辛卜生公式：3333)(11ffdxxf例例 399529. 25 . 111dxx401848. 2577350. 05 . 1577350. 05 . 15 . 111dxx288246. 25 . 11215 . 112125 . 111dxx395742. 25 . 25 . 145 . 0315 . 111dxx解：运用高斯解：运用高斯勒让德公式计算积分勒让德公式计算积分三点辛卜生公式：

40、三点辛卜生公式：两点梯形公式两点梯形公式20sinxdxI),1(4tx例例 用二点高斯用二点高斯-勒让德公式计算积分勒让德公式计算积分解解 作变量代换作变量代换则 tabbx22a3333)(11ffdxxf5773503. 0it20sinxdxI),1(4tx11204) 1(sin4sindttxdxI4) 1(sin)(ttf5773503. 0it例例解解 作变量代换作变量代换则则,因为节点因为节点得得3333)(11ffdxxf 32589. 04/)5773503. 01sin(0tf 94541. 04/)5773503. 01sin(1tf94541. 0)94541. 0

41、32589. 0(4)()(104tftfI20sinxdxI计算结果比用复合梯形公式计算结果比用复合梯形公式7个节点计算的结果个节点计算的结果还要好还要好 。所以所以,由二点高斯公式由二点高斯公式One widely used example of this is Gauss quadrature which enables exact integration of cubics with only two function evaluations (in contrast Simpsons Rule, which is also exact for cubics, requires thr

42、ee function evaluations). The Gauss Quadrature accurate to order 2M 1 may be determined using the same approach required for the two-point scheme. This may be derived by comparing the Taylor Series expansion for the integral with that)()()(0100010 xxfAxxfAdxxfxxxxthat for the points x0+x and x0+x:)(

43、)()(0100010 xxfAxxfAdxxfxxxx10AA2/10 xAAthat for the points x0+x and x0+x:=Gauss quadrature has the formula=)()()(0100010 xxfAxxfAdxxfxxxx2/10 xAA2, 10 xx3333)(11ffdxxf20sinxdxI20sinxdxI94541. 0)94541. 032589. 0(4)()(104tftfI4/)5773503. 01sin(4/)5773503. 01sin(4/Using three function evaluations: sym

44、metry suggests for the interval x0- x to x0+x the points should be evaluated at x0 and x0. The Taylor Series expansion of the integral)515(95)0(98)515(95)(11fffdxxf1, 00 xx)(157501)()6(ffR)515(95)0(98)515(95)(11fffdxxf三点高斯三点高斯- -勒让德公式为勒让德公式为 积分区间为无限的积分积分区间为无限的积分变量代换法变量代换法 x= -ln(t)dx= -dt/t6.8 Examp

45、le of numerical integrationC Trapezium rulenx = 1DO i=1,20Value = TrapeziumRule(x0,x1,nx)WRITE(6,*)nx,Value,Value - Exactnx = 2*nxENDDOREAL*8 FUNCTION TrapeziumRule(x0,x1,nx)C=parametersINTEGER*4 nxREAL*8 x0,x1C=functionsREAL*8 fC=local variablesINTEGER*4 iREAL*8 dx,xa,xb,fa,fb,Sumdx = (x1 - x0)/DFLOAT(nx)Sum = 0.0DO i=0,nx-1xa = x0 + DFLOAT(i)*dxxb = x0 + DFLOAT(i+1)*dxfa = f(xa) fb = f(xb)Sum = Sum + fa + fbENDDOSum = Sum * dx / 2.0Tra