c5-13 (779497), страница 2
Текст из файла (страница 2)
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)0− 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.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).10 − 65.13 Rational Chebyshev Approximation207}Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).xs=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.
Then e is supposed to become an approximation to the equal-rippledeviation.for (j=1;j<=mm+1;j++) {u[i][j]=power;power *= xs[i];}power = -bb[i];for (j=mm+2;j<=ncof;j++) {power *= xs[i];u[i][j]=power;}}dsvdcmp(u,npt,ncof,w,v);Singular Value Decomposition.In especially singular or difficult cases, one might here edit the singular values w[1..ncof],replacing small values by zero. Note that dsvbksb works with one-based arrays, so wemust subtract 1 when we pass it the zero-based array coff.dsvbksb(u,w,v,npt,ncof,bb,coff-1);devmax=sum=0.0;for (j=1;j<=npt;j++) {Tabulate the deviations and revise the weights.ee[j]=ratval(xs[j],coff,mm,kk)-fs[j];wt[j]=fabs(ee[j]);Use weighting to emphasize most deviant points.sum += wt[j];if (wt[j] > devmax) devmax=wt[j];}e=sum/npt;Update e to be the mean absolute deviation.if (devmax <= *dev) {Save only the best coefficient set found.for (j=0;j<ncof;j++) cof[j]=coff[j];*dev=devmax;}printf(" ratlsq iteration= %2d max error= %10.3e\n",it,devmax);}free_dvector(xs,1,npt);free_dvector(wt,1,npt);free_dvector(w,1,ncof);free_dmatrix(v,1,ncof,1,ncof);free_dmatrix(u,1,npt,1,ncof);free_dvector(fs,1,npt);free_dvector(ee,1,npt);free_dvector(coff,0,ncof-1);free_dvector(bb,1,npt);208Chapter 5.Evaluation of FunctionsFigure 5.13.1 shows the discrepancies for the first five iterations of ratlsq when it isapplied to find the m = k = 4 rational fit to the function f (x) = cos x/(1 + ex ) in theinterval (0, π).
One sees that after the first iteration, the results are virtually as good as theminimax solution. The iterations do not converge in the order that the figure suggests: Infact, it is the second iteration that is best (has smallest maximum deviation). The routineratlsq accordingly returns the best of its iterations, not necessarily the last one; there is noadvantage in doing more than five iterations.Ralston, A.
and Wilf, H.S. 1960, Mathematical Methods for Digital Computers (New York: Wiley),Chapter 13. [1]5.14 Evaluation of Functions by PathIntegrationIn computer programming, the technique of choice is not necessarily the mostefficient, or elegant, or fastest executing one. Instead, it may be the one that is quickto implement, general, and easy to check.One sometimes needs only a few, or a few thousand, evaluations of a specialfunction, perhaps a complex valued function of a complex variable, that has manydifferent parameters, or asymptotic regimes, or both. Use of the usual tricks (series,continued fractions, rational function approximations, recurrence relations, and soforth) may result in a patchwork program with tests and branches to differentformulas.
While such a program may be highly efficient in execution, it is often notthe shortest way to the answer from a standing start.A different technique of considerable generality is direct integration of afunction’s defining differential equation – an ab initio integration for each desiredfunction value — along a path in the complex plane if necessary. While this may atfirst seem like swatting a fly with a golden brick, it turns out that when you alreadyhave the brick, and the fly is asleep right under it, all you have to do is let it fall!As a specific example, let us consider the complex hypergeometric function 2 F1 (a, b, c; z), which is defined as the analytic continuation of the so-calledhypergeometric series,a(a + 1)b(b + 1) z 2ab z++···c 1!c(c + 1)2!a(a + 1) . .
. (a + j − 1)b(b + 1) . . . (b + j − 1) z j+···+c(c + 1) . . . (c + j − 1)j!(5.14.1)The series converges only within the unit circle |z| < 1 (see [1]), but one’s interestin the function is often not confined to this region.The hypergeometric function 2 F1 is a solution (in fact the solution that is regularat the origin) of the hypergeometric differential equation, which we can write as2 F1 (a, b, c; z)=1+z(1 − z)F 00 = abF − [c − (a + b + 1)z]F 0(5.14.2)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).CITED REFERENCES AND FURTHER READING:.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
