Atkinson K. An Introduction to Numerical Analysis
Подождите немного. Документ загружается.
...
----- -
·--
-
.
·-
--
·-
-·--··
..
NUMERICAL DIFFERENTIATION 319
The
actual numerical derivative computed
is
(5.7 .20)
For
its error, substitute (5.7.19) into (5.7.20), obtaining
-h
2
£ -
2£
+ (
=
--t<4>(~)
+
2.
1·
0
12 h
2
(5.7.21)
For
the term involving { t:;
},
assume these errors are random within some
interval -
E
:$
t:
:$E.
Then
(5.7.22)
and the last bound would be attainable in many situations. An example
of
such
errors would
be
rounding errors, with E a bound
on
their magnitude.
The
error bound in (5. 7.22) will initially get smaller as h decreases,
but
for h
sufficiently close to zero, the error will begin to increase again. There is
an
optimal value of h, call it h*, to minimize
the
right side of (5.7.22),
and
presumably there is a similar value for the actual error
/"(x
1
) -
D~
2
lf(x
1
).
Example Let
f(x)
=
-cos(x),
and compute
/"(0)
using the numerical ap-
proximation (5.7.17).
In
Table 5.25,
we
give the errors in (1) Dflf(O), computed
exactly, and (2)
Jj~
2
>J(O),
computed using 8-digit rounded decimal arithmetic.
In
Table 5.25 Example of
D~
2
>J(O)
and
D~
1
>j(O)
h f"(O) -
D~
2
>J(O)
Ratio
f"(O) -
bJ?>J(O)
.5
2.07E-
2
2.07E-
2
.25
5.20E-
3 3.98
5.20E-
3
.125
1.30E-
3 3.99
1.30E-
3
.0625
3.25E-
4 4.00
3.25E-
4
.03125
8.14E-
5 4.00
8.45E-
5
.015625
2.03E-
5 4.00
2.56E-
6
.0078125
5.09E-
6 4.00
-7.94E-
5
.00390625
1.27E-
6
4.00
-7.94E-
5
.001953125
3.18E-
7
4.00
-1.39E-
3
320
NUMERICAL
INTEGRATION
this last case,
(5.7.23)
This
bound
is minimized
at
h* = .0022, which is consistent with the errors
f"(O) - bf>J(O) given in the table. For the exactly computed
Dh
2
'f(O), note that
the errors decrease by four whenever h is halved, consistent with the error
formula (5.7.18).
Discussion
of
the Literature
Even though the topic
of
numerical integration is one
of
the
oldest in numerical
analysis
and
there is a very large literature, new papers continue to appear
at
a
fairly high rate. Many
of
these results give methods for special classes
of
problems, for example, oscillatory integrals,
and
others
are
a response to changes
in computers, for example, the use of vector pipeline architectures. The best
survey
of
numerical integration
is
the large
and
detailed work
of
Davis
and
Rabinowitz (1984).
It
contains a comprehensive survey
of
most quadrature
methods, a very extensive bibliography, a set
of
computer programs,
and
a
bibliography
of
published quadrature programs.
It
also contains the article
"On
the practical evaluation of integrals" by Abramowitz, which gives some excellent
suggestions
on
analytic approaches to quadrature. Other important texts in
numerical integration are Engels (1980), Krylov (1962),
and
Stroud (1971).
For
a
history of the classical numerical integration methods, see Goldstine (1977).
For
reasons
of
space, we have had to omit some important ideas. Chief among
these are (1) Clenshaw-Curtis quadrature,
and
(2) multivariable quadrature.
The
former
is
based
on
integrating a Chebyshev expansion
of
the integrand; em-
pirically the method has proved excellent for a wide variety
of
integrals.
The
original method is presented in Clenshaw and Curtis (1960); a current account
of
the method is given in Piessens et
al.
(1983, pp.
28-39).
The
area
of
multivariable
quadrature is
an
active area
of
research,
and
the
texts
of
Engels (1980) and
Stroud (1971) are the best introductions to the area. Because
of
the widespread
use
of
multivariable quadrature in the finite element method for solving partial
differential equations, texts on the finite element method will often contain
integration formulas for triangular and rectangular regions.
Automatic numerical integration was a very active area of research in the
1960s
and
1970s, when it was felt that most numerical integrations could be done
in
this way. Recently, there has been a return
to
a greater use of nonautomatic
quadrature especially adapted to the integral
at
hand.
An
excellent discussion
of
the
relative advantages
and
disadvantages
of
automatic quadrature is given in
Lyness (1983).
The
most powerful and flexible
of
the current automatic programs
are probably those given in QUADPACK, which is discussed and illustrated
in
Piessens et al. (1983). Versions
of
QUADPACK
are included in the IMSL
and
NAG
libraries.
-----
---
----·---·
--------------
-·
·----·
----
·--
-
·-
..\
BIBLIOGRAPHY 321
For
microcomputers and hand computation, Simpson's rule
is
still popular
because of its simplicity. Nonetheless, serious consideration should be given to
Gaussian quadrature because of its much greater accuracy. The nodes and
weights are readily available, in Abramowitz
and
Stegun (1964) and Stroud and
Secrest (1966), and programs for their calculation are also available.
Numerical differentiation
is
an ill-posed problem in the sense of Section 1.6.
Numerical differentiation procedures that account for this have been developed
in the past ten to fifteen years. In particular, see Anderssen and Bloomfield
(1974a), (1974b), Cullum (1971), Wahba (1980),
and
Woltring (1986).
Bibliography
Abramowitz, M., and
I.
Stegun, Eds. (1964). Handbook
of
Mathematical Func-
tions. National Bureau of Standards, U.S. Government Printing Office,
Washington,
D.C.
Anderssen, R., and P. Bloomfield (1974a), Numerical differentiation procedures
for non-exact data,
Numer. Math.,
22,
157-182.
Anderssen, R., and
P.
Bloomfield (1974b). A time series approach to numerical
differentiation,
Technometrics
16,
69-75.
Atkinson,
K.
(1976). A Survey
of
Numerical Methods for the Solution
of
Fredholm
Integral Equations
of
the Second Kind. Society for Industrial and Applied
Mathematics, Philadelphia.
Atkinson,
K.
(1982). Numerical integration on the sphere, J. Austr. Math. Soc.
(Ser.
B)
23, 332-347.
Bauer, F., H. Rutishauser,
and
E.
Stiefel (1963). New aspects in numerical
quadrature. In
Experimental Arithmetic, High Speed Computing, and
Mathematics, pp. 199-218. Amer. Math. Soc., Providence, R.I.
de Boor,
C. (1971). CADRE: An algorithm for numerical quadrature.
In
Mathematical Software, pp. 201-209. Academic Press, New York.
Clenshaw, C., and
A.
Curtis (1960). A method for numerical integration on an
automatic computer,
Numer. Math.
2,
197-205.
Cryer,
C.
(1982). Numerical Functional Analysis. Oxford Univ. Press (Clarendon),
Oxford, England.
Cullum, J. (1971). Numerical differentiation and regularization,
SIAM
J. Numer.
Anal.
8,
254-265.
Davis,
P.
(1963). Interpolation and Approximation.
Ginn
(Blaisdell), Boston.
Davis, P.,
and
P.
Rabinowitz (1984). Methods
of
Numerical Integration, 2nd ed.
Academic Press, New York.
Dixon,
V.
(1974). Numerical quadrature: A survey
of
the available algorithms.
In
Software for Numerical Mathematics, D. Evans, Ed., pp. 105-137. Academic
Press, London.
322 NUMERICAL INTEGRATION
Donaldson,
J.,
and
D. Elliott (1972). A unified
approach
to
quadrature
rules with
asymptotic estimates of their remainders,
SIAM
J. Numer. Anal.
9,
573-602.
de Doncker, E.,
and
R. Piessens (1976). Algorithm 32:
Automatic
computation
of
integrals with singular integrand, over a finite
or
an infinite interval,
Computing 17, 265-279.
Engels, H. (1980).
Numerical Quadrature and
C~bature.
Academic Press, New
York.
Goldstine, H. (1977). A History
of
Numerical Analysis. Springer-Verlag, New
York.
Golub,
G.,
and
J. Welsch (1969). Calculation
of
Gauss
quadrature
rules, Math.
Comput.
23, 221-230.
de
Hoog, F.,
and
R. Weiss (1973). Asymptotic expansions for
product
integra-
tion,
Math. Comput. 27, 295-306.
Iri, M.,
S.
Moriguti,
andY.
Takasawa (1970).
On
a numerical integration formula
(in Japanese), J. Res. lost.
Math. Sci., 91, 82. Kyoto Univ., Kyoto, Japan.
Isaacson, E.,
and
H. Keller (1966). Analysis
of
Numerical Methods. Wiley, New
York.
Kronrod,
A.
(1965). Nodes and Weights
of
Quadrature Formulas, Consultants
Bureau,
New
York.
Krylov,
V.
(1962). Approximate Calculation
of
Integrals. Macmillan, New York.
Lyness, J. (1983). When
not
to use an automatic
quadrature
routine,
SIAM
Rev.
25,
63-88.
Lyness, J.,
and
C.
Moler (1967). Numerical differentiation
of
analytic functions,
SIAM
J. Num. Anal.
4,
202-210.
Lyness, J.,
and
J. Kaganove (1976). Comments
on
the
nature
of
automatic
quadratic
routines,
ACM
Trans. Math. Softw. 2,
65-81.
Lyness, J.,
and
K.
Puri (1973).
The
Euler-MacLaurin expansion for the simplex,
Math. Comput. 27, 273-293.
Lyness, J.,
and
B.
Ninham
(1967). Numerical
quadrature
and
asymptotic expan-
sions, Math. Comput. 21, 162-178.
Patterson,
T.
(1968).
The
optimum
addition
of
points
to.
quadrature
formulae,
Math. Comput. 22, 847-856. (Includes a microfiche enclosure; a microfiche
correction is in 23, 1969.)
Patterson, T. (1973) Algorithm 468: Algorithm for numerical integration over a
finite interval.
Commun.
ACM
16, 694-699.
Piessens, R., E. deDoncker-Kapenga,
C.
Oberhuber,
and
D.
Kahaner
(1983).
QUADPACK: A Subroutine Package
for.
Automatic Integration. Springer-
Verlag,
New
York.
Ralston, A. (1965). A First Course in Numerical Analysis. McGraw-Hill, New
York.
;
\
----
__
....;
PROBLEMS 323
Robinson,
I.
(1979).
A comparison of numerical integration programs,
J.
Com-
put. Appl. Math.
5,
207-223.
Robinson,
1.,
and
E.
de Doncker (1981). Automatic computation of improper
integrals over a bounded or unbounded planar region,
Computing
27,
253-284.
Stenger, F. (1981). Numerical methods based on Whittaker cardinal or sine
functions,
SIAM
Rev.
23,
165-224.
Stroud,
A.
(1971). Approximate Calculation
of
Multiple Integrals. Prentice-Hall,
Englewood
Cliffs,
N.J.
Stroud,
A., and D. Secrest (1966). Gaussian Quadrature Formulas. Prentice-Hall,
Englewood
Cliffs,
N.J.
Wahba, G. (1980). Ill-posed problems: Numerical and statistical methods for
mildly, moderately, and severely ill-posed problems with noisy data, Tech.
Rep.
#
595,
Statistics Dept., Univ. of Wisconsin, Madison. (Prepared for
the Proc. Int. Symp. on Ill-posed Problems, Newark, Del., 1979.)
Woltring, H.
(1986).
A Fortran package for generalized, cross-validatory spline
smoothing and differentiation,
Adv. Eng. Softw.
8,
104-113.
Problems
1.
Write a program to evaluate
I=
J:f(x)
dx
using the trapezoidal rule with
n subdivisions, calling the result ln. Use the program to calculate the
following integrals with
n =
2,
4,
8,
16,
...
, 512.
(d)
(c)
l
z,.
dx
o
2+cos(x)
(e)
fa""
ex
cos(4x)
dx
Analyze empirically the rate of convergence of
In
to I by calculating
the ratios of (5.4.27):
2.
Repeat Problem 1 using Simpson's rule.
3. Apply the corrected trapezoidal rule (5.1.12) to the integrals in Problem 1.
Compare the results with those of Problem 2 for Simpson's rule.
324
NUMERICAL INTEGRATION
4. As another approach
to
the corrected trapezoidal rule (5.1.12), use the cubic
Hermite interpolation polynomial to
f(x)
to obtain
( )
2
h
b-a
b-a
£
J(x)
dx
~
(-
2
-)[/(a)
+/(b)]
-
12
[/'(b)-
f'(a)]
Use the error formula for Hermite interpolation to obtain an error formula
for the preceding approximation. Generalize these results to
(5.1.12), with n
subdivisions of [a,
b].
5. (a) Assume that
f(x)
is
continuous and that
f'(x)
is
integrable on
[0,
1].
Show that the error in the trapezoidal rule for calculating
/JJ(x)
dx
has the form
tj-1
+
tj
K(t} =
2
-
t
j = 1,
...
, n
This contrasts with (5.1.22) in which
/'(x)
is continuous and
f"(x)
is
integrable.
(b) Apply the result to
f(x)
=
xa
and to
f(x)
=
xa
log(x), 0
<a
< 1.
This gives an order
of
convergence, although it is less than the true
order.
[See
Problem 6 and (5.4.23).]
6. ·Using the program of Problem 1, determine empirically the rate of conver-
gence of the trapezoidal rule applied to
JJxa
ln(x)
dx, 0
~
a~
1,
with a
range
of
values of a, say a = .25, .5, .75, 1.0.
7. Derive the composite form of Boole's rule, which is given as the n = 4 entry
in Table
5.8. Develop error formulas analogous to those given in (5.1.17)
and (5.1.18) for Simpson's rule.
8. Repeat
Problem 1 using the composite Boole's rule obtained in Problem
7.
9. Let
P2(x)
be the quadratic polynomial interpolating
/(x)
at x =
0,
h, 2h.
Use this to derive a numerical integration formula
Ih
for 1 =
Mhf(x)
dx.
Use a Taylor series expansion of
/(x)
to show
PROBLEMS 325
10.
For
the midpoint integration formula (5.2.17), derive the Peano kernel error
formula
{
1
2
-t
K(t)
=
~
2(h-
t)2
h
0
~
t
~
2
h
- < t
<h.
2-
-
Use this to derive the error term in (5.2.17).
11.
Considerfunctions defined in the following way. Let n >
0,
h =
(b-
a)jn,
ti
= a +
jh
for j =
0,
1,
...
, n. Let
f(x)
be
linear on each subinterval
[ti_
1
,t),
for j =
1,
...
,
n.
Show that the set
of
all such f, for all n
~
1,
is
dense in C[a,
b].
12.
Let
w(x)
be an integrable function on [a, b],
jblw(x)!dx
<
oo
a
and let
b n
j
w(x)f(x)
dx
~
L
wj,nf(xj,n)
a
j=l
be a sequence of numerical integration rules. Generalize Theorem
5.2
to
this case.
1~.
Assume the numerical integration rule
is convergent for all continuous functions. Consider the effect of errors in
the function values. Assume
we
use[;~
f(x;),
with
What
is
the effect on the numerical integration
of
these errors in the
function values?
14.
Apply Gauss-Legendre quadrature
to
the integrals
in
Problem
1.
Compare
the results with
those for the trapezoidal and Simpson methods.
15.
Use Gauss-Legendre quadrature
to
evaluate t:
..
4
dxj(1
+ x
2
)
with the
n =
2,
4,
6, 8 node-point formulas. Compare with the results in Table
-5.9,
obtained using the Newton-Cotes formula.
326 NUMERICAL INTEGRATION
16.
Prove that the Gauss-Legendre nodes and weights on [
-1,
1]
are symmetri-
cally distributed about the origin
x =
0.
17. Derive the two-point Gaussian quadrature formula for
I(!)=
{t{x)log(~)
dx
in which the weight function
is
w(x)
=
log(l/x).
Hint:
See
Problem 20(a)
of
Chapter
4.
Also, use the analogue of (5.3.7) to compute the weights, not
formula (5.3.11).
18.
Derive the one- and two-poini Gaussian quadrature formulas for
1 n
I=
1
xf(x)
dx,;,
L
wJ(x)
0
j=l
with weight function
w(x)
= x.
19.
For
the integral
I=
f~
1
Vl-
x
2
f(x)
dx with weight
w(x)
=VI-x
2
,
find explicit formulas
for
the nodes and weights of the Gaussian quadrature
formula. Also
give
the error formula. Hint:
See
Problem
24
of Chapter
4.
20. Using the column en.
4
of Table 5.13, produce the fourth-order error
formula analogous to the second-order formula (5.3.40) for
en.
2
•
Compare
it with the fourth-order Simpson error formula
(5.1.17).
21. The weights in the Kronrod formula (5.3.48) can
be
calculated as the
solution of four simultaneous linear equations. Find that system and then
solve it to verify the values given following (5.3.48).
Hint: Use the approach
leading to (5.3.7).
22. Compare the seven-point Gauss-Legendre formula with the seven-point
Kronrod formula of (5.3.48).
Use each
of
them on a variety of integrals,
and then compare their respective errors.
23. (a) Derive the relation (5.4.7) for the Bernoulli polynomials
B/x)
and
Bernoulli numbers
Bj. Show that Bj = 0 for all
odd
integers j
~
3.
(b) Derive the identities
Bj(x)
=
jBj_
1
(x)
j
~
4 and even
j
~
3 and odd
These can be used to give a general proof of the Euler-MacLaurin
formula (5.4.9).
PROBLEMS 327
24.
Using the Euler-MacLaurin summation formula (5.4.17), obtain an esti-
mate of
t(
1
),
accurate to three decimal places. The zeta function
t(
p)
is
defined in (5.4.22).
25. Obtain the asymptotic error formula for the trapezoidal rule applied to
f~
':;xj(x) dx. Use the estimate from Problem
24.
26. Consider the following table of approximate integrals
In
produced using
Simpson's rule. Predict the order of convergence of
In
to
I:
n
In
2 .28451779686
4 .28559254576
8 .28570248748
16
.28571317731
32
.28571418363
64
.28571427643
That is, if
I-
In
= cjnP, then what
is
p?
Does this appear to be a valid
form for the
error for these data? Predict a value of c and the error in I
64
•
How large should n be chosen if
In
is
to be in error by less than
w-u?
27. Assume that the error in an integration formula has the asymptotic expan-
sion
Generalize the Richardson extrapolation process of Section 5.4 to obtain
formulas for C
1
and C
2
.
Assume that three values
In,
I
2
n,
and I
4
n have
been computed, and use these to compute C
1
,
C
2
,
and an estimate of I,
with an error of order 1 j n
2../n.
28. For the trapezoidal rule (denoted by IY>) fofevaluating
I=
J:J(x) dx,
we
have the asymptotic error formula
and for the midpoint formula
I~M>,
we
have
provided
f is sufficiently differentiable on [a,
b].
Using these results, obtain
a new numerical integration formula
I:
combining
IY>
and
I~M>,
with a
higher order of convergence. Write out the weights to the new formula
in-
328
NUMERICAL
INTEGRATION
29.
Obtain
an
asymptotic error formula for Simpson's rule, comparable
to
the
Euler-MacLaurin formula
(5.4.9) for the trapezoidal rule. Use (5.4.9),
(5.4.33), and (5.4.36),
as
in (5.4.37).
30. Show
that
the formula (5.4.40) for
I~
2
l
is
the composite Boole's rule. See
Problem
7.
31. Implement the algorithm Romberg of Section 5.4, and then apply it to the
integrals of
Problem
1.
Compare the results with those for the trapezoidal
·
and
Simpson rules.
32. Consider evaluating
I =
JJxa
dx by using the adaptive Simpson's rule that
was described following
(5.5.3). To see that the error test (5.5.5) may fail,
consider the case of
-'-1
< a <
0,
with the integrand arbitrarily set to zero
when
x =
0.
Show that for sufficiently small t:, the test (5.5.5) will never be
satisfied for the subinterval
[0,
h]. Note that for
t:
specified on
[0,
1],
the
error tolerance for
[0,
h]
will be
t:h.
33. Use Simpson's rule with even spacing to integrate
I=
JJ
log(x)
dx.
For
x = 0, set the integrand to zero. Compare the results with those in Tables
5.19
and
5.20, for integral I
4
'.
34. Use
an
adaptive integration program [e.g.,
DQAGP
from Piessens et al.
(1983)]
to
calculate the integrals in Problem
1.
Compare the results with
those
of
Problems
1,
2,
14, and
31.
35. Decrease the singular behavior of the integrand in
I=
f!(x)log(x)
dx
by using the change of variable x = t ', r >
0.
Analyze the smoothness
of
the resulting integrand. Also explore the empirical behavior of the
trapezoidal and
Simpson rules for various
r.
36. Apply the
IMT
method, described following (5.6.5), to
the
calculation of
I=
f
0
00
J(x)
dx, using some change of variable from
[0,
oo)
to a finite
interval. Apply this to the calculation of the integrals following
(5.6.11).
37. Use Gauss-Laguerre quadrature with n = 2,
4,
6,
and 8 node points to
evaluate the following integrals.
Use (5.6.11) to
put
the integrals in proper
form.
1
00
2
.;;
(a)
0
e-x
dx = 2
(c)
1
oo
sin
(x)
'IT
---dx
=-
0 X 2
oo
xdx
1
(b) 1
o
(1
+ x
2
)
2
2