Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 55
Текст из файла (страница 55)
Forexample, in some cases you might want to integrate from (0, 1/2) to (3/2, 1/2),and go from there to any point with Re z > 1 — with either sign of Im z. (Ifyou are, for example, finding roots of a function by an iterative method, you donot want the integration for nearby values to take different paths around a branchpoint. If it does, your root-finder will see discontinuous function values, and willlikely not converge correctly!)In any case, be aware that a loss of numerical accuracy can result if you integratethrough a region of large function value on your way to a final answer where thefunction value is small. (For the hypergeometric function, a particular case of this iswhen a and b are both large and positive, with c and x >∼ 1.) In such cases, you’llneed to find a better dog-leg path.The general technique of evaluating a function by integrating its differentialequation in the complex plane can also be applied to other special functions.
For5.14 Evaluation of Functions by Path Integration211example, the complex Bessel function, Airy function, Coulomb wave function, andWeber function are all special cases of the confluent hypergeometric function, with adifferential equation similar to the one used above (see, e.g., [1] §13.6, for a table ofspecial cases). The confluent hypergeometric function has no singularities at finite z:That makes it easy to integrate.
However, its essential singularity at infinity meansthat it can have, along some paths and for some parameters, highly oscillatory orexponentially decreasing behavior: That makes it hard to integrate. Some case bycase judgment (or experimentation) is therefore required.CITED REFERENCES AND FURTHER READING:Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York). [1]Chapter 6.Special Functions6.0 IntroductionThere is nothing particularly special about a special function, except thatsome person in authority or textbook writer (not the same thing!) has decided tobestow the moniker.
Special functions are sometimes called higher transcendentalfunctions (higher than what?) or functions of mathematical physics (but they occur inother fields also) or functions that satisfy certain frequently occurring second-orderdifferential equations (but not all special functions do). One might simply call them“useful functions” and let it go at that; it is surely only a matter of taste whichfunctions we have chosen to include in this chapter.Good commercially available program libraries, such as NAG or IMSL, containroutines for a number of special functions. These routines are intended for users whowill have no idea what goes on inside them. Such state of the art “black boxes” areoften very messy things, full of branches to completely different methods dependingon the value of the calling arguments. Black boxes have, or should have, carefulcontrol of accuracy, to some stated uniform precision in all regimes.We will not be quite so fastidious in our examples, in part because we wantto illustrate techniques from Chapter 5, and in part because we want you tounderstand what goes on in the routines presented.
Some of our routines have anaccuracy parameter that can be made as small as desired, while others (especiallythose involving polynomial fits) give only a certain accuracy, one that we believeserviceable (typically six significant figures or more). We do not certify that theroutines are perfect black boxes. We do hope that, if you ever encounter troublein a routine, you will be able to diagnose and correct the problem on the basis ofthe information that we have given.In short, the special function routines of this chapter are meant to be used —we use them all the time — but we also want you to be prepared to understandtheir inner workings.CITED REFERENCES AND FURTHER READING:Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York) [full of useful numerical approximations to a great varietyof functions].IMSL Sfun/Library Users Manual (IMSL Inc., 2500 CityWest Boulevard, Houston TX 77042).NAG Fortran Library (Numerical Algorithms Group, 256 Banbury Road, Oxford OX27DE, U.K.),Chapter S.2122136.1 Gamma, Beta, and Related FunctionsHart, J.F., et al.
1968, Computer Approximations (New York: Wiley).Hastings, C. 1955, Approximations for Digital Computers (Princeton: Princeton University Press).Luke, Y.L. 1975, Mathematical Functions and Their Approximations (New York: Academic Press).6.1 Gamma Function, Beta Function, Factorials,Binomial CoefficientsThe gamma function is defined by the integral&∞Γ(z) =tz−1 e−t dt(6.1.1)0When the argument z is an integer, the gamma function is just the familiar factorialfunction, but offset by one,n! = Γ(n + 1)(6.1.2)The gamma function satisfies the recurrence relationΓ(z + 1) = zΓ(z)(6.1.3)If the function is known for arguments z > 1 or, more generally, in the half complexplane Re(z) > 1 it can be obtained for z < 1 or Re (z) < 1 by the reflection formulaΓ(1 − z) =πzπ=Γ(z) sin(πz)Γ(1 + z) sin(πz)(6.1.4)Notice that Γ(z) has a pole at z = 0, and at all negative integer values of z.There are a variety of methods in use for calculating the function Γ(z)numerically, but none is quite as neat as the approximation derived by Lanczos [1].This scheme is entirely specific to the gamma function, seemingly plucked fromthin air.
We will not attempt to derive the approximation, but only state theresulting formula: For certain integer choices of γ and N , and for certain coefficientsc1 , c2 , . . . , cN , the gamma function is given byΓ(z + 1) = (z + γ + 12 )z+ 2 e−(z+γ+ 2 )√c2cNc1++···++× 2π c0 +z +1 z+2z +N11(6.1.5)(z > 0)You can see that this is a sort of take-off on Stirling’s approximation, but with aseries of corrections that take into account the first few poles in the left complexplane. The constant c0 is very nearly equal to 1. The error term is parametrized by .For γ = 5, N = 6, and a certain set of c’s, the error is smaller than || < 2 × 10−10.Impressed? If not, then perhaps you will be impressed by the fact that (with thesesame parameters) the formula (6.1.5) and bound on apply for the complex gammafunction, everywhere in the half complex plane Re z > 0.214Chapter 6.Special FunctionsIt is better to implement ln Γ(x) than Γ(x), since the latter will overflow manycomputers’ floating-point representation at quite modest values of x.
Often thegamma function is used in calculations where the large values of Γ(x) are divided byother large numbers, with the result being a perfectly ordinary value. Such operationswould normally be coded as subtraction of logarithms. With (6.1.5) in hand, we cancompute the logarithm of the gamma function with two calls to a logarithm and 25or so arithmetic operations. This makes it not much more difficult than other built-infunctions that we take for granted, such as sin x or ex :#include <math.h>float gammln(float xx)Returns the value ln[Γ(xx)] for xx > 0.{Internal arithmetic will be done in double precision, a nicety that you can omit if five-figureaccuracy is good enough.double x,y,tmp,ser;static double cof[6]={76.18009172947146,-86.50532032941677,24.01409824083091,-1.231739572450155,0.1208650973866179e-2,-0.5395239384953e-5};int j;y=x=xx;tmp=x+5.5;tmp -= (x+0.5)*log(tmp);ser=1.000000000190015;for (j=0;j<=5;j++) ser += cof[j]/++y;return -tmp+log(2.5066282746310005*ser/x);}How shall we write a routine for the factorial function n!? Generally thefactorial function will be called for small integer values (for large values it willoverflow anyway!), and in most applications the same integer value will be called formany times.
It is a profligate waste of computer time to call exp(gammln(n+1.0))for each required factorial. Better to go back to basics, holding gammln in reservefor unlikely calls:#include <math.h>float factrl(int n)Returns the value n! as a floating-point number.{float gammln(float xx);void nrerror(char error_text[]);static int ntop=4;static float a[33]={1.0,1.0,2.0,6.0,24.0};int j;Fill in table only as required.if (n < 0) nrerror("Negative factorial in routine factrl");if (n > 32) return exp(gammln(n+1.0));Larger value than size of table is required.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
