Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 59
Текст из файла (страница 59)
The latterof these is straightforward, since its upward recurrence is always stable:float bessy(int n, float x)Returns the Bessel function Yn (x) for positive x and n ≥ 2.{float bessy0(float x);float bessy1(float x);void nrerror(char error_text[]);int j;float by,bym,byp,tox;if (n < 2) nrerror("Index n less than 2 in bessy");tox=2.0/x;by=bessy1(x);Starting values for the recurrence.bym=bessy0(x);for (j=1;j<n;j++) {Recurrence (6.5.7).byp=j*tox*by-bym;bym=by;by=byp;}return by;}The cost of this algorithm is the call to bessy1 and bessy0 (which generate acall to each of bessj1 and bessj0), plus O(n) operations in the recurrence.As for Jn (x), things are a bit more complicated.
We can start the recurrenceupward on n from J0 and J1 , but it will remain stable only while n does not exceedx. This is, however, just fine for calls with large x and small n, a case whichoccurs frequently in practice.The harder case to provide for is that with x < n. The best thing to do hereis to use Miller’s algorithm (see discussion preceding equation 5.5.16), applyingthe recurrence downward from some arbitrary starting value and making use of theupward-unstable nature of the recurrence to put us onto the correct solution. Whenwe finally arrive at J0 or J1 we are able to normalize the solution with the sum(5.5.16) accumulated along the way.The only subtlety is in deciding at how large an n we need start the downwardrecurrence so as to obtain a desired accuracy by the time we reach the n that wereally want. If you play with the asymptotic forms (6.5.3) and (6.5.5), you shouldbe able to convince yourself that the answer is to start larger than the desired n by6.5 Bessel Functions of Integer Order235an additive amount of order [constant × n]1/2 , where the square root of the constantis, very roughly, the number of significant figures of accuracy.The above considerations lead to the following function.#include <math.h>#define ACC 40.0#define BIGNO 1.0e10#define BIGNI 1.0e-10Make larger to increase accuracy.float bessj(int n, float x)Returns the Bessel function Jn (x) for any real x and n ≥ 2.{float bessj0(float x);float bessj1(float x);void nrerror(char error_text[]);int j,jsum,m;float ax,bj,bjm,bjp,sum,tox,ans;if (n < 2) nrerror("Index n less than 2 in bessj");ax=fabs(x);if (ax == 0.0)return 0.0;else if (ax > (float) n) {Upwards recurrence from J0 and J1 .tox=2.0/ax;bjm=bessj0(ax);bj=bessj1(ax);for (j=1;j<n;j++) {bjp=j*tox*bj-bjm;bjm=bj;bj=bjp;}ans=bj;} else {Downwards recurrence from an even m here comtox=2.0/ax;puted.m=2*((n+(int) sqrt(ACC*n))/2);jsum=0;jsum will alternate between 0 and 1; when it isbjp=ans=sum=0.0;1, we accumulate in sum the even terms inbj=1.0;(5.5.16).for (j=m;j>0;j--) {The downward recurrence.bjm=j*tox*bj-bjp;bjp=bj;bj=bjm;if (fabs(bj) > BIGNO) {Renormalize to prevent overflows.bj *= BIGNI;bjp *= BIGNI;ans *= BIGNI;sum *= BIGNI;}if (jsum) sum += bj;Accumulate the sum.jsum=!jsum;Change 0 to 1 or vice versa.if (j == n) ans=bjp;Save the unnormalized answer.}sum=2.0*sum-bj;Compute (5.5.16)ans /= sum;and use it to normalize the answer.}return x < 0.0 && (n & 1) ? -ans : ans;}236Chapter 6.Special FunctionsCITED 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), Chapter 9.Hart, J.F., et al. 1968, Computer Approximations (New York: Wiley), §6.8, p. 141. [1]6.6 Modified Bessel Functions of Integer OrderThe modified Bessel functions In (x) and Kn (x) are equivalent to the usualBessel functions Jn and Yn evaluated for purely imaginary arguments.
In detail,the relationship isIn (x) = (−i)n Jn (ix)Kn (x) =π n+1i[Jn (ix) + iYn (ix)]2(6.6.1)The particular choice of prefactor and of the linear combination of Jn and Yn to formKn are simply choices that make the functions real-valued for real arguments x.For small arguments x n, both In (x) and Kn (x) become, asymptotically,simple powers of their argumentIn (x) ≈1 # x $nn! 2n≥0K0 (x) ≈ − ln(x)Kn (x) ≈(6.6.2)(n − 1)! # x $−n22n>0These expressions are virtually identical to those for Jn (x) and Yn (x) in this region,except for the factor of −2/π difference between Yn (x) and Kn (x).
In the regionx n, however, the modified functions have quite different behavior than theBessel functions,1exp(x)In (x) ≈ √2πxπexp(−x)Kn (x) ≈ √2πx(6.6.3)The modified functions evidently have exponential rather than sinusoidalbehavior for large arguments (see Figure 6.6.1). The smoothness of the modifiedBessel functions, once the exponential factor is removed, makes a simple polynomialapproximation of a few terms quite suitable for the functions I0 , I1 , K0 , and K1 .The following routines, based on polynomial coefficients given by Abramowitz andStegun [1], evaluate these four functions, and will provide the basis for upwardrecursion for n > 1 when x > n.2376.6 Modified Bessel Functions of Integer Order4modified Bessel functions3I02K0K1K2I1I21I30012x3Figure 6.6.1.
Modified Bessel functions I0 (x) through I3 (x), K0 (x) through K2 (x).#include <math.h>float bessi0(float x)Returns the modified Bessel function I0 (x) for any real x.{float ax,ans;double y;Accumulate polynomials in double precision.if ((ax=fabs(x)) < 3.75) {Polynomial fit.y=x/3.75;y*=y;ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492+y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));} else {y=3.75/ax;ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1+y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2+y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1+y*0.392377e-2))))))));}return ans;}#include <math.h>float bessk0(float x)Returns the modified Bessel function K0 (x) for positive real x.{float bessi0(float x);double y,ans;Accumulate polynomials in double precision.4238Chapter 6.Special Functionsif (x <= 2.0) {Polynomial fit.y=x*x/4.0;ans=(-log(x/2.0)*bessi0(x))+(-0.57721566+y*(0.42278420+y*(0.23069756+y*(0.3488590e-1+y*(0.262698e-2+y*(0.10750e-3+y*0.74e-5))))));} else {y=2.0/x;ans=(exp(-x)/sqrt(x))*(1.25331414+y*(-0.7832358e-1+y*(0.2189568e-1+y*(-0.1062446e-1+y*(0.587872e-2+y*(-0.251540e-2+y*0.53208e-3))))));}return ans;}#include <math.h>float bessi1(float x)Returns the modified Bessel function I1 (x) for any real x.{float ax,ans;double y;Accumulate polynomials in double precision.if ((ax=fabs(x)) < 3.75) {Polynomial fit.y=x/3.75;y*=y;ans=ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934+y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));} else {y=3.75/ax;ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1-y*0.420059e-2));ans=0.39894228+y*(-0.3988024e-1+y*(-0.362018e-2+y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));ans *= (exp(ax)/sqrt(ax));}return x < 0.0 ? -ans : ans;}#include <math.h>float bessk1(float x)Returns the modified Bessel function K1 (x) for positive real x.{float bessi1(float x);double y,ans;Accumulate polynomials in double precision.if (x <= 2.0) {Polynomial fit.y=x*x/4.0;ans=(log(x/2.0)*bessi1(x))+(1.0/x)*(1.0+y*(0.15443144+y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1+y*(-0.110404e-2+y*(-0.4686e-4)))))));} else {y=2.0/x;ans=(exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619+y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2+y*(0.325614e-2+y*(-0.68245e-3)))))));}return ans;}6.6 Modified Bessel Functions of Integer Order239The recurrence relation for In (x) and Kn (x) is the same as that for Jn (x)and Yn (x) provided that ix is substituted for x.
This has the effect of changinga sign in the relation,2nIn+1 (x) = −In (x) + In−1 (x)x 2nKn (x) + Kn−1 (x)Kn+1 (x) = +x(6.6.4)These relations are always unstable for upward recurrence. For Kn , itself growing,this presents no problem. For In , however, the strategy of downward recursion istherefore required once again, and the starting point for the recursion may be chosenin the same manner as for the routine bessj. The only fundamental difference isthat the normalization formula for In (x) has an alternating minus sign in successiveterms, which again arises from the substitution of ix for x in the formula usedpreviously for Jn1 = I0 (x) − 2I2 (x) + 2I4 (x) − 2I6 (x) + · · ·(6.6.5)In fact, we prefer simply to normalize with a call to bessi0.With this simple modification, the recursion routines bessj and bessy becomethe new routines bessi and bessk:float bessk(int n, float x)Returns the modified Bessel function Kn (x) for positive x and n ≥ 2.{float bessk0(float x);float bessk1(float x);void nrerror(char error_text[]);int j;float bk,bkm,bkp,tox;if (n < 2) nrerror("Index n less than 2 in bessk");tox=2.0/x;bkm=bessk0(x);Upward recurrence for all x...bk=bessk1(x);for (j=1;j<n;j++) {...and here it is.bkp=bkm+j*tox*bk;bkm=bk;bk=bkp;}return bk;}#include <math.h>#define ACC 40.0#define BIGNO 1.0e10#define BIGNI 1.0e-10Make larger to increase accuracy.float bessi(int n, float x)Returns the modified Bessel function In (x) for any real x and n ≥ 2.{float bessi0(float x);void nrerror(char error_text[]);240Chapter 6.Special Functionsint j;float bi,bim,bip,tox,ans;if (n < 2) nrerror("Index n less than 2 in bessi");if (x == 0.0)return 0.0;else {tox=2.0/fabs(x);bip=ans=0.0;bi=1.0;for (j=2*(n+(int) sqrt(ACC*n));j>0;j--) {Downward recurrence from evenbim=bip+j*tox*bi;m.bip=bi;bi=bim;if (fabs(bi) > BIGNO) {Renormalize to prevent overflows.ans *= BIGNI;bi *= BIGNI;bip *= BIGNI;}if (j == n) ans=bip;}ans *= bessi0(x)/bi;Normalize with bessi0.return x < 0.0 && (n & 1) ? -ans : ans;}}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), §9.8. [1]Carrier, G.F., Krook, M. and Pearson, C.E. 1966, Functions of a Complex Variable (New York:McGraw-Hill), pp.
220ff.6.7 Bessel Functions of Fractional Order, AiryFunctions, Spherical Bessel FunctionsMany algorithms have been proposed for computing Bessel functions of fractional ordernumerically. Most of them are, in fact, not very good in practice. The routines given here arerather complicated, but they can be recommended wholeheartedly.Ordinary Bessel FunctionsThe basic idea is Steed’s method, which was originally developed [1] for Coulomb wavefunctions. The method calculates Jν , Jν , Yν , and Yν simultaneously, and so involves fourrelations among these functions. Three of the relations come from two continued fractions,one of which is complex. The fourth is provided by the Wronskian relation2W ≡ Jν Yν − Yν Jν =(6.7.1)πxThe first continued fraction, CF1, is defined byfν ≡JννJν+1= −JνxJν1ν1= −···x2(ν + 1)/x − 2(ν + 2)/x −(6.7.2)6.7 Bessel Functions of Fractional Order241You can easily derive it from the three-term recurrence relation for Bessel functions: Start withequation (6.5.6) and use equation (5.5.18).
Forward evaluation of the continued fraction byone of the methods of §5.2 is essentially equivalent to backward recurrence of the recurrencerelation.(The rate of convergence of CF1 is determined by the position of the turning pointxtp = ν(ν + 1) ≈ ν, beyond which the Bessel functions become oscillatory. If x <∼ xtp,convergence is very rapid. If x >x,theneachiterationofthecontinuedfractioneffectivelytp∼increases ν by one until x <∼ xtp ; thereafter rapid convergence sets in. Thus the numberof iterations of CF1 is of order x for large x. In the routine bessjy we set the maximumallowed number of iterations to 10,000. For larger x, you can use the usual asymptoticexpressions for Bessel functions.One can show that the sign of Jν is the same as the sign of the denominator of CF1once it has converged.The complex continued fraction CF2 is defined byp + iq ≡Jν + iYνi (1/2)2 − ν 2 (3/2)2 − ν 21+i+···=−Jν + iYν2xx 2(x + i) + 2(x + 2i) +(6.7.3)(We sketch the derivation of CF2 in the analogous case of modified Bessel functions in thenext subsection.) This continued fraction converges rapidly for x >∼ xtp, while convergencefails as x → 0.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
