Atkinson K. An Introduction to Numerical Analysis
Подождите немного. Документ загружается.
AUTOMATIC
NUMERICAL
INTEGRATION
299
12.
If
iRk-
Rk-d
> t:, then go to step 3
13.
Since
IRk-
Rk_Ji
.:::;;
t:, accept
int
=
Rk-l
and n;turn.
There are many variants of Romberg integration.
For
example, other ways
of
increasing the number of nodes have been studied.
For
a very complete survey
of
the literature on Romberg integration, see Davis
and
Rabinowitz (1984, pp.
434-446). They also give a Fortran program for Romberg integration.
5.5 Automatic Numerical Integration
An automatic numerical integration program calculates
an
approximate integral
to within an accuracy specified by the user of the program. The user does not
need to specify either the method or the number
of
nodes to be used. There are
some excellent automatic integration programs,
and
many people use them. Such
a program saves you the time of writing your own program, and for many people,
it avoids having to understand the needed numerical integration theory. Nonethe-
less, it is almost always possible to improve
upon
an automatic program,
although
it
usually
~equires
a good knowledge
of
the numerical integration
needed for your particular problem. When doing only a small number
of
numerical integrations, automatic integration
is
often
a good way to save time.
But for problems involving many integrations, it is probably better to invest the
time to find a less expensive numerical integration procedure.
An automatic numerical integration program functions as a
"black
box,"
without the user being able to see the intermediate steps
of
the computation.
Because
of
this, the most important characteristic
of
such a program is that it be
reliable: The approximate integral that is returned
by
the program and that the
program says satisfies the user's error tolerance must, in fact, be that accurate. In
theory, no such algorithm exists,
as
we
explain in
the
next paragraph. But for the
type of integrands that one usually considers in practice, there are programs that
have a high order of reliability. This reliability will
be
improved if the user reads
the program description, to see the restrictions
and
assumptions of the program.
To understand the theoretical impossibility
of
a:
perfectly reliable automatic
integration program, note that the program will evaluate the integrand
j(x)
at
only a finite number
of
points, say x
1
,
••.
,
xn-
Then
there are an infinity of
continuous functions
/(x)
for which
/-(x;) = f(x;)
i = 1,
...
, n
and
In
fact, there are an infinity of such functions
/(x)
that
are infinitely differentia-
300 NUMERICAL INTEGRATION
ble.
For
practical problems, it
is
unlikely that a well-constructed automatic
integration program
will
be unreliable, but it
is
possible. An automatic integra-
tion program can be made more reliable by increasing the stringency of its error
tests,
but
this also makes the program
less
efficient. Generally there
is
a tradeoff
between reliability and efficiency. For a further discussion of the questions of
reliability and efficiency of automatic quadrature programs, see Lyness and
Kaganove (1976).
Adaptive
quadrature
Automatic programs can
be
divided into
(1)
those using a
global rule, such
as
Gaussian quadrature or the trapezoidal rule with even
spacing, and (2) those using an
adaptive strategy, in which the integration rule
varies its placement of node points and even its definition to reflect the varying
local behavior of the integrand. Global strategies use the type of error estimation
that we have discussed in previous sections. We now discuss the concept and
practice of
an
adaptive strategy.
Many integrands vary in their smoothness
or
differentiability at different
points of the interval of integration [a,
b].
For example, with
I=
frx dx
0
the integrand has infinite slope at x =
0,
but the function is
well
behaved at
points
x near 1. Most numerical methods use a uniform grid of node points, that
is, the density of node points
is
about equal throughout the integration interval.
This includes composite Newton-Cotes formulas,. Gaussian quadrature, and
Romberg integration. When the integrand
is
badly behaved at some point a in
the interval
[a, b
],
many node points must be placed near a to compensate for
this. But this forces many more node points than necessary to be used at all other
parts of
[a, b]. Adaptive integration attempts to place node points according to
the behavior of the integrand, with the density
of
node points being greater near
points of
bad
behavior.
We now explain the basic concept of adaptive integration using a simplified
adaptive Simpson's rule. To see more precisely why variable spacing
is
necessary,
consider Simpson's rule with such a spacing of the nodes:
nj2
nj2
( )
x2j
•
x2j-
x2j-2
I(I)
=
~
J
f(x)
dx
=In(!)=
~
6
(/2j-2
+
4/2j-l
+
/2)
;=1
Xzj-2
;=1
with x
2
j_
1
=
(x
2
j_
2
+ x
2
)/2.
Using (5.1.15),
1
n/2
I(/)
-In(!)
= -
2880
~
(xzj-
Xzj-2)
5
J<
4
)(
~)
;~1
(5.5.1)
(5.5.2)
with x
2
j_
2
~
~j
~
x
2
j.
Clearly, you want to choose x
2
j-
x
2
j_
2
according to the
size
of
J<
4
)(~),
which is unknown in general.
If
j<
4
)(x)
varies greatly in magni-
tude, you
do
not
want even spacing of the node points.
---·--··--
---
·-----~-----
-··-·--·'-·'~
AUTOMATIC NUMERICAL INTEGRATION
301
As notation, introduce
h[
(a+{))
·]
~~~~
= 3
f(a)
+
4f
-
2
-
+
f(P)
/(2)
=
J(l)
+
J(l)
a,p
a,y
y,p
{)+a
y=--
2
To describe the adaptive algorithm for computing
{1-a
h=--
2
(5.5.3)
we
use a recursive definition. Suppose that
t:
> 0 is given, and that
we
want to
find an approximate integral
f for which
(5
.5
.4)
Begin by setting
a=
a,
{1
=b.
Compute~~~~
and
I!~~-
If
,
1
~2>
_
J<I>
I
:$
t:
a,p
a,P
(5.5.5)
then accept
~~~1
as
the adaptive integral approXimation to Ia,p· Otherwise let
t:
==
t:/2, and set the adaptive integral for
Ia,p
equal to the sum of the adaptive
integrals for
/a,y
and
I.r.P•
y
=(a+
{1)j2, each to be computed with an error
tolerance of
t:.
In an actual implementation as a computer program, many extra limitations
are included as safeguards; and the error estimation
is
usually much more
sophisticated. All function evaluations are handled carefully in order to ensure
that the integrand
is
never evaluated twice at the same point. This requires a
clever stacking procedure for those values of
f(x)
that must be temporarily
stored because they will be needed again later in the computation. There are
many small modifications that can be made to improve the performance of
the
program, but generally a great deal of experience and empirical investigation is
first necessary. For that and other reasons, it is recommended that standard
well-tested adaptive procedures be used [e.g., de Boor (1971), Piessens et al.
(1983)]. This
is
discussed further at the end of the section.
Table
5.18
Adaptive Simpson's example (5.5.6)
[a,,B]
[<2)
I-
J<2l
I-
JCI>
11(2)-
[(ill
t:
[0.0,
.0625]
.010258
1.6E-
4
4.5E-
4
2.9E-
4
.0003125
[.0625,
.0125]
.019046
1.2E-
7
l.lE-6
l.OE-
6
.0003125
[.125,
.25]
.053871
4.5E-
7
3.6E-
6
4.0E-
6
.000625
[.25,.5] .152368
9.3E-
7
l.lE-5
l.OE-
5
.00125
[.5,
1.0]
.430962
2.4E-
6
3.0E-
5
2.8E-
5
.0025
i
I
I
I
i
·-·-
----·---
--------··-
----
-------
J
i
I
!
;
302
NUMERICAL
INTEGRATION
Example Consider using the preceding simpleminded adaptive Simpson proce-
dure to evaluate
I=
[IX
dx
0
(5.5.6)
with t: = .005 on
[0,
1).
The final intervals [a,
fJJ
and integrals
I~~~
are given in
Table
5.18. The column labeled t: gives the error tolerance used in the test (5.5.5),
which estimates the error
in
I~~1-
The error estimated for
I~~~
on
[0,
.0625]
was
inaccurate,
but
it
was
accurate for the remaining subintervals. The value used to
estimate
Ia.fJ
is actually
I~~1,
and it
is
sufficiently accurate on all subintervals.
The total integral, obtained by summing
all
I~~1.
is
i = .666505
and the calculated bound
is
I-f=
1.6£-4
I I -
il
~
3.3£ - 4
obtained
by
summing the column labeled
jJC
2
> -
JC
1
>j.
Note that the error is
concentrated on the first subinterval,
as
could have been predicted from the
behavior of the integrand near x = 0. For an example where the test (5.5.5)
is
not
adequate, see
Problem
32.
Some automatic integration programs One
of
the better known automatic
integration programs
is
the adaptive program CADRE (Cautious Adaptive
Romberg Extrapolation), given in de Boor
(1971).
It
includes a means of
recognizing algebraic singularities at the endpoints of the integration interval.
The asymptotic error formulas
of
Lyness and Ninham (1967), given in (5.4.21) in
a special case, are used to produce a more rapidly convergent integration method,
again based
on
repeated Richardson extrapolation. The routine CADRE has
been found empirically to be both quite reliable and efficient.
Table 5.19
Integration examples for CADRE
Desired Error
10-2 10-5
10-s
Integral
Error
N Error N Error
N
Jl
A=
2.49E-
6
76
A =
1.40E-
7
73
A=
4.60E-
11
225
P =
5.30E-
4 P =
4.45E-
6
P =
2.48E-
9
/2
A=
1.18E-
5 9
A=
3.96E-
7
17
P = 2.73E
-10
129
P =
3.27E-
3
P =
3.56E-
6
P =
2.81E-
9
/3
A=
1.03E-
4 17
A=
3.23E-
8
33
A=
1.98E-
9
65
P =
2.98E-
3
P =
4.43E-
8 P =
2.86E-
9
/4
A=
6.57E-
5
105
A=
6.45E-
8
209
A=
4.80E-
9
281
P =
4.98E-
3
P =
9.22E-
6
P =
1.55E-
8
Is
A=
2.77E-
5
226
A=
7.41E-
8
418
A=
5.89E-
9
562
P =
3.02E-
3
P =
l.OOE-
5
P =
1.11E-
8
/6
A=
8.49E-
6
955
A=
2.37E-
8
1171
A=
4.30E-
11
2577
P =
8.48E-
3
P =
1.67E-
5
P =
2.07E-
8
/7
A=
4.54E-
4
98
A=
7.72E-
7
418
A=**
1506
P =
1.30E-
3
P =
8.02E-
6
p = **
AUTOMATIC NUMERICAL INTEGRATION
303
A more recently developed package
is
QUAD PACK, some of whose programs
are general purpose, while others deal with special classes of integrals. The
package
was
a collaborative effort, and a complete description of it
is
given in
Piessens et al. (1983). The package
is
well
tested and appears to be an excellent
collection of programs.
We illustrate the preceding by calculating numerical approximations to the
following integrals:
1
1 4 dx 1 .
I1
=
2
=
-[tan-
1
(10) +
tan-
1
(6)]
o 1 + 256(x - .375) 4
1
1 2
I
2
=
xVx
dx
=-
0 7
1
1 2
I
3
=
IX
dx
=-
0 3
I
4
= f log (
x)
dx =
~
1
Is=
flog
lx-
.71
dx
= .3log (.3) + .?log (.7) - 1
!
10000 dx
I1
=
--
=
206
-9
v1xT
From QUADPACK,
we
chose DQAGP.
It
too contains ways to recognize'
algebraic singularities at the endpoints and to compensate for their presence. To
improve performance, it allows the user to specify points interior to the integra-
tion interval
at
which the integrand is singular.
We used both CADRE and DQAGP to calculate the preceding integrals, with
error tolerances of
10-
2
,
10-
5
,
and
10-
8
•
The results are shown in Tables 5.19
and 5.20. To more fairly compare DQAGP and CADRE,
we
applied CADRE to
two integrals in both
Is
and I
7
,
to have the singularities occur at endpoints.
For·
example,
we
used CADRE for each of the integrals
in
J
o dx
1
10000 dx
I-
--+
-
7
-
-91-x
o
IX
(5.5.7)
In the tables, P denotes the error bound predicted
by
the program and A denotes
the actual absolute error in the calculated answer. Column
N gives the number of
integrand evaluations. At all points at which
the integrand was undefined, it was
arbitrarily set to zero. The examples were computed in double precision on a
Prime
850,
with a unit round of r
46
= 1.4 x
10-
14
•
In Table 5.19, CADRE failed for I
7
with the tolerance
10-
8
,
even though
(5.5.7)
was
used. Otherwise, it performed quite well. When the decomposition
I
i
I
I
I
J
304
NUMERICAL INTEGRATION
Table 5.20 Integration examples for DQAGP
Desired Error
10
2
10
5
w-s
Integral Error N Error
N
Error
N
II
A=
2.88E-
9
105
A=
5.40E-
13
147
A=
5.40E-
13
147
P =
2.96E-
3 P =
5.21E-
10
P = 5.21E
-10
I2
A =
1.17E-
11
21
A=
1.17E-
11
21
A=
1.17E-
11
21
P =
7.46E-
9
P =
7.46E-
9
P =
7.46E-
9
IJ
A=
4.79E-
6
21
A=
4.62E
-13
189
A=
4.62E-
13
189
P =
4.95E-
3
P =
4.77E-
14
P =
4.77E-
14
I4
A=
5.97E-
13
231
A=
5.97E
-13
231
A=
5.97E-
13
231
P = 7.15E
-14
P =
7.15E-
14
P =
7.15E-
14
Is
A=
8.67E-
13
462
A=
8.67E-
13
462
A=
8.67E-
13
462
P
=USE-
13
P =
1.15E-
13
P =
1.15E-
13
I6
A=
l.OOE-
3
525
A=
6.33E-
14
861
A =
5.33E-
14
1239
P =
4.36E-
3 P =
8.13E-
6 P =
7.12E-
9
I1
A=
1.67E-
10
462
A=
1.67E-
10
462
A =
1.67E-
10
462
P = 1.16E
-10
P =
1.16E-
10 P =
1.16E-
10
(5.5.7) is
not
used and CADRE
is
called only once for the single interval
[-
9,
10000], it fails for
all
three error tolerances.
In Table
5.20, the predicted error
is
in some cases smaller than the actual
error. This difficulty appears to be due to working at the limits of the machine
arithmetic precision, and in all cases the final error was well within the limits
requested.
In comparing the two programs,
DQAGP and
CADRE
are both quite reliable
and efficient. Also, both programs perform relatively poorly for the highly
oscillatory integral /
6
,
showing
that
/
6
should be evaluated using a program
designed for oscillatory integrals (such as
DQAWO in QUADPACK, for Fourier
coefficient calculations). From the tables,
DQAGP
is
somewhat more able to deal
with difficult integrals, while remaining about equally efficient compared to
CADRE. Much more detailed examples for
CADRE
are given in Robinson
(1979).
Automatic quadrature programs can be easily misused in large calculations,
resulting
in
erroneous results and great inefficiency.
For
comments on the use of
such programs in large calculations and suggestions for choosing when to use
them, see Lyness (1983). The following are from his concluding remarks.
The Automatic Quadrature Rule (AQR)
is
an impressive and practical item
of numerical software. Its main advantage for the
user is that it is conveni-
ent.
He
can take it from the library shelf, plug
it
in, and feel confident that
it will work. For this convenience, there
is
in general
.a
modest charge in
CPU
time, this surcharge being
.a
factor of about
3.
The Rule Evalua;ion
Quadrature Routine (REQR) [non-automatic quadrature rule] does not
carry this surcharge, but to
code and check
out
an
REQR
might take a
SINGULAR INTEGRALS 305
couple
of
hours
of
the user's time. So unless the expected
CPU
time is high,
many user's willingly pay the surcharge in
order
to save themselves time
and
trouble.
However there are
certain-usually
large
scale-problems
for which the
AQR
is
not
designed and in which its uncritical use can lead to
CPU
time
surcharges by factors of
100
or
more. . . . These are characterized by the
circumstances that a large number of separate quadratures are involved,
and
that
the results of these quadratures are subsequently used as input to
some other numerical process.
In
order to recognize this situation, it is
necessary to examine the subsequent numerical process to see whether it
requires a smooth input function. . . .
For
some
of
these problems, an
REQR
is quite suitable while
an
AQR
may lead to a numerical disaster.
5.6 Singular Integrals
We
discuss the approximate evaluation
of
integrals for which methods
of
the type
discussed
in
Sections
5.1
through 5.4
do
not perform well: these methods include
the composite Newton-Cotes rules (e.g., the trapezoidal rule),
Gauss-Legendre
quadrature,
and
Romberg integration. The integrals discussed here lead to poorly
convergent numerical integrals when evaluated using the latter integration rules,
for a variety
of
reasons. We discuss (1) integrals whose integrands contain a
singularity
in
the interval
of
integration
(a,
b),
and
(2) integrals with an infinite
interval
of
integration. Adaptive integration methods can be used for these
integrals,
but
it
is usually possible to obtain more rapidly convergent approxima-
tions
by
carefully examining the nature
of
the singular behavior
and
then
compensating for it.
Change of the variable of integration We illustrate the importance
of
this idea
with several examples.
For
1
=
ibf(x)
dx
o
IX
(5.6.1)
with
f(x)
a function with several continuous derivatives, let x = u
2
,
0
::s;
u
::s;
lb.
Then
.
This integral has a smooth integrand and
standard
techniques
can
be applied
to
it.
Similarly,
using the change
of
variable u =
..;r-::x.
The
right-hand integrand has
an
infinite
number
of
continuous derivatives
on
[0,
1], whereas the derivative
of
the
first integrand was singular
at
x = 1.
i
-i
I
-····--·-----
-----
·-----------------------·
_,._
--
------~--------
----~
'
'
!
I
____
_j
306 NUMERICAL INTEGRATION
For
an infinite interval of integration, the change of variable technique
is
also
useful. Suppose
p>1
(5.6.2)
with Limitx
...
oo
f(x)
existing. Also assume
f(x)
is
smooth on
[1,
oo).
Then use
the change of variable
Then
-a
dx =
--du
ul+a
for some
a>
0
(5.6.3)
Maximize the smoothness of the new integrand
at
u = 0 by picking a to produce
a large value
for the exponent
(p-
1)a-
1.
For
example, with
1
=
joof(x)
dx·
1 x/X
the change
of
variable x =
1ju
4
leads to
If
we assume a behavior at x =
oo
of
then
and
(5.6.4) has a smooth integrand at u =
0.
An
interesting idea has been given in Iri et al. (1970) to deal with endpoint
singularities in the integral
Define
1/J(t),;,
exp(~)
1-
t
b-
a
1
<p(t)
=a+
--j
1/J(u)
du
y
-1
(5.6.5)
(5.6.6)
-1:;;;
t:;;; 1
(5.6.7)
SINGULAR
INTEGRALS 307
where c
is
a positive constant and
y = r
1/J(u)
du
-1
As t varies from
-1
to
1,
<p{t) varies from a to
b.
Using x = <p{t) as a change
of variable in
(5.6.5),
we
obtain
I=
t /(!J!(t))!J!'(t)
dt
-1
(5.6.8)
The function qJ'(t) =
((b-
a)jy)I/J(t)
is
infinitely differentiable on [
-1,
1], and
"it
and
all of derivatives are zero
at
t = ±
1.
In (5.6.8), the integrand and all of
derivatives will vanish
at
t = ± 1 for virtually all functions
f(x)
of interest.
Using the error formula (5.4.9) for the trapezoidal rule
on
[
-1,
1], it can be seen
that the trapezoidal rule will converge very rapidly when applied to
(5.6.~).
We
will call this method the
IMT
method.
This method has been implemented in de Doncker
and
Piessens (1976), and
in
the general comparisons
of
Robinson (1979), it is rated as an extremely reliable
and quite efficient way of handling integrals
(5.6.4) that have endpoint singulari-
ties.
De
Doncker and Piessens (1976) also treat integrals over
[0,
oo)
by first using
the change
of
variable x =
(1
+
u)/(1
- u),
-1
~
u < 1, followed by the change
of
variable u = <p(t).
Example Use the preceding method (5.6.5)-(5.6.8) with the trapezoidal rule, to
evaluate
1
1
.;;
I=
...j-Inx
dx
=-
~
.8862269
0 2
(5.6.9)
Note that the integrand has singular behavior
at
both
endpoints, although it is
different in the two cases. The constant in
(5.6.6) is c =
4,
and the evaluation of
(5.6.7) is taken from Robinson and de Doncker (1981). The results are shown in
Table 5.21. The column labeled nodes
gives
the
number
of
nodes interior to
[0,
1].
Table 5.21 Example of the
IMT
method
Nodes Error
2
-6.54E-
2
4
5.82E-
3
8
-1.30E-
4·
16
7.42E-
6
32
1.17E-
8
64
1.18E-
12
i
1
I
_______
_\
308 NUMERICAL INTEGRATION
Gaussian quadrature In Section
5.3,
we
developed a general theory for Gauss-
ian quadrature formulas
n~1
that have degree of precision 2n -
1.
The construction of the nodes and weights,
and the form of the error, are given
in Theorem 5.3. For our work in this section
we
note
that
(1) the interval
(a,
b)
is
allowed to be infinite, and (2)
w(
x)
can
have singularities on
(a, b), provided it
is
nonnegative and satisfies the assump-
tions (4.3.8) and (4.3.9) of Section 4.3. For rapid convergence,
we
would also
expect that
f(x)
would need to be a smooth function, as was illustrated with
Gauss-Legendre quadrature
in Section 5.3.
The weights and nodes for a wide variety
of
weight functions w(x) and
intervals
(a,
b)
are known. The tables of Stroud and Secrest (1966) include the
integrals
fIn
(
~
)t(x)
dx (5.6.10)
and
others. The constant a >
-1.
There are additional books containing tables
for integrals other than those
in
(5.6.10). In addition, the paper by Golub and
Welsch
(1969) describes a procedure for constructing the nodes and weights in
(5.6.10), based on solving a matrix eigenvalue problem. A program
is
given, and
it includes most of the more popular weighted integrals to which Gaussian
quadrature is applied. For an additional discussion of Gaussian quadrature, with
references
to
the literature (including tables
and
programs), see Davis and
Rabinowitz (1984, pp. 95-132, 222-229).
Example We illustrate the use of Gaussian quadrature for evaluating integrals
I=
loco
g(x)
dx
We use Gauss-Laguerre quadrature, in which
w(x)
=e-x.
Then write I
as
(5.6.11)
We give results for three integrals:
oo
xdx
.,
2
1
<1>
= r ___ = _
lo
ex-
1 6
xdx
1
J(2)-
(
00
=-
-
lo
(1
+ x
2
)
5
8
oo
dx
.,
1
<3>
= r
·
lo
1 + x
2
2