Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 58
Текст из файла (страница 58)
A certain statistic F , essentially the ratio of the observeddispersion of the first sample to that of the second one, is calculated. (For furtherdetails, see Chapter 14.) The probability that F would be as large as it is if thefirst sample’s underlying distribution actually has smaller variance than the second’sis denoted Q(F |ν1, ν2 ), where ν1 and ν2 are the number of degrees of freedomin the first and second samples, respectively.
In other words, Q(F |ν1, ν2) is thesignificance level at which the hypothesis “1 has smaller variance than 2” can berejected. A small numerical value implies a very significant rejection, in turnimplying high confidence in the hypothesis “1 has variance greater or equal to 2.”Q(F |ν1, ν2 ) has the limiting valuesQ(0|ν1 , ν2) = 1Q(∞|ν1, ν2) = 0(6.4.10)Its relation to the incomplete beta function Ix (a, b) as evaluated by betai above isQ(F |ν1, ν2 ) = Iν2ν2 +ν1 Fν2 ν1,2 2(6.4.11)Cumulative Binomial Probability DistributionSuppose an event occurs with probability p per trial.
Then the probability P ofits occurring k or more times in n trials is termed a cumulative binomial probability,and is related to the incomplete beta function Ix (a, b) as follows:P ≡n nj=kjpj (1 − p)n−j = Ip (k, n − k + 1)(6.4.12)For n larger than a dozen or so, betai is a much better way to evaluate the sum in(6.4.12) than would be the straightforward sum with concurrent computation of thebinomial coefficients. (For n smaller than a dozen, either method is acceptable.)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), Chapters 6 and 26.Pearson, E., and Johnson, N. 1968, Tables of the Incomplete Beta Function (Cambridge: Cambridge University Press).230Chapter 6.Special Functions6.5 Bessel Functions of Integer OrderThis section and the next one present practical algorithms for computing variouskinds of Bessel functions of integer order. In §6.7 we deal with fractional order. Infact, the more complicated routines for fractional order work fine for integer ordertoo.
For integer order, however, the routines in this section (and §6.6) are simplerand faster. Their only drawback is that they are limited by the precision of theunderlying rational approximations. For full double precision, it is best to work withthe routines for fractional order in §6.7.For any real ν, the Bessel function Jν (x) can be defined by the seriesrepresentation ν ∞(− 14 x2 )k1x(6.5.1)Jν (x) =2k!Γ(ν + k + 1)k=0The series converges for all x, but it is not computationally very useful for x 1.For ν not an integer the Bessel function Yν (x) is given byYν (x) =Jν (x) cos(νπ) − J−ν (x)sin(νπ)(6.5.2)The right-hand side goes to the correct limiting value Yn (x) as ν goes to someinteger n, but this is also not computationally useful.For arguments x < ν, both Bessel functions look qualitatively like simplepower laws, with the asymptotic forms for 0 < x ν ν11xJν (x) ∼ν ≥0Γ(ν + 1) 22Y0 (x) ∼ ln(x)(6.5.3)π −νΓ(ν) 1Yν (x) ∼ −xν>0π2For x > ν, both Bessel functions look qualitatively like sine or cosine waves whoseamplitude decays as x−1/2 .
The asymptotic forms for x ν are)112cos x − νπ − πJν (x) ∼πx24(6.5.4))112Yν (x) ∼sin x − νπ − ππx24In the transition region where x ∼ ν, the typical amplitudes of the Bessel functionsare on the orderJν (ν) ∼21/30.44731∼ 1/3ν32/3 Γ( 23 ) ν 1/321/30.77481Yν (ν) ∼ − 1/6 2 1/3 ∼ − 1/3ν3 Γ( 3 ) ν(6.5.5)2316.5 Bessel Functions of Integer Order1J0J1Bessel functions.5J2J30− .5Y0Y1Y2−1− 1.5−20Figure 6.5.1.24x6810Bessel functions J0 (x) through J3 (x) and Y0 (x) through Y2 (x).which holds asymptotically for large ν. Figure 6.5.1 plots the first few Besselfunctions of each kind.The Bessel functions satisfy the recurrence relationsJn+1 (x) =2nJn (x) − Jn−1 (x)x(6.5.6)Yn+1 (x) =2nYn (x) − Yn−1 (x)x(6.5.7)andAs already mentioned in §5.5, only the second of these (6.5.7) is stable in thedirection of increasing n for x < n.
The reason that (6.5.6) is unstable in thedirection of increasing n is simply that it is the same recurrence as (6.5.7): A smallamount of “polluting” Yn introduced by roundoff error will quickly come to swampthe desired Jn , according to equation (6.5.3).A practical strategy for computing the Bessel functions of integer order dividesinto two tasks: first, how to compute J0 , J1 , Y0 , and Y1 , and second, how to use therecurrence relations stably to find other J’s and Y ’s.
We treat the first task first:For x between zero and some arbitrary value (we will use the value 8),approximate J0 (x) and J1 (x) by rational functions in x. Likewise approximate byrational functions the “regular part” of Y0 (x) and Y1 (x), defined as221J1 (x) ln(x) −(6.5.8)andY1 (x) −Y0 (x) − J0 (x) ln(x)ππxFor 8 < x < ∞, use the approximating forms (n = 0, 1)) 288Jn (x) =Pncos(Xn ) − Qnsin(Xn )πxxx(6.5.9)232Chapter 6.)Yn (x) =Special Functions 288Pnsin(Xn ) + Qncos(Xn )πxxx(6.5.10)2n + 1π4(6.5.11)whereXn ≡ x −and where P0 , P1 , Q0 , and Q1 are each polynomials in their arguments, for 0 <8/x < 1.
The P ’s are even polynomials, the Q’s odd.Coefficients of the various rational functions and polynomials are given byHart [1], for various levels of desired accuracy. A straightforward implementation is#include <math.h>float bessj0(float x)Returns the Bessel function J0 (x) for any real x.{float ax,z;double xx,y,ans,ans1,ans2;Accumulate polynomials in double precision.if ((ax=fabs(x)) < 8.0) {Direct rational function fit.y=x*x;ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7+y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));ans2=57568490411.0+y*(1029532985.0+y*(9494680.718+y*(59272.64853+y*(267.8532712+y*1.0))));ans=ans1/ans2;} else {Fitting function (6.5.9).z=8.0/ax;y=z*z;xx=ax-0.785398164;ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4+y*(-0.2073370639e-5+y*0.2093887211e-6)));ans2 = -0.1562499995e-1+y*(0.1430488765e-3+y*(-0.6911147651e-5+y*(0.7621095161e-6-y*0.934935152e-7)));ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);}return ans;}#include <math.h>float bessy0(float x)Returns the Bessel function Y0 (x) for positive x.{float bessj0(float x);float z;double xx,y,ans,ans1,ans2;Accumulate polynomials in double precision.if (x < 8.0) {Rational function approximation of (6.5.8).y=x*x;ans1 = -2957821389.0+y*(7062834065.0+y*(-512359803.6+y*(10879881.29+y*(-86327.92757+y*228.4622733))));ans2=40076544269.0+y*(745249964.8+y*(7189466.438+y*(47447.26470+y*(226.1030244+y*1.0))));ans=(ans1/ans2)+0.636619772*bessj0(x)*log(x);} else {Fitting function (6.5.10).6.5 Bessel Functions of Integer Order233z=8.0/x;y=z*z;xx=x-0.785398164;ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4+y*(-0.2073370639e-5+y*0.2093887211e-6)));ans2 = -0.1562499995e-1+y*(0.1430488765e-3+y*(-0.6911147651e-5+y*(0.7621095161e-6+y*(-0.934945152e-7))));ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);}return ans;}#include <math.h>float bessj1(float x)Returns the Bessel function J1 (x) for any real x.{float ax,z;double xx,y,ans,ans1,ans2;Accumulate polynomials in double precision.if ((ax=fabs(x)) < 8.0) {Direct rational approximation.y=x*x;ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1+y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));ans2=144725228442.0+y*(2300535178.0+y*(18583304.74+y*(99447.43394+y*(376.9991397+y*1.0))));ans=ans1/ans2;} else {Fitting function (6.5.9).z=8.0/ax;y=z*z;xx=ax-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/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);if (x < 0.0) ans = -ans;}return ans;}#include <math.h>float bessy1(float x)Returns the Bessel function Y1 (x) for positive x.{float bessj1(float x);float z;double xx,y,ans,ans1,ans2;Accumulate polynomials in double precision.if (x < 8.0) {Rational function approximation of (6.5.8).y=x*x;ans1=x*(-0.4900604943e13+y*(0.1275274390e13+y*(-0.5153438139e11+y*(0.7349264551e9+y*(-0.4237922726e7+y*0.8511937935e4)))));ans2=0.2499580570e14+y*(0.4244419664e12+y*(0.3733650367e10+y*(0.2245904002e8+y*(0.1020426050e6+y*(0.3549632885e3+y)))));234Chapter 6.Special Functionsans=(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;}We now turn to the second task, namely how to use the recurrence formulas(6.5.6) and (6.5.7) to get the Bessel functions Jn (x) and Yn (x) for n ≥ 2.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
