Thompson - Computing for Scientists and Engineers (523188), страница 25
Текст из файла (страница 25)
Here we could estimate afirst central derivative at x + h (in the figure at x = 0.1 by the slope of C+C+),then a similar derivative at x - h (at x = -0.1 by the slope of line C-C-).FIGURE 4.5 Schemes for numerical second derivatives, illustrated for the working function(4.1) and Figure 4.1 near x = 0. Consecutive central derivatives (CCD) are obtained from centraldifference first derivatives along C+C+ and C-C-. The 3CD estimate (4.35) and the 5CD estimate(4.37) use the central 3 points and all 5 points, respectively.4.4HOW TO APPROXIMATE DERIVATIVES127This method, consecutive central derivatives, or CCD, can be analyzed algebraically, as follows.
First, write the formula (4.26) for central derivatives of the function which is already the first derivative, thus(4.28)Next, use the central-difference estimates for the first derivative, (4.26), in this expression, to obtain(4.29)For optimum accuracy, minimizing subtractive cancellation, it is better to computethe second-derivative estimate this way rather than from further algebraic simplification of the expression.Exercise 4.13(a) For the working function (4.2) calculate its central-difference derivative estimate (4.29) at x = 0 for the values of h in Table 4.6. In Program 4.2, Working Function, this second-derivative estimate is variable yCCD from the function Der_2CCD.(b) Write down the first few terms of the Taylor expansion of a functiony (x + nh) with n = -2, 0, and 2, in terms of the function y (x) and its derivatives at x.
Thus show that the first additional term beyond the second derivativeitself in the central-difference derivative estimate (4.29) is y(4) (x) h2 / 2 .(c) Show that this neglected term is only a fair estimate of the error in the numerical derivatives in Table 4.6, except for the smallest h value. (You can obtain the value of the fourth derivative at x = 0 from Table 4.3.) For this example and values of h, why is the estimate of the error inaccurate?(d) Use the program Working Function (or a program of your own devising)to show that if one wants an accuracy of about 6 significant figures in the secondderivative at x = 0, then the stepsize required is about h = 10-4. Show thatthis agrees with the error estimate derived in (b). nTABLE 4.6 Central-difference second derivatives for the working function, evaluated at x = 0.The exact second derivative is 6.128NUMERICAL DERIVATIVES AND INTEGRALSCompare Table 4.6 for the second derivative with Table 4.5 for the central-difference first derivative, to note that about the same value of h is required for thesame accuracy in first and second derivatives.
This is not a usual behavior. By inspecting the derivative values in Table 4.3 you can see that, although the increasingpowers of h (which is < 1) bring down the error estimates, the steadily increasingderivatives increase them. Our choice of the wiggly sixth-order polynomial, definedby (4.1) and shown in Figures 4.1 and 4.2, is responsible for this pathological behavior. By contrast, the derivatives at x = 0 for the exponential function are allequal to unity, and the convergence is therefore much more rapid, as we discover inSection 4.5.Better algorithms for second derivativesThe algorithm for numerical second derivatives, (4.29), can be significantly improved without increased complexity, as we now show. The method that we use alsoshows how optimized formulas can be produced and how their errors can be estimated from Taylor series. We derive two formulas, one using three points and theother using five points, all having the x value of interest as the central value.Suppose that we go n steps, each of size h, to the left or right of the point x atwhich we want a numerical estimate of the second derivative.
Let us use Taylor’stheorem (3.6) for up to five points surrounding x, as in Figure 4.5. We then have(4.30)(4.3 1)where the coefficients(4.32)contain the desired (and undesired) derivatives. For second derivatives, we need toeliminate the term v1 and as many higher derivatives as possible. All odd-order derivatives (all vi with i odd) can be eliminated by summing values that are symmetricabout the midpoint. Thus(4.33)(4.34)From these two formulas, depending on how many points are available, we obtainvarious estimates for the second derivative that is hiding in v2 .4.4HOW TO APPROXIMATE DERIVATIVES129Exercise 4.14By using (4.33) show that an estimate of the second derivative that uses threepoints is y 3CD, given by(4.35)with an error,3CD,estimated as(4.36)which is the first omitted term of the Taylor expansion.
nNotice that this formula is predicted to be a factor of 6 more accurate than our CCDformula (4.29).If we may use all five points centered on x, as in Figure 4.5, then the estimateof the second derivative can be significantly improved.Exercise 4.15Use both (4.33) and (4.34) to eliminate the fourth-derivative term between them,and thence to show that an estimate of the second derivative that uses five pointsis y5CD, given by(4.37)with an error,5CD,estimated as(4.38)which is the first omitted term of the Taylor expansion. nFormula (4.37) predicts many improvements over the two other estimates of the second derivative, in that the error depends upon the sixth derivative, which is usuallysmaller than the fourth derivative (except for our pathological working function!),and scales as h 4 rather than h,2 so it will become relatively much more accurate as his decreased.Now that we have done much intricate algebra, let us try out these secondderivative formulas on our working function.
The example is again for the derivative at x = 0.130NUMERICAL DERIVATIVES AND INTEGRALSExercise 4.16(a) For the working function (4.1) calculate its central-difference derivative estimate (4.29) at x = 0 for the values of h in Table 4.7. In Program 4.2, Working Function, this second-derivative estimate is variable yCCD from the function Der_2CCD.(b) Show that if you want the second derivative of the working function (4.1) atx = 0 accurate to about 6 significant figures, then you have to chooseh 2 x 1O-4 for the three-point second derivative, but only h 1 x 10-2 forthe five-point derivative. nTABLE 4.7 Central-difference second derivatives and error estimates for the working functionevaluated at x = 0.
The exact second derivative is 6.We now have considerable experience with estimating first and second derivatives numerically for the very badly behaved working function (Figures 4.1 and4.2). In the next section we explore the derivatives of some functions that are betterbehaved.4.5 PROJECT 4A: COMPUTING DERIVATIVES NUMERICALLYAlthough in Section 4.4 we acquired considerable experience with numerical firstand second derivatives of the badly behaved working function (Section 4.1), it isnow interesting to explore the numerics with functions that are more regular. In thisproject we explore two common functions that are paradigms for many other functions and data, the exponential function and the cosine function. These characterizesteadily increasing (or decreasing, depending on your point of view) and oscillatoryfunctions, respectively.Derivatives of the exponential functionThe exponential function y = exis straightforward to work with, especially atthe origin x = 0, because there its derivatives have value unity.4.5 COMPUTING DERIVATIVES NUMERICALLY131Exercise 4.17(a) Modify the program Working Function by replacing the sections forpolynomial evaluation (including the Horner-polynomial function) by the exponential function, y = eX.
Keep the five functions for estimating derivatives.(Since all the derivatives of the exponential are just the function itself, a singleuse of the function exp (x) will produce all the derivatives you need.) Addprogram statements to compute the difference between these derivatives andthose obtained from the five formulas for first and second derivatives that are already programmed. Test this new program, called Exponential FunctionNumerical Derivatives.(b) Run this program with x = 0 and successively smaller values of stepsize h,say h = 0.2, 0.1, 0.01, 0.001, calculating all five numerical derivatives andtheir errors for each h. Make a log-log plot of the absolute value of the errorversus h, as in Figure 4.6.(c) To the extent that the graph for the error in each method versus stepsize h isa straight line on a log-log plot, argue that the error has a power-law dependenceon stepsize h, and that the overall scale of the error estimates the factor precedingthe power of h.
Do the powers estimated for the five derivative methods fromyour graph agree with those predicted from our Taylor series analyses? nFIGURE 4.6 Errors in derivatives of the exponential function at x = 0 as a function of stepsize h for various derivative estimators.
First-derivative estimators are forward (F) from (4.25) andcentral (C) from (4.26). Second-derivative estimators are consecutive central derivatives (CCD)from (4.29), three-point derivatives (3CD) from (4.35), and five-point derivatives (5CD) from(4.37).132NUMERICAL DERIVATIVES AND INTEGRALSNow that you have a working program for estimating numerical derivatives, byreplacing the exponential function in Exponential Function Numerical Derivatives by the function of your choice, you can explore its application to variousother functions, such as the cosine function in the next subsection. Also the program may be used directly with data (if they are equally spaced) to estimate slopesand higher-order derivatives.Differentiating the cosine functionThe cosine function, y = cos X, is interesting to investigate because it characterizesoscillatory behavior, such as we obtain from the numerical solution of differentialequations in Chapters 7 and 8, as well as from the Fourier expansions in Chapters 9 and 10.
Indeed, an understanding of numerical derivatives provides an essential foundation for developing numerical solutions of differential equations.FIGURE 4.7 Errors in derivatives of the cosine function at x = 0 as a function of stepsize hfor the first-derivative forward (F) from (4.25), and for second derivatives - consecutive centralderivatives (CCD) from (4.29) three-point derivatives (3CD) from (4.35) and five-point derivatives(5CD) from (4.37). The errors in the central-difference first-derivative estimate (4.26) are zero.The procedure for the second part of Project 4A is similar to that in the first part,and it can be done almost effortlessly.Exercise 4.18(a) Modify either the program Working Function or (more readily) the program Exponential Function by writing sections for evaluating the cosinefunction, y = cos x, for its value, for its first derivative -sin x, and for its second derivative -cos x. Keep the five functions for estimating derivatives.