NUMAL Section 3.6.2

BEGIN SECTION : 3.6.2 (December, 1978)

AUTHORS:        M. BAKKER, P.J. HARINGHUIZEN AND C.G. VAN DER LAAN.

CONTRIBUTORS:   M. BAKKER, P.J. HARINGHUIZEN,
                C.G. VAN DER LAAN AND M. VOORINTHOLT.

INSTITUTE:      MATHEMATICAL CENTRE AND RIJKSUNIVERSITEIT GRONINGEN.

RECEIVED:       780601.

BRIEF DESCRIPTION:

    THIS SECTION CONTAINS FIVE PROCEDURES FOR CALCULATING ZEROS OF
    ORTHOGONAL POLYNOMIALS WHICH ARE GIVEN BY THE COEFFICIENTS OF
    THEIR RECURRENCE RELATION:
    ALLZERORTPOL: CALCULATES ALL ZEROS,
    LUPZERORTPOL: CALCULATES A NUMBER OF ADJACENT UPPER OR LOWER ZEROS,
    SELZERORTPOL: CALCULATES A NUMBER OF ADJACENT ZEROS. IT IS
                  EFFICIENT TO USE ALLZERORTPOL IF MORE THAN 50
                  PERCENT OF EXTREME ZEROS OR MORE THAN 25 PERCENT OF
                  SELECTED ZEROS ARE WANTED.
    ALLJACZER   : CALCULATES THE ZEROS OF THE N-TH JACOBIAN POLYNOMIAL.
    ALLLAGZER   : CALCULATES THE ZEROS OF THE N-TH LAGUERRE POLYNOMIAL.

KEYWORDS:

    ZEROS,
    ORTHOGONAL POLYNOMIALS,
    CHRISTOFFEL ABSCISSAS.

REFERENCES:

    ABRAMOWITZ, M. AND I.A. STEGUN (1964):
    HANDBOOK OF MATHEMATICAL FUNCTIONS.
    DOVER PUBLICATIONS INC.

    GOLUB, G.H. AND J.H. WELSCH (1969):
    CALCULATION OF GAUSS QUADRATURE RULES.
    MATH. COMP. VOL. 23, P.221-230.

    LANCZOS, C. (1957):
    APPLIED ANALYSIS.
    PRENTICE HALL.

    STOER, J. (1972):
    EINFUEHRUNG IN DIE NUMERISCHE MATHEMATIK 1.
    HEIDELBERG TASCHENBUECHER 105, SPRINGER.

    WILKINSON,J AND REINSCH,C. :
    HANDBOOK OF AUTOMATIC COMPUTATION. VOL. 2.
    LINEAR ALGEBRA
    HEIDELBERG (1971).


SUBSECTION:         ALLZERORTPOL.

CALLING SEQUENCE:

    THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:

    "PROCEDURE" ALLZERORTPOL (N, B, C, ZER, EM);
    "VALUE" N; "INTEGER" N; "ARRAY" B, C, ZER, EM;
    "CODE" 31362;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:      <ARITHMETIC EXPRESSION>;
            ENTRY:  THE DEGREE OF THE ORTHOGONAL POLYNOMIAL OF WHICH
                    THE ZEROS ARE TO BE CALCULATED;
    B, C:   <ARRAY IDENTIFIER>;
            "ARRAY" B, C [0:N-1];
            ENTRY:  THE ELEMENTS B[I] AND C[I], I = 0, 1, ... , N-1,
                    CONTAIN THE COEFFICIENTS OF THE RECURRENCE RELATION
                    P[I+1](X) = (X - B[I]) * P[I](X) - C[I] * P[I-1](X)
                    I = 0, 1, ... , N-1; ASSUMED IS C[0]=0, WHILE THE
                    CONTENTS OF THE ARRAYS ARE PRESERVED;
    ZER:    <ARRAY IDENTIFIER>;
            "ARRAY" ZER[1:N];
            EXIT:   THE ZEROS OF THE N-TH DEGREE ORTHOGONAL POLYNOMIAL;
                    (B MAY BE USED FOR ZER, BUT THEN THIS RECURRENCE
                     COEFFICIENTS ARE OVERWRITTEN BY THE ZEROS AND
                     THE ORIGINAL CONTENTS OF B ARE NOT PRESERVED.)
    EM:     <ARRAY IDENTIFIER>;
            "ARRAY" EM[0:5];
            ENTRY:  EM[0]:  THE MACHINE PRECISION;
                    EM[2]:  THE RELATIVE TOLERANCE OF THE ZEROS;
                    EM[4]:  THE MAXIMUM ALLOWED NUMBER OF ITERATIONS,
                            E.G. 5 * N;
            EXIT:   EM[1]: THE MAXIMUM OF ABS(B[0])+1,C[I]+ABS(B[I])+1,
                            (I=1,...N-2),C[N-1]+ABS(B[N-1]);
                    EM[3]:  INFORMATION CONCERNING THE PROCESS USED;
                            I.E. THE MAXIMUM ABSOLUTE VALUE OF THE
                            CODIAGONAL ELEMENTS NEGLECTED,
                            (SEE ALSO SECTION 3.3.1.1.1.);
                    EM[5]:  THE NUMBER OF ITERATIONS PERFORMED.

PROCEDURES USED:

    QRIVALSYMTRI    =  CP34160,
    DUPVEC          =  CP31030.

METHOD AND PERFORMANCE : SEE SELZERORTPOL (THIS SECTION).

EXAMPLE OF USE:

    AS A FORMAL TEST OF THE PROCEDURE WE CALCULATE THE ZEROS OF
    THE CHEBYSHEV POLYNOMIAL (OF THE FIRST KIND) OF THE THIRD DEGREE.
    THE RECURRENCE COEFFICIENTS ARE:
    B[I] = 0, I = 0, 1, ....;
    C[0] = 0, C[1] = .5, C[I] = .25 , I = 2, 3, .....
    (IT IS RECOMMENDED TO STORE THE ELEMENTS OF THE ARRAYS B AND C IN
     REVERSED ORDER IF THESE ELEMENTS ARE STRONGLY INCREASING).

"BEGIN" "ARRAY" B, C[0:3], ZER[1:3], EM[0:5];
    EM[0]:= EM[2]:= "-14; EM[4]:=15;
    B[2]:=B[1]:=B[0]:=0;
    C[0]:= 0; C[1]:= .5; C[2]:= .25;
    ALLZERORTPOL (3, B, C, ZER, EM);
    OUTPUT(61,"(""("THE THREE ZEROS:")",/,3(/ZD5B,N),2/,
              "("EM[1]:")",5BD.2D"+2D, /,
              "("EM[3]:")",5BD.2D"+2D, /,"("EM[5]:")",5ZD")",
              1,ZER[1],2,ZER[2],3,ZER[3],EM[1],EM[3],EM[5])
"END"

THE THREE ZEROS:

 1     -8.6602540378444"-001
 2     +8.6602540378444"-001
 3     -1.0000000000000"-014

EM[1]:     1.50"+00
EM[3]:     7.07"-15
EM[5]:     1


SUBSECTION:         LUPZERORTPOL.

CALLING SEQUENCE:

    THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:

    "PROCEDURE" LUPZERORTPOL (N, M, B, C, ZER, EM);
    "VALUE" N, M; "INTEGER" N, M; "ARRAY" B, C, ZER, EM;
    "CODE" 31363;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:      <ARITHMETIC EXPRESSION>;
            ENTRY:  THE DEGREE OF THE ORTHOGONAL POLYNOMIAL OF WHICH
                    THE ZEROS ARE TO BE CALCULATED;
    M:      <ARITHMETIC EXPRESSION>;
            ENTRY:  THE NUMBER OF ZEROS TO BE CALCULATED;
    B, C:   <ARRAY IDENTIFIER>;
            "ARRAY" B, C [0:N-1];
            ENTRY:  THE ELEMENTS B[I] AND C[I], I = 0, 1, ... , N-1,
                    CONTAIN THE COEFFICIENTS OF THE RECURRENCE RELATION
                    P[I+1](X) = (X - B[I]) * P[I](X) - C[I] * P[I-1](X)
                    I = 0, 1, ... , N-1;
                    ASSUMED IS C[0]=0, WHILE THE CONTENTS
                    OF THE ARRAYS ARE NOT PRESERVED;
    ZER:    <ARRAY IDENTIFIER>;
            "ARRAY" ZER[1:M];
            EXIT:   THE M LOWEST ZEROS ARE DELIVERED;
                    IF HOWEVER THE ARRAY B[0:N-1] CONTAINED THE OPPOSIT
                    VALUES OF THE CORRESPONDING RECURRENCE COEFFICIENTS
                    THEN THE OPPOSITE VALUES OF THE M UPPER ZEROS
                    ARE DELIVERED.
                    IN EITHER CASE, ZER[I]<ZER[I+1], I = 1,...,M-1;
    EM:     <ARRAY IDENTIFIER>;
            "ARRAY" EM[0:6];
            ENTRY:  EM[0]:  THE MACHINE PRECISION.
                    EM[2]:  THE RELATIVE TOLERANCE OF THE ZEROS;
                    EM[4]:  THE MAXIMUM ALLOWED NUMBER OF ITERATIONS,
                            E.G. 15 * M;
                    EM[6]:  IF ALL ZEROS ARE KNOWN TO BE POSITIVE
                            THEN 1 ELSE 0;
            EXIT:   EM[1]:  THE MAXIMUM OF ABS(B[0]) + 1 ,
                            C[I]+ABS(B[I])+1,(I=1,...N-2),
                            C[N-1]+ABS(B[N-1]);
                    EM[3]:  INFORMATION CONCERNING THE PROCESS USED,
                            I.E. THE MAXIMUM ABSOLUTE VALUE OF THE
                            THEORETICAL ERRORS OF THE ZEROS
                            (SEE WILKINSON AND REINSCH, 1971, P.263);
                    EM[5]:  THE NUMBER OF ITERATIONS PERFORMED.

PROCEDURES USED:

    DUPVEC          =  CP31030,
    INFNRMVEC       =  CP31061.

METHOD AND PERFORMANCE : SEE SELZERORTPOL (THIS SECTION).

EXAMPLE OF USE:

    WE CALCULATE THE TWO LOWER AND THE TWO UPPER ZEROS OF THE
    LAGUERRE POLYNOMIAL OF THE THIRD DEGREE.
    THE RECURRENCE COEFFICIENTS ARE OBTAINED FROM [1], P.782:
    B[I] = - A2I / A3I = 2I + 1;
    C[I] = A4I / (A3I * A3(I-1)) * A1(I-1) = I * I, I = 0, 1, 2.
    (IT IS RECOMMENDED TO STORE THE ELEMENTS OF THE ARRAYS B AND C IN
     REVERSED ORDER IF THESE ELEMENTS ARE STRONGLY DECREASING).

    "BEGIN" "ARRAY" B, C[0:3], ZER[1:2], EM[0:6];
    "INTEGER" I;
    EM[0]:= EM[2]:= "-14; EM[4]:= 45; EM[6]:= 1;
    "FOR" I:= 0, 1, 2 "DO"
    "BEGIN" B[I]:= 2 * I + 1; C[I]:= I * I "END";
    LUPZERORTPOL (3, 2, B, C, ZER, EM);
    OUTPUT(61,"(""("THE TWO LOWER ZEROS:")",/,2(/ZD5B,N),2/,
              "("EM[1]:")",5BD.2D"+2D, /,
              "("EM[3]:")",5BD.2D"+2D, /,"("EM[5]:")",5ZD")",
              1,ZER[1],2,ZER[2],EM[1],EM[3],EM[5]);
    EM[6]:= 0;
    "FOR" I:= 0, 1, 2 "DO"
    "BEGIN" B[I]:= - 2 * I - 1; C[I]:= I * I "END";
    LUPZERORTPOL (3, 2, B, C, ZER, EM);
    OUTPUT(61,"("3/,"("THE TWO UPPER ZEROS:")",/,2(/ZD5B,N),2/,
              "("EM[1]:")",5BD.2D"+2D, /,
              "("EM[3]:")",5BD.2D"+2D, /,"("EM[5]:")",5ZD")",
              1,-ZER[1],2,-ZER[2],EM[1],EM[3],EM[5]);
    "END"

    THE TWO LOWER ZEROS:

     1     +4.1577455678348"-001
     2     +2.2942803602791"+000

    EM[1]:     9.00"+00
    EM[3]:     5.72"-16
    EM[5]:    12

    THE TWO UPPER ZEROS:

     1     +6.2899450829375"+000
     2     +2.2942803602791"+000

    EM[1]:     9.00"+00
    EM[3]:     4.70"-20
    EM[5]:    14


SUBSECTION:         SELZERORTPOL.

CALLING SEQUENCE:

    THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:

    "PROCEDURE" SELZERORTPOL (N, N1, N2, B, C, ZER, EM);
    "VALUE" N, N1, N2; "INTEGER" N, N1, N2; "ARRAY" B, C, ZER, EM;
    "CODE" 31364;

    THE MEANING OF THE FORMAL PARAMETERS IS :
    N,N1,N2:<ARITHMETIC EXPRESSION>;
            ENTRY:  N IS THE DEGREE OF THE ORTHOGONAL POLYNOMIAL OF
                    WHICH THE N1-ST UP TO AND INCLUDING N2-ND ZEROS ARE
                    TO BE CALCULATED(ZER[N1]>=ZER[N2]);
    B, C:   <ARRAY IDENTIFIER>;
            "ARRAY" B,C [0 : N-1];
            ENTRY:  THE ELEMENTS B[I] AND C[I], I = 0, 1, ... , N-1,
                    CONTAIN THE COEFFICIENTS OF THE RECURRENCE RELATION
                    P[I+1](X) = (X - B[I]) * P[I](X) - C[I] * P[I-1](X)
                    I = 0, 1, ... , N-1,
                    ASSUMED IS C[0]=0, WHILE THE CONTENTS
                    OF THE ARRAYS IS PRESERVED;
    ZER:    <ARRAY IDENTIFIER>;
            "ARRAY" ZER [N1:N2];
            EXIT:   THE N2-N1+1 CALCULATED ZEROS IN DECREASING ORDER.
    EM:     <ARRAY IDENTIFIER>;
            "ARRAY" EM[0:5];
            ENTRY:  EM[0]:  THE MACHINE PRECISION.
                    EM[2]:  THE RELATIVE TOLERANCE OF THE ZEROS;
            EXIT:   EM[1]:  THE MAXIMUM OF ABS(B[0]) + 1,
                            C[I]+ABS(B[I])+1 (I=1,...N-2) AND
                            C[N-1]+ABS(B[N-1]);
                    EM[5]:  THE NUMBER OF ITERATIONS PERFORMED.

PROCEDURES USED:

   VALSYMTRI       =  CP34151.

METHOD AND PERFORMANCE:

    THE ZEROS OF AN ORTHOGONAL POLYNOMIAL ARE THE EIGENVALUES OF A
    SYMMETRIC TRIDIAGONAL MATRIX (SEE GOLUB AND WELSCH (1969),
    LANCZOS (1957),P.375,376, STOER (1972),P.120).
    THE ORTHOGONAL POLYNOMIAL IS DEFINED BY A LINEAR THREE-TERM
    HOMOGENEOUS RECURRENCE RELATION.

EXAMPLE OF USE :

    WE CALCULATE THE THIRD ZERO OF THE LEGENDRE POLYNOMIAL OF THE
    FOURTH DEGREE. THE RECURRENCE COEFFICIENTS ARE OBTAINED FROM
    ABRAMOWITZ AND STEGUN (1964),P.782:
    B[I] = 0, I = 0, 1, ....;
    C[I] = A4I / (A3I * A3(I-1)) * A1(I-1) = I * I / ( 4 * I * I - 1),
            I = 0, 1, .....
    (IT IS RECOMMENDED TO STORE THE ELEMENTS OF THE ARRAYS B AND C IN
     REVERSED ORDER IF THESE ELEMENTS ARE STRONGLY DECREASING).

    "BEGIN" "ARRAY" B, C[0:4], ZER[3:3], EM[0:5];
    "INTEGER" I;
    EM[0]:= EM[2]:= "-14;
    "FOR" I:= 0, 1, 2, 3 "DO"
    "BEGIN" B[I]:= 0; C[I]:= I * I / (4 * I * I - 1) "END";
    SELZERORTPOL (4, 3, 3, B, C, ZER, EM);
    OUTPUT(61,"(""("THE THIRD ZERO:")",2/,ZD5B,N,2/,
              "("EM[1]:")",5BD.2D"+2D, /,
              "("EM[5]:")",5ZD")",3,ZER[3],EM[1],EM[5])
    "END"

    THE THIRD ZERO:

     3     -3.3998104358486"-001

    EM[1]:     1.33"+00
    EM[5]:    12


SUBSECTION: ALL JAC ZER.

CALLING SEQUENCE:

    THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:

    "PROCEDURE" ALL JAC ZER(N, ALFA, BETA, ZER);
    "VALUE" N, ALFA, BETA;
    "INTEGER" N; "REAL" ALFA, BETA; "ARRAY" ZER;
    "CODE" 31370;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:      <ARITHMETIC EXPRESSION>;
            THE UPPER BOUND OF THE ARRAY ZER; N >= 1;
    ALFA,BETA:
            <ARITHMETIC EXPRESSION>;
            THE PARAMETERS OF THE JACOBI POLYNOMIAL,
            SEE ABRAMOWITZ AND STEGUN (1964);
            ALFA, BETA > - 1;
    ZER:    <ARRAY IDENTIFIER>;
            "ARRAY" ZER[1 : N];
            EXIT: ZER[1], ..., ZER[N] ARE THE ZEROS OF THE N-TH
                  JACOBI POLYNOMIAL WITH PARAMETERS ALFA AND BETA.

PROCEDURES USED:

    ALL ZER ORT POL = CP 31362.

REQUIRED CENTRAL MEMORY:

    IF ALFA = BETA THEN TWO AUXILIARY ARRAYS OF N//2 REALS ARE
    USED, OTHERWISE TWO AUXILIARY ARRAYS OF N REALS ARE DECLARED.

METHOD AND PERFORMANCE:

    THE JACOBI POLYNOMIALS ARE A SPECIAL CASE OF ORTHOGONAL
    POLYNOMIALS (SEE ABRAMOWITZ AND STEGUN (1964)); ALL JAC ZER
    COMPUTES THE COEFFICIENTS OF THE THREE-TERM RECURRENCE RELATION
    AND CALLS THE PROCEDURE ALL ZER ORT POL TO COMPUTE THE ZEROS;
    IF ALFA=BETA, THE POLYNOMIALS ARE ODD OR EVEN, HENCE ONLY THE
    POSITIVE ZEROS ARE CALCULATED; THIS IS DONE BY MEANS OF THE
    FORMULAS

          P(2*M, ALFA, ALFA, X) = C(M)*P(M, ALFA, -0.5, 2*X*X - 1),

      P(2*M - 1, ALFA, ALFA, X) = D(M)*P(M, ALFA, +0.5, 2*X*X - 1)*X

    (SEE ABRAMOWITZ AND STEGUN (1964), FORMULAS 22.5.20 - 22.5.27).

EXAMPLE OF USE:

    THE PROGRAM

    "BEGIN" "ARRAY" X[1:3];
      ALL JAC ZER(3,-.5,-.5,X);
      OUTPUT(61,"("3(4B-D.13D"-ZD)")",X[1],X[2],X[3])
    "END"

    DELIVERS THE FOLOWING RESULTS:

    -8.6602540378444"-1    0.0000000000000" 0     8.6602540378444"-1


SUBSECTION: ALL LAG ZER.

CALLING SEQUENCE:

    THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:

    "PROCEDURE" ALL LAG ZER(N, ALFA, ZER);
    "VALUE" N, ALFA;
    "INTEGER" N; "REAL" ALFA; "ARRAY" ZER;
    "CODE" 31371;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:      <ARITHMETIC EXPRESSION>;
            THE UPPER BOUND OF THE ARRAY ZER; N >= 1;
    ALFA:   <ARITHMETIC EXPRESSION>;
            THE PARAMETER OF THE LAGUERRE POLYNOMIAL,
            SEE ABRAMOWITZ AND STEGUN (1964); ALFA > -1;
    ZER:    <ARRAY IDENTIFIER>;
            "ARRAY" ZER[1 : N];
            EXIT: ZER[1], ..., ZER[N] ARE THE ZEROS OF THE N-TH
                  LAGUERRE POLYNOMIAL WITH PARAMETER ALFA.

PROCEDURES USED:

    ALL ZER ORT POL = CP 31362.

REQUIRED CENTRAL MEMORY:

    TWO AUXILIARY ARRAYS OF N REALS ARE USED.

METHOD AND PERFORMANCE:

    THE LAGUERRE POLYNOMIALS ARE A SPECIAL CASE OF ORTHOGONAL
    POLYNOMIALS (SEE ABRAMOWITZ AND STEGUN (1964)); ALL LAG ZER
    COMPUTES THE COEFFICIENTS OF THE THREE-TERM RECURRENCE RELATION
    AND CALLS THE PROCEDURE ALL ZER ORT POL TO COMPUTE THE ZEROS.

EXAMPLE OF USE:

    "BEGIN" "ARRAY" X[1:3];
      ALL LAG ZER(3,-.5,X);
      OUTPUT(61,"("3(4B-D.13D"-ZD)")",X[1],X[2],X[3])
    "END"

    DELIVERS THE FOLOWING RESULTS:

    5.5253437422633" 0.    1.7844927485432" 0     1.9016350919350" -1

SOURCE TEXT(S) :
"CODE" 31362;
"CODE" 31363;

"CODE" 31364;

"CODE" 31370;

"CODE" 31371;