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;