c6-5 (779504)
Текст из файла
230Chapter 6.Special Functions6.5 Bessel Functions of Integer Orderk=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 ν ν11Jν (x) ∼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 ν arer112cos x − νπ − πJν (x) ∼πx24(6.5.4)r112Yν (x) ∼sin x − νπ − ππx24In the transition region where x ∼ ν, the typical amplitudes of the Bessel functionsare on the orderJν (ν) ∼21/3132/3 Γ( 23 )1/3ν 1/3Yν (ν) ∼ −2∼131/6 Γ( 23 ) ν 1/30.4473ν 1/30.7748∼ − 1/3ν(6.5.5)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).This 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 ν X∞(− 14 x2 )k1x(6.5.1)Jν (x) =2k!Γ(ν + k + 1)2316.5 Bessel Functions of Integer Order1J0J1J2J30− .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)r 288Jn (x) =Pncos(Xn ) − Qnsin(Xn )πxxx(6.5.9)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).Bessel functions.5232Chapter 6.rYn (x) =Special Functions 288Pnsin(Xn ) + Qncos(Xn )πxxx(6.5.10)where2n + 1π4(6.5.11)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).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).Xn ≡ x −6.5 Bessel Functions of Integer Order233}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)))));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).z=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);234Chapter 6.Special Functions}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.
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 bySample 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.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
