Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 53
Текст из файла (страница 53)
(We have omitted some technicalassumptions in this theorem. See Ralston [1] for a precise statement.) We thus learn that thesituation with rational functions is quite analogous to that for minimax polynomials: In §5.8we saw that the error term of an nth order approximation, with n + 1 Chebyshev coefficients,was generally dominated by the first neglected Chebyshev term, namely Tn+1 , which itselfhas n + 2 extrema of equal magnitude and alternating sign.
So, here, the number of rationalcoefficients, m + k + 1, plays the same role of the number of polynomial coefficients, n + 1.A different way to see why r(x) should have m + k + 2 extrema is to note that R(x)can be made exactly equal to f (x) at any m + k + 1 points xi . Multiplying equation (5.13.1)by its denominator gives the equationskp0 + p1 xi + · · · + pm xmi = f (xi )(1 + q1 xi + · · · + qk xi )i = 1, 2, . . .
, m + k + 1(5.13.3)This is a set of m + k + 1 linear equations for the unknown p’s and q’s, which can besolved by standard methods (e.g., LU decomposition). If we choose the xi ’s to all be inthe interval (a, b), then there will generically be an extremum between each chosen xi andxi+1 , plus also extrema where the function goes out of the interval at a and b, for a totalof m + k + 2 extrema. For arbitrary xi ’s, the extrema will not have the same magnitude.The theorem says that, for one particular choice of xi ’s, the magnitudes can be beaten downto the identical, minimal, value of r.Instead of making f (xi ) and R(xi ) equal at the points xi , one can instead force theresidual r(xi ) to any desired values yi by solving the linear equationskp0 + p1 xi + · · · + pm xmi = [f (xi ) − yi ](1 + q1 xi + · · · + qk xi )i = 1, 2, . .
. , m + k + 1(5.13.4)5.13 Rational Chebyshev Approximation205In fact, if the xi ’s are chosen to be the extrema (not the zeros) of the minimax solution,then the equations satisfied will bekp0 + p1 xi + · · · + pm xmi = [f (xi ) ± r](1 + q1 xi + · · · + qk xi )i = 1, 2, . . .
, m + k + 2(5.13.5)where the ± alternates for the alternating extrema. Notice that equation (5.13.5) is satisfied atm + k + 2 extrema, while equation (5.13.4) was satisfied only at m + k + 1 arbitrary points.How can this be? The answer is that r in equation (5.13.5) is an additional unknown, so thatthe number of both equations and unknowns is m + k + 2. True, the set is mildly nonlinear(in r), but in general it is still perfectly soluble by methods that we will develop in Chapter 9.We thus see that, given only the locations of the extrema of the minimax rationalfunction, we can solve for its coefficients and maximum deviation. Additional theorems,leading up to the so-called Remes algorithms [1], tell how to converge to these locations byan iterative process.
For example, here is a (slightly simplified) statement of Remes’ SecondAlgorithm: (1) Find an initial rational function with m + k + 2 extrema xi (not having equaldeviation). (2) Solve equation (5.13.5) for new rational coefficients and r. (3) Evaluate theresulting R(x) to find its actual extrema (which will not be the same as the guessed values).(4) Replace each guessed value with the nearest actual extremum of the same sign. (5) Goback to step 2 and iterate to convergence.
Under a broad set of assumptions, this method willconverge. Ralston [1] fills in the necessary details, including how to find the initial set of xi ’s.Up to this point, our discussion has been textbook-standard. We now reveal ourselvesas heretics. We don’t much like the elegant Remes algorithm. Its two nested iterations (onr in the nonlinear set 5.13.5, and on the new sets of xi ’s) are finicky and require a lot ofspecial logic for degenerate cases. Even more heretical, we doubt that compulsive searchingfor the exactly best, equal deviation, approximation is worth the effort — except perhaps forthose few people in the world whose business it is to find optimal approximations that getbuilt into compilers and microchips.When we use rational function approximation, the goal is usually much more pragmatic:Inside some inner loop we are evaluating some function a zillion times, and we want tospeed up its evaluation.
Almost never do we need this function to the last bit of machineaccuracy. Suppose (heresy!) we use an approximation whose error has m + k + 2 extremawhose deviations differ by a factor of 2. The theorems on which the Remes algorithmsare based guarantee that the perfect minimax solution will have extrema somewhere withinthis factor of 2 range – forcing down the higher extrema will cause the lower ones to rise,until all are equal.
So our “sloppy” approximation is in fact within a fraction of a leastsignificant bit of the minimax one.That is good enough for us, especially when we have available a very robust methodfor finding the so-called “sloppy” approximation. Such a method is the least-squares solutionof overdetermined linear equations by singular value decomposition (§2.6 and §15.4). Weproceed as follows: First, solve (in the least-squares sense) equation (5.13.3), not just form + k + 1 values of xi , but for a significantly larger number of xi ’s, spaced approximatelylike the zeros of a high-order Chebyshev polynomial. This gives an initial guess for R(x).Second, tabulate the resulting deviations, find the mean absolute deviation, call it r, and thensolve (again in the least-squares sense) equation (5.13.5) with r fixed and the ± chosen to bethe sign of the observed deviation at each point xi .
Third, repeat the second step a few times.You can spot some Remes orthodoxy lurking in our algorithm: The equations we solveare trying to bring the deviations not to zero, but rather to plus-or-minus some consistentvalue. However, we dispense with keeping track of actual extrema; and we solve only linearequations at each stage. One additional trick is to solve a weighted least-squares problem,where the weights are chosen to beat down the largest deviations fastest.Here is a program implementing these ideas. Notice that the only calls to the function fnoccur in the initial filling of the table fs.
You could easily modify the code to do this fillingoutside of the routine. It is not even necessary that your abscissas xs be exactly the onesthat we use, though the quality of the fit will deteriorate if you do not have several abscissasbetween each extremum of the (underlying) minimax solution. Notice that the rationalcoefficients are output in a format suitable for evaluation by the routine ratval in §5.3.206Chapter 5.Evaluation of Functionsm=k=4f(x) = cos(x)/(1 + e x )0<x<π2 × 10 − 6R(x) − f (x)10 − 60− 1 × 10 − 6− 2 × 10 − 60.511.522.53xFigure 5.13.1.
Solid curves show deviations r(x) for five successive iterations of the routine ratlsqfor an arbitrary test problem. The algorithm does not converge to exactly the minimax solution (shownas the dotted curve). But, after one iteration, the discrepancy is a small fraction of the last significantbit of accuracy.#include <stdio.h>#include <math.h>#include "nrutil.h"#define NPFAC 8#define MAXIT 5#define PIO2 (3.141592653589793/2.0)#define BIG 1.0e30void ratlsq(double (*fn)(double), double a, double b, int mm, int kk,double cof[], double *dev)Returns in cof[0..mm+kk] the coefficients of a rational function approximation to the functionfn in the interval (a, b).
Input quantities mm and kk specify the order of the numerator anddenominator, respectively. The maximum absolute deviation of the approximation (insofar asis known) is returned as dev.{double ratval(double x, double cof[], int mm, int kk);void dsvbksb(double **u, double w[], double **v, int m, int n, double b[],double x[]);void dsvdcmp(double **a, int m, int n, double w[], double **v);These are double versions of svdcmp, svbksb.int i,it,j,ncof,npt;double devmax,e,hth,power,sum,*bb,*coff,*ee,*fs,**u,**v,*w,*wt,*xs;ncof=mm+kk+1;npt=NPFAC*ncof;bb=dvector(1,npt);coff=dvector(0,ncof-1);ee=dvector(1,npt);fs=dvector(1,npt);u=dmatrix(1,npt,1,ncof);v=dmatrix(1,ncof,1,ncof);w=dvector(1,ncof);wt=dvector(1,npt);Number of points where function is evaluated,i.e., fineness of the mesh.5.13 Rational Chebyshev Approximation207xs=dvector(1,npt);*dev=BIG;for (i=1;i<=npt;i++) {Fill arrays with mesh abscissas and function valif (i < npt/2) {ues.hth=PIO2*(i-1)/(npt-1.0);At each end, use formula that minimizes roundxs[i]=a+(b-a)*DSQR(sin(hth));off sensitivity.} else {hth=PIO2*(npt-i)/(npt-1.0);xs[i]=b-(b-a)*DSQR(sin(hth));}fs[i]=(*fn)(xs[i]);wt[i]=1.0;In later iterations we will adjust these weights toee[i]=1.0;combat the largest deviations.}e=0.0;for (it=1;it<=MAXIT;it++) {Loop over iterations.for (i=1;i<=npt;i++) {Set up the “design matrix” for the least-squarespower=wt[i];fit.bb[i]=power*(fs[i]+SIGN(e,ee[i]));Key idea here: Fit to fn(x) + e where the deviation is positive, to fn(x) − e whereit is negative.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
