0 - Semantic Scholar
CS 81
.
..
CALCULATION OF GAUSS QUADRATURE RULES BY GENE H. GOLUB AND JOHN H. WELSCH
TECHNICAL REPORT NO. CS 81 NOVEMBER 3, 1967
COMPUTER SC IENCE DEPARTMENT School of Humanities and Sciences STANFORD UNIVERSITY
CALCULATION OF GAUSS QUADRATURE RULES*
Gene H. Golub and John H. Welsch
Introduction Most numerical integration techniques consist of approximating the integrand by a polynomial in a region or regions and then integrating the polynomial exactly.
Often a complicated integrand can be factored
into a non-negative=."weight" function and another function better approximated by a polynomial, thus
= s; cu(t)f(t)dtm f wif(ti) o i=l Hopefully, the quadrature rule
corresponding to the weight
fw
function cu(t) is available in tabulated form, but more likely it is not0
We present here two algorithms for generating the Gaussian quadra-
ture rule defined by the weight function when: 4 .
the three term recurrence relation is known for the orthogonal polynomials generated by cu(t), and
b) the moments of the weight function are known or can be calculated.
*The work of the first author was in part supported by the Office of Naval Research and the National Science Foundation; the work of the second author was‘in part supported by the Atomic Energy Commission. 1
1.
Definitions and Preliminaries Let w(x) 10 be a fixed weight function defined on [a, b] . For
dx),
it is possible to define a sequence-. of polynomials
which are orthonormal with respect to U(X) exact degree n
P,(X), P,(X), * *
and in which p,(x) is of
so that
J'; w(x)pm(x)pn(x)dx =
when m = n when m # n
l
The polynomial --. P,(X) = kn
n" Cxuti) 9 i=l
kn ’ 0,
has n real roots
a < tl < t2
~jY
I.
"j+l-
'j+2 .
j = 1, 2, . . . . N-l .
1
It is shown in [l] that
xPj(x)--. = pj-lpjml (x> +
+
ajPj(x)
for j = 1, 2,
l
*
‘,
BjPj+l(x)
(4.1)
N
where F a. = 'D(Fj-1 - Fj-2 D 3 j-l j-2'
0 -1
j = 1, 2, . . . . N
-1
j = 1, 2, . . . . N-l O
= 1, Do = II,)
Note that the tri-diagonal matrix so generated is symmetric. 1 In [7], Rutishauser gives explicit formulas in terms of determinants for the Gauss-Crout decomposition of a matrix.
We may use these relation-
ships to evalute the coefficients N 1 a.3 J j=l'
N-l CB . I
.
J j=l
Let R denote the Cholesky decomposition of M so that
T M =RR
and R is an upper triangular matrix whose elements are ..
r
ii
= ( mii - % rEi)' i=l
and i-l m ij - k=l c
r
for i and j
(4.2)
i < j,
between 1 and n .
--.
Then, from the formulas of Rutishauser
cl. =
3
'j,j+L
-
rj-l,
j
3 = 1, 2, . . . . N
rj-l,j-l
'j,j
7
(4.3) rtj+l,cj+l Bj
3 = 1, 2, . . . . N-l
=
J
'j,j
with r 0'0 = l'
rO,l =
There are other means for evaluating
0 .
(a! ), (pj] but it is the j
opinion of the authors that the above method will lead to the most accurate formulas.
50
Description of Computational Procedures In the following section there are three ALGOL 60 procedures for
performing the algorithms presented above.
We have tried to keep the
identifiers as close to the notation of the equations as possible without
14
.
sacrificing storage or efficiency.
The weights and abscissas of the
quadrature rule are the result of the procedure GAUSSQUADRULE which must be supplied with the recurrence relation by either procedure .. GENORTHOPOLY or procedure CLASSICORTHOPOLY . The former requires the moments of the weight function and the latter the name of the particular orthogonal polynomial.
A short description of each procedure follows.,
CLASSICORTHOPOLY produces ~0 and the three term recurrence relationship (a., b cj) J j' polynomials: KING = 1,
for six well-known kinds of orthogonal
Legendre Polynomials
P,(x)
on
[-1.0, +l.Ol,
44 = 1.0 . KIND = 2,
Chebyshev Polynomials of the first kind T,(x) on [-1.0, +1.01, CD(x) = (1-x2)-$
KIND=3,
l
Chebyshev Polynomials of the second kind V,(x) on c-1.0, +1.0], w(x) = (l-x2)+& .
*
KIND = 4, Jacobi Polynomials P(a'p)(x) n 4x1 = (l-x)a(l+x)@ IKIND = 5,
for cx > -1 and @ > -1 .
( 1(x) Laguerre Polynomials L," w ( x ) = ewxxa
KIND = 6, Hermite
on [-1.0, +l.O]'
on LO, +4,
for a > -1 .
Polynomials
2 H,(x) on [-a, +a], m(x) = e-x o
Notice that this procedure requires a real procedure to evaluate the
gamma function P(X) .
15
GENORTHOPOLY uses the 2N+l moments of the weight function which are supplied in JfUXl through Md=!@N] of formula (4.1).
to compute the CX~'s and p.'s J First, The Cholesky decomposition (formula 4.2) of the ._
moment matrix is placed in the upper right triangular part of the array R’
then the formulas (4.3) are used to compute the
which are placed in the arrays
cx 's and p.'s J 3 A and B respectively.
GAUSSQUADRULE has two modes of operation controlled by the Boolean parameter
SYMM which indicates whether the tri-diagonal matrix is
symmetric or not.
When the recursion relation is produced by GENORTHOPOLY,
SYMM is true; when--_ produced by CLASSICQRTHQPOLY, SYMM is false. If SYMM is false, the matrix is symmetricized using the formulas (2.2). The diagonal elements ai 'i
are stored in A[I] and the off diagonal elements
are stored in B[I] . Beginning at label SETUP,
are done:
the Ll
several calculations and initializations
norm of the tri-diagonal matrix and the relative zero
tolerance are computed; the first component of each eigenvector W[I] and the Q-R iteration are initialized.
LAMBDA is a variable subtracted off
the diagonal elements to accelerate convergence of the Q-R iteration and control to some extent in what order the eigenvalues (abscissas) are found.
It begins with a value outside and to the right of the interval
containing the abscissas
(=NORM)
and moves to the left as the abscissas
are found; thus the abscissas will be in ascending order in the array T (just to be sure an exchange sort is used at label S@RT ). The maximum (EIGMAX) of the lower
of the eigenvalues
( LAMBDA1 and LAMBDA2 )
2 X 2 submatrix is compared to the maximum (RHO) from
the last iteration. If they are clase,
16
LAMBDA is replaced by EIGMAX .
’. onu[i-11; ger;orthopoly(n, mu, a, b); gaussquadrule(n, a, b, c, muzero, true, t, w);
26
nutstring(l,
‘
abscissas:'
outstring(l, c weights: 9
); out.array (1, t>;
); outarray: (1, w);
end;
27
w 0 0
r--!-lr-lri .+ -1. + + 0 0 0 0
8. References [l] N, Ahiezer and M, Krein, Some Questions -in the - Theory of Moments, Translations of MathematixMono.graph, Vol, 2, American Mathematical Society, Providence, 1962, Article VI. [23 P, Concus, D0 Cassatt, G. Jaehnig and E. Melky, "Tables for the evaluation of
s
co xBewx f(x)dx
0
by Gauss-Laguerre Quadrature," Math, Corny., 17 (196~)~ ppO 1~45-256~
[3] P, Davis and I0 Polonsky, "Numerical interpolation,differentiatlon, and integration," Handbook of Mathematical Functions, (edited 'by M0 Abramowitz --. and I. Stegunr National Bureau of Standards, Applied Mathematics Series 55, U, S, Government Printing Office, Washington, D, C,, ppe 875-924, f4] P., Davis and P. Rabinowitz, Numerical Integration, Blaisdell Publishing Company, Waltham, Mass.,-
[5] J0 Francis, "The Q-R transformation. A unitary analogue to the L-R transformation," Comput. J., 4 (1961, 1962), ppe 265-271, 3~2,~ 3450
[6] U0 Hochstrasser, "Orthogonal polynomials, "Handbook of Mathematical Functions, (edited by M. Abramowitz and I. 'm ztional Bureau of Standards, Applied Mathematics Series 55, U. S, Government Printing Office, Washington, D. C., ppO 771-802. H0 Rutishauser, "Solution of eigenvalue problems with the L-R transformation," Further Contributions & the Solution of Simultaneous Linear Equations and the - Determination of Eigenvalues, Et,ional Bureau of Standards, Applied Mathematics Series -9, U. So Government Printi.ng Office, Washington, De C., ppO 47-81, [8] H, Rutishauser, "On a modification of the Q-D algorithm with Graeffe!type convergence," zm' 13 O-963), PP~ 493-496.
[9] A. Stroud and D. Secrest, Gaussian Quadr,ature Formulas, Prentice-Hall, Englewood Cliffs, N. J., lr [lOI H0 Wilf, Mathematics -for the Physical Sciences, John Wiley and Sons, Inc, New York, 1962.
29
.
..
CALCULATION OF GAUSS QUADRATURE RULES BY GENE H. GOLUB AND JOHN H. WELSCH
TECHNICAL REPORT NO. CS 81 NOVEMBER 3, 1967
COMPUTER SC IENCE DEPARTMENT School of Humanities and Sciences STANFORD UNIVERSITY
CALCULATION OF GAUSS QUADRATURE RULES*
Gene H. Golub and John H. Welsch
Introduction Most numerical integration techniques consist of approximating the integrand by a polynomial in a region or regions and then integrating the polynomial exactly.
Often a complicated integrand can be factored
into a non-negative=."weight" function and another function better approximated by a polynomial, thus
= s; cu(t)f(t)dtm f wif(ti) o i=l Hopefully, the quadrature rule
corresponding to the weight
fw
function cu(t) is available in tabulated form, but more likely it is not0
We present here two algorithms for generating the Gaussian quadra-
ture rule defined by the weight function when: 4 .
the three term recurrence relation is known for the orthogonal polynomials generated by cu(t), and
b) the moments of the weight function are known or can be calculated.
*The work of the first author was in part supported by the Office of Naval Research and the National Science Foundation; the work of the second author was‘in part supported by the Atomic Energy Commission. 1
1.
Definitions and Preliminaries Let w(x) 10 be a fixed weight function defined on [a, b] . For
dx),
it is possible to define a sequence-. of polynomials
which are orthonormal with respect to U(X) exact degree n
P,(X), P,(X), * *
and in which p,(x) is of
so that
J'; w(x)pm(x)pn(x)dx =
when m = n when m # n
l
The polynomial --. P,(X) = kn
n" Cxuti) 9 i=l
kn ’ 0,
has n real roots
a < tl < t2
~jY
I.
"j+l-
'j+2 .
j = 1, 2, . . . . N-l .
1
It is shown in [l] that
xPj(x)--. = pj-lpjml (x> +
+
ajPj(x)
for j = 1, 2,
l
*
‘,
BjPj+l(x)
(4.1)
N
where F a. = 'D(Fj-1 - Fj-2 D 3 j-l j-2'
0 -1
j = 1, 2, . . . . N
-1
j = 1, 2, . . . . N-l O
= 1, Do = II,)
Note that the tri-diagonal matrix so generated is symmetric. 1 In [7], Rutishauser gives explicit formulas in terms of determinants for the Gauss-Crout decomposition of a matrix.
We may use these relation-
ships to evalute the coefficients N 1 a.3 J j=l'
N-l CB . I
.
J j=l
Let R denote the Cholesky decomposition of M so that
T M =RR
and R is an upper triangular matrix whose elements are ..
r
ii
= ( mii - % rEi)' i=l
and i-l m ij - k=l c
r
for i and j
(4.2)
i < j,
between 1 and n .
--.
Then, from the formulas of Rutishauser
cl. =
3
'j,j+L
-
rj-l,
j
3 = 1, 2, . . . . N
rj-l,j-l
'j,j
7
(4.3) rtj+l,cj+l Bj
3 = 1, 2, . . . . N-l
=
J
'j,j
with r 0'0 = l'
rO,l =
There are other means for evaluating
0 .
(a! ), (pj] but it is the j
opinion of the authors that the above method will lead to the most accurate formulas.
50
Description of Computational Procedures In the following section there are three ALGOL 60 procedures for
performing the algorithms presented above.
We have tried to keep the
identifiers as close to the notation of the equations as possible without
14
.
sacrificing storage or efficiency.
The weights and abscissas of the
quadrature rule are the result of the procedure GAUSSQUADRULE which must be supplied with the recurrence relation by either procedure .. GENORTHOPOLY or procedure CLASSICORTHOPOLY . The former requires the moments of the weight function and the latter the name of the particular orthogonal polynomial.
A short description of each procedure follows.,
CLASSICORTHOPOLY produces ~0 and the three term recurrence relationship (a., b cj) J j' polynomials: KING = 1,
for six well-known kinds of orthogonal
Legendre Polynomials
P,(x)
on
[-1.0, +l.Ol,
44 = 1.0 . KIND = 2,
Chebyshev Polynomials of the first kind T,(x) on [-1.0, +1.01, CD(x) = (1-x2)-$
KIND=3,
l
Chebyshev Polynomials of the second kind V,(x) on c-1.0, +1.0], w(x) = (l-x2)+& .
*
KIND = 4, Jacobi Polynomials P(a'p)(x) n 4x1 = (l-x)a(l+x)@ IKIND = 5,
for cx > -1 and @ > -1 .
( 1(x) Laguerre Polynomials L," w ( x ) = ewxxa
KIND = 6, Hermite
on [-1.0, +l.O]'
on LO, +4,
for a > -1 .
Polynomials
2 H,(x) on [-a, +a], m(x) = e-x o
Notice that this procedure requires a real procedure to evaluate the
gamma function P(X) .
15
GENORTHOPOLY uses the 2N+l moments of the weight function which are supplied in JfUXl through Md=!@N] of formula (4.1).
to compute the CX~'s and p.'s J First, The Cholesky decomposition (formula 4.2) of the ._
moment matrix is placed in the upper right triangular part of the array R’
then the formulas (4.3) are used to compute the
which are placed in the arrays
cx 's and p.'s J 3 A and B respectively.
GAUSSQUADRULE has two modes of operation controlled by the Boolean parameter
SYMM which indicates whether the tri-diagonal matrix is
symmetric or not.
When the recursion relation is produced by GENORTHOPOLY,
SYMM is true; when--_ produced by CLASSICQRTHQPOLY, SYMM is false. If SYMM is false, the matrix is symmetricized using the formulas (2.2). The diagonal elements ai 'i
are stored in A[I] and the off diagonal elements
are stored in B[I] . Beginning at label SETUP,
are done:
the Ll
several calculations and initializations
norm of the tri-diagonal matrix and the relative zero
tolerance are computed; the first component of each eigenvector W[I] and the Q-R iteration are initialized.
LAMBDA is a variable subtracted off
the diagonal elements to accelerate convergence of the Q-R iteration and control to some extent in what order the eigenvalues (abscissas) are found.
It begins with a value outside and to the right of the interval
containing the abscissas
(=NORM)
and moves to the left as the abscissas
are found; thus the abscissas will be in ascending order in the array T (just to be sure an exchange sort is used at label S@RT ). The maximum (EIGMAX) of the lower
of the eigenvalues
( LAMBDA1 and LAMBDA2 )
2 X 2 submatrix is compared to the maximum (RHO) from
the last iteration. If they are clase,
16
LAMBDA is replaced by EIGMAX .
’. onu[i-11; ger;orthopoly(n, mu, a, b); gaussquadrule(n, a, b, c, muzero, true, t, w);
26
nutstring(l,
‘
abscissas:'
outstring(l, c weights: 9
); out.array (1, t>;
); outarray: (1, w);
end;
27
w 0 0
r--!-lr-lri .+ -1. + + 0 0 0 0
8. References [l] N, Ahiezer and M, Krein, Some Questions -in the - Theory of Moments, Translations of MathematixMono.graph, Vol, 2, American Mathematical Society, Providence, 1962, Article VI. [23 P, Concus, D0 Cassatt, G. Jaehnig and E. Melky, "Tables for the evaluation of
s
co xBewx f(x)dx
0
by Gauss-Laguerre Quadrature," Math, Corny., 17 (196~)~ ppO 1~45-256~
[3] P, Davis and I0 Polonsky, "Numerical interpolation,differentiatlon, and integration," Handbook of Mathematical Functions, (edited 'by M0 Abramowitz --. and I. Stegunr National Bureau of Standards, Applied Mathematics Series 55, U, S, Government Printing Office, Washington, D, C,, ppe 875-924, f4] P., Davis and P. Rabinowitz, Numerical Integration, Blaisdell Publishing Company, Waltham, Mass.,-
[5] J0 Francis, "The Q-R transformation. A unitary analogue to the L-R transformation," Comput. J., 4 (1961, 1962), ppe 265-271, 3~2,~ 3450
[6] U0 Hochstrasser, "Orthogonal polynomials, "Handbook of Mathematical Functions, (edited by M. Abramowitz and I. 'm ztional Bureau of Standards, Applied Mathematics Series 55, U. S, Government Printing Office, Washington, D. C., ppO 771-802. H0 Rutishauser, "Solution of eigenvalue problems with the L-R transformation," Further Contributions & the Solution of Simultaneous Linear Equations and the - Determination of Eigenvalues, Et,ional Bureau of Standards, Applied Mathematics Series -9, U. So Government Printi.ng Office, Washington, De C., ppO 47-81, [8] H, Rutishauser, "On a modification of the Q-D algorithm with Graeffe!type convergence," zm' 13 O-963), PP~ 493-496.
[9] A. Stroud and D. Secrest, Gaussian Quadr,ature Formulas, Prentice-Hall, Englewood Cliffs, N. J., lr [lOI H0 Wilf, Mathematics -for the Physical Sciences, John Wiley and Sons, Inc, New York, 1962.
29
