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、的，在插值型求積公式中，插值節點是事先固定的，有時還進一步限定是等距的。現在進一步考慮，有時還進一步限定是等距的?，F在進一步考慮，在節點個數一定的情況下，是否可以在在節點個數一定的情況下，是否可以在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