Numerical Dierentiation Hector D. Ceniceros 1 Numerical Dierentiation Suppose f is a dierentiable function. We would like to approximate f (x0 ) given values of f at x0 and at neighboring points x1 , x2 , ..., xn . We could approximate f by its interpolating polynomial Pn (x) at those points and use f (x0 ) Pn (x0 ). There are several other possibilities. For example, we can approximate f (x0 ) by the derivative of the cubic spline of f evaluated at x0 , or by the derivative of the Least Squares Chebyshev expansion of f , n f (x0 ) = aj Tj (x0 ), j=1 etc. If we use polynomial interpolation then (1) f (x) = Pn (x) + 1 f (n+1) ((x)) wn (x), (n + 1)! where (2) wn (x) = (x x0 )(x x1 ) (x xn ). Thus, f (x0 ) = Pn (x0 ) + d 1 wn (x) f (n+1) ((x)) + f (n+1) ((x)) wn (x) (n + 1)! dx . x=x0 These are lecture notes for Math 104 A. These notes and all course materials are protected by United States Federal Copyright Law, the California Civil Code. The UC Policy 102.23 expressly prohibits students (and all other persons) from recording lectures or discussions and from distributing or selling lectures notes and all other course materials without the prior written permission of the instructor. 1 But wn (x0 ) = 0 and wn (x0 ) = (x0 x1 ) (x0 xn ), thus (3) f (x0 ) = Pn (x0 ) + 1 f (n+1) (0 )(x0 x1 ) (x0 xn ), (n + 1)! where 0 = (x0 ). Example 1. Take n = 1 and x1 = x0 + h (h > 0). In Newton's form (4) P1 (x) = f (x0 ) + f (x0 + h) f (x0 ) (x x0 ), h 1 and P1 (x0 ) = h [f (x0 +h)f (x0 )]. We obtain the so-called Forward Dierence Formula for approximating f (x0 ) + Dh f (x0 ) := (5) f (x0 + h) f (x0 ) . h From (3) the error in this approximation is + f (x0 ) Dh f (x0 ) = (6) 1 1 f (0 )(x0 x1 ) = f (0 )h, 2! 2 where x0 < 0 < x0 + h. Example 2. Take again n = 1 but now x1 = x0 h. Then P1 (x0 ) = 1 [f (x0 ) f (x0 h)] and we get the so-called Backward Dierence Formula h for approximating f (x0 ) Dh f (x0 ) := (7) f (x0 ) f (x0 h) . h Its error is (8) 1 f (x0 ) Dh f (x0 ) = f (0 )h, 2 where now x0 h < 0 < x0 . Example 3. Let n=2 and x1 = x0 h, x2 = x0 + h. P2 in Newton's form is P2 (x) = f [x1 ] + f [x1 , x0 ](x x1 ) + f [x1 , x0 , x2 ](x x1 )(x x0 ). 2 Let us obtain the nite dierence table: x0 h f (x0 h) f (x0 )f (x0 h) h x0 f (x0 ) f (x0 +h)f (x0 ) h f (x0 +h)2f (x0 )+f (x0 +h) 2h2 x0 + h f (x0 + h) Therefore, P2 (x0 ) = f (x0 ) f (x0 h) f (x0 + h) 2f (x0 ) + f (x0 h) + h h 2h2 and thus (9) P2 (x0 ) = f (x0 + h) f (x0 h) . 2h This denes the Centered Dierence Formula to approximate f (x0 ) 0 Dh f (x0 ) := (10) f (x0 + h) f (x0 h) . 2h Its error is (11) 0 f (x0 ) Dh f (x0 ) = 1 1 f (0 )(x0 x1 )(x0 x2 ) = f (0 )h2 , 3! 6 for some 0 between x0 h and x0 + h. Example 4. Let n = 2 and x1 = x0 + h, x2 = x0 + 2h. The table of nite dierences is x0 f (x0 ) f (x0 +h)f (x0 ) h x0 + h f (x0 + h) f (x0 +2h)f (x0 +h) h f (x0 +2h)2f (x0 +h)+f (x0 ) 2h2 x0 + 2h f (x0 + 2h) and P2 (x0 ) = f (x0 + h) f (x0 ) f (x0 + 2h) 2f (x0 + h) + f (x0 ) h h 2h2 3 thus P2 (x0 ) = (12) f (x0 + 2h) + 4f (x0 + h) 3f (x0 ) 2h If we use this sided dierence to approximate f (x0 ) the error is f (x0 ) P2 (x0 ) = (13) 1 1 f (0 )(x0 x1 )(x0 x2 ) = h2 f (0 ), 3! 3 where x0 < 0 < x0 + 2h, and the error is (asymptotically) twice as large as that in the Centered Finite Dierence Formula. 2 The eect of round-o errors In numerical dierentiation we take dierences of values which could be very close to each other, for small h. As we know, this leads to loss of accuracy because of the nite precision oating point arithmetic. Consider for example the centered dierence formula. For simplicity let us suppose that h has an exact oating point representation and that we make no rounding error when doing the division by h. That is, suppose that the the only source of round-o error is in the computation of the dierence f (x0 + h) f (x0 h). Then f (x0 + h) and f (x0 h) are replaced by f (x0 + h)(1 + + ) and f (x0 h)(1 + ), respectively with |+ | and | | . Then f (x0 + h) f (x0 h) f (x0 + h)(1 + + ) f (x0 h)(1 + ) = + rh , 2h 2h where rh = f (x0 + h)+ f (x0 h) . 2h Clearly, |rh | (|f (x0 + h)| + |f (x0 h)|) 2h |f (x0 )| h . The approximation error or truncation error for the centered nite dierence approximation is 1 6 f (0 )h2 . Thus, the total error E(h) can be approximately bounded by 1 2 h M3 + |f (x0 )| h . The minimum error occurs at h0 such that E (h0 ) = 0, 6 i.e. (14) h0 = 3 |f (x0 )| M3 4 1 3 c 1 3 2 and E(h0 ) = O( 3 ). We do not get machine precision! Higher order nite dierences exacerbate the problem of digit cancellation. When f can be extended to an analytic function in the complex plane, Cauchy Integral Theorem can be used to evaluate the derivative: (15) f (z0 ) = 1 2i C f (z) dz, (z z0 )2 where C is a simple closed contour around z0 and f is analytic on and inside C. Parametrizing C as a circle of radius r centered at x0 we get (16) f (x0 ) = 1 2r 2 f (x0 + rei )ei d 0 The integrand is periodic and smooth so it can be approximated with spectral accuracy with the composite trapezoidal rule. Another approach to obtain nite dierence formulas to approximate derivatives is through Taylor explasions. For example, (17) (18) 1 1 f (x0 )h2 + f (3) (x0 )h3 + . . . , 2! 3! 1 1 2 f (x0 h) = f (x0 ) f (x0 )h + f (x0 )h f (3) (x0 )h3 + . . . . 2! 3! f (x0 + h) = f (x0 ) + f (x0 )h + Then subtracting (17) from (18) we have f (x0 + h) f (x0 h) = 2f (x0 )h + 2 (3) f (x0 )h3 + and therefore 3! (19) f (x0 + h) f (x0 h) = f (x0 ) + c2 h2 + c4 h4 + 2h Similarly if we add (17) and (18) we obtain f (x0 + h) + f (x0 h) = 2f (x0 ) + 2 f (x0 )h2 + 4! f (4) (x0 )h4 + and consequently (20) f (x0 ) = f (x0 + h) 2f (x0 ) + f (x0 h) + C2 h2 + 2 h The nite dierence (21) 2 Dh f (x0 ) = f (x0 + h) 2f (x0 ) + f (x0 h) h2 2 is thus a second order approximation to f (x0 ), i.e., f (x0 ) Dh f (x0 ) = 2 O(h ). 5 2.1 Richardson's extrapolation From (19) we know that, asymptotically 0 Dh f (x0 ) = f (x0 ) + c2 h2 + c4 h4 + (22) We can apply Richardson extrapolation once to obtain a fourth order approximation. Evaluating (22) at h/2 we get (23) 1 1 0 Dh/2 f (x0 ) = f (x0 ) + c2 h2 + c4 h4 + 4 16 and multiplying this equation by 4, subtracting (22) to the result and dividing by 3 we get (24) ext Dh f (x0 ) := 0 0 4Dh/2 f (x0 ) Dh f (x0 ) 3 = f (x0 ) + c4 h4 + ext The method Dh f (x0 ) has order of convergence 4 for about twice the amount 0 of work of that Dh f (x0 ). Round-o errors are still O( /h) and the minimum total error will be when O(h4 ) is O( /h), i.e. when h = 1/5 . The minimum ext error is thus O( 4/5 ) for Dh f (x0 ), about 1014 in double precision with h = O(103 ). 6 \f