Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 49
Текст из файла (страница 49)
Are there no bettermethods?Yes, there are. All, however, involve exploration of the function’s behaviorover scales comparable to xc , plus some assumption of smoothness, or analyticity,so that the high-order terms in a Taylor expansion like equation (5.7.2) have somemeaning. Such methods also involve multiple evaluations of the function f, so theirincreased accuracy must be weighed against increased cost.The general idea of “Richardson’s deferred approach to the limit” is particularlyattractive.
For numerical integrals, that idea leads to so-called Romberg integration(for review, see §4.3). For derivatives, one seeks to extrapolate, to h → 0, the resultof finite-difference calculations with smaller and smaller finite values of h. By theuse of Neville’s algorithm (§3.1), one uses each new finite-difference calculation toproduce both an extrapolation of higher order, and also extrapolations of previous,lower, orders but with smaller scales h.
Ridders [2] has given a nice implementationof this idea; the following program, dfridr, is based on his algorithm, modified byan improved termination criterion. Input to the routine is a function f (called func),a position x, and a largest stepsize h (more analogous to what we have called xcabove than to what we have called h). Output is the returned value of the derivative,and an estimate of its error, err.#include <math.h>#include "nrutil.h"#define CON 1.4#define CON2 (CON*CON)#define BIG 1.0e30#define NTAB 10#define SAFE 2.0Stepsize is decreased by CON at each iteration.Sets maximum size of tableau.Return when error is SAFE worse than the best sofar.float dfridr(float (*func)(float), float x, float h, float *err)Returns the derivative of a function func at a point x by Ridders’ method of polynomialextrapolation.
The value h is input as an estimated initial stepsize; it need not be small, butrather should be an increment in x over which func changes substantially. An estimate of theerror in the derivative is returned as err.{int i,j;float errt,fac,hh,**a,ans;if (h == 0.0) nrerror("h must be nonzero in dfridr.");a=matrix(1,NTAB,1,NTAB);hh=h;a[1][1]=((*func)(x+hh)-(*func)(x-hh))/(2.0*hh);*err=BIG;for (i=2;i<=NTAB;i++) {Successive columns in the Neville tableau will go to smaller stepsizes and higher orders ofextrapolation.hh /= CON;a[1][i]=((*func)(x+hh)-(*func)(x-hh))/(2.0*hh);Try new, smaller stepfac=CON2;size.for (j=2;j<=i;j++) {Compute extrapolations of various orders, requiringa[j][i]=(a[j-1][i]*fac-a[j-1][i-1])/(fac-1.0);no new function evalfac=CON2*fac;uations.errt=FMAX(fabs(a[j][i]-a[j-1][i]),fabs(a[j][i]-a[j-1][i-1]));5.7 Numerical Derivatives189The error strategy is to compare each new extrapolation to one order lower, bothat the present stepsize and the previous one.if (errt <= *err) {If error is decreased, save the improved answer.*err=errt;ans=a[j][i];}}if (fabs(a[i][i]-a[i-1][i-1]) >= SAFE*(*err)) break;If higher order is worse by a significant factor SAFE, then quit early.}free_matrix(a,1,NTAB,1,NTAB);return ans;}In dfridr, the number of evaluations of func is typically 6 to 12, but is allowedto be as great as 2×NTAB.
As a function of input h, it is typical for the accuracyto get better as h is made larger, until a sudden point is reached where nonsensicalextrapolation produces early return with a large error. You should therefore choosea fairly large value for h, but monitor the returned value err, decreasing h if it isnot small. For functions whose characteristic x scale is of order unity, we typicallytake h to be a few tenths.Besides Ridders’ method, there are other possible techniques. If your functionis fairly smooth, and you know that you will want to evaluate its derivative manytimes at arbitrary points in some interval, then it makes sense to construct aChebyshev polynomial approximation to the function in that interval, and to evaluatethe derivative directly from the resulting Chebyshev coefficients.
This method isdescribed in §§5.8–5.9, following.Another technique applies when the function consists of data that is tabulatedat equally spaced intervals, and perhaps also noisy. One might then want, at eachpoint, to least-squares fit a polynomial of some degree M , using an additionalnumber nL of points to the left and some number nR of points to the right of eachdesired x value. The estimated derivative is then the derivative of the resultingfitted polynomial. A very efficient way to do this construction is via Savitzky-Golaysmoothing filters, which will be discussed later, in §14.8. There we will give aroutine for getting filter coefficients that not only construct the fitting polynomial but,in the accumulation of a single sum of data points times filter coefficients, evaluateit as well. In fact, the routine given, savgol, has an argument ld that determineswhich derivative of the fitted polynomial is evaluated.
For the first derivative, theappropriate setting is ld=1, and the value of the derivative is the accumulated sumdivided by the sampling interval h.CITED REFERENCES AND FURTHER READING:Dennis, J.E., and Schnabel, R.B. 1983, Numerical Methods for Unconstrained Optimization andNonlinear Equations (Englewood Cliffs, NJ: Prentice-Hall), §§5.4–5.6. [1]Ridders, C.J.F.
1982, Advances in Engineering Software, vol. 4, no. 2, pp. 75–76. [2]190Chapter 5.Evaluation of Functions5.8 Chebyshev ApproximationThe Chebyshev polynomial of degree n is denoted Tn (x), and is given bythe explicit formulaTn (x) = cos(n arccos x)(5.8.1)This may look trigonometric at first glance (and there is in fact a close relationbetween the Chebyshev polynomials and the discrete Fourier transform); however(5.8.1) can be combined with trigonometric identities to yield explicit expressionsfor Tn (x) (see Figure 5.8.1),T0 (x) = 1T1 (x) = xT2 (x) = 2x2 − 1T3 (x) = 4x3 − 3x(5.8.2)T4 (x) = 8x − 8x + 142···Tn+1 (x) = 2xTn (x) − Tn−1 (x)n ≥ 1.(There also exist inverse formulas for the powers of x in terms of the Tn ’s — seeequations 5.11.2-5.11.3.)The Chebyshev polynomials are orthogonal in the interval [−1, 1] over a weight(1 − x2 )−1/2 .
In particular,&1−1Ti (x)Tj (x)√dx =1 − x2.0π/2πi = ji=j= 0i=j=0(5.8.3)The polynomial Tn (x) has n zeros in the interval [−1, 1], and they are locatedat the pointsx = cosπ(k − 12 )nk = 1, 2, . . . , n(5.8.4)In this same interval there are n + 1 extrema (maxima and minima), located atx = cosπknk = 0, 1, . . . , n(5.8.5)At all of the maxima Tn (x) = 1, while at all of the minima Tn (x) = −1;it is precisely this property that makes the Chebyshev polynomials so useful inpolynomial approximation of functions.1915.8 Chebyshev Approximation1T0T1T2Chebyshev polynomials.5T30T6−.5T5T4−1−1−.8−.6−.4−.20x.2.4.6.81Figure 5.8.1. Chebyshev polynomials T0 (x) through T6 (x). Note that Tj has j roots in the interval(−1, 1) and that all the polynomials are bounded between ±1.The Chebyshev polynomials satisfy a discrete orthogonality relation as well asthe continuous one (5.8.3): If xk (k = 1, . .
. , m) are the m zeros of Tm (x) givenby (5.8.4), and if i, j < m, then.0Ti (xk )Tj (xk ) = m/2k=1mi = ji=j= 0i=j=0m(5.8.6)It is not too difficult to combine equations (5.8.1), (5.8.4), and (5.8.6) to provethe following theorem: If f(x) is an arbitrary function in the interval [−1, 1], andif N coefficients cj , j = 0, . . . , N − 1, are defined byN2 f(xk )Tj (xk )Nk=1Nπj(k − 12 )π(k − 12 )2 cos=f cosNNNcj =(5.8.7)k=1then the approximation formulaf(x) ≈ N−1k=0!1ck Tk (x) − c02(5.8.8)192Chapter 5.Evaluation of Functionsis exact for x equal to all of the N zeros of TN (x).For a fixed N , equation (5.8.8) is a polynomial in x which approximates thefunction f(x) in the interval [−1, 1] (where all the zeros of TN (x) are located).
Whyis this particular approximating polynomial better than any other one, exact on someother set of N points? The answer is not that (5.8.8) is necessarily more accuratethan some other approximating polynomial of the same order N (for some specifieddefinition of “accurate”), but rather that (5.8.8) can be truncated to a polynomial oflower degree m N in a very graceful way, one that does yield the “most accurate”approximation of degree m (in a sense that can be made precise). Suppose N isso large that (5.8.8) is virtually a perfect approximation of f(x).
Now considerthe truncated approximation! m−11ck Tk (x) − c0(5.8.9)f(x) ≈2k=0with the same cj ’s, computed from (5.8.7). Since the Tk (x)’s are all boundedbetween ±1, the difference between (5.8.9) and (5.8.8) can be no larger than thesum of the neglected ck ’s (k = m, . . . , N − 1). In fact, if the ck ’s are rapidlydecreasing (which is the typical case), then the error is dominated by cm Tm (x),an oscillatory function with m + 1 equal extrema distributed smoothly over theinterval [−1, 1].
This smooth spreading out of the error is a very important property:The Chebyshev approximation (5.8.9) is very nearly the same polynomial as thatholy grail of approximating polynomials the minimax polynomial, which (among allpolynomials of the same degree) has the smallest maximum deviation from the truefunction f(x).
The minimax polynomial is very difficult to find; the Chebyshevapproximating polynomial is almost identical and is very easy to compute!So, given some (perhaps difficult) means of computing the function f(x), wenow need algorithms for implementing (5.8.7) and (after inspection of the resultingck ’s and choice of a truncating value m) evaluating (5.8.9). The latter equation thenbecomes an easy way of computing f(x) for all subsequent time.The first of these tasks is straightforward. A generalization of equation (5.8.7)that is here implemented is to allow the range of approximation to be between twoarbitrary limits a and b, instead of just −1 to 1.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
