c6-5 (779504), страница 2
Текст из файла (страница 2)
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).ans=(ans1/ans2)+0.636619772*(bessj1(x)*log(x)-1.0/x);} else {Fitting function (6.5.10).z=8.0/x;y=z*z;xx=x-2.356194491;ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4+y*(0.2457520174e-5+y*(-0.240337019e-6))));ans2=0.04687499995+y*(-0.2002690873e-3+y*(0.8449199096e-5+y*(-0.88228987e-6+y*0.105787412e-6)));ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);}return ans;6.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.Make 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;}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).#include <math.h>#define ACC 40.0#define BIGNO 1.0e10#define BIGNI 1.0e-10236Chapter 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]The 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) ≈(n − 1)! x −n22(6.6.2)n>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πKn (x) ≈ √exp(−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.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).6.6 Modified Bessel Functions of Integer Order.















