Numerical Differentiation - Karen A. Kopecky

1 LECTURENOTESECO613/614 FALL2007 KARENA. KOPECKYN umerical DifferentiationThe simplest way to compute a function s derivatives numerically is to use finite differ-ence approximations. Suppose we are interested in computing the first and second deriva-tives of a smooth functionf:R R. The definition of a derivative,f (x) = limh 0f(x+h) f(x)h,suggests a natural approximation. Take a small numberh, (more on how small latter) andf (x) f(x+h) f(x)h.(1)This is the easiest and most intuitive finite difference formula and it is called the for-ward difference. The forward difference is the most widely used way to compute numericalderivatives but often it is not the best choice as we will see. In order to compare to alterna-tive approximations we need to derive an error bound for the forward difference. This canbe done by taking a Taylor expansion off,f(x+h) =f(x) +hf (x) +h22f (x) +h36f (x) + (2)A little manipulation yieldsf (x) =f(x+h) f(x)h+O(h).

2 (3)We say that this approximation is first-order accurate since the dominate term in the trun-cation error isO(h).An alternative formula to the forward difference is to use a two-sided difference or centerdifference. The center difference formula can be derived by taking the two second-orderTaylor seriesf(x+h) =f(x) +hf (x) +h22f (x) +O(h3),(4)1andf(x h) =f(x) hf (x) +h22f (x) +O(h3),(5)and subtracting the second series from the the first and dividing by2h, obtainingf (x) =f(x+h) f(x h)2h+O(h2).(6)Notice that the center difference approximation is second-order accurate since the dom-inate term in its truncation error isO(h2). Thus the center difference is more accurate thanthe forward difference due to its smaller truncation error. For one-dimensional case it isalso just as costly to compute (they both require two function evaluations) and thereforeshould almost always be chosen over the forward difference approximation.

3 Unfortunatelythe usual tradeoff between accuracy and efficiency comes into play though as we increasethe dimensionality off. Iff:Rn Rmthen computing the Jacobian matrix using theforward difference requiresm(n+1)function evaluations, while using the center differencerequires2mnfunction evaluations. So center differences take approximately twice as longto compute whennis approximations can also be derived. Suppose we take the pointsh,x+h,andx+ hand write our approximate tof asf (x) af(x) +bf(x+h) +cf(x+ h). using the Taylor expansions off(x+h)andf(x+ h)aroundxand applying the meanvalue theorem, there existsx1 [x, x+h]andx2 [x, x+ h]such thataf(x) +bf(x+h) +cf(x+ h) = (a+b+c)f(x) +h(b+c )f (x)+h22(b+c 2)f (x) +h36[bf (x1) +c 3f (x2)]We now can see how to choosea,b, andcto make the right-hand-side approximatelyf (x).They must satisfya+b+c= 0,b+c = 1/h,andb+c 2= yieldsa= 2 1 (1 )h,b= 2 (1 )h,andc=1 (1 ) note the in fact the center difference is the special case of the three point approxi-mation when = 1.

4 Second notice that the error from this approximation isO(h2). Notnecessarily an improvement over the center difference but more costly to compute (we nowhave to do three function evaluations instead of two). Why would we want to use it?One situation is when we want to approximate the derivative at the boundary of thedomain. In this situation the center difference is not an option since it requires evaluatingthe function outside of the domain. We could choose to use the forward or equivalent of theforward but withh <0(called the backward difference) at these endpoints. But what if wewant to obtain a more accurate approximation? In this case it may payoff to use the threepoint approximation. Often is set to2yieldingf (x) =12h[ 3f(x) + 4f(x+h) f(x+ 2h)] +O(h2).This will generate a second-order accurate approximate to the derivative at either endpointby settinghgreater than or less can find finite difference approximations for second derivatives and other higherorder derivatives using a similar approach.

5 For example, a centered finite difference ap-proximation to the second derivative can be derived by adding together the two Taylorseriesf(x+h) =f(x) +hf (x) +h22f (x) +h36f (x) +O(h4),(7)andf(x h) =f(x) hf (x) +h22f (x) h36f (x) +O(h4),(8)and solving forf (x). This yieldsf (x) =f(x+h) 2f(x) +f(x h)h2+O(h2).3In addition we can derive general second-order accurate approximations tof using weightedsums offevaluated at various points, only now we would need4points instead of3. Inan analogous way to thef case we can use these (usually more costly) general formulas tohandle special situations like approximating the second derivative at the boundary of to ChoosehNotice that the truncation error that arises from the finite difference approximations isincreasing inhsuggesting thathshould be as small as possible to reduce this error. Butalso notice that the computation of the finite differences involves the computer taking dif-ferences and sums of floating-point numbers and dividing the result by a small number,h.

6 This is equivalent to multiplying the result by a large number. The multiplication willmagnify any round-off errors in the numerator. The smallerhis the larger is this magni-fication. Hence in facthshould be chosen to minimize total error which includes both theerror that arises from truncation (and is decreasing inh) and the round-off error that arisesfrom doing floating-point arithmetic with floating-point numbers (and is increasing inh).Let us first look at the case of the forward difference. We want to derive an upperbound on the roundoff error incurred when computing the approximation. Let f(x)be thenumerical representation off(x), and assume that the only error arising in the computationoffis round-off error. Denote machine epsilon byuandLf= supx [x,x+h]|f(x)|so that| f(x+h) f(x+h)| uLf,(9)and| f(x) f(x)| uLf,(10)and the|fl( f(x+h) f(x)) ( f(x+h) f(x))| 2uLf,(11)Then the roundoff error is bounded by f(x+h) f(x)h f(x+h) f(x)h 4uLfh.

7 (12)4 Notice that ifhis too small relative tou, the error can easily be huge. If you pickh=uyourerror is4Lf. The total roundoff error increases as we decrease the we need an upper bound on the truncation error. We can show, using the meanvalue theorem, that there existsc [x, x+h]such thatf (x) =f(x+h) f(x)h f (c)2h.(13)LetM= supx [x,x+h]|f (x)|.Then the truncation error is bounded by f(x+h) f(x)h f (x) Mh2.(14)Notice that the truncation error decreases as we decreaseh. Smaller values ofhimproveour seen above, choosinghrequires a trade-off between more round-off error, on theone hand, and more truncation error, on the other. Thus the optimalhis the one whichminimizes the sum of these two errors. Combining the roundoff and truncation error weget that the total error incurred in the computation off (x) using the forward difference isbounded by f(x+h) f(x)h f (x) 4uLfh+Mh2(15)Minimizing this bound with respect tohleads to an optimal bandwidth ofh =(8 LfuM)1/2.

8 (16)In practice, we will usually not be able to computeLforMbut a good approximation is toset the optimal bandwidth toh? LfM u |x| uThis rule works well as long asxis not to small. Smallxwill lead us to bandwidth s thatare too small for floating point operations. Thus the general rule of thumb is to seth?= max(|x|,1) double precision computing, u 10 8. Finally we can compute the upper bound onthe total absolute error incurred in the computation,4uLfM1/2(8uLf)1/2+M(8uLf)1/2 2M1/2= 2 2 LfMuand the relative error is2 2 LfMu|f (x)| LfM|f (x)| u using the forward difference with the optimal bandwidth, you shouldn t expect yourderivatives to accurate to more than about 8 significant , instead of using the forward difference, we use the center difference formula we havea different optimal bandwidth. The derivation is identical to that for the forward a similar argument to before we can derive an upper bound on the roundoff error f(x+h) f(x h)2h f(x+h) f(x h)2h 4uLf2h.

9 (17)By using the center difference formula, we do not modify the order of the rounding should not be surprising since the motivation for the center difference formula is toreduce errors due to curvature whereas the rounding error is due to the computer we need an upper bound on the truncation error. Again we can show, using themean value theorem, that there existsc [x, x+h]such thatf (x) =f(x+h) f(x h)2h f (c)6h2.(18)LetS= supx [x h,x+h]|f (x)|.Then the truncation error is bounded by f(x+h) f(x h)2h f (x) Sh26.(19)The total error in the computation off (x)is now bounded by: f(x+h) f(x h)h f (x) 2uLfh+Sh26.(20)Minimizing with respect tohleads to the optimal bandwidth:h =(12uLfS)1/3.(21)6 Since, in practice, we cannot computeLforSnotice that a good approximation toh ish (LfS)1/3u1/3 |x|u1 in the case of the forward difference it is best to rule outhbecoming too small so therule of thumb is to seth = max(|x|,1)u1 double precision,u1/3is approximately10 we can compute the upper boundon the total error incurred in the computation,4(SL2f12)1/3u2/3(22)which is on the order of10 11in double precision.

10 So we can expect are center differenceapproximation to be accurate to approximately the is one further problem that can arise. If the floating point value ofx+hdoes notequalx+hbutx+h+ethen we are actually using the approximationf(x+h+e) f(x)h=f(x+h+e) f(x+h)eeh+f(x+h) f(x)hNotice that this is approximatelyf (x+h)eh+f (x) (1 +eh)f (x).If the rounding erroreis on the order of magnitude ofuandhis on the order of magnitudeof uwe have now introduced an relative error on the order of uinto the calculation. Toavoid this problem it is recommended that you definehas follows. For one-sided (forwardor backward) difference first define has h= u max(|x|,1)then setxhasxh=x+ hand finally seth=xh two-sided (or center-difference) definehash=u1/3max(|x|,1)then setxh1asxh1 =x+handxh0 =x hand finallyhh=xh1 the roundoff error. For the second derivative approximation, thesame reasoning can be used to derive a reasonable value forh.