NUMAL Section 6.10.1
BEGIN SECTION : 6.10.1 (December, 1978)
AUTHORS: M.BAKKER AND N.M.TEMME.
CONTRIBUTOR: R.MONTIJN.
INSTITUTE: MATHEMATICAL CENTRE.
RECEIVED: 781101.
BRIEF DESCRIPTION:
THIS SECTION CONTAINS THE PROCEDURES:
BESS JAPLUSN:
THIS PROCEDURE CALCULATES THE BESSEL FUNCTIONS
OF THE FIRST KIND OF ORDER A+K (0<=K<=N, 0<=A<1) AND
ASSIGNS THEM TO AN ARRAY. THE ARGUMENT MUST BE NON-NEGATIVE.
BESS YA01:
THIS PROCEDURE CALCULATES THE BESSEL FUNCTIONS
OF THE SECOND KIND (ALSO CALLED NEUMANN'S FUNCTIONS)
OF ORDER A AND A+1 AND ARGUMENT X>0.
BESS YAPLUSN:
THIS PROCEDURE GENERATES AN ARRAY OF BESSEL FUNCTIONS OF THE
SECOND KIND OF ORDER A+N, N=0, 1, 2, ..., NMAX, AND
ARGUMENT X>0.
THE BESSEL FUNCTIONS OF THE SECOND KIND CORRESPOND TO THE
FUNCTION DEFINED IN FORMULA 9.1.2 OF REFERENCE [1].
BESS PQA01:
THIS PROCEDURE IS AN AUXILIARY PROCEDURE FOR THE COMPUTATION OF
THE BESSEL FUNCTIONS FOR LARGE VALUES OF THEIR ARGUMENT.
BESS ZEROS:
THIS PROCEDURE CALCULATES THE FIRST N ZEROS OF A BESSEL
FUNCTION OF THE FIRST OR THE SECOND KIND OR ITS DERIVATIVE.
START:
THIS IS AN AUXILIARY PROCEDURE WHICH COMPUTES A STARTING VALUE
OF AN ALGORITHM USED IN SEVERAL BESSEL FUNCTION PROCEDURES.
KEYWORDS:
BESSEL FUNCTION, BESSEL FUNCTION OF THE SECOND KIND, NEUMANN'S
FUNCTION, ZEROS OF BESSEL FUNCTIONS.
REFERENCES:
[1]. ABRAMOWITZ, M., AND STEGUN, I. (EDS),
HANDBOOK OF MATHEMATICAL FUNCTIONS WITH FORMULAS, GRAPHS AND
MATHEMATICAL TABLES.
APPL. MATH. SER. 55, U.S. GOVT. PRINTING OFFICE,
WASHINGTON, D.C. , 1974.
[2]. GAUTSCHI, W., COMPUTATIONAL ASPECTS OF
THREE TERM RECURRENCE RELATIONS.
SIAM REVIEW, VOLUME 9(1967), NUMBER 1, P.24 FF.
[3]. TEMME, N.M. ON THE NUMERICAL EVALUATION OF THE
ORDINARY BESSEL FUNCTION OF THE SECOND KIND.
J. COMP. PHYS., 21, P. 343 FF, 1976.
[4]. WATSON, G.N.
A TREATISE ON THE THEORY OF BESSEL FUNCTIONS.
CAMBRIDGE UNIV. PRESS, LONDON AND NEW YORK, 1945.
[5]. TEMME, N.M., SPECIALE FUNCTIES, IN:
COLLOQUIUM NUMERIEKE PROGRAMMATUUR,
J.C.P. BUS (RED.), MC SYLLABUS 29.1B,
MATHEMATICAL CENTRE, AMSTERDAM, 1976.
[6]. TEMME, N.M., AN ALGOLRITHM WITH ALGOL 60 IMPLEMENTATION
FOR THE CALCULATION OF THE ZEROS OF A BESSEL FUNCTION,
REPORT TW 179 MATHEMATICAL CENTRE, AMSTERDAM, 1978.
SUBSECTION: BESS JAPLUSN.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" BESS JAPLUSN(A, X, N, JA);
"VALUE" A, X, N;
"INTEGER" N; "REAL" A, X; "ARRAY" JA;
"CODE" 35180;
THE MEANING OF THE FORMAL PARAMETERS IS:
A: < ARITHMETIC EXPRESSION > ;
THE NONINTEGER PART OF THE ORDER; 0 <= A < 1;
X: < ARITHMETIC EXPRESSION >;
THE ARGUMENT VALUE; X > = 0;
N: < ARITHMETIC EXPRESSION >;
THE UPPER BOUND OF THE INDICES OF THE ARRAY JA;
JA: < ARRAY IDENTIFIER >;
"ARRAY" JA[0:N];
EXIT: JA[K] IS ASSIGNED THE VALUE OF THE BESSEL
FUNCTION OF THE FIRST KIND J[K+A](X),
0 < = K < = N.
PROCEDURES USED:
BESS J = CP 35162,
SPHER BESS J = CP 35150,
GAMMA = CP 35061,
START = CP 35185.
REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED.
METHOD AND PERFORMANCE:
IN ALL THE CASES THE BESSEL FUNCTIONS ARE COMPUTED
ACCORDING TO THE MILLER METHOD DISCRIBED IN [2, P.46-52].
THE STARTING VALUE IS COMPUTED BY THE PROCEDURE START.
RUNNING TIME: ROUGHLY PROPORTIONAL TO THE MAXIMUM OF
X AND N.
EXAMPLE OF USE:
"BEGIN" "INTEGER" N; "REAL" A, X; "ARRAY" JA[0:2];
X:= 2; A:= .78; N:= 2;
BESS JAPLUSN(A, X, N, JA);
OUTPUT(61, "("/, "("X=")"D, 3B"("A=")".DD, 3B"("N=")"D,
/, 3(3B-.14D"-ZD)")", X, A, N, JA[0], JA[1], JA[2])
"END"
RESULTS:
X=2 A= .78 N=2
.57306126928364"0 .41529475124424" 0 .16616338793111" 0
SUBSECTION: BESS YA01.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" BESS YA01(A, X, YA, YA1);
"VALUE" A, X; "REAL" A, X, YA, YA1;
"CODE" 35181;
THE MEANING OF THE FORMAL PARAMETERS IS:
A: <ARITHMETIC EXPRESSION>;
THE ORDER;
X: <ARITHMETIC EXPRESSION>;
THE ARGUMENT. THIS ARGUMENT SHOULD SATISFY X>0;
YA: <VARIABLE>;
EXIT: THE NEUMANN FUNCTION OF ORDER A
AND ARGUMENT X;
YA1: <VARIABLE>;
EXIT: THE NEUMANN FUNCTION OF ORDER A+1.
PROCEDURES USED:
RECIP GAMMA = CP 35060;
BESS PQA01 = CP 35183;
SINH = CP 35111.
REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED.
METHOD AND PERFORMANCE:
FOR 0<X<3 THE BESSEL FUNCTIONS ARE COMPUTED BY USING TAYLOR
SERIES. THE METHOD IS DESCRIBED IN REFERENCE [3].
FOR X>=3 THE PROCEDURE CALLS FOR THE PROCEDURE BESS PQA01
(SEE SUBSECTION BESS PQA01).
THE RELATIVE ACCURACY IS ABOUT "-13, EXCEPT FOR LARGE VALUES OF X;
IN THAT CASE THE ACCURACY MAINLY DEPENDS ON THE ACCURACY OF THE
FUNCTIONS SIN(X) AND COS(X).
EXAMPLE OF USE:
THE PROGRAM:
"BEGIN" "REAL" P, Q;
BESS YA01(0, 1, P, Q);
OUTPUT(61, "("2(N)")", P, Q)
"END"
YIELDS THE FOLLOWING RESULTS
+8.8256964215677"-002 -7.8121282130028"-001.
SUBSECTION: BESS YAPLUSN.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" BESS YAPLUSN(A, X, NMAX, YAN); "VALUE" A, X, NMAX;
"REAL" A, X; "INTEGER" NMAX; "ARRAY" YAN;
"CODE" 35182;
THE MEANING OF THE FORMAL PARAMETERS IS:
A: <ARITHMETIC EXPRESSION>;
THE ORDER;
X: <ARITHMETIC EXPRESSION>;
THE ARGUMENT; THIS ARGUMENT SHOULD SATISFY X>0;
NMAX: <ARITHMETIC EXPRESSION>;
THE UPPER BOUND OF THE INDICES OF THE ARRAY YAN;
YAN: <ARRAY IDENTIFIER>;
"ARRAY" YAN[0:NMAX]; NMAX>=0;
EXIT: THE VALUES OF THE BESSEL FUNCTIONS OF
THE SECOND KIND OF ORDER A+K, FOR THE ARGUMENT X
ARE ASSIGNED TO YAN[K],0<=K<=NMAX.
PROCEDURES USED: BESS YA01 = CP 35181.
REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED.
METHOD AND PERFORMANCE:
THE RECURRENCE RELATION
YAN[N+1]= -YAN[N-1] + 2*(N+A)*YAN[N]/X
IS USED. THE INITIAL VALUES ARE OBTAINED FROM THE
PROCEDURE BESS YA01. THE RECURRENCE RELATION IS NUMERICALLY
STABLE IN THE FORWARD DIRECTION (IF A >= 0).
EXAMPLE OF USE:
THE PROGRAM:
"BEGIN" "ARRAY" YAN[0:2];
BESS YAPLUSN(0, 1, 2, YAN);
OUTPUT(61, "("3(N)")", YAN[0], YAN[1], YAN[2])
"END"
YIELDS THE FOLLOWING RESULTS
+8.8256964215677"-002 -7.8121282130028"-001 -1.6506826068163"+000.
SUBSECTION: BESS PQA01.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" BESS PQA01(A, X, PA, QA, PA1, QA1); "VALUE" X, A;
"REAL" X, A, PA, QA, PA1, QA1;
"CODE" 35183;
THE MEANING OF THE FORMAL PARAMETERS IS:
A: <ARITHMETIC EXPRESSION>;
THE ORDER;
X: <ARITHMETIC EXPRESSION>;
THE ARGUMENT. THIS ARGUMENT SHOULD SATISFY X>0;
PA: <VARIABLE>;
EXIT: THIS FUNCTION CORRESPONDS TO THE FUNCTION
P(X, A) DEFINED ON P. 205 OF REFERENCE [4].
SEE ALSO REFERENCE [1], FORMULA 9.2.6;
QA: <VARIABLE>;
EXIT: THIS FUNCTION CORRESPONDS TO THE FUNCTION
Q(X, A) DEFINED ON P.205 OF REFERENCE [4].
SEE ALSO REFERENCE [1], FORMULA 9.2.6;
PA1: <VARIABLE>;
EXIT: THE FUNCTION P(X, A+1);
QA1: <VARIABLE>;
EXIT: THE FUNCTION Q(X, A+1).
PROCEDURES USED:
BESS JAPLUSN = CP35180,
BESS YA01 = CP35181.
REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED.
METHOD AND PERFORMANCE:
X < 3 :
PA, QA, PA1, QA1 ARE COMPUTED FROM THE RELATIONS
PA = B * (YA0 * S + JA0 * C),
QA = B * (YA0 * C - JA0 * S),
PA1 = B * (JA1 * S - YA1 * C),
QA1 = B * (JA1 * C + YA1 * S),
WHERE
B = SQRT(HALFPI * X),
C = COS(X - HALFPI * (A + .5)),
S = SIN(X - HALFPI * (A + .5)),
HALFPI = 1.57079 63267 9489,
YA0 = Y[A](X),
YA1 = Y[A + 1](X),
JA0 = J[A](X),
JA1 = J[A + 1](X);
X >= 3:
THE METHOD IS DESCRIBED IN REFERENCE [3]. IT DEPENDS ON USING
A MILLER ALGORITHM FOR CONFLUENT HYPERGEOMETRIC FUNCTIONS.
THE ACCURACY IS ABOUT "-13 AND IS BETTER FOR LARGE X.
THE FUNCTIONS PA AND QA CAN ALSO BE USED FOR THE COMPUTATION
OF THE BESSEL FUNCTION J OF THE FIRST KIND.
SEE REFERENCE[1], FORMULA 9.2.5.
EXAMPLE OF USE:
FROM SOME PROPERTIES OF THE BESSEL FUNCTIONS IT CAN BE PROVED
THAT PA*PA1+QA*QA1=1, WHATEVER X AND A. IN THE FOLLOWING PROGRAM
WE VERIFY THIS RELATION.
"BEGIN" "REAL" A, X, P, Q, R, S;
"FOR" X:= 1, 3, 5, 10, 15, 20, 50 "DO"
"BEGIN" BESS PQA01(0, X, P, Q, R, S);
OUTPUT(61, "("BB, D.2D"+3D")", ABS(P*R+Q*S-1))
"END"
"END"
THIS PROGRAM GIVES THE FOLLOWING RESULTS:
1.42"-014 7.11"-015 7.11"-015 7.11"-015 1.42"-014 0.00"+000
2.13"-014.
SUBSECTION: BESS ZEROS.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"PROCEDURE" BESS ZEROS(A,N,Z,D);
"VALUE" A,N,D; "REAL" A;
"INTEGER" N,D; "ARRAY" Z;
"CODE" 35184;
THE MEANING OF THE FORMAL PARAMETERS IS:
A: <ARITHMETIC EXPRESSION>;
THE ORDER OF THE BESSEL FUNCTION, A>=0.
N: <ARITHMETIC EXPRESSION>;
THE NUMBER OF ZEROS TO BE EVALUATED, N>=1.
Z: <ARRAY IDENTIFIER>;
"ARRAY" Z[1:N];
EXIT: Z[J] IS THE J-TH ZERO OF THE
SELECTED BESSEL FUNCTON;
D: <ARITHMETIC EXPRESSION>;
THE CHOICE OF D DETERMINES THE TYPE OF THE
BESSEL FUNCTION OF WHICH THE ZEROS ARE COMPUTED:
IF D=1 THEN JA ,
IF D=2 THEN YA ,
IF D=3 THEN JA-PRIME,
IF D=4 THEN YA-PRIME.
PROCEDURES USED: BESS PQA01 = CP 35183.
REQUIRED CENTRAL MEMMORY: NO AUXILIARY ARRAYS ARE USED.
RUNNING TIME: DEPENDS ON THE VALUES OF A AND N AND ON
THE MUMBER OF ITERATIONS IN THE ALGORITHM.
FROM TESTS IT FOLLOWS THAT FOR EACH ZERO AT MOST 3
EVALUATIONS OF THE PROCEDURE BESS PQA01 ARE NEEDED.
METHOD AND PERFORMANCE:
A FIRST APPROXIMATION OF THE ZEROS OF THE SELECTED BESSEL
FUNCTION IS CALCULATED BY MEANS OF THE ASYMPTOTIC EXPANTIONS
( SEE THE FORMULAS 9.5.12, 9.5.13 ( FOR A < 3 ) AND 9.5.22,
9.5.24( FOR A >= 3 ) OF REF [1] ). THIS VALUE IS CORRECTED BY THE
USE OF A FOURTH ORDER NEWTON-RAPHSON METHOD AS DISCRIBED ON P. 179
OF REF [6]. MORE DETAILS CAN BE FOUND IN REF [7].
A RELATIVE PRECISION OF 13 DIGITS IS PERSUED.
THE COMPUTATION OF A ZERO IS TERMINATED IF THIS ACCURRACY
IS ACHIEVED OR IF MORE THAN 5 ITERATIONS ARE NEEDED.
THE PROCEDURE DOES NOT CHECK ON THE RANGE OF THE PARAMETERS
A,N AND D.
EXAMPLE OF USE:
THE PROGRAM
"BEGIN" "REAL" A; "INTEGER" N,D; "ARRAY" Z[1:2];
A:=3.14; N:= 2; D:= 2;
BESS ZEROS(A,N,Z,D);
OUTPUT(61,"("N,/,N")",Z[1],Z[2])
"END"
PRINTS THE FIRST TWO ZEROS OF THE BESSEL FUNCTION Y OF
THE ORDER 3.14; THE RESULT IS:
+4.6847847078799"+000
+8.2765898338392"+000
SUBSECTION: START.
CALLING SEQUENCE:
THE HEADING OF THE PROCEDURE READS:
"INTEGER" "PROCEDURE" START(X,N,T);
"VALUE" X,N,T; "REAL" X;
"INTEGER" N,T;
"CODE" 35185;
START:= A STARTING VALUE FOR THE MILLER ALGORITHM
FOR COMPUTING AN ARRAY OF BESSEL FUNCTIONS;
THE MEANING OF THE FORMAL PARAMETERS IS:
X: <ARITHMETIC EXPRESSION>;
THE ARGUMENT OF THE BESSEL FUNCTIONS, X > 0;
N: <ARITHMETIC EXPRESSION>;
THE NUMBER OF BESSEL FUNCTIONS TO BE COMPUTED, N >= 0;
T: <ARITHMETIC EXPRESSION>;
THE TYPE OF BESSEL FUNCTION IN QUESTION,
T = 0 CORRESPONDS TO ORDINARY BESSEL FUNCTIONS,
T = 1 CORRESPONDS TO MODIFIED BESSEL FUNCTIONS.
PROCEDURES USED: NONE.
REQUIRED CENTRAL MEMORY: NO AUXILIARY ARRAYS ARE DECLARED.
METHOD AND PERFORMANCE:
THE PROCEDURE IS CALLED IN THE FOLLOWING PROCEDURES:
BESS J CODE 35162
NON EXP BESS I CODE 35177
BESS JAPLUSN CODE 35180
BESS KAPLUSN CODE 35192
NON EXP BESS IAPLUSN CODE 35193
SPHER BESS J CODE 35150
NON EXP SPHER BESS I CODE 35154.
IN THESE PROCEDURES AN ARRAY OF BESSEL FUNCTIONS IS GENERATED
BY USING MILLER 'S ALGORITHM (SEE REF[5]). FOR STARTING THIS
ALGORITHM ONE NEEDS AN INTEGER NU WHICH CAN BE COMPUTED BY USING
GAUTSCHI 'S ESTIMATES OF THE ERROR ( SEE REF[5,FORMULA (5.11)] ).
WE COMPUTE THIS STARTING VALUE NU BY USING ASYMPTOTIC APPROXIMA-
TIONS OF THE BESSEL FUNCTIONS, AS GIVEN IN REF[1, FORMULA 9.3.7,
9.3.8, 9.7.7, AND 9.7.8]. GAUTSCHI USED DIFFERENT FORMULAS, BUT
THOSE USED HERE GIVE FOR LARGE X AND N MORE REALISTIC ESTIMATES.
THE PERSUED ACCURACY IN THE ABOVE MENTIONED PROCEDURES IS ABOUT
"-14 . FOR OBTAINING AN ACCURACY OF "-D THE NUMBERS 36 AND 18
APPEARING IN THE FOURTH AND SIXTH LINE OF THE SOURCE TEXT OF START
SHOULD BE REPLACED BY (D+1)* LN(10) AND .5*(D+1)* LN(10),
RESPECTIVELY. FOR MODIFIED BESSEL FUNCTIONS THE ACCURRACY IS IN A
RELATIVE SENSE; FOR ORDINARY BESSEL FUNCTIONS THE ACCURRACY IS
ABSOLUTE IF THE ORDER OF THE BESSEL FUNCTION IS SMALLER THAN X,
OTHERWISE IT IS RELATIVE.
RUNNING TIME: NEGLECTABLE IF COMPARED WITH THE TIME NEEDED
FOR THE BESSEL FUNCTION PROCEDURES.
EXAMPLE OF USE: SEE THE ABOVE MENTIONED PROCEDURES.
SOURCE TEXT(S):
"CODE" 35180;
"CODE" 35181;
"CODE" 35182;
"CODE" 35183;
"CODE" 35184;
"CODE" 35185;