Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 54

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 54)

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.CITED REFERENCES AND FURTHER READING: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 = abF − [c − (a + b + 1)z]F (5.14.2)5.14 Evaluation of Functions by Path Integration209Here prime denotes d/dz. One can see that the equation has regular singular pointsat z = 0, 1, and ∞. Since the desired solution is regular at z = 0, the values 1 and∞ will in general be branch points.

If we want 2 F1 to be a single valued function,we must have a branch cut connecting these two points. A conventional position forthis cut is along the positive real axis from 1 to ∞, though we may wish to keepopen the possibility of altering this choice for some applications.Our golden brick consists of a collection of routines for the integration of setsof ordinary differential equations, which we will develop in detail later, in Chapter16. For now, we need only a high-level, “black-box” routine that integrates sucha set from initial conditions at one value of a (real) independent variable to finalconditions at some other value of the independent variable, while automaticallyadjusting its internal stepsize to maintain some specified accuracy. That routine iscalled odeint and, in one particular invocation, calculates its individual steps witha sophisticated Bulirsch-Stoer technique.Suppose that we know values for F and its derivative F at some value z0 , andthat we want to find F at some other point z1 in the complex plane.

The straight-linepath connecting these two points is parametrized byz(s) = z0 + s(z1 − z0 )(5.14.3)with s a real parameter. The differential equation (5.14.2) can now be written asa set of two first-order equations,dF= (z1 − z0 )F dsdF abF − [c − (a + b + 1)z]F = (z1 − z0 )dsz(1 − z)(5.14.4)to be integrated from s = 0 to s = 1.

Here F and F are to be viewed as twoindependent complex variables. The fact that prime means d/dz can be ignored; itwill emerge as a consequence of the first equation in (5.14.4). Moreover, the real andimaginary parts of equation (5.14.4) define a set of four real differential equations,with independent variable s. The complex arithmetic on the right-hand side can beviewed as mere shorthand for how the four components are to be coupled. It isprecisely this point of view that gets passed to the routine odeint, since it knowsnothing of either complex functions or complex independent variables.It remains only to decide where to start, and what path to take in the complexplane, to get to an arbitrary point z.

This is where consideration of the function’ssingularities, and the adopted branch cut, enter. Figure 5.14.1 shows the strategythat we adopt. For |z| ≤ 1/2, the series in equation (5.14.1) will in general convergerapidly, and it makes sense to use it directly. Otherwise, we integrate along a straightline path from one of the starting points (±1/2, 0) or (0, ±1/2). The former choicesare natural for 0 < Re(z) < 1 and Re(z) < 0, respectively.

The latter choices areused for Re(z) > 1, above and below the branch cut; the purpose of starting awayfrom the real axis in these cases is to avoid passing too close to the singularity atz = 1 (see Figure 5.14.1). The location of the branch cut is defined by the fact thatour adopted strategy never integrates across the real axis for Re (z) > 1.An implementation of this algorithm is given in §6.12 as the routine hypgeo.210Chapter 5.Evaluation of FunctionsImuse power seriesbranch cut01ReFigure 5.14.1. Complex plane showing the singular points of the hypergeometric function, its branchcut, and some integration paths from the circle |z| = 1/2 (where the power series converges rapidly)to other points in the plane.A number of variants on the procedure described thus far are possible, and easyto program. If successively called values of z are close together (with identical valuesof a, b, and c), then you can save the state vector (F, F ) and the corresponding valueof z on each call, and use these as starting values for the next call.

The incrementalintegration may then take only one or two steps. Avoid integrating across the branchcut unintentionally: the function value will be “correct,” but not the one you want.Alternatively, you may wish to integrate to some position z by a dog-leg paththat does cross the real axis Re z > 1, as a means of moving the branch cut.

Характеристики

Список файлов книги