Some new methods for performing numerical quadrature of an integrable function over a finite interval are described. Each method provides a sequence of approximations of increasing order to the value of the integral. Each approximation makes use of all previously computed values of the integrand. The points at which new values of the integrand are computed are selected in such a way that the order of the approximation is maximized. The methods are compared with the quadrature methods of Clenshaw and Curtis, Gauss, Patterson, and Romberg using several examples

Feagin, T.

English

NASA Technical Reports Server

- X-582-73-176
"-a_
- PREPRINT
THE-DEVELOPMENT OF ACCURATE
AND EFFICIENT METHODS
OF NUMERICAL QUADRATURE-
(NASA-TM-X-662 77) THE DEVELOPMIENT OF
ACCURATE AND EFFICIENT METHODS OF
NUMERICAL QUADRATURE (NASA) 1T p HC
CSCL 12A
N73-27518
Unclas
G3/19 09187
Reproduced by
NATIONAL TECHNICAL
INFORMATION SERVICE
US Department of Commerce
Springfield, VA. 22151
JUNE 1973
-' -  -
- GODDARD SPACE FLIGHT CENTER
.--i . GREENBELT,-MARYLAND - -I-
.~ ~ ~ I
.
(
I
bP
;
https://ntrs.nasa.gov/search.jsp?R=19730018791 2020-03-23T00:09:40+00:00Z
THE DEVELOPMENT OF ACCURATE AND EFFICIENT
METHODS OF NUMERICAL QUADRATURE*
By
T. Feagin
June 1973
*This research was conducted while the author held a NRC Resident Research Associateship at
the Goddard Space Flight Center, Greenbelt, Maryland.
I
PAEE LBANK NOT FILMS l
CONTENTS
ABSTRACT . .........................
INTRODUCTION .......................
DESCRIPTION AND DERIVATION OF THE METHODS
RESULTS, COMPARISONS AND CONCLUSIONS.....
ACKNOWLEDGEMENTS ................... ..
REFERENCES .... .. .. ..0 .. . 0........
.0....0....
TABLES
Table
I
II
III
Method I ..... ....0 0 ......................
Method II . 0 .......................
Errors Obtained Using Certain Quadrature Methods ..........
iii
Page
v
1
5
10
13
13
Page
6
7
12
. .. .. 0.
. .. .. .. ..
PIiRCEDING PAGE BLANK NOT FILMED
THE DEVELOPMENT OF ACCURATE AND EFFICIENT
METHODS OF NUMERICAL QUADRATURE*
By
T. Feagin
ABSTRACT
Some new methods for performing the numerical quadrature of an
integrable function over a finite interval are described. Each method provides
a sequence of approximations of increasing order to the value of the integral.
Each approximation makes use of all previously computed values of the inte-
grand. The points at which new values of the integrand are computed are se-
lected in such a way that the order of the approximation is maximized. The
methods are compared with the quadrature methods of Clenshaw and Curtis,
Gauss, Patterson, and Romberg using several examples.
v
THE DEVELOPMENT OF ACCURATE AND EFFICIENT
METHODS OF NUMERICAL QUADRATURE*
INTRODUCTION
One of the most basic problems in numerical mathematics concerns the
numerical computation of the definite integral of an integrable function defined
over some finite interval. Most numerical methods that have been devised for
performing such quadratures approximate the value of the integral by a linear
combination of function values, that is,
b n
I= f(x) dx= wi f(x i ) + En,
where f(x) represents the function defined on the interval (a, b); wi denotes the
weight of the value of the function at the ith point or abscissa, xi; and En
designates the error committed when the approximation is made.
There have generally been two basic approaches to the problem of numeri-
cal quadratures, adaptive and non-adaptive. If the abscissae are selected at
some intermediate stage of the process depending upon the behavior of the
integrand, then the pl cu.ss is said to be adaptive. If the process evolves inde-
pendently of the behavior of the integrand, the process is said to be non-adaptive.
*This research was conducted while the author held a NRC Resident Research Associateship at
Goddard Space Flight Cr.iv,', Greenbelt, Maryland.
1
The methods developed here may be used readily in either approach, because an
estimate of the error is available at each stage. The error estimate could be
used to make a decision regarding the subdivision of the interval into smaller
parts and the formulae could then be used on the smaller intervals.
If the behavior of the integrand is sufficiently well understood at the outset,
it may be feasible to determine a value of n for which the error En is accep-
table. In such a context, the Gaussian. formulae generally provide the most ac-
curate and economical results. The Gaussian weights and abscissae are
tabulated, for example, by Stroud and Secrest (1966).
On the other hand, if little is known about the behavior of the integrand, it
may be necessary to compute many approximate values of the integral before
an acceptable result is found. Furthermore, it is usually necessary to determine
whether or not a given result is acceptable. It is therefore recommended that
any practical procedure for performing quadratures be able to provide a sequence
of results of increasing accuracy, providing simultaneously reliable error esti-
mates which may be used to decide when to terminate the procedure. The
methods presented here have the above characteristics and also appear to be
very efficient with regard to the number of integrand values required in order
to achieve a given accuracy.
One method of numerical quadratures that has the ability to provide a
sequence of approximations of increasing order would be to use a sequence of
Gaussian formulae and to estimate the error as the difference between the
2
is evaluated include all old points, thereby causing the algorithm to be quite
economical. It should be remarked that the additional effort required in order
to obtain the coefficients of the series is only necessary if the error estimate
or an indefinite integral is desired.
It is important to realize that no estimate of the error En is completely
foolproof. There is always some function which will cause any estimate to be a
gross underestimate or a gross overestimate of the true error. Usually the
effort expended in order to improve the estimate of E is equivalent to improv-
n
ing the approximate value of the integral itself. The most essential feature of
a good quadrature method is therefore the ability to obtain a sequence of accu-
rate results with as few function evaluations as is possible.
Patterson (1968) has derived a set of quadrature formulae which have
several desirable features. The method is comprised of a sequence of quad-
rature formulae (beginning with the midpoint rule, the Gauss three-point rule,
and Kronrod's seven-point rule) each of which contains all previous abscissae
and the new abscissae and the weights of each formulae are selected so as to
optimize the degree of the formula. The new set of n-point formulae are of
degree (3n + 1)/2. These formulae may be differenced to provide an estimate of
the error, or perhaps a series of Chebyshev polynomials could be fitted to the
values of the integrand and integrated term-by-term and the last few terms of
the series could be used as an estimate of the error.
4
results obtained for the different orders. The drawback with such a method is
that the Gaussian formulae have few abscissae in common for different values
of n, thereby causing the total number of function evaluations required to obtain
an acceptable result to grow rapidly with the order of the approximation.
Another method which might be used to generate a sequence of results of
increasing order is based upon the use of Newton-Cotes formulae. In this case,
the formulae do have many abscissae in common for different values of n. The
problem here is that the higher-order Newton-Cotes have undesirable features
with regard to round-off errors and convergence properties (Hildebrand, 1956).
Romberg (1955) has developed a method which is also based upon equally-spaced
points that avoids these difficulties. However, the Romberg method achieves
high-order results very slowly as the number of function evaluations increases.
Bauer, Rutishauser, and Stiefel (1963) have succeeded in reducing somewhat the
number of function values required by modifying the Romberg algorithm, but the
modification causes the method to be more susceptible to round-off errors.
Clenshaw and Curtis (1960) recommend a method of numerical quadrature
based upon the development of the integrand in a Chebyshev series of polynomials.
The quadrature is accomplished by integrating the series term-by-term. The
error is estimated by examining the last few terms of the integrated series. If
the coefficients of the series are converging rapidly enough, the last few terms
generally provide a conservative estimate of the error. This method is
particularly advantageous if an indefinite integral is desired. If the number of
terms in the series is increased from N to 2N-1, the points at which the function
3
DESCRIPTION AND DERIVATION OF THE METHODS
The methods presented here are founded upon the same basic philosophy as
that adopted by Patterson (1968). Such methods are characterized by the initial
formula or rule to which all subsequent abscissae are added and by the way in
which these points are added. In the same way that Patterson's method has the
midpoint rule as its foundation, the two methods developed here are based upon
the trapezoid rule and the Gauss two-point rule respectively.
The first method is based upon the trapezoid rule and is comprised of a
sequence of quadrature formulae with 2 + 1 abscissae for k = 0, 1, 2, . ..
The weights and abscissae of the first five members of the sequence are given
in Table I, where the interval (a, b) has been chosen to be (-1, 1) for simplifica-
tion. The error term of each formulae of the sequence is also tabulated.
The second method is based upon the Gauss two-point rule and is com-
prised of a sequence of quadrature formulae with 3 ' 2 k -1 abscissae for
k = 0, 1, 2, .... The weights, abscissae, and the error term for each of the
first four members of this sequence are given in Table II.
The formulae may be derived in the following way. For each formulae of
a given method, the n abscissae may be characterized as the roots of an nth
degree polynomial, Gn (x), where
n
Gn(x) = b x j .
j =0
5
Table I
Method I
Two-point Formula of Degree One (Trapezoid Rule) E. = -6.67 x 10-1f(2)(
Xi Wi
1.0000 0000 0000 1.0000 0000 0000
Three-point Formula of Degree Three (Simpson's Rule) En = -1.11 x 10-2f(4)(
X.i  W
0.0000 0000 0000 1.3333 3333 3333
1.0000 0000 0000 0.3333 3333 3333
Five-point Formula of Degree 7 (Lobatto5-ptRule) En = -3.60 x 10-7f(8)(8)
Xi  Wi
0.0000 0000 0000 0.7111 1111 1111
0.6546 5367 0708 0.5444 4444 4444
1.0000 0000 0000 0.1000 0000 0000
Nine-point Formula (Degree 13) E, = -6.16 x 10-16 f (14)(6)
- X i  W i
0.0000 0000 0000 0.3437 6208 7210
0.3409 8226 5911 0.3342 3373 9816
0.6546 5367 0708 0.2839 7877 8048
0.8904 0552 7513 0.1792 6269 9553
1.0000 0000 0000 0.0306 4373 8977
Seventeen-point Formula (Degree 25) En = -2.35 x 10- 35 f( 26 )(6)
X XWi
0.0000 0000 0000 0.1588 8791 8790
0.1643 3743 9165 0.1726 1252 6200
0.3409 8226 5911 0.1761 6511 9745
0.5089 0120 1269 0.1572 3704 4683
0.6546 5367 0708 0.1356 0287 5607
0.7824 1604 1120 0.1196 2516 8318
0.8904 0552 7513 0.0939 8151 8229
0.9661 9986 3094 0.0560 0038 0783
1.0000 0000 0000 0.0093 3140 7070
6
Table II
Method II
Two-point Formula of Degree Three (2-point Gauss) En = 7.41 x 10- 3 f(4 )(f)
' X i  W.
0.5773 5026 9190 1.0000 0000 0000
Five-point Formula of Degree Seven (Kronrod) En = -9.00 x 10-8f(8)(4)
Xi  Wi
0.0000 0000 0000 0.6222 2222 2222
0.5773 5026 9190 0.4909 0909 0909
0.9258 2009 9773 0.1979 7979 7980
Eleven-point Formula of Degree 17 En = 5.16 x 10-22 f(l 8 )(6)
Xi W
0.0000 0000 0000 0.3102 6774 3075
0.3047 6216 9153 0.2938 2188 8625
0.5773 5026 9190 0.2466 4959 0846
0.7900 8322 3401 0.1757 8249 0943
0.9258 2009 9773 0.0961 C769 3081
0.9878 5368 2853 0.0325 0446 4966
Twenty-three-point Formula of Degree 35 E = -1.03 x 10-5 3 f(36) ()
fXi Wi
0.0000 0000 0000 0.1551 3596 2892
0.1544 4378 2428 0.1530 6173 0468
0.3047 6216 9153 0.1469 0821 8576
0.4469 6796 4050 0.1368 8253 6133
0.5773 5026 9190 0.1233 3089 7738
0.6926 1884 7335 0.1067 5329 7040
0.7900 8322 3401 0.0878 6743 0877
0.8679 4828 6028 0.0677 8753 4352
0.9258 2009 9773 0.0481 9474 7432
0.9650 2643 1039 0.0306 0569 3330
0.9878 5368 2853 0.0155 9062 3275
0.9978 9629 3295 0.0054 4930 9333
7
Similarly, let Qp (x) be a polynomial of degree p with zeroes at the p abscissae
that are to be added, where
p
Q(x):21 axi.
j=o
Then, the next formula in the sequence has n + p abscissae located at the
zeroes of G+ p (x) where
n+p
Gn+ p (x) = cj xi j ,
j=o
and where
r
c. = L a. b.
j=q
i = 0, 1, ... , n + p
where q = max (0, i-n) and r = min (p, i). In order to eliminate any ambiguity,
let us assume that ao, b0, and co are equal to one. In order to optimize the order
or degree of the new formula, it is required that the conditions
Gn+p ( x) x k dx = 0 k = O, 1, ... , p-l
be satisfied (Patterson, 1968). These conditions comprise a system of linear
equations in the p unknowns { aj )P . For the two methods presented here,j--1 I
8
1i
-1I
n + p is always an odd number and for such cases the system of linear equations
is
n+p
0 = cj/(j +k+l) k = 0 , p-l.
j=o
j+k even
So long as the coefficients, bj, are rational numbers, then so are the aj and the
c . Therefore, if the abscissae of the initial formula are the zeroes of a poly-
nomial with rational coefficients, all subsequently added abscissae will also be
the zeroes of a polynomial with rational coefficients. This fact permits the use
of rational arithmetic in the derivation of the methods presented here in order to
avoid the propagation of numerical errors from one stage of the process to the
next. The only errors committed at a given stage are those arising from the
determination of the new abscissae as the zeroes of these polynomials and from
the computation of the corresponding weights of the quadrature formula. The
secant method is used to obtain the zeroes of the polynomials in such a way that
each of the zeroes is determined to lie within an interval less than 10 -14 in
length. Gaussian elimination is then used in order to determine the weights from
the linear equations,
n
, w i xi = (1 -[-])/(j + 1) j = O, 1, ... , n - 1.
i=1
9
RESULTS, COMPARISONS, AND CONCLUSIONS
The two methods presented in Tables I and II have been applied to a few
simple quadrature problems for which the true value of the integral is known.
The results are compared with those obtained using the following well-known
methods: Romberg's method, the Clenshaw-Curtis method, Patterson's method,
and a method composed of Gauss formulae. The composite Gauss method con-
sists of Gauss formulae of degrees one, three, seven, and fifteen. The total
number of function evaluations required at each stage are therefore one, three,
nine, and twenty-three respectively.
The following set of problems are considered:
T 8 dx
= 7-, (1)
1 x2 + 2x+ 5
I x1/ 2 dx = 2/3, (2)
1 x3 /2 dx = 2/5, (3)
dx n 2, (4)
x1
and
l dx - 7/2. (5)
-1 + x 2
10
The results of applying the above methods to this set of problems are
summarized in Table III.
The first three methods shown in the table use the function values at the
endpoints of the interval of integration and the last three methods do not. For
this reason, let us examine these first three methods as a group. In problems
1, 4, and 5, the method which results in the smallest error is Method I, followed
by the Clenshaw-Curtis method, and Romberg's method. In examples 2 and 3,
which possess singular derivatives at the origin (first and second derivative,
respectively), the Clenshaw-Curtis method produces the best approximation,
followed by Method I and the Romberg method in that order.
For the last three methods in the table, the comparison is more difficult
because the different methods produce results for different values of n, the total
number of function evaluations. In any case, if we consider the number of sig-
nificant figures gained per function evaluation as a measure of the efficiency
of a method, some conclusions can be drawn. In each of the five problems, the
Patterson method and Method II provide results that are almost identical using
the above measure of efficiency. The results using the composite Gauss method
are somewhat inferior.
Comparison of all six methods tends to indicate that the last three methods
are more efficient, especially on problems 2 and 3. In all cases, the composite
Gauss method performed as well or slightly better than the Clenshaw-Curtis
method. For each problem, the best methods are Patterson's method and
11
Table III
Errors Obtained Using Certain Quadrature Methods
Problem Romberg's Clenshaw-Curtis Method I from Patterson's Composite Gauss Method II from
n o . n nNo. n Method n Method n Table I n Method n Method (1-3-7-15) Table II
3 8.26 x 10- 3  3 8.26 x 10- 3  3 8.26 x 10- 3  1 5.84 x 10- 2 1 5.84 x 10 - 2 2 5.94 x 10 - 3
5 5.25 x 10-4 5 2.36 x 10- 4  5 2.40 x 10-5 3 5.25 x 10 - 4  3 5.25 x 10 - 4  5 6.04 x 10 - 6
9 6.87 x 10 -6 9 3.31 x 10 - 9  9 6.91 x 10-1 °  7 1.84 x 10-8 9 2.66 x 10 - 9  11 6.32 x 10 - 13
11 8 10 7
17 1.17x10 17 - 17 - 15 - 23 - 23
3 2.86 x 10-2 3 2.86 x 10- 2  3 2.86 x 10 - 2  1 4.04 x 10 - 2 1 4.04 x 10 - 2  2 7.22 x 10- 3
5 8.91 x 10 - 3  5 2.07 x 10- 3  5 4.53 x 10-3 3 2.51 x 10- 3  3 2.51 x 10 - 3  5 3.41 x 10 - 4
9 3.06 x 10-3 9 2.24 x 10- 4  9 7.61 x 10-4 7 1.42 x 10 - 4  9 2.46 x 10- 4  11 2.48 x 10-5
17 1.07 x 10-3 17 2.69 x 10 - 5 17 1.27 x 10- 4  15 7.05 x 10-6 23 2.77 x 10 - 5  23 2.13 x 10 - 6
3 2.37 x 10- 3  3 2.37 x 10- 3  3 2.37 x 10 - 3  1 4.64 x 10 - 2  1 4.64 x 10- 2  2 1.22 x 10- 3
5 3.03 x 10 -4  5 1.25 x 10-5 5 9.21 x 10-5 3 1.88 x 10 - 4  3 1.88 x 10 - 4  5 2.76 x 10 -6
9 4.96 x 1 0 -s 9 7.79 x 10-7 9 4.53 x 10-6 7 1.03 x 10-6 9 3.60 x 10-6 11 7.71 x 10-8
17 8.62 x 10 - 6 17 2.73 x 10-8 17 2.27 x 10 - 7  15 1.78 x 10 - 9  23 9.29 x 10 - 8  23 1.63 x 10 - 9
3 1.30 x 10- 3  3 1.30 x 10-3 3 1.30 x 10 - 3  1 2.65 x 10 - 2 1 2.65 < 10-2 2 8.39 x 10- 4
5 2.74 x 10 - 5  5 9.93 x 10 - 5  5 9.68 x 10 - 7  3 2.55 x 10- 5 3 2.55 x 10-5 5 2.18 x 10 - 7
4
9 2.97 x 10 - 7  9 6.40 x 10 - l °0  9 6.44 x 10- 12  7 2.39 x 10 - °  9 1.99 x 10 l  11
17 1.36 x 10-9 17 - 17 - 15 - 23 - 23 -
3 9.59 x 10- 2  3 9.59 x 10 - 2  3 9.59 x 10 - 2  1 4.29 x 10-' 1 4.29 x1 10-2
5 1.08 x 10- 2  5 6.98 x 10-3 5 2.54 x 10- 3  3 1.25 x 10-2 3 1.25 x 10- 2  5 9.99 x 10-45
9 4.38 x 10- 4  9 9.77 x 10-6 9 1.06 x 10-6 7 3.35 x 10-5 9 1.11 x 10-5 11 1.37 x 10-7
17 5.17 x 10-6 17 5.30 x 10- °' 0  17 4.88 x 10-11 15 1.87 x 10 - °10 23 8.46 x 10 - 12 23
-~~~~~~~~~~~ ~ ~ ~   2 -ni-e nerris hn1x1.3 n oa ubro ucineautos
tŽ4
- Indicates an error less than 1 x 10-13 n - Total number of function evaluations.
Method II, each giving about the same over-all efficiency. Following these for
examples 2 and 3 are the composite Gauss method, the Clenshaw-Curtis method,
Method I, and Romberg's method, in that order. For examples 1, 4, and 5,
Method I is the third most efficient method, followed by the composite Gauss,
the Clenshaw-Curtis method, and Romberg's method, in that order.
These preliminary results seem to indicate that the methods of the Patter-
son type are certainly worthy of further development and study.
ACKNOWLEDGEMENTS
This research was accomplished while the author held a National Research
Council Resident Research Associateship at NASA's Goddard Space Flight
Center, Greenbelt, Maryland.
REFERENCES
Bauer, F. L., Rutishauser, M., and Stiefel, E., 1963, "New Aspects in Numerical
Quadrature," Proc. Symp. Applied Math., Vol. XV, Experimental Arithmetic,
High Speed Computing, and Mathematics, Amer. Math. Soc.
Clenshaw, C. W., and Curtis, A. R., 1960, "A method for numerical integration
on an automatic computer," Numer. Math. 2, 197-205.
Fox, L., and Parker, I. B., 1968, Chebyshev Polynomials in Numerical Analysis,
Oxford Univ. Press, London, pp. 34-36.
Hildebrand, F. B., 1956, Introduction to Numerical Analysis, McGraw-Hill,
New York.
13
Patterson, T. N. L., 1968, "The Optimum Addition of Points to Quadrature
Formulae," Math. Comp. 22, pp. 847-856.
Romberg, W., 1955, "Verein fachte numerische Integration," Det. Kong. Norske
Videnskaber Selskab Forhandlinger, Band 28, Nr. 7, Trondheim.
Stroud, A. M., and Secrest, D., 1966, Gaussian Quadrature Formulas, Prentice-
Hall, Inc., Englewood Cliffs, N.J., pp. 100-151.
NASA-GSFC
14


The development of accurate and efficient methods of numerical quadrature

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server