Numerical Analysis I Numerical Differentiation and Integration - PowerPoint PPT Presentation

Numerical Analysis I Numerical Differentiation and Integration Instructor: Wei-Cheng Wang 1 Department of Mathematics National TsingHua University Fall 2010 1These slides are based on Prof. Tsung-Ming Huang(NTNU)’s original slides Wei-Cheng Wang (NTHU) Numerical Diff. & Integ. Fall 2010 1 / 66

Outline 1 Numerical Differentiation 2 Richardson Extrapolation Method 3 Elements of Numerical Integration 4 Composite Numerical Integration 5 Gaussian Quadrature Wei-Cheng Wang (NTHU) Numerical Diff. & Integ. Fall 2010 2 / 66

Numerical Differentiation f ( x 0 + h ) − f ( x 0 ) f ′ ( x 0 ) = lim . h h → 0 Question How accurate is f ( x 0 + h ) − f ( x 0 ) ? h Suppose a given function f has continuous first derivative and f ′′ exists. From Taylor’s theorem f ( x + h ) = f ( x ) + f ′ ( x ) h + 1 2 f ′′ ( ξ ) h 2 , where ξ is between x and x + h , one has f ′ ( x ) = f ( x + h ) − f ( x ) − h 2 f ′′ ( ξ ) = f ( x + h ) − f ( x ) + O ( h ) . h h Wei-Cheng Wang (NTHU) Numerical Diff. & Integ. Fall 2010 3 / 66

Hence it is reasonable to use the approximation f ′ ( x ) ≈ f ( x + h ) − f ( x ) h which is called forward finite difference, and the error involved is | e | = h 2 | f ′′ ( ξ ) | ≤ h t ∈ ( x,x + h ) | f ′′ ( t ) | . max 2 Similarly one can derive the backward finite difference approximation f ′ ( x ) ≈ f ( x ) − f ( x − h ) (1) h which has the same order of truncation error as the forward finite difference scheme. Wei-Cheng Wang (NTHU) Numerical Diff. & Integ. Fall 2010 4 / 66

The forward difference is an O ( h ) scheme. An O ( h 2 ) scheme can also be derived from the Taylor’s theorem f ( x ) + f ′ ( x ) h + 1 2 f ′′ ( x ) h 2 + 1 6 f ′′′ ( ξ 1 ) h 3 f ( x + h ) = f ( x ) − f ′ ( x ) h + 1 2 f ′′ ( x ) h 2 − 1 6 f ′′′ ( ξ 2 ) h 3 , f ( x − h ) = where ξ 1 is between x and x + h and ξ 2 is between x and x − h . Hence f ( x + h ) − f ( x − h ) = 2 f ′ ( x ) h + 1 6[ f ′′′ ( ξ 1 ) + f ′′′ ( ξ 2 )] h 3 and f ′ ( x ) = f ( x + h ) − f ( x − h ) − 1 12[ f ′′′ ( ξ 1 ) + f ′′′ ( ξ 2 )] h 2 2 h Let z ∈ [ x − h,x + h ] f ′′′ ( z ) z ∈ [ x − h,x + h ] f ′′′ ( z ) . M = max and m = min Wei-Cheng Wang (NTHU) Numerical Diff. & Integ. Fall 2010 5 / 66

If f ′′′ is continuous on [ x − h, x + h ], then by the intermediate value theorem, there exists ξ ∈ [ x − h, x + h ] such that f ′′′ ( ξ ) = 1 2[ f ′′′ ( ξ 1 ) + f ′′′ ( ξ 2 )] . Hence f ′ ( x ) = f ( x + h ) − f ( x − h ) − 1 6 f ′′′ ( ξ ) h 2 = f ( x + h ) − f ( x − h ) + O ( h 2 ) . 2 h 2 h This is called center difference approximation and the truncation error is | e | = h 2 6 f ′′′ ( ξ ) Similarly, we can derive an O ( h 2 ) scheme from Taylor’s theorem for f ′′ ( x ) f ′′ ( x ) = f ( x + h ) − 2 f ( x ) + f ( x − h ) − 1 12 f (4) ( ξ ) h 2 , h 2 where ξ is between x − h and x + h . Wei-Cheng Wang (NTHU) Numerical Diff. & Integ. Fall 2010 6 / 66

Polynomial Interpolation Method Suppose that ( x 0 , f ( x 0 )), ( x 1 , f ( x 1 )) · · · , ( x n , f ( x n )) have been given, we apply the Lagrange polynomial interpolation scheme to derive n � P ( x ) = f ( x i ) L i ( x ) , i =0 where n x − x j � L i ( x ) = . x i − x j j =0 ,j � = i Since f ( x ) can be written as n 1 � ( n + 1)! f ( n +1) ( ξ x ) w ( x ) , f ( x ) = f ( x i ) L ( x ) + i =0 where n � w ( x ) = ( x − x j ) , j =0 Wei-Cheng Wang (NTHU) Numerical Diff. & Integ. Fall 2010 7 / 66

we have, 1 f ′ ( x ) � nf ( x i ) L ′ ( n + 1)! f ( n +1) ( ξ x ) w ′ ( x ) = i ( x ) + i =0 ( n + 1)! w ( x ) d 1 dxf ( n +1) ( ξ x ) . + Note that n n w ′ ( x ) = � � ( x − x i ) . j =0 i =0 ,i � = j Hence a reasonable approximation for the first derivative of f is n f ′ ( x ) ≈ � f ( x i ) L ′ i ( x ) . i =0 When x = x k for some 0 ≤ k ≤ n , n � w ′ ( x k ) = w ( x k ) = 0 and ( x k − x i ) . i =0 ,i � = k Wei-Cheng Wang (NTHU) Numerical Diff. & Integ. Fall 2010 8 / 66

Hence n n 1 � ( n + 1)! f ( n +1) ( ξ x ) � f ′ ( x k ) = f ( x i ) L ′ i ( x k ) + ( x k − x i ) , (2) i =0 i =0 ,i � = k which is called an ( n + 1)-point formula to approximate f ′ ( x ). • Three Point Formulas Since ( x − x 1 )( x − x 2 ) L 0 ( x ) = ( x 0 − x 1 )( x 0 − x 2 ) we have 2 x − x 1 − x 2 L ′ 0 ( x ) = ( x 0 − x 1 )( x 0 − x 2 ) . Similarly, 2 x − x 0 − x 2 2 x − x 0 − x 1 L ′ ( x 1 − x 0 )( x 1 − x 2 ) and L ′ 1 ( x ) = 2 ( x ) = ( x 2 − x 0 )( x 2 − x 1 ) . Wei-Cheng Wang (NTHU) Numerical Diff. & Integ. Fall 2010 9 / 66

Hence � 2 x j − x 1 − x 2 � � 2 x j − x 0 − x 2 � f ′ ( x j ) = f ( x 0 ) + f ( x 1 ) ( x 0 − x 1 )( x 0 − x 2 ) ( x 1 − x 0 )( x 1 − x 2 ) 2 � 2 x j − x 0 − x 1 � + 1 � 6 f (3) ( ξ j ) + f ( x 2 ) ( x j − x k ) , ( x 2 − x 0 )( x 2 − x 1 ) k =0 ,k � = j for each j = 0 , 1 , 2. Assume that x 1 = x 0 + h and x 2 = x 0 + 2 h, for some h � = 0 . Then + h 2 1 � − 3 2 f ( x 0 ) + 2 f ( x 1 ) − 1 � f ′ ( x 0 ) 3 f (3) ( ξ 0 ) , = 2 f ( x 2 ) h − h 2 1 � − 1 2 f ( x 0 ) + 1 � f ′ ( x 1 ) 6 f (3) ( ξ 1 ) , = 2 f ( x 2 ) h + h 2 1 � 1 2 f ( x 0 ) − 2 f ( x 1 ) + 3 � f ′ ( x 2 ) 3 f (3) ( ξ 2 ) . = 2 f ( x 2 ) h Wei-Cheng Wang (NTHU) Numerical Diff. & Integ. Fall 2010 10 / 66

That is + h 2 � � 1 − 3 2 f ( x 0 ) + 2 f ( x 0 + h ) − 1 f ′ ( x 0 ) 3 f (3) ( ξ 0 ) , = 2 f ( x 0 + 2 h ) h − h 2 1 � − 1 2 f ( x 0 ) + 1 � f ′ ( x 0 + h ) 6 f (3) ( ξ 1 ) , = 2 f ( x 0 + 2 h ) (3) h + h 2 1 � 1 2 f ( x 0 ) − 2 f ( x 0 + h ) + 3 � 3 f (3) ( ξ 2 ) . f ′ ( x 0 + 2 h ) = 2 f ( x 0 + 2 h ) (4) h Using the variable substitution x 0 for x 0 + h and x 0 + 2 h in (3) and (4), respectively, we have 2 h [ − 3 f ( x 0 ) + 4 f ( x 0 + h ) − f ( x 0 + 2 h )] + h 2 1 3 f (3) ( ξ 0 ) , (5) f ′ ( x 0 ) = 2 h [ − f ( x 0 − h ) + f ( x 0 + h )] − h 2 1 6 f (3) ( ξ 1 ) , f ′ ( x 0 ) = 2 h [ f ( x 0 − 2 h ) − 4 f ( x 0 − h ) + 3 f ( x 0 )] + h 2 1 3 f (3) ( ξ 2 ) . f ′ ( x 0 ) = (6) Note that (6) can be obtained from (5) by replacing h with − h . Wei-Cheng Wang (NTHU) Numerical Diff. & Integ. Fall 2010 11 / 66

• Five-point Formulas 1 f ′ ( x 0 ) = 12 h [ f ( x 0 − 2 h ) − 8 f ( x 0 − h ) + 8 f ( x 0 + h ) − f ( x 0 + 2 h )] h 4 30 f (5) ( ξ ) , + where ξ ∈ ( x 0 − 2 h, x 0 + 2 h ) and 1 f ′ ( x 0 ) = 12 h [ − 25 f ( x 0 ) + 48 f ( x 0 + h ) − 36 f ( x 0 + 2 h ) 16 f ( x 0 + 3 h ) − 3 f ( x 0 + 4 h )] + h 4 5 f (5) ( ξ ) , + where ξ ∈ ( x 0 , x 0 + 4 h ). Wei-Cheng Wang (NTHU) Numerical Diff. & Integ. Fall 2010 12 / 66

Round-off Error Consider 2 h [ − f ( x 0 − h ) + f ( x 0 + h )] − h 2 1 6 f (3) ( ξ 1 ) , f ′ ( x 0 ) = where h 2 6 f (3) ( ξ 1 ) is called truncation error. Let ˜ f ( x 0 + h ) and ˜ f ( x 0 − h ) be the computed values of f ( x 0 + h ) and f ( x 0 − h ), respectively. Then f ( x 0 + h ) = ˜ f ( x 0 + h ) + e ( x 0 + h ) and f ( x 0 − h ) = ˜ f ( x 0 − h ) + e ( x 0 − h ) . Therefore, the total error in the approximation f ( x 0 + h ) − ˜ ˜ − h 2 f ( x 0 − h ) = e ( x 0 + h ) − e ( x 0 − h ) f ′ ( x 0 ) − 6 f (3) ( ξ 1 ) 2 h 2 h is due in part to round-off error and in part to truncation error. Wei-Cheng Wang (NTHU) Numerical Diff. & Integ. Fall 2010 13 / 66

Assume that | e ( x 0 ± h ) | ≤ ε and | f (3) ( ξ 1 ) | ≤ M. Then � � f ( x 0 + h ) − ˜ ˜ h + h 2 f ( x 0 − h ) � ≤ ε � � � f ′ ( x 0 ) − 6 M ≡ e ( h ) . � � 2 h � � � 3 Note that e ( h ) attains its minimum at h = 3 ε/M . In double precision arithmetics, for example, ε ≈ | f ( x 0 ± h ) | × 10 − 16 . The √ 3 Mε 2 ) = O (10 − 10 ). minimum is O ( Wei-Cheng Wang (NTHU) Numerical Diff. & Integ. Fall 2010 14 / 66

Richardson’s Extrapolation Suppose ∀ h � = 0 we have a formula N 1 ( h ) that approximates an unknown value M M − N 1 ( h ) = K 1 h + K 2 h 2 + K 3 h 3 + · · · , (7) for some unknown constants K 1 , K 2 , K 3 , . . . . If K 1 � = 0, then the truncation error is O ( h ). For example, h 2 − f (4) ( x ) = − f ′′ ( x ) h − f ′′′ ( x ) f ′ ( x ) − f ( x + h ) − f ( x ) h 3 − · · · . h 2! 3! 4! Goal Find an easy way to produce formulas with a higher-order truncation error. Replacing h in (7) by h/ 2, we have h 2 h 3 � h � h M = N 1 + K 1 2 + K 2 4 + K 3 8 + · · · . (8) 2 Wei-Cheng Wang (NTHU) Numerical Diff. & Integ. Fall 2010 15 / 66

Subtracting (7) with twice (8), we get M = N 2 ( h ) − K 2 2 h 2 − 3 K 3 4 h 3 − · · · , (9) where � h � � h � � � h � � N 2 ( h ) = 2 N 1 − N 1 ( h ) = N 1 + N 1 − N 1 ( h ) , 2 2 2 which is an O ( h 2 ) approximation formula. Replacing h in (9) by h/ 2, we get � h � − K 2 8 h 2 − 3 K 3 32 h 3 − · · · . M = N 2 (10) 2 Subtracting (9) from 4 times (10) gives � h � − N 2 ( h ) + 3 K 3 8 h 3 + · · · , 3 M = 4 N 2 2 which implies that � � h � + N 2 ( h/ 2) − N 2 ( h ) � + K 3 8 h 3 + · · · ≡ N 3 ( h ) + K 3 8 h 3 + · · · M = N 2 2 3 Wei-Cheng Wang (NTHU) Numerical Diff. & Integ. Fall 2010 16 / 66