Quadpack: A Subroutine Package for Automatic Integration

310 Pages·1983·7.69 MB·English

by  Robert Piessens

SUMMARY OF THE CALLING SEQUENCES OF THE QUADPACK ROUTINES
CALL QNG(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER)
CALL QAG(F,A,B,EPSABS,EPSREL,KEY,RESULT,ABSERR,NEVAL,IER)
CALL QAGS(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER)
CALL QAGP(F,A,B,NPTS2,POINTS,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER)
CALL QAGI(F,BOUND,INF,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER)
CALL QAWO(F,A,B,OMEGA, INTEGR,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER)
CALL QAWF(F,A,OMEGA,INTEGR,EPSABS,RESULT,ABSERR,NEVAL,IER)
CALL QAWS(F,A,B,ALFA,BETA, INTEGR,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER)
CALL QAWC(F,A,B,C,EPSABS,EPSREL,RESULT,ABSERR,NEVAL,IER)
PARAMETERS
ON ENTRY
F FUNCTION SUBPROGRAM DEFINING THE INTEGRAND FUNCTION OR
A FACTOR OF THE INTEGRAND FUNCTION
A,B LIMITS OF THE INTEGRATION INTERVAL
EPSABS ABSOLUTE ACCURACY REQUESTED
EPSREL RELATIVE ACCURACY REQUESTED
KEY KEY FOR CHOICE OF LOCAL INTEGRATION RULE IN QAG (SEE P. 137)
NPTS2 DETERMINE THE USER PROVIDED BREAK POINTS IN QAGP (SEE P. 156)
POINTS
BOUND DETERMINE THE INTEGRATION INTERVAL IN QAGI (SEE P. 167)
INF
OMEGA PARAMETERS OF THE WEIGHT FUNCTION IN QAWO AND QAWF
INTEGR (SEE P. 176 AND P. 180)
ALFA PARAMETERS OF THE WEIGHT FUNCTION IN QAWS (SEE P. 192)
BETA
INTEGR
C PARAMETER IN THE WEIGHT FUNCTION IN QAWC (SEE P. 203)
ON RETURN
RESULT APPROXIMATION TO THE INTEGRAL
ABSERR ESTIMATE OF THE ABSOLUTE ERROR
NEVAL NUMBER OF INTEGRAND EVALUATIONS
IER PARAMETER PROVIDING INFORMATION ON THE TERMINATION OF THE
ROUTINE
Springer Series in
Computational
Mathematics
1
Editorial Board
R. Graham, Murray Hill
J. Stoer, WOrzburg
R. Varga, Cleveland
R. Piessens E. de Doncker-Kapenga
w.
C. Oberhuber D. K Kahaner
QUADPACK
A Subroutine Package
for Automatic Integration
With 26 Figures
Springer-Verlag
Berlin Heidelberg New York Tokyo 1983
Robert Piessens
Computer Science Department, University of Leuven,
Celestijnenlaan 200 A, 3030 Heverlee, Belgium
Elise de Doncker-Kapenga
Computer Science Department, Western Michigan University,
Kalamazoo, Mi 49008, USA
Christoph W. Oberhuber
Department of Applied and Numerical Mathematics,
Technical University Vienna, Gusshausstr. 27-29,
1040 Wien, Austria
David K. Kahaner
National Bureau of Standards, Washington, DC 20234, USA
AMS Subject Classifications (1980): 65 D 30
ISBN-13:978-3-540-12553-2 e-ISBN-13:978-3-642-61786-7
DOl: 10.1007/978-3-642-61786-7
Library of Congress Cataloging in Publication Data. Main entry under title:
QUAD PACK : a subroutine package for automatic integration. (Springer series in
computational mathematics; 1) Includes bibliographical references. 1. QUADPACK
(Computer programs) 2. Numerical integration-Computer programs. I. Piessens,
R. II. Series. QA299.3.Q36. 1983. 515.4'028'5425. 83-10359
This work is subject to copyright. All rights are reserved, whether the whole or
part of the material is concerned, specifically those of translation, reprinting,
re-use of illustrations, broadcasting, reproduction by photocopying machine or
similar means, and storage in data banks. Under § 54 of the German Copyright
Law where copies are made for other than private use, a fee is payable to
"Verwertungsgesellschaft Wort", Munich.
© by Springer-Verlag Berlin Heidelberg 1983
Reprint of the original edition 1983
214113140-543210
Acknowledgement
We owe the financial support to the "Nationaal Fonds voor Wetenschap
pelijk Onderzoek" and the "Onderzoeksfonds van de Katholieke Universi
teit van Leuven".
We thank Philip Rabinowitz, James Lyness, Ann Haegemans and May
Branders for their comments and suggestions regarding either the
manuscript or the routines and their results.
We thank Brian Ford and the NAG staff, and Bob Huddleston, for their
faith in our product at an early stage, and their assistance with the
implementation of the package as a part of the NAG and of the SLATEC
libraries.
We thank everyone involved in testing the routines, in the first place
Ian Robinson and Tom Patterson.
We thank our analyst-programmers Marc De Meue, Liliane De Roose and
Anita Ceulemans our typists Denise Brams, Lieve Swinnen and Bea
Stroobants, and for drawing the figures Rudy de Doncker.
Table of Contents
I. INTRODUCTION ••......•...•..••............•.........•...•........ 1
II. THEORETICAL BACKGROUND
2.1. Automatic integration with QUADPACK ..............•......... 9
2.2. Integration methods ....................................... 11
2.2.1. Quadrature sum ..................................... 11
2.2.2. Construction of Gauss-type rules with preassigned
abscissae ................•............•............ 15
2.2.3. Modified Clenshaw-Curtis integration •.............. 28
2.2.4. Numerical quadrature and extrapolation ............. 39
III. ALGORITHM DESCRIPTIONS
3.1. QUADPACK contents ....................•......•......•...... 56
3.2. Prototype of algorithm description ....•••.......•......... 59
3.3. Algor i thm schemes •..........••....•....................... 60
3.4. Heuristics used in the algorithms ......•.................. 67
IV. GUIDELINES FOR THE USE OF QUADPACK
4.1. General remarks ........••..............•...•.............. 75
4.2. Decision tree for finite-range integration .•.............. 79
4.3. Decision tree for infinite-range integration .............. 80
4.4. Numerical examples .........................•.............. 81
4.5. Sample programs illustrating the use of the QUADPACK
integrators ...............................•.............. 100
V. SPECIAL APPLICATIONS OF QUADPACK
5.1. Two-dimensional integration .............................. 112
5.2. Hankel transform ......................................... 118
5.3. Numerical inversion of the Laplace transform ............. 121
VI. IMPLEMENTATION NOTES AND ROUTINE LISTINGS
6.1. Implementation notes ........................•...•........ 128
VIII
6.2. Routine listings ••••••••••••••••••••••••••••••••••••••••• 129
QNG •••••••••••••••••••••••••••••••••••••••••••••••• 130
QAG •••••••••••••••••••••••••••••••••••••••••••••••• 137
QAGE ••••••••••••••••••••••••••••••••••••••••••••••• 140
GAGS ••••••••••••••••••••••••••••••••••••••••••••••• 147
QAGP ••••••••••••••••••••••••••••••••••••••••••••••• 156
QAGI ••••••••••••••••••••••••••••••••••••••••••••••• 167
QAWO ••••••••••••••••••••••••••••••••••••••••••••••• 176
QAWF ••••••••••••••••••••••••••••••••••••••••••••••• 180
QAWFE •••••••••••••••••••••••••••••••••••••••••••••• 184
QAWS ••••••••••••••••••••••••••••••••••••••••••••••• 192
QAWSE •••••••••••••••••••••••••••••••••••••••••••••• 195
QAWC ••••••••••••••••••••••••••••••••••••••••••••••• 203
QAWCE •••••••••••••••••••••••••••••••••••••••••••••• 206
QFOUR •••••••••••••••••••••••••••••••••••••••••••••• 213
QK15 ••••••••••••••••••••••••••••••••••••••••••••••• 225
QK21 ••••••••••••••••••••••••••••••••••••••••••••••• 229
QK31 ••••••••••••••••••••••••••••••••••••••••••••••• 233
QK41 ••••••••••••••••••••••••••••••••••••••••••••••• 237
QK51 ••••••••••••••••••••••••••••••••••••••••••••••• 241
QK61 ••••••••••••••••••••••••••••••••••••••••••••••• 245
QK15I •••••••••••••••••••••••••••••••••••••••••••••• 250
QK15W •••••••••••••••••••••••••••••••••••••••••••••• 254
QEXT ••••••••••••••••••••••••••••••••••••••••••••••• 258
QSORT •••••••••••••••••••••••••••••••••••••••••••••• 262
QC250 •••••••••••••••••••••••••••••••••••••••••••••• 265
QC25S •••••••••••••••••••••••••••••••••••••••••••••• 273
QC25C •••••••••••••••••••••••••••••••••••••••••••••• 281
QMOMO •••••••••••••••••••••••••••••••••••••••••••••• 285
QCHEB •••••••••••••••••••••••••••••••••••••••••••••• 288
QWGTO •••••••••••••••••••••••••••••••••••••••••••••• 291
QWGTS •••••••••••••••••••••••••••••••••••••••••••••• 292
QWGTC •••••••••••••••••••••••••••••••••••••••••••••• 293
QMACO •••••••••••••••••••••••••••••••••••••••••••••• 294
REFERENCES •••••••••••••••••••••••••••••••••••••••••••••••••••• 295
I. Introduction
1.1. Overview of Numerical Quadrature
The numerical evaluation of integrals is one of the oldest problems
in mathematics. One can trace its roots back at least to Archimedes.
The task is to compute the value of the definite integral of a given
function. This is the area under a curve in one dimension or a volume
in several dimensions. In addition to being a problem of great practi
cal interest it has also lead to the development of mathematics of much
beauty and insight. Many portions of approximation theory are directly
applicable to integration and results from areas as diverse as orthogo
nal polynomials, Fourier series and number theory have had important
implications for the evaluation of integrals. We denote the problem
addressed here as numerical integration or numerical quadrature. Over
the years analysts and engineers have contributed to a growing body of
theorems, algorithms and lately, programs, for the solution of this
specific problem.
Much effort has been devoted to techniques for the analytic evalua
tion of integrals. However, most routine integrals in practical scien
tific work are incapable of being evaluated in closed form. Even if an
expression can be derived for the value of an integral, often this
reveals itself only after inordinate amounts of error prone algebraic
manipulation. Recently some computer procedures have been developed
which can perform analytic integration when it is possible. Unfor
tunately, these programs are not at all portable and consequently are
only available in a very few installations. Furthermore, successful
analytic integration of rather simple functions can result in lengthy
2
formulae. Such a formula can involve the evaluation of complicated
functions Which can only be computed approximately. Thus many expres
sions for integrals which appear in books or tables as "exact" require
approximations before they can be evaluated numerically. One is lead
then, naturally, to seek methods which attempt to approximate integrals
directly for general integrands rather than by performing problem
specific algebraic manipulation.
An approximate integration sum, also known as an approximate quadra
ture or quadrature evaluation sum, is an expression which is a linear
combination of values of the integrand and which serves as an approxi
mation to the definite integral of this integrand. A Riemann sum is a
simple example. A quadrature or integration rule or formula is any
rule Which yields an integral approximation and may thus include non
linear combinations of integrand or other values. Quadrature rules, or
sums, usually occur in families depending upon a parameter such as the
spacing between the integrand eval~ation points or the number of these
points. Historically, some quadrature rules have been considered which
utilize not only integrand, but derivative evaluations as well. While
these have interesting mathematical properties they do not have much
practical utility for solving quadrature problems. Henceforth we will
be concerned only with rules Which use integrand evaluations, both
linearly and nonlinearly.
Quadrature rules are the building blocks upon which quadrature algo
rithms can be built. We think of an algorithm as a concrete implemen
tation of one or more rules along with a strategy, albeit sometimes
primitive, for evaluating the results. A quadrature rule, by itself,
is not especially useful. Any approximation to an integral involves
errors and unless some assessment of these errors is available no sig
nificance can be attached to the approximation. A Riemann sum without
further specialization, for example, requires far too many points for
practical work. Often during the derivation of the rule a formula will
also be found for the error in the approximation in terms of certain
characteristics of the integrand. The most common of these expresses
the error (of the integral estimate) in terms of a point value of some
derivative of the integrand or the integral of a derivative times an
appropriate kernel function. Alternatively, errors are often given in
terms of a complex path integral, an infinite or asymptotic series, or
a norm in an appropriate function space. In a practical sense, though,
many of these error formulae are not realistically computable. It is
3
sometimes possible to produce an error bound by taking upper bounds of
the various indeterminates that appear in the error expressions.
Unfortunately, such bounds can be wildly pessimistic. The determina-
tion of practical error estimates is as important as the derivation of
new formulae and a subject that lends itself more to engineering and
experiment than to traditional mathematical analysis.
As computer simulations have grown in popularity so have the needs
for numerical methods for the evaluation of integrals. Every computer
center now has some programs which address themselves to these prob
lems. From the point of view of scientists and engineers the diversity
of needs is very great. In principle it is possible to approximate any
integral by a simple average of integrand values. Such general but
primitive methods are too slow by orders of magnitude for any realistic
problem except those with very specific kinds of integrands. Hence
special methods have been developed which make various assumptions
about the integral which is to be approximated. Examples of explicit
assumptions are : the integrand has an algebraic singularity at left or
right end-point, the interval of integration extends to plus infinity,
integrand has a multiplicative factor of sin( wx), etc. An implicit
assumption is that practical integrands are composed of piecewise ana
lytic functions. Programs for "general" integrals usually assume that
Taylor series or other expansions provide good approximations, at least
on small enough intervals. Anyone studying this subject even as a user
or consumer of programs should be aware that the number of potential
references is staggering. The basic theory is presented in the text
books of Krylov (1962), Davis and Rabinowitz (1975), and Engels (1980).
Both new results and programs are continually being pUblished in jour
nals like Mathematics of Computation, Numerische Mathematik, The Com
puter Journal, ACM Transactions on Mathematical Software, SIAM Journal
on Numerical Analysis, Computing, BIT and many others.
Many of the programs published in the course of the past 20 years
are automatic. That is, they are meant to be called with the problem
specifications such as the integrand, the interval of integration, the
level of accuracy desired, etc., given by the user. Upon return, they
provide an estimate of the integral and an error estimate which hope
fully satisfy the accuracy requirement determined by the call sequence.
Papers which are concerned with automatic integration are listed in the
bibliography of de Doncker and Piessens (1976). Important surveys and
comparative studies were published by Kahaner (1971), Dixon (1974) and
Robinson (1979).