Download Report

NUMERICAL INTEGRATION: ANOTHER APPROACH
We look for numerical integration formulas
Z
1
f (x) dx ≈
−1
n
X
wj f (xj )
j=1
which are to be exact for polynomials of as large a degree as
possible. There are no restrictions placed on the nodes {xj } nor
the weights {wj } in working towards that goal. The motivation is
that if it is exact for high degree polynomials, then it is expected
to be accurate when integrating functions that can be well
approximated by polynomials.
The case n = 1. We want a formula
Z 1
w1 f (x1 ) ≈
f (x) dx
−1
The weight w1 and the node x1 are to be so chosen that the
formula is exact for polynomials of as large a degree as possible.
To do this we substitute f (x) = 1 and f (x) = x. The first choice
leads to
Z
1
w1 · 1 =
1 dx,
w1 = 2
−1
The choice f (x) = x leads to
Z 1
w1 x1 =
x dx = 0,
x1 = 0
−1
The desired formula is
Z
1
f (x) dx ≈ 2f (0)
−1
It is called the midpoint rule and was introduced in the problems
of Section 5.1.
The case n = 2. We want a formula
Z
1
w1 f (x1 ) + w2 f (x2 ) ≈
f (x) dx
−1
The weights w1 , w2 and the nodes x1 , x2 are to be so chosen that
the formula is exact for polynomials of as large a degree as
possible. We substitute and force equality for
f (x) = 1, x, x 2 , x 3
This leads to the system
Z
1
w1 + w2 =
1 dx = 2
−1
Z 1
w1 x1 + w2 x2 =
w1 x12 + w2 x22 =
w1 x13 + w2 x23 =
x dx = 0
−1
Z 1
−1
1
Z
−1
x 2 dx =
2
3
x 3 dx = 0
The solution is given by
w1 = w2 = 1,
x1 =
−1
√
3
,
x2 =
√1
3
This yields the formula
Z 1
√ √ f (x) dx ≈ f −1/ 3 + f 1/ 3
(1)
−1
We say it has degree of precision equal to 3 since it integrates
exactly all polynomials of degree ≤ 3. We can verify directly that
it does not integrate exactly f (x) = x 4 .
Z
1
2
x 4 dx =
5
−1
√ √ f −/ 3 + f 1/ 3 = 2/9
Thus (1) has degree of precision exactly 3.
EXAMPLE
Integrate
Z
1
−1
dx
.
= log 2 = 0.69314718
3+x
The formula (1) yields
1
1
+
= 0.69230769
3 + x1 3 + x2
Error = .000839
THE GENERAL CASE
We want to find the weights {wi } and nodes {xi } so as to have
Z 1
n
X
f (x) dx ≈
wj f (xj )
−1
j=1
be exact for a polynomials f (x) of as large a degree as possible. As
unknowns, there are n weights wi and n nodes xi . Thus it makes
sense to initially impose 2n conditions so as to obtain 2n equations
for the 2n unknowns. We require the quadrature formula to be
exact for f (x) = x i , i = 0, 1, ..., 2n − 1. Then we obtain the
system of equations
Z 1
i
i
i
w1 x1 + w2 x2 + · · · + wn xn =
x i dx
−1
for i = 0, 1, 2, ..., 2n − 1. For the right sides,

Z 1
 2 , i = 0, 2, ..., 2n − 2
i +1
x i dx =

−1
0,
i = 1, 3, ..., 2n − 1
The system of equations
w1 x1i
+ ··· +
wn xni
Z
1
=
x i dx,
i = 0, ..., 2n − 1
−1
has a solution, and the solution is unique except for re-ordering the
unknowns. The resulting numerical integration rule is called
Gaussian quadrature.
For n ≥ 3, it is difficult to solve the system directly. Fortunately,
the nodes and weights have other properties which enable them to
be found more easily by other methods. There are programs to
produce them; and most subroutine libraries have either a program
to produce them or tables of them for commonly used cases.
SYMMETRY OF FORMULA
The nodes and weights possess symmetry properties. In particular,
xi = −xn−i ,
wi = wn−i ,
i = 1, 2, ..., n
A table of these nodes and weights for n = 2, ..., 8 is given in the
text in Table 5.7.
In addition, it can be shown that all weights satisfy
wi > 0
for all n > 0. This is a very desirable property from a practical
point of view, since this property leads to less accumulation of
roundoff errors in implementing the numerical integration formula.
Moreover, it permits us to develop a useful error formula, which we
illustrate later.
CHANGE OF INTERVAL OF INTEGRATION
Integrals on other finite intervals [a, b] can be converted to
integrals over [−1, 1], as follows:
Z
Z b
b−a 1
b + a + t(b − a)
dt
F (x) dx =
F
2
2
a
−1
based on the change of integration variables
x=
b + a + t(b − a)
,
2
−1 ≤ t ≤ 1
EXAMPLE Over the interval [0, π], use
x = (1 + t) π2
Then
Z
π
F (x) dx =
0
π
2
Z
1
−1
F (1 + t) π2 dt
EXAMPLE
Consider again the integrals used as examples in Section 5.1:
Z 1
Z 4
Z 2π
dx
dx
−x 2
(1)
(2)
(3)
e dx,
I =
I =
,
I =
2
1
+
x
2
+
cos x
0
0
0
n
2
3
4
5
6
7
10
15
20
I − I (1)
2.29E − 4
9.55E − 6
−3.35E − 7
6.05E − 9
−7.77E − 11
8.60E − 13
∗
∗
∗
I − I (2)
−2.33E − 2
−3.49E − 2
−1.90E − 3
1.70E − 3
2.74E − 4
−6.45E − 5
1.27E − 6
7.40E − 10
∗
I − I (3)
8.23E − 1
−4.30E − 1
1.77E − 1
−8.12E − 2
3.55E − 2
−1.58E − 2
1.37E − 3
−2.33E − 5
3.96E − 7
Compare these results with those of Section 5.1.
AN ERROR FORMULA
For a Gaussian quadrature formula, the error
Z
1
f (x) dx −
En (f ) =
−1
n
X
wj f (xj )
j=1
is equal to
En (f ) = en
f (2n) (cn )
,
(2n)!
en =
22n+1 (n!)4
2
(2n + 1) [(2n)!]
≈
π
4n
for some a ≤ cn ≤ b.
To help in understanding the implications of this error formula,
introduce
(k) f (x)
Mk = max
−1≤x≤1
k!
With many integrands f (x), this sequence {Mk } is bounded or
even decreases to zero. For example,
f (x) = cos x
f (x) = (2 + x)
⇒
−1
Mk ≤ 1/k!
⇒
Mk ≤ 1
From our error formula,
|En (f )| ≤ en M2n ,
en ≈ π/4n
(2)
Under an assumption of uniform boundedness for {Mk }, we have
the error decreases by a factor of at least 4 with each increase of n
to n + 1. Compare this to the convergence of the trapezoidal and
Simpson rules for such functions, to help explain the very rapid
convergence of Gaussian quadrature.
A SECOND ERROR FORMULA
Let f (x) be continuous for a ≤ x ≤ b; let n ≥ 1. Then, for the
Gaussian numerical integration formula
Z
I ≡
b
f (x) dx ≈
a
n
X
wj f (xj ) ≡ In
j=1
on [a, b], the error in In satisfies
|I (f ) − In (f )| ≤ 2 (b − a) ρ2n−1 (f )
(3)
Here ρ2n−1 (f ) is the minimax error of degree 2n − 1 for f (x) on
[a, b]:
ρm (f ) = min
max |f (x) − p(x)| , m ≥ 0
deg(p)≤m a≤x≤b
EXAMPLE
2
Let f (x) = e −x . Then the minimax errors ρm (f ) are given in the
following table.
m
1
2
3
4
5
ρm (f )
5.30E − 2
1.79E − 2
6.63E − 4
4.63E − 4
1.62E − 5
m
6
7
8
9
10
ρm (f )
7.82E − 6
4.62E − 7
9.64E − 8
8.05E − 9
9.16E − 10
Using this table, apply (3) to
Z
I =
1
2
e −x dx
0
For n = 3, (3) implies
2
.
|I − I3 | ≤ 2ρ5 e −x = 3.24 × 10−5
The actual error is 9.55E − 6.
INTEGRATING A NON-SMOOTH INTEGRAND
Z
Use Gaussian quadrature to evaluate I =
1√
0
n
2
4
8
16
32
64
I − In
−7.22E − 3
−1.16E − 3
−1.69E − 4
−2.30E − 5
−3.00E − 6
−3.84E − 7
x dx =
2
3
Ratio
6.2
6.9
7.4
7.6
7.8
The column labeled Ratio is defined by
I − I1n
2
I − In
c
It is consistent with I − In ≈ 3 , which can be proven theoretically.
n
c
In comparison for the trapezoidal and Simpson rules, I − In ≈ 1.5 .
n
WEIGHTED GAUSSIAN QUADRATURE
Consider needing to evaluate integrals such as
Z 1
Z 1
1
f (x) log x dx,
x 3 f (x) dx
0
0
How do we proceed? Consider numerical integration formulas
Z b
n
X
w (x)f (x) dx ≈
wj f (xj )
a
j=1
in which f (x) is considered a “nice” function (one with several
continuous derivatives). The function w (x) is allowed to be
singular, but must be integrable. We assume here that [a, b] is a
finite interval. The function w (x) is called a “weight function”,
and it is implicitly absorbed into the definition of the quadrature
weights {wi }. We again determine the nodes {xi } and weights
{wi } so as to make the integration formula exact for f (x) a
polynomial of as large a degree as possible.
The resulting numerical integration formula
Z
b
w (x)f (x) dx ≈
a
n
X
wj f (xj )
j=1
is called a Gaussian quadrature formula with weight function w (x).
We determine the nodes {xi } and weights {wi } by requiring
exactness in the above formula for
f (x) = x i ,
i = 0, 1, 2, ..., 2n − 1
To make the derivation more understandable, consider the
particular case
Z 1
n
X
1
x 3 f (x) dx ≈
wj f (xj )
0
j=1
We follow the same pattern as used earlier.
The case n = 1. We want a formula
Z 1
1
w1 f (x1 ) ≈
x 3 f (x) dx
0
The weight w1 and the node x1 are to be so chosen that the
formula is exact for polynomials of as large a degree as possible.
Choosing f (x) = 1, we have
Z 1
1
w1 =
x 3 dx = 34
0
Choosing f (x) = x, we have
Z
w1 x1 =
1
1
x 3 x dx =
0
x1 =
Thus
Z
1
0
has degree of precision 1.
1
3
7
4
7
x 3 f (x) dx ≈ 43 f
4
7
The case n = 2. We want a formula
Z
1
w1 f (x1 ) + w2 f (x2 ) ≈
1
x 3 f (x) dx
0
The weights w1 , w2 and the nodes x1 , x2 are to be so chosen that
the formula is exact for polynomials of as large a degree as
possible. We determine them by requiring equality for
f (x) = 1, x, x 2 , x 3
This leads to the system
Z
1
1
x 3 dx =
w1 + w2 =
0
Z
1
w1 x1 + w2 x2 =
3
4
1
x x 3 dx =
0
w1 x12 + w2 x22 =
Z
w1 x13 + w2 x23 =
Z
1
1
3
10
1
3
13
x 2 x 3 dx =
0
0
1
3
7
x 3 x 3 dx =
The solution is
7
3 √
35,
−
13 65
3
3 √
w1 = −
35,
8 392
x1 =
7
3 √
35
+
13 65
3
3 √
w2 = +
35
8 392
x2 =
Numerically,
x1 = .2654117024, x2 = .8115113746
w1 = .3297238792, w2 = .4202761208
The formula
Z
1
1
x 3 f (x) dx ≈ w1 f (x1 ) + w2 f (x2 )
0
has degree of precision 3.
(4)
EXAMPLE
Consider evaluating the integral
Z 1
1
x 3 cos x dx
(5)
0
In applying (4), we take f (x) = cos x. Then
w1 f (x1 ) + w2 f (x2 ) = 0.6074977951
The true answer is
Z
1
1
.
x 3 cos x dx = 0.6076257393
0
.
and our numerical answer is in error by E2 = .000128. This is quite
a good answer involving very little computational effort (once the
formula has been determined). In contrast, the trapezoidal and
Simpson rules applied to (5) would converge very slowly because
the first derivative of the integrand is singular at the origin.
CHANGE OF INTEGRATION VARIABLE
As a side note to the preceding example, we observe that the
change of variables x = t 3 transforms the integral (5) to
Z 1
3
t 3 cos t 3 dt
0
and both the trapezoidal and Simpson rules will perform better
with this formula, although still not as good as our weighted
Gaussian quadrature.
A change of the integration variable can often improve the
performance of a standard method, usually by increasing the
differentiability of the integrand.
EXAMPLE Using x = t r for some r > 1, we have
Z 1
Z 1
g (x) log x dx = r
t r −1 g (t r ) log t dt
0
0
The new integrand is generally smoother than the original one.