NumericalDifferentiation andIntegration

1 CHAPTER 11 Numerical Differentiationand IntegrationDifferentiation and integration are basic mathematical operations with a widerange of applications in many areas of science. It is therefore important to havegood methods to compute and manipulate derivatives and integrals. You proba-bly learnt the basic rules of differentiation and integration in school symbolicmethods suitable for pencil-and-paper calculations. These are important, andmost derivatives can be computed this way. Integration however, is different,and most integrals cannot be determined with symbolic methods like the onesyou learnt in complication is the fact that in practical applications a function isonly known at a few points. For example, we may measure the position of acar every minute via a GPS (Global Positioning System) unit, and we want tocompute its speed.

2 If the position is known as a continuous function of time, wecan find the speed by differentiating this function. But when the position is onlyknown at isolated times, this is not possible. The same applies to solution, both when it comes to integrals that cannot be determinedby the usual methods, and functions that are only known at isolated points, isto use approximate methods of differentiation and integration. In our context,these are going to be numerical methods. We are going to present a number ofmethods for doing numerical integration and differentiation, but more impor-tantly, we are going to present a general strategy for deriving such methods. Inthis way you will not only have a number of methods available to you, but youwill also be able to develop new methods, tailored to special situations that youmay use the same general strategy for deriving both numerical integrationand numerical differentiation methods.

3 The basic idea is to evaluate a func-tion at a few points, find the polynomial that interpolates the function at thesepoints, and use the derivative or integral of the polynomial as an approximationto the function. This technique also allows us to keep track of the so-calledtrun-cation error, the mathematical error committed by integrating or differentiatingthe polynomial instead of the function itself. However, when it comes to round-off error, we have to treat differentiation and integration differently: Numericalintegration is veryinsensitiveto round-off errors, while numerical differentia-tion behaves in the opposite way; it is verysensitiveto round-off simple method for numerical differentiationWe start by studying numerical differentiation. We first introduce the simplestmethod, derive its error, and its sensitivity to round-off errors.

4 The procedureused here for deriving the method and analysing the error is used over again inlater sections to derive and analyse additional us first make it clear what numerical differentiation (Numerical differentiation).Letfbe a given function that isonly known at a number of isolated points. The problem of numerical differ-entiation is to compute an approximation to the derivativef offby suitablecombinations of the known values typical example is thatfis given by a computer program (more specifi-cally a function, procedure or method, depending on you choice of program-ming language), and you can call the program with a floating-point argumentxand receive back a floating-point approximation off(x). The challenge is tocompute an approximation tof (a) for some real numberawhen the only aidwe have at our disposal is the program to compute values basic ideaSince we are going to compute derivatives, we must be clear about they are de-fined.

5 Recall thatf (a) is defined byf (a)=limh 0f(a+h) f(a)h.( )In the following we will assume that this limit exists; , thatfis ( ) we immediately have a natural approximation tof (a); we simply228pick a positivehand use the approximationf (a) f(a+h) f(a)h.( )Note that this corresponds to approximatingfby the straight linep1that inter-polatesfataanda h, and then usingp 1(a) as an approximation tof (a).Observation derivative offatacan be approximated byf (a) f(a+h) f(a) a practical situation, the numberawould be given, and we would have tolocate the two nearest valuesa1anda2to the left and right ofasuch thatf(a1)andf(a2) can be found. Then we would use the approximationf (a) f(a2) f(a1)a2 later sections, we will derive several formulas like ( ). Which formulato use for a specific example, and exactly how to use it, will have to be decidedin each us test the approximation ( ) for the functionf(x)=sinxata= (using 64-bit floating-point numbers).

6 In this case we havef (x)=cosxsof (a)= This makes it is easy to check the accuracy. We trywith a few values ofhand findh(f(a+h) f(a))/h E1(f;a,h)10 10 210 10 310 10 410 10 510 10 610 10 7whereE1(f;a,h)=f(a) (f(a+h) f(a))/h. In other words, the approximationseems to improve with decreasingh, as expected. More precisely, whenhisreduced by a factor of 10, the error is reduced by the same truncation errorWhenever we use approximations, it is important to try and keep track of theerror, if at all possible. To analyse the error in numerical differentiation, Taylorpolynomials with remainders are useful. To analyse the error in the approxima-tion above, we do a Taylor expansion off(a+h). We havef(a+h)=f(a)+h f (a)+h22f ( h),where hlies in the interval (a,a+h). If we rearrange this formula , we obtainf (a) f(a+h) f(a)h= h2f ( h).

7 ( )This is often referred to as thetruncation errorof the approximation, and is areasonable error formula , but it would be nice to get rid of h. We first takeabsolute values in ( ), f (a) f(a+h) f(a)h =h2 f ( h) .Recall from the Extreme value theorem that if a function is continuous, then itsmaximum always exists on any closed and bounded interval. In our setting here,it is natural to let the closed and bounded interval be [a,a+h]. This leads to thefollowing thatfhas continuous derivatives up to order two neara. If the derivativef (a)is approximated byf(a+h) f(a)h,then thetruncation erroris bounded byE(f;a,h)= f (a) f(a+h) f(a)h h2maxx [a,a+h] f (x) .( )Let us check that the error formula ( ) confirms the numerical values inexample We havef (x)= sinx, so the right-hand side in ( ) becomesE(sin; ,h)=h2sin h,230where h ( , +h).

8 Forh= we therefore have that the error must lie inthe interval[ sin , sin ]=[ 10 2, 10 2],and the right end of the interval is the maximum value of the right-hand sidein ( ). Whenhis reduced by a factor of 10, the factorh/2 is reduced by thesame factor. This means that hwill approach so sin hwill approach thelower value sin Forh=10 n, the error will therefore tend to10 nsin 10 , which is in complete agreement with example is true in general. Iff is continuous, then hwill approachawhenhgoes to zero. But even whenh>0, the error in using the approximationf ( h) f (a) is small. This is the case since it is usually only necessary to know themagnitude of the error, , it is sufficient to know the error with one or twocorrect truncation error is approximately given by f (a) f(a+h) f(a)h h2 f (a).

9 Round-off errorSo far, we have just considered the mathematical error committed whenf (a)is approximated by(f(a+h) f(a))/h. But what about the round-off error? Infact, when we compute this approximation we have to perform the one criticaloperationf(a+h) f(a) subtraction of two almost equal numbers whichwe know from chapter 5 may lead to large round-off errors. Let us continueexample and see what happens if we use smaller values that we estimated the derivative off(x)=sinxata= that the correct value with ten digits isf ( ) If we checkvalues ofhfrom 10 7and smaller we findh(f(a+h) f(a))/h E(f;a,h)10 70. 10 810 10 1010 10 1010 10 610 10 410 10 210 10 110 10 1231 This shows very clearly that something quite dramatic happens, and when wecome toh=10 17, the derivative is computed as (a) is the floating-point number closest tof(a), we know from lemma the relative error will be bounded by 5 2 53since floating-point numbersare represented in binary ( =2) with 53 bits for the significand (m=53).

10 Wetherefore have| | 5 2 53 6 10 16. In practice, the real upper bound on isusually smaller, and in the following we will denote this upper bound by . Thismeans that a definite upper bound on is 6 10 maximum relative error when a real number is repre-sented by a floating-point number is denoted by .There is a handy way to express the relative error inf(a). If we denote thecomputed value off(a) byf(a), we will havef(a)=f(a)(1+ )which corresponds to the relative error being| |.Observation thatf(a)is computed with 64-bit floating-pointnumbers and that no underflow or overflow occurs. Then the computed valuef(a)satisfiesf(a)=f(a)(1+ )( )where| | , and depends on computation off(a+h) is of course also affected by round-off error, sowe havef(a)=f(a)(1+ 1),f(a+h)=f(a+h)(1+ 2)( )where| i| fori=1, 2.