Atkinson K. An Introduction to Numerical Analysis
Подождите немного. Документ загружается.
NEWTON-COTES INTEGRATION FORMULAS
269
Combining these results, the weights w
1
.
n must vary in sign
as
n becomes
sufficiently large. For example, using
n =
8,
Such formulas can cause loss-of-significance errors, although it
is
unlikely to be a
serious problem until
n
is
larger. But because of this problem, people have
generally avoided using Newton-Cotes formulas for
n
~
8,
even in forming
composite formulas.
The most serious problem of the Newton-Cotes method
(5.2.2)
is
that it may
not converge for perfectly well-behaved integrands,
as
in
(5.2.9).
The midpoint rule There are additional Newton-Cotes formulas in which one
or both
of
the endpoints of integration are deleted from the interpolation (and
integration) node points. The best known of these is also the simplest, the
midpoint rule.
It
is
based on interpolation of the integrand
f(x)
by
the constant
f((a
+
b)/2);
and the resulting integration formula
is
f
h
(a+b)
(b-a)
3
f(x)dx=(b-a)f
--
+ /"(7J)
a 2
24
some
7J
E [ a,
b]
(5.2.17)
For
its composite form, define
x
1
=a+(j-!)h
j=1,2,
...
,n
the midpoints of the intervals
[a+
(j-
1)h,
a+
jh].
Then
jhf(x)
dx
=In(/)+
En(/)
a
(5.2.18)
some
7J
E [
a,
b]
(5.2.19)
The proof
of
these results
is
left as Problem 10.
These integration formulas in which one or both endpoints are missing are
called
open
Newton-Cotes formulas, and the previous formulas are called closed
formulas. The open formulas of higher order were used classically in deriving
numerical formulas for the solution of ordinary differential equations.
i
I
i
i
I
I
i
I
_j
270 NUMERICAL INTEGRATION
5~3
Gaussian
Quadrature
The composite trapezoidal and Simpson rules are based on using a low-order
polynomial approximation of the integrand
f(x)
on subintervals of decreasing
size.
In
this section,
we
investigate a class of methods that use polynomial
approximations of
f(x)
of increasing degree. The resulting integration formulas
are extremely accurately in most cases, and they should be considered seriously
by anyone faced with many integrals to evaluate.
For
greater generality,
we
will
consider formulas
n b
I,(!)=
L
wj_,f(xj,,)
= j
w(x)f(x)
dx
=I(/)
j=1
a
(5.3.1)
The weight function
w(x)
is
assumed to be nonnegative and integrable on [a, b],
and it is to also satisfy the hypotheses (4.3.8) and (4.3.9) of Section 4.3. The
nodes {
xj.,}
and weights { wj,,} are to be chosen so that
/,(/)
equals I
(f)
exactly for polynomials
f(x)
of
as
large a degree as possible.
It
is
hoped that this
will result in a formula
/,(/)
that
is
nearly exact for integrands
f(x)
that are
well approximable by polynomials. In Section 5.2, the Newton-Cotes formulas
have an increasing degree of precision as
n increased, but nonetheless they do not
converge for many well-behaved integrands. The difficulty with the Newton-Cotes
formulas is
that the nodes
{xi.,}
must be evenly spaced.
By
omitting this
restriction, we
will
be able to obtain new formulas
/,(/)
that converge for all
f E
C[a,
b].
To
obtain some intuition for the determination of
/,(/),
consider the special
case
1 n
J
J(x)
dx = L
wJ(x)
-1
j=l
(5.3.2)
where
w(x)
= 1 and the explicit dependence of {
wi}
and {xi} on n has been
dropped.
The
weights {
wi}
and nodes {xi} are
to
be determined to make the
error
(5.3.3)
equal zero for as high a degree polynomial
f(x)
as
possible. To derive equations
for the nodes and weights,
we
first note that
Thus
E,(f)
= 0 for every polynomial of degree
:s:
m if and only if
i=O,l,
...
,m
(5.3.5)
GAUSSIAN QUADRATURE
271
Case
1.
n = 1. Since there are two parameters, w
1
and
x
1
,
we
consider
requiring
This gives
J
l .
1dx-
w
1
= 0
-1
This implies w
1
= 2 and x
1
=
0.
Thus the formula (5.3.2) becomes
f
f(x)
dx
~
2/(0)
-1
the midpoint rule.
Case 2. n = 2. There are four parameters, w
1
, w
2
, x
1
, x
2
, and thus
we
put
four constraints on these parameters:
i=0,1,2,3
or
These equations lead to the unique formula
1 (
{3)
({3)
f_J(x)
dx
~1
-3
+ f 3
(5.3.6)
which has degree of precision three. Compare this with Simpson's rule (5.1.13),
which uses three nodes to attain the same degree
of
precision.
Case 3.
For
a general n there are 2n free parameters {X;} and {
W;
},
and we
would guess
that
there
is
a formula (5.3.2) that uses n nodes and gives a degree of
precision of
2n -
1.
The equations to be solved are
i =
0,1,
...
,2n
-1
or
n
{0
L
wjxj=
_2_
j=l
i + 1
i =
1,3,
...
,2n-
1
i =
0,2,
...
,2n-
2
(5.3.7)
I
I
---
---------
------------------------
-----------
-
______
__\
272 NUMERICAL INTEGRATION
These are nonlinear equations, and their solvability
is
not at
all
obvious. Because
of the difficulty in working with this nonlinear system,
we
use
another approach
to the theory for
(5.3.2), one that
is
.somewhat circuitous.
Let
{qJ,(x)ln;;::
0}
be the orthogonal polynomials on
(a,
b)
with respect to
the weight function
w(x);;::
0.
Denote the zeros of rp,(x) by
a < x
1
< · · · <
x,
< b
(5.3.8)
Also, recall the notation from (4.4.18)-(4.4.20):
rp,(x) =
A,x"
+ · · ·
(5.3.9)
Theorem 5.3 For each n
;;::
1,
there
is
a unique numerical integration formula
(5.3.1) of degree of precision
2n
-
1.
Assuming
f(x)
is
2n
times
continuously differentiable on
[a, b], the formula for
I,(f)
and its
error
is
given by
for some
a<
11
<b.
The nodes
{xJ
are the zeros ot rp,(x), and
the weights {
wj}
are given by
j = 1,
...
, n
(5.3.11)
Proof The proof
is
divided into three parts. We first obtain a formula with
degree
of
precision
2n
-
1,
using the nodes (5.3.8). We then show that it
is
unique. Finally,
we
sketch the derivation of the error formula and the
wei~t~.
(a) Construction of
the
formula. Hermite interpolation
is
used as the
vehicle for the construction (see Section
3.6
to review the notation and
results). For the nodes in
(5.3.8), the Hermite polynomial interpolating
f(x)
and
f'(x)
is
II
II
H,(x)
=
2:
fix)h;{x) +
2:
f'(x)hj(x)
(5.3.12)
j=l j=l
with
hix)
and
hix)
defined in (3.6.2)
of
Section 3.6. The interpolation
error is given by
tf,(x)
=J(x)
-H,.(x)
= [o/,.(x)]
2
/[x
1
,x
1
,
...
,x,.,x,.,x]
(
~~~(~)]
2
J<
2
">(E)
~
E (a,
b]
(5.3.13)
2n
!
·----------·
---
-···
·-----------·--·
----····.
------···
------
--
...
..J
GAUSSIAN QUADRATURE
273
with
o/n(x) =
{x-
xJ
· · ·
{x-
xJ
Note that
{5.3.14)
since both
cpn(x)
and
o/n(x)
are of degree n and have the same zeros.
Using (5.3.12), if
f(x)
is
continuously differentiable, then
Jhw(x)f(x)
dx
= Jbw(x)Hn(x)
dx
+ jbw(x)tffn(x)
dx
a a a
(5.3.15)
The degree of precision is at least 2n - 1, since tfn(x) = 0 if
f(x)
is
a
polynomial of degree <
2n, from (5.3.13). Also from (5.3.13),
En(x
2
n)
=
tw(x)tfn(x)
dx
=
jhw(x)[
o/n(x)]
2
dx
> 0 (5.3.16)
a a
Thus the degree of precision of
In(/)
is exactly
2n
-
1.
To derive a simpler formula for In(/),
n b n b
In(/)=
L
/(xi)
j
w(x)hi(x)
dx
+ L
f'(x)
J
w(x)h/x)
dx
j=l
a
j=l
a
{5.3.17)
we
show that all of the integrals in the second sum are zero. Recall that
from
(3.6.2),
h/x)
=
(x-
x)[l/x)]
2
l/x)
=
1/Jn(x
)'
CJln(x)
(x-
x)tfln(xi)
(x-
x)cp~(xi)
The last step
uses
(5.3.14). Thus
(5.3.18)
Since degree(/i)
=
n-
1, and since
CJln(x)
is orthogonal
to
all polynomi-
als
of
degree<
n,
we
have
j = 1,
...
, n
(5.3.19)
274 NUMERICAL INTEGRATION
The integration formula (5.3.15) becomes
b n
1
w(x)f(x)
dx = L
wJ(x)
+En(!)
a
j=1
~-=
1bw(x)hj(x)dx
j=
1,
...
,n
a
(5.3.20)
(b) Uniqueness of formula (5.3.19). Suppose that
we
have a numerical
integration formula
(5.3.21)
that has degree of precision
~
2n -
1.
Construct the Hermite interpola-
tion formula to
f(x)
at
the nodes z
1
,
•••
,
zn.
Then for any pqlynomial
f(x)
of
degree~
2n-
1,
n n
f(x)=
L
f(z)h/x)+
L
f'(z)hj(x)
deg
(!)
~
2n-
1 {5.3.22)
j=l
j=l
where
h/x)
and
hix)
are defined using {zj}· Multiply (5.3.22) by
w(x), use the assumption on the degree of precision of (5.3.21), and
integrate to get
n n b n b
L
vjf(zj)
= L
f(z)
1
w(x)hix)
dx + L
f'(zj)
J
w(x)hix)
dx
j=l
j=l
a
j=1
a
(5.3.23)
for any polynomial
f(x)
of
degree~
2n
-
1.
Let
f(x}
= hi(x). Use the properties (3.6.3)
of
hi(x)
to obtain from
(5.3.23) that
1
b -
0 =
w(x)hi(x)
dx
a
i = 1,
...
, n
As
before in (5.3.18),
we
can write
w(x)
=
(x-
z
1
)
· • ·
(x-
zJ
Then (5.3.24) becomes
jbw(x)wn(x)l;(x)
dx = 0
a
i=1,2,
...
,n
(5.3.24)
;
-
_j
(iAUSSIAN
QUADRATURE
275
Since all polynomials of degree
.:;;;
n - 1 can be written as a combination
of /
1
(x),
...
, /n(x),
we
have that wn(x)
is
orthogonal to every polynomial
of
degree
.:;;;
n -
1.
Using the uniqueness of the orthogonal polynomials
[from Theorem
4.2], wn(x) must be a constant multiple of
<p
11
(x).
Thus
they must have the same zeros, and
Z;
=X;
i=1,
...
,n
To complete the proof of the uniqueness,
we
must show that
W;
=
u;,
where
u;
is
the weight in (5.3.21) and
w;
in (5.3.10). Use (5.3.23) with
(5.3.24) and
f(x)
=
h;(x).
The result will follow immediately, since
h
;(
x)
is now constructed using { x
;}.
(c) The error formula. We begin by deducing some further useful
properties about the weights {
w;}
in (5.3.10). Using the definition (3.6.2)
of
h;(x),
W;
=
jh
w(x)h;(x)
dx
=
jh
w(x)[1-
2/f(x;)(x-
X;)]
[l;(x)]
2
dx
a a
=
jhw(x)[l;(x)Y
dx-
2/f(x;)
jhw(x)(x-
x;)[l;(x)]
2
dx
a a
The last integral
is
zero from (5.3.19), since
h;(x)
=
(x-
x;)[/;(x)]
2
•
Thus
i =
1,2,
...
, n (5.3.25)
and all the weights are positive, for all n.
To
construct
W;,
begin by substituting
f(x)
=
l;(x)
into (5.3.20), and
note that
En(f)
=
0,
since
degree(/;)=
n - 1. Then using
l;(x)
= 8ij•
we
have
i = 1,
...
, n (5.3.26)
To further simplify the formula, the Christoffel-Darboux identity (Theo-
rem
4.6) can be used, followed by much manipulation, to give the
formula
(5.3.11). For the details, see Isaacson and Keller (1966, pp.
333-334).
For
the integration error, if
f(x)
is
2n
times continuously differentia-
ble on
[a,
b], then
En{!)=
Jh
w_(x)tfn(x)
dx
a
=
f[xl,
X1,
...
,
xn,
xn,
E]
jb
w(x)[
tltn(x )]
2
dx
some E E [a, b]
a
i
I
I
I
I
I
-.
------
~
--·
··--··
·····-·-·.
--
-·---
··-··
·--·-
...
n H -
•••
I
276
NUMERICAL INTEGRATION
·the last step using the integral mean value theorem. Using (5.3.14) in the
last integral, and replacing the divided difference by a derivative,
we
have
E
(f)=
J<2n){TI)
Jhw(x)
[cpn(x)}2
dx
n {2n)! a
A~
whlch gives the error formula in (5.3.10).
•
Gauss-Legendre
quadrature
For
w(x)
= 1,
the
Gaussian formula
on
[
-1,
1]
is
given
by
(5.3.27)
with the nodes equal to the zeros
of
the degree n Legendre polynomial Pn(x) on
[
-1,
1].
The weights are
and
-2
w.=
~--~~~--~~
'
(n
+
1)P;(x;)Pn+l(x;)
i = 1,2,
...
, n
2
2n+l(n!)4 /(2n)(
Tl)
/(2n)(
Tl)
E
(!)-
-e
--
n - (2n + 1)[(2n)!]
2
• (2n)! = n (2n)! '
Table
5.10 Gauss-Legendre nodes
and weights
n
X;
w.·
I
2 ±.5773502692 1.0
3
±.7745966692
.5555555556
0.0
.8888888889
4
± .8611363116
.3478546451
±.3399810436
.6521451549
5
±.9061798459
.2369268851
±.5384693101
.4786286705
0.0
.5688888889
6
±.9324695142
.1713244924
±.6612093865
.3607615730
±.2386191861
.4679139346
7 ±.9491079123
.1294849662
± .7415311856
.2797053915
±.4058451514
.3818300505
0.0
.4179591837
8 ±.9602898565
.1012285363
±.7966664774
.2223810345
±
.5255324099
.3137066459
±.1834346425
.3626837834
(5.3.28)
(5.3.29)
GAUSSIAN
QUADRATURE
277
Table
5.11
Gaussian quadrature
for (5.1.11)
n
I,
!-!,
2
-12.33621046570
2.66E-
1
3
-12.12742045017
5.71E-
2
4 -12.07018949029
-1.57E-
4
5
-12.07032853589
-1.78E-
5
6
-12.07034633110
1.47E-
8
7
-12.07034631753
1.14E-
9
8
-12.07034631639
-4.25E
-13
for some
-1
<
7J
<
1.
For integrals on other finite intervals with weight function
w(x)
=
1,
use the following linear change of variables:
l
b ( b - a )
JI
( a + b + x ( b -
a)
)
J(t)
dt
=
--
f dx
u 2
-1
2
(5.3.30)
reduc_ing
the integral to the standard interval [
-1,
1].
For
convenience,
we
include Table 5.10, which gives the nodes and weights for
formula
(5.3.27) with small values of n.
For
larger values of n, see the very
complete tables in Stroud and Secrest
(1966), which go up to n = 512.
Example Evaluate the integral (5.1.11),
I=
['ex
cos
(x)
dx =
-12.0703463164
0
which
was
used previously
in
Section
5.1
as an example for the trapezoidal rule
(see Table
5.1) and Simpson's rule (see Table 5.3).
The
results given in Table 5.11
show the marked superiority of Gaussian quadrature.
A general error result We
give
a useful result in trying to explain the excellent
convergence
of
Gaussian quadrature. In the next subsection,
we
consider in more
detail the error in Gauss-Legendre quadrature.
Theorem 5.4 Assume [a,
bJ
is
finite. Then the error in Gaussian quadrature,
b n
En(!)=
J
w(x)J(x)
dx-
L
wjf(x)
u
j=l
satisfies
n
~
1 (5.3.31)
with p
2
n_
1
{f)
the minimax error from (4.2.1).
i
I
_J
278
NUMERICAL
INTEGRATION
Proof
E,(p)
= 0 for any polynomial
p(x)
of
degree~
2n-
1.
Also. the error
function
En
satisfies
for all
F, G E
C[a,
b]. Let
p(x)
=
qin-
1
(x),
the minimax approxima-
tion of
degree~
2n-
1 to
f(x)
on
[a,
b]. Then
From
(5.3.25), all
wj
>
0.
Also, since
p(x)
= 1 is
of
degree
0,
n h
[wj=jw(x)dx
j=l
"
This completes the proof of (5.3.31).
•
From the results in Sections
4.6
and 4.7, the speed
of
convergence to zero of
PmU)
increases with the smoothness of the integrand. From (5.3.31), the same
is
true
of
Gaussian quadrature. In contrast, the composite trapezoidal rule
will
usually not converge faster than order h
2
[in particular, if
/'(b)
-
/'(a)
=f.
OJ,
regardless
of
the smoothness of
f(x).
Gaussian quadrature takes advantage
of
additional smoothness in the integrand, in contrast to most composite rules.
Example Consider using Gauss-Legendre quadrature to integrate
Table
5.12 contains error bounds based on (5.3.31),
Table 5.12 Gaussian quadrature
of
(5.3.32)
n
E,(f)
(5.3.33)
1
-3.20E-
2
1.06E-l
2
2.29E-
4
1.33E-
3
3
9.55E-
6
3.24E-
5
4
-3.35E-
7
9.24E-
7
5
6.05E-
9
1.61E-
8
(5.3.32)
(5.3.33)