NUMAL Section 4.2.3.1

BEGIN SECTION : 4.2.3.1 (November, 1978)

AUTHOR:         C.G. VAN DER LAAN.

CONTRIBUTORS:   C.G. VAN DER LAAN AND M. VOORINTHOLT.

INSTITUTE:      REKENCENTRUM RIJKSUNIVERSITEIT GRONINGEN.

RECEIVED:       780601.

BRIEF DESCRIPTION:

    THIS SECTION CONTAINS THE PROCEDURES GSSWTS, GSSWTSSYM AND RECCOF.
    RECCOF CALCULATES FROM A GIVEN WEIGHT FUNCTION ON [-1,1] THE
    RECURRENCE COEFFICIENTS OF THE CORRESPONDING ORTHOGONAL
    POLYNOMIALS; GSSWTS AND GSSWTSSYM CALCULATE FROM THE RECURRENCE
    COEFFICIENTS THE GAUSSIAN WEIGHTS OF THE CORRESPONDING WEIGHT
    FUNCTION.

KEYWORDS:

    RECURRENCE COEFFICIENTS ORTHOGONAL POLYNOMIALS,
    GAUSSIAN WEIGHTS,
    GAUSSIAN QUADRATURE.


SUBSECTION: RECCOF.

CALLING SEQUENCE:

    THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:

    "PROCEDURE" RECCOF(N,M,X,WX,B,C,L,SYM);
    "VALUE" N,M,SYM; "INTEGER" N,M; "BOOLEAN" SYM;
    "REAL" X,WX; "ARRAY" B,C,L;
    "CODE" 31254;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:       <ARITHMETIC EXPRESSION>;
             ENTRY: UPPER BOUND FOR THE INDICES OF THE ARRAYS B, C, L
                    (N>=0);
    M:       <ARITHMETIC EXPRESSION>;
             ENTRY: THE NUMBER OF POINTS USED IN THE GAUSS-CHEBYSHEV
                    QUADRATURE RULE FOR CALCULATING THE APPROXIMATION
                    OF THE INTEGRAL REPRESENTATIONS OF B[K],C[K]

    SYM:     <BOOLEAN EXPRESSION>;
             ENTRY: "IF" SYM
                        "THEN" WEIGHT FUNCTION ON [-1,1] IS EVEN
                        "ELSE" WEIGHT FUNCTION ON [-1,1] IS NOT EVEN;
    X,WX:    <ARITHMETIC EXPRESSION>;
             ENTRY: JENSEN VARIABLES WITH WX AN EXPRESSION OF X
                    DENOTING THE WEIGHT FUNCTION ON [-1,1];
    B,C,L:   <ARRAY IDENTIFIER>;
             "ARRAY" B,C,L[0:N];
             EXIT:  THE APPROXIMATE RECURRENCE COEFFICIENTS FOR
                    P[K+1](X) = (X-B[K])*P[K](X) - C[K]*P[K-1](X),
                                                     K=0,1,2,...N,
                    AND THE APPROXIMATE SQUARE LENGTHS
                            X = +1
                    L[K] = INTEGRAL ( W(X) * P[K](X) **  2 ) DX
                            X = -1

PROCEDURES USED:    ORTPOL = CP31044.

RUNNING TIME:       PROPORTIONAL TO M*N**2.

METHOD AND PERFORMANCE:

    THE RECURRENCE COEFFICIENTS ARE REPRESENTED BY
                       X = +1
             B[K] = ( INTEGRAL ( W(X) * X * P[K](X) ** 2 ) DX ) / L[K],
                       X = -1

             C[K] = L[K] / L[K-1]],
    WHERE P[K](X) IS THE K-TH ORTHOGONAL POLYNOMIAL.
    THE INTEGRALS ARE APPROXIMATED BY THE M-POINTS GAUSS-CHEBYSHEV
    RULE AS
          X = +1                       M
        INTEGRAL(F(X))DX := PI / M * SUM SIN(THETA[J])*F(COS(THETA[J]))
          X = -1                      J=1
    WITH THETA[J] = (2*J-1) * PI / (2*M)      (SEE GAUTSCHI, 1968A).
    THE VALUE OF M IS TO BE SUPPLIED BY THE USER.

REFERENCES:

    GAUTSCHI, W. (1968A):
    CONTRUCTION OF GAUSS-CHRISTOFFEL FORMULAS.
    MATH. COMP., 22,P.251-270.

    GAUTSCHI, W. (1968B):
    GAUSSIAN QUADRATURE FORMULAS.
    COMM. ACM. CALGO 331.

EXAMPLE OF USE:

    THE FOLLOWING PROGRAM DELIVERS AN APPROXIMATION
    FOR THE RECURSION COEFFICIENTS C[1] AND C[2], OF THE CHEBYSHEV
    POLYNOMIALS OF THE SECOND KIND;

    "BEGIN"
        "REAL" X; "ARRAY" B,C,L[0:2];
        RECCOF(2,200,X,SQRT(1 - X**2),B,C,L,"TRUE");
        OUTPUT(61,"("2/,2(3B,-ZD.3D)")",C[1],C[2]);
    "END";

    RESULTS:

         0.250     0.250


SUBSECTION: GSSWTS.

CALLING SEQUENCE:

    THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:

    "PROCEDURE" GSSWTS(N,ZER,B,C,W);
    "VALUE" N; "INTEGER" N; "ARRAY" ZER,B,C,W;
    "CODE" 31253;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:       <ARITHMETIC EXPRESSION>;
             ENTRY: THE NUMBER OF WEIGHTS TO BE COMPUTED; UPPER
                    BOUND FOR THE ARRAYS Z AND W (N>=1);
    ZER:     <ARRAY IDENTIFIER>;
             "ARRAY" ZER[1:N];
             ENTRY: THE ZEROS OF THE N-TH DEGREE ORTHOGONAL POLYNOMIAL;
    B,C:     <ARRAY IDENTIFIER>;
             "ARRAY" B[0:N-1], C[1:N-1];
             ENTRY: THE RECURRENCE COEFFICIENTS;
    W:       <ARRAY IDENTIFIER>;
             "ARRAY" W[1:N];
             EXIT : THE GAUSSIAN WEIGHTS DIVIDED BY THE
                    INTEGRAL OVER THE WEIGHT FUNCTION.

PROCEDURES USED:    ALLORTPOL = CP 31045.

METHOD AND PERFORMANCE:

    THE GAUSSIAN WEIGHTS OF AN N-POINTS RULE DIVIDED BY THE INTEGRAL
    OF THE WEIGHT FUNCTION MAY BE REPRESENTED AS
            W[K] = 1/(...((P[N-1](Z)**2/C[N-1]+P[N-2](Z)**2)/C[N-2]+...
                                 ...+P[1](Z)**2)/C[1]+1)  , K=1,2,...N
    WITH    Z= K-TH ZERO OF P[N](X).  (GAUTSCHI, 1970).

    ALLZERORTPOL AND GSSWTS MAY BE USED TO GENERATE GAUSSIAN
    QUADRATURE RULES PROVIDED THE RECURRENCE COEFFICIENTS AND THE
    INTEGRAL OF THE WEIGHT FUNCTION ARE KNOWN.
    FOR EXAMPLE THE GAUSS-LAGUERRE QUADRATURE RULE APPLIED TO F
    MAY BE OBTAINED BY THE CALLS

             "FOR" K:=1 "STEP" 1 "UNTIL" N-1 "DO"
             "BEGIN"
                B[K]:=2*K+ALPHA+1;
                C[K]:=K*(K+ALPHA)
             "END";
             B[0]:=ALPHA+1;
             ALLZERORTPOL(N,ZER,B,C);
             GSSWTS(N,ZER,B,C,W);
             GAUSSRULE:=0;
             "FOR" K:=1 "STEP" 1 "UNTIL" N "DO"
             GAUSSRULE:=GAUSSRULE+W[K]*F(ZER[K]);
             GAUSSRULE:=GAUSSRULE*GAMMA(ALPHA+1)

    GAUSSRULE CONTAINS THE VALUE OF THE APPOXIMATING GAUSS
    QUADRATURE RULE AND ZER[1:N],W[1:N] CONTAIN THE GAUSSIAN
    ABSCISSA AND WEIGHTS.

    IN THE FOLLOWING TABEL WE SUMMARIZE CLASSICAL QUADRATURE RULES

           :  WEIGHT   :    RECURRENCE COEFFICIENTS    :   INTEGRAL
 GAUSSIAN  : FUNCTION  :---------------:---------------:      OF
QUADRATURE :   W(X)    :      B[K]     :     C[K]      :WEIGHT FUNCTION
-----------:-----------:---------------:---------------:---------------
           :           :               :               :
LEGENDRE   :     1     :       0       :K**2*(4*K**2-1):      2
           :           :               :               :
 CHEBYSHEV : 1/SQRT(1- :       0       :  1/2  ,  K=1  :      PI
(1-ST KIND): X**2)     :               :  1/4  ,  K>1  :
           :           :               :               :
 CHEBYSHEV : SQRT(1-   :       0       :      1/4      :    PI/2
(2-ND KIND): X**2)     :               :               :
           :           :               :               :
JACOBI     : (1-X)**   : -(ALPHA**2-   : 4*(1+ALPHA)*  : 2**(ALPHA+
           : ALPHA*(1+ : BETA**2)/((2* : (1+BETA)/     : BETA+1)*
           : X)**BETA  : K+ALPHA+BETA)*: ((ALPHA+BETA+ : GAMMA(ALPHA+
           :           : (2*K+ALPHA+   : 2)**2*(ALPHA+ : 1)*GAMMA(BETA+
           :           : BETA+2))      : BETA+3)) ,K=1 : 1)/GAMMA(ALPHA
           :           :               :               : +BETA+2)
           :           :               : 4*K*(K+ALPHA)*:
           :           :               : (K+BETA)*(K+  :
           :           :               : ALPHA+BETA)/  :
           :           :               : ((2*K+ALPHA+  :
           :           :               : BETA)**2*     :
           :           :               : ((2*K+ALPHA+  :
           :           :               : BETA)**2-1))  :
           :           :               :          ,K>1 :
           :           :               :               :
           :           :               :(ALPHA,BETA>-1):
           :           :               :               :
LAGUERRE   : EXP(-X)*  : 2*K+ALPHA+1   : K*(K+ALPHA)   : GAMMA(ALPHA+1)
           : X**ALPHA  :               :               :
           :           :               :               :
HERMITE    : EXP(-X**2):       0       :      K/2      :   SQRT(PI)

(THE INTEGRATION INTERVALS ARE:  [-INFINITY,INFINITY]  FOR HERMITE;
                                         [0,INFINITY]  FOR LAGUERRE;
                                        [-1,1]         FOR THE OTHERS.)
FOR NON-CLASSICAL WEIGHT FUNCTIONS ON A FINITE INTERVAL THE RECURSION
COEFFICIENTS (AND THE SQUARE LENGTHS OF THE CORRESPONDING ORTHOGONAL-
POLYNOMIALS) CAN BE OBTAINED BY THE PROCEDURE RECCOF (THIS SECTION).

REFERENCES:

    GAUTSCHI, W. (1970):
    GENERATION OF GAUSSIAN QUADRATURE RULES AND
    ORTHOGONAL POLYNOMIALS.
    IN:  COLLOQUIUM APPROXIMATIETHEORIE,
         MC SYLLABUS 14.

EXAMPLE OF USE: SEE SUBSECTION GSSWTSSYM.


SUBSECTION: GSSWTSSYM.

CALLING SEQUENCE:

    THE DECLARATION OF THE PROCEDURE IN THE CALLING PROGRAM READS:

    "PROCEDURE" GSSWTSSYM(N,ZER,C,W);
    "VALUE" N; "INTEGER" N; "ARRAY" ZER,C,W;
    "CODE" 31252;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:       <ARITHMETIC EXPRESSION>;
             ENTRY: THE WEIGHTS OF AN N-POINTS GAUSS RULE ARE
                    TO BE CALCULATED (BECAUSE OF SYMMETRY ONLY
                    (N+1)//2 OF THE VALUES ARE DELIVERED);
    ZER:     <ARRAY IDENTIFIER>;
             "ARRAY" ZER[1:N//2]
             ENTRY: THE NEGATIVE ZEROS OF THE N-TH DEGREE ORTHOGONAL
                    POLYNOMIAL (ZER[I] < ZER[I+1], I=1,2,...,N//2-1);
             (IF N IS ODD THEN 0 IS ALSO A ZERO.)
    C:       <ARRAY IDENTIFIER>;
             "ARRAY" C[1:N-1];
             ENTRY: THE RECURRENCE COEFFICIENTS;
    W:       <ARRAY IDENTIFIER>;
             "ARRAY" W[1:(N+1)//2];
             EXIT : PART OF THE GAUSSIAN WEIGHTS DIVIDED BY THE
                    INTEGRAL OF THE WEIGHT FUNCTION.
             (NOTE THAT W[N+1-K]=W[K]  , K=1,2,...,(N+1)//2.)

PROCEDURES USED:    ALLORTPOLSYM = CP 31049.

METHOD AND PERFORMANCE: SEE SUBSECTION GSSWTS; THIS PROCEDURE IS
                        SUPPLIED FOR STORAGE ECONOMICAL REASONS.

REFERENCES:  SEE SUBSECTION GSSWTS.

EXAMPLE OF USE:

    THE FOLLOWING PROGRAM DELIVERS THE GAUSSIAN WEIGHTS
    FOR THE 5-POINTS GAUSS-CHEBYSHEV QUADRATURE RULE BY MEANS OF
    THE PROCEDURE GSSWTSSYM (C[1]=0.5; C[K]=0.25, K=2,3,...;
    ZER[I] = COS((2*(N-I) - 1) / (2*N) * PI), I=1,2,...,N//2.

    "BEGIN"
        "REAL" PI; "INTEGER" I;
        "ARRAY" ZER[1:2], W[1:3], C[1:4];
        PI:=4*ARCTAN(1);
        C[1]:=.5;
        "FOR" I:=2 "STEP" 1 "UNTIL" 4 "DO"
        C[I]:=.25;
        ZER[1]:=COS(.9*PI);
        ZER[2]:=COS(.7*PI);
        GSSWTSSYM(5,ZER,C,W);
        OUTPUT(61,"("2/,5(3B,-ZD.3D)")",W[1]*PI,W[2]*PI,W[3]*PI,
                                                W[2]*PI,W[1]*PI);
    "END";

    RESULTS:
         0.628     0.628     0.628     0.628     0.628

SOURCE TEXT(S):
"CODE" 31254;
"CODE" 31253;

"CODE" 31252;