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;