Chapter 4 Numerical Differentiation and Integration Per-Olof - PowerPoint PPT Presentation

Chapter 4 Numerical Differentiation and Integration Per-Olof Persson persson@berkeley.edu Department of Mathematics University of California, Berkeley Math 128A Numerical Analysis

Numerical Differentiation Forward and Backward Differences Inspired by the definition of derivative: f ( x 0 + h ) − f ( x 0 ) f ′ ( x 0 ) = lim , h h → 0 choose a small h and approximate f ′ ( x 0 ) ≈ f ( x 0 + h ) − f ( x 0 ) h The error term for the linear Lagrange polynomial gives: f ′ ( x 0 ) = f ( x 0 + h ) − f ( x 0 ) − h 2 f ′′ ( ξ ) h Also known as the forward-difference formula if h > 0 and the backward-difference formula if h < 0

General Derivative Approximations Differentiation of Lagrange Polynomials Differentiate n f ( x k ) L k ( x ) + ( x − x 0 ) · · · ( x − x n ) � f ( n +1) ( ξ ( x )) f ( x ) = ( n + 1)! k =0 to get n k ( x j ) + f ( n +1) ( ξ ( x j )) � � f ′ ( x j ) = f ( x k ) L ′ ( x j − x k ) ( n + 1)! k =0 k � = j This is the ( n + 1) -point formula for approximating f ′ ( x j ) .

Commonly Used Formulas Using equally spaced points with h = x j +1 − x j , we have the three-point formulas 2 h [ − 3 f ( x 0 ) + 4 f ( x 0 + h ) − f ( x 0 + 2 h )] + h 2 f ′ ( x 0 ) = 1 3 f (3) ( ξ 0 ) 2 h [ − f ( x 0 − h ) + f ( x 0 + h )] − h 2 f ′ ( x 0 ) = 1 6 f (3) ( ξ 1 ) 2 h [ f ( x 0 − 2 h ) − 4 f ( x 0 − h ) + 3 f ( x 0 )] + h 2 f ′ ( x 0 ) = 1 3 f (3) ( ξ 2 ) h 2 [ f ( x 0 − h ) − 2 f ( x 0 ) + f ( x 0 + h )] − h 2 f ′′ ( x 0 ) = 1 12 f (4) ( ξ ) and the five-point formula 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) ( ξ )

Optimal h Consider the three-point central difference formula: 2 h [ f ( x 0 + h ) − f ( x 0 − h )] − h 2 f ′ ( x 0 ) = 1 6 f (3) ( ξ 1 ) Suppose that round-off errors ε are introduced when computing f . Then the approximation error is � � f ( x 0 + h ) − ˜ ˜ h + h 2 f ( x 0 − h ) � ≤ ε � � � f ′ ( x 0 ) − 6 M = e ( h ) � � 2 h � � where ˜ f is the computed function and | f (3) ( x ) | ≤ M Sum of truncation error h 2 M/ 6 and round-off error ε/h � Minimize e ( h ) to find the optimal h = 3 3 ε/M

Richardson’s Extrapolation Suppose N ( h ) approximates an unknown M with error M − N ( h ) = K 1 h + K 2 h 2 + K 3 h 3 + · · · then an O ( h j ) approximation is given for j = 2 , 3 , . . . by � h � + N j − 1 ( h/ 2) − N j − 1 ( h ) N j ( h ) = N j − 1 2 j − 1 − 1 2 The results can be written in a table: O ( h 2 ) O ( h 3 ) O ( h 4 ) O ( h ) 1: N 1 ( h ) ≡ N ( h ) 2: N 1 ( h 2 ) ≡ N ( h 2 ) 3: N 2 ( h ) 4: N 1 ( h 4 ) ≡ N ( h 5: N 2 ( h 4 ) 2 ) 6: N 3 ( h ) 7: N 1 ( h 8 ) ≡ N ( h 8: N 2 ( h 9: N 3 ( h 8 ) 4 ) 2 ) 10: N 4 ( h )

Richardson’s Extrapolation If some error terms are zero, different and more efficient formulas can be derived Example: If M − N ( h ) = K 2 h 2 + K 4 h 4 + · · · then an O ( h 2 j ) approximation is given for j = 2 , 3 , . . . by � h � + N j − 1 ( h/ 2) − N j − 1 ( h ) N j ( h ) = N j − 1 4 j − 1 − 1 2

Numerical Quadrature Integration of Lagrange Interpolating Polynomials Select { x 0 , . . . , x n } in [ a, b ] and integrate the Lagrange polynomial P n ( x ) = � n i =0 f ( x i ) L i ( x ) and its truncation error term over [ a, b ] to obtain � b n � f ( x ) dx = a i f ( x i ) + E ( f ) a i =0 with � b a i = L i ( x ) dx a and � b n 1 � ( x − x i ) f ( n +1) ( ξ ( x )) dx E ( f ) = ( n + 1)! a i =0

Trapezoidal and Simpson’s Rules The Trapezoidal Rule Linear Lagrange polynomial with x 0 = a , x 1 = b , h = b − a , gives � b 2[ f ( x 0 ) + f ( x 1 )] − h 3 f ( x ) dx = h 12 f ′′ ( ξ ) a Simpson’s Rule Second Lagrange polynomial with x 0 = a , x 2 = b , x 1 = a + h , h = ( b − a ) / 2 gives � x 2 3[ f ( x 0 ) + 4 f ( x 1 ) + f ( x 2 )] − h 5 dx = h 90 f (4) ( ξ ) x 0 Definition The degree of accuracy , or precision , of a quadrature formula is the largest positive integer n such that the formula is exact for x k , for each k = 0 , 1 , . . . , n .

The Newton-Cotes Formulas The Closed Newton-Cotes Formulas Use nodes x i = x 0 + ih , x 0 = a , x n = b , h = ( b − a ) /n : � b n � f ( x ) dx ≈ a i f ( x i ) a i =0 � x n � x n ( x − x j ) � a i = L i ( x ) dx = ( x i − x j ) dx x 0 x 0 j � = i n = 1 gives the Trapezoidal rule, n = 2 gives Simpson’s rule. The Open Newton-Cotes Formulas Use nodes x i = x 0 + ih , x 0 = a + h , x n = b − h , h = ( b − a ) / ( n + 2) . Setting n = 0 gives the Midpoint rule: � x 1 f ( x ) dx = 2 hf ( x 0 ) + h 3 3 f ′′ ( ξ ) x − 1

Composite Rules Theorem Let f ∈ C 2 [ a, b ] , h = ( b − a ) /n , x j = a + jh , µ ∈ ( a, b ) . The Composite Trapezoidal rule for n subintervals is   � b n − 1 f ( x ) dx = h  − b − a � 12 h 2 f ′′ ( µ )  f ( a ) + 2 f ( x j ) + f ( b ) 2 a j =1 Theorem Let f ∈ C 4 [ a, b ] , n even, h = ( b − a ) /n , x j = a + jh , µ ∈ ( a, b ) . The Composite Simpson’s rule for n subintervals is   � b ( n/ 2) − 1 n/ 2 f ( x ) dx = h � �  f ( a ) + 2 f ( x 2 j ) + 4 f ( x 2 j − 1 ) + f ( b )  3 a j =1 j =1 − b − a 180 h 4 f (4) ( µ )

Error Estimation The error term in Simpson’s rule requires knowledge of f (4) : � b f ( x ) dx = S ( a, b ) − h 5 90 f (4) ( µ ) a Instead, apply it again with step size h/ 2 : � b � h 5 � � � a + b � � a, a + b − 1 f (4) (˜ f ( x ) dx = S + S , b µ ) 2 2 16 90 a The assumption f (4) ( µ ) ≈ f (4) (˜ µ ) gives the error estimate � � b �� � a, a + b � � a + b � � f ( x ) dx − S − S , b � � � 2 2 � a � � � �� ≈ 1 � a, a + b � � a + b � � � S ( a, b ) − S − S , b � � 15 2 2 �

Adaptive Quadrature � b To compute a f ( x ) dx within a tolerance ε > 0 , first apply Simpson’s rule with h = ( b − a ) / 2 and with h/ 2 If � � � � a + b �� a, a + b � � � S ( a, b ) − S − S , b � < 15 ε � � 2 2 then the integral is sufficiently accurate If not, apply the technique to [ a, ( a + b ) / 2] and [( a + b ) / 2 , b ] , and compute the integral within a tolerance of ε/ 2 Repeat until each portion is within the required tolerance

Gaussian Quadrature Basic idea: Calculate both nodes x 1 , . . . , x n and coefficients c 1 , . . . , c n such that � b n � f ( x ) dx ≈ c i f ( x i ) a i =1 Since there are 2 n parameters, we might expect a degree of precision of 2 n − 1 Example: n = 2 gives the rule √ � √ � 1 � � � − 3 3 f ( x ) dx ≈ f + f 3 3 − 1 with degree of precision 3

Legendre Polynomials The Legendre polynomials P n ( x ) have the properties For each n , P n ( x ) is a monic polynomial of degree n (leading 1 coefficient 1) � 1 − 1 P ( x ) P n ( x ) dx = 0 when P ( x ) is a polynomial of degree 2 less than n The roots of P n ( x ) are distinct, in the interval ( − 1 , 1) , and symmetric with respect to the origin. Examples: P 0 ( x ) = 1 , P 1 ( x ) = x P 2 ( x ) = x 2 − 1 P 3 ( x ) = x 3 − 3 5 x 3 P 4 ( x ) = x 4 − 6 7 x 2 + 3 35

Gaussian Quadrature Theorem Suppose x 1 , . . . , x n are roots of P n ( x ) and � 1 n x − x j � c i = dx x i − x j − 1 j � = i If P ( x ) is any polynomial of degree less than 2 n , then � 1 n � P ( x ) dx = c i P ( x i ) − 1 i =1

Computing Gaussian Quadrature Coefficients MATLAB Implementation function [x, c] = gaussquad(n) % Compute Gaussian quadrature points and coefficients. P = zeros(n+1,n+1); P([1,2],1) = 1; for k = 1:n − 1 P(k+2,1:k+2) = ((2*k+1)*[P(k+1,1:k+1) 0] − ... k*[0 0 P(k,1:k)]) / (k+1); end x = sort(roots(P(n+1,1:n+1))); A = zeros(n,n); for i = 1:n A(i,:) = polyval(P(i,1:i),x)'; end c = A \ [2; zeros(n − 1,1)];

Arbitrary Intervals � b Transform integrals a f ( x ) dx into integrals over [ − 1 , 1] by a change of variables: t = 2 x − a − b ⇔ x = 1 2[( b − a ) t + a + b ] b − a Gaussian quadrature then gives � ( b − a ) � b � 1 � ( b − a ) t + ( b + a ) f ( x ) dx = f dt 2 2 a − 1

Double Integrals Consider the double integral �� f ( x, y ) dA, R = { ( x, y ) | a ≤ x ≤ b, c ≤ y ≤ d } R Partition [ a, b ] and [ c, d ] into even number of subintervals n, m Step sizes h = ( b − a ) /n and k = ( d − c ) /m Write the integral as an iterated integral � b �� d � �� f ( x, y ) dA = f ( x, y ) dy dx R a c and use any quadrature rule in an iterated manner.

Composite Simpson’s Rule Double Integration The Composite Simpson’s rule gives � b �� d n m � dx = hk � � f ( x, y ) dy w i,j f ( x i , y j ) + E 9 a c i =0 j =0 where x i = a + ih , y j = c + jk , w i,j are the products of the nested Composite Simpson’s rule coefficients (see below), and the error is h 4 ∂ 4 f µ ) + k 4 ∂ 4 f E = − ( d − c )( b − a ) � � ∂x 4 (¯ η, ¯ ∂y 4 (ˆ η, ˆ µ ) 180 1 4 2 4 1 d 4 16 8 16 4 1 4 2 4 1 c a b

Non-Rectangular Regions The same technique can be applied to double integrals of the form � b � d ( x ) f ( x, y ) dy dx a c ( x ) The step size for x is still h = ( b − a ) /n , but for y it varies with x : k ( x ) = d ( x ) − c ( x ) m