<< . .

. 54
( : 95)



. . >>

where x and x are two points in (a, b). Then
¯
s s s
δj ¤ δj u(xj ) ¤ uM
um δj . (9.10)
j=0 j=0 j=0
9.2 Interpolatory Quadratures 375

Let σs = s δj u(xj ) and consider the continuous function U (x) = u(x) s δj .
j=0 j=0
Thanks to (9.10), U (¯) ¤ σs ¤ U (x). Applying the mean-value theorem, there
¯
x
exists a point · between a and b such that U (·) = σs , which is (9.9). A similar
3
proof can be carried out if the coe¬cients δj are negative.

The composite midpoint formula is implemented in Program 71. Through-
out this chapter, we shall denote by a and b the end points of the integration
interval and by m the number of quadrature subintervals. The variable fun
contains the expression of the function f , while the output variable int
contains the value of the approximate integral.

Program 71 - midpntc : Midpoint composite formula

function int = midpntc(a,b,m,fun)
h=(b-a)/m; x=[a+h/2:h:b]; dim = max(size(x)); y=eval(fun);
if size(y)==1, y=diag(ones(dim))*y; end; int=h*sum(y);




9.2.2 The Trapezoidal Formula
This formula is obtained by replacing f with Π1 f , its Lagrange interpolat-
ing polynomial of degree 1, relative to the nodes x0 = a and x1 = b (see
Figure 9.2, left). The resulting quadrature, having nodes x0 = a, x1 = b
and weights ±0 = ±1 = (b ’ a)/2, is

b’a
I1 (f ) = [f (a) + f (b)] . (9.11)
2

If f ∈ C 2 ([a, b]), the quadrature error is given by

h3
E1 (f ) = ’ f (ξ), h = b ’ a (9.12)
12

where ξ is a point within the integration interval.

f (x) f (x)




x
x
a+b
a = x0 b = x1 a = x0 = x1 b = x2
2

FIGURE 9.2. Trapezoidal formula (left) and Cavalieri-Simpson formula (right)
376 9. Numerical Integration

Indeed, from the expression of the interpolation error (8.7) one gets

b b
1
(f (x) ’ Π1 f (x))dx = ’ f (ξ(x))(x ’ a)(b ’ x)dx.
E1 (f ) =
2
a a


Since ω2 (x) = (x ’ a)(x ’ b) < 0 in (a, b), the mean-value theorem yields

b

E1 (f ) = (1/2)f (ξ) ω2 (x)dx = ’f (ξ)(b ’ a)3 /12,
a


for some ξ ∈ (a, b), which is (9.12). The trapezoidal quadrature therefore
has degree of exactness equal to 1, as is the case with the midpoint rule.
To obtain the composite trapezoidal formula, we proceed as in the case
where n = 0, by replacing f over [a, b] with its composite Lagrange polyno-
mial of degree 1 on m subintervals, with m ≥ 1. Introduce the quadrature
nodes xk = a + kH, for k = 0, . . . , m and H = (b ’ a)/m, getting

m’1
H
m ≥ 1.
I1,m (f ) = (f (xk ) + f (xk+1 )) , (9.13)
2
k=0


Each term in (9.13) is counted twice, except the ¬rst and the last one, so
that the formula can be written as

1 1
I1,m (f ) = H f (x0 ) + f (x1 ) + . . . + f (xm’1 ) + f (xm ) . (9.14)
2 2

As was done for (9.8), it can be shown that the quadrature error associated
with (9.14) is

b’a 2
E1,m (f ) = ’ H f (ξ),
12

provided that f ∈ C 2 ([a, b]), where ξ ∈ (a, b). The degree of exactness is
again equal to 1.

The composite trapezoidal rule is implemented in Program 72.

Program 72 - trapezc : Composite trapezoidal formula

function int = trapezc(a,b,m,fun)
h=(b-a)/m; x=[a:h:b]; dim = max(size(x)); y=eval(fun);
if size(y)==1, y=diag(ones(dim))*y; end;
int=h*(0.5*y(1)+sum(y(2:m))+0.5*y(m+1));
9.2 Interpolatory Quadratures 377

9.2.3 The Cavalieri-Simpson Formula
The Cavalieri-Simpson formula can be obtained by replacing f over [a, b]
with its interpolating polynomial of degree 2 at the nodes x0 = a, x1 =
(a + b)/2 and x2 = b (see Figure 9.2, right). The weights are given by
±0 = ±2 = (b ’ a)/6 and ±1 = 4(b ’ a)/6, and the resulting formula reads

b’a a+b
I2 (f ) = f (a) + 4f + f (b) . (9.15)
6 2
It can be shown that the quadrature error is

b’a
h5 (4)
E2 (f ) = ’ (9.16)
f (ξ), h =
90 2
provided that f ∈ C 4 ([a, b]), and where ξ lies within (a, b). From (9.16) it
turns out that (9.15) has degree of exactness equal to 3.
Replacing f with its composite polynomial of degree 2 over [a, b] yields
the composite formula corresponding to (9.15). Introducing the quadrature
nodes xk = a + kH/2, for k = 0, . . . , 2m and letting H = (b ’ a)/m, with
m ≥ 1 gives
m’1 m’1
H
I2,m = f (x0 ) + 2 f (x2r ) + 4 f (x2s+1 ) + f (x2m ) . (9.17)
6 r=1 s=0

The quadrature error associated with (9.17) is
b’a
E2,m (f ) = ’ (H/2)4 f (4) (ξ),
180
provided that f ∈ C 4 ([a, b]) and where ξ ∈ (a, b); the degree of exactness
of the formula is 3.
The composite Cavalieri-Simpson quadrature is implemented in Program
73.

Program 73 - simpsonc : Composite Cavalieri-Simpson formula
function int = simpsonc(a,b,m,fun)
h=(b-a)/m; x=[a:h/2:b]; dim = max(size(x)); y=eval(fun);
if size(y)==1, y=diag(ones(dim))*y; end;
int=(h/6)*(y(1)+2*sum(y(3:2:2*m-1))+4*sum(y(2:2:2*m))+y(2*m+1));


Example 9.1 Let us employ the midpoint, trapezoidal and Cavalieri-Simpson
composite formulae to compute the integral

3(e’2π ’ 1) ’ 10πe’2π
’x
’0.122122.
xe cos(2x)dx = (9.18)
25
0
378 9. Numerical Integration

Table 9.1 shows in even columns the behavior of the absolute value of the er-
ror when halving H (thus, doubling m), while in odd columns the ratio Rm =
|Em |/|E2m | between two consecutive errors is given. As predicted by the previous
theoretical analysis, Rm tends to 4 for the midpoint and trapezoidal rules and

to 16 for the Cavalieri-Simpson formula.


|E0,m | Rm |E1,m | Rm |E2,m | Rm
m
1 0.9751 1.589e-01 7.030e-01
2 1.037 0.9406 0.5670 0.2804 0.5021 1.400
3.139 · 10’3
4 0.1221 8.489 0.2348 2.415 159.96
2.980 · 10’2 5.635 · 10’2 1.085 · 10’3
8 4.097 4.167 2.892
6.748 · 10’3 1.327 · 10’2 7.381 · 10’5
16 4.417 4.245 14.704
1.639 · 10’3 3.263 · 10’3 4.682 · 10’6
32 4.118 4.068 15.765
4.066 · 10’4 8.123 · 10’4 2.936 · 10’7
64 4.030 4.017 15.946
1.014 · 10’4 2.028 · 10’4 1.836 · 10’8
128 4.008 4.004 15.987
2.535 · 10’5 5.070 · 10’5 1.148 · 10’9
256 4.002 4.001 15.997
TABLE 9.1. Absolute error for midpoint, trapezoidal and Cavalieri-Simpson com-
posite formulae in the approximate evaluation of integral (9.18)




9.3 Newton-Cotes Formulae
These formulae are based on Lagrange interpolation with equally spaced
nodes in [a, b]. For a ¬xed n ≥ 0, let us denote the quadrature nodes
by xk = x0 + kh, k = 0, . . . , n. The midpoint, trapezoidal and Simpson
formulae are special instances of the Newton-Cotes formulae, taking n = 0,
n = 1 and n = 2 respectively. In the general case, we de¬ne:
b’a
(n ≥ 1);
- closed formulae, those where x0 = a, xn = b and h =
n
b’a
(n ≥ 0).
- open formulae, those where x0 = a+h, xn = b’h and h =
n+2
A signi¬cant property of the Newton-Cotes formulae is that the quadra-
ture weights ±i depend explicitly only on n and h, but not on the integration
interval [a, b]. To check this property in the case of closed formulae, let us
introduce the change of variable x = Ψ(t) = x0 + th. Noting that Ψ(0) = a,
Ψ(n) = b and xk = a + kh, we get
x ’ xk a + th ’ (a + kh) t’k
= = .
xi ’ xk a + ih ’ (a + kh) i’k
Therefore, if n ≥ 1
n
t’k
0 ¤ i ¤ n.
li (x) = = •i (t),
i’k
k=0,k=i
9.3 Newton-Cotes Formulae 379

The following expression for the quadrature weights is obtained
b n n

±i = li (x)dx = •i (t)hdt = h •i (t)dt,
a 0 0

from which we get the formula
n
n
In (f ) = h wi f (xi ), wi = •i (t)dt.
i=0 0

Open formulae can be interpreted in a similar manner. Actually, using again
the mapping x = Ψ(t), we get x0 = a + h, xn = b ’ h and xk = a + h(k + 1)
for k = 1, . . . , n ’ 1. Letting, for sake of coherence, x’1 = a, xn+1 = b and
n+1
proceeding as in the case of closed formulae, we get ±i = h ’1 •i (t)dt,
and thus
n+1
n
In (f ) = h wi f (xi ), wi = •i (t)dt.
i=0 ’1

In the special case where n = 0, since l0 (x) = •0 (t) = 1, we get w0 = 2.
The coe¬cients wi do not depend on a, b, h and f , but only depend on n,
and can therefore be tabulated a priori. In the case of closed formulae, the
polynomials •i and •n’i , for i = 0, . . . , n ’ 1, have by symmetry the same
integral, so that also the corresponding weights wi and wn’i are equal for
i = 0, . . . , n ’ 1. In the case of open formulae, the weights wi and wn’i are
equal for i = 0, . . . , n. For this reason, we show in Table 9.2 only the ¬rst
half of the weights.
Notice the presence of negative weights in open formulae for n ≥ 2. This can
be a source of numerical instability, in particular due to rounding errors.

n 1 2 3 4 5 6 n 0 1 2 3 4 5
1 1 3 14 95 41 3 8 55 66 4277
w0 w0 2
2 3 8 45 288 140 2 3 24 20 1440
4 9 64 375 216
’4 ’ 84 ’ 3171
5
w1 0 w1 0 0
3 8 45 288 140 3 24 20 1440
24 250 27 156 3934
w2 0 0 0 w2 0 0 0 0
45 288 140 20 1440
272
w3 0 0 0 0 0 140

TABLE 9.2. Weights of closed (left) and open Newton-Cotes formulae (right)

Besides its degree of exactness, a quadrature formula can also be quali¬ed
by its order of in¬nitesimal with respect to the integration stepsize h, which
is de¬ned as the maximum integer p such that |I(f ) ’ In (f )| = O(hp ).
Regarding this, the following result holds
380 9. Numerical Integration

Theorem 9.2 For any Newton-Cotes formula corresponding to an even
value of n, the following error characterization holds

Mn
hn+3 f (n+2) (ξ),
En (f ) = (9.19)
(n + 2)!

provided f ∈ C n+2 ([a, b]), where ξ ∈ (a, b) and
± n



 t πn+1 (t)dt < 0 for closed formulae,



0
Mn =
 n+1



 t πn+1 (t)dt > 0 for open formulae,


’1

n
having de¬ned πn+1 (t) = i=0 (t ’ i). From (9.19), it turns out that the
degree of exactness is equal to n + 1 and the order of in¬nitesimal is n + 3.

Similarly, for odd values of n, the following error characterization holds

Kn
hn+2 f (n+1) (·),
En (f ) = (9.20)
(n + 1)!

provided f ∈ C n+1 ([a, b]), where · ∈ (a, b) and
± n



 πn+1 (t)dt < 0 for closed formulae,



0
Kn =
 n+1



 πn+1 (t)dt > 0 for open formulae.


’1


The degree of exactness is thus equal to n and the order of in¬nitesimal is
n + 2.

Proof. We give a proof in the particular case of closed formulae with n even,
referring to [IK66], pp. 308-314, for a complete demonstration of the theorem.
Thanks to (8.20), we have

b

En (f ) = I(f ) ’ In (f ) = f [x0 , . . . , xn , x]ωn+1 (x)dx. (9.21)
a

x
Set W (x) = a ωn+1 (t)dt. Clearly, W (a) = 0; moreover, ωn+1 (t) is an odd func-
tion with respect to the midpoint (a + b)/2 so that W (b) = 0. Integrating by
9.3 Newton-Cotes Formulae 381

parts (9.21) we get
b b
d
f [x0 , . . . , xn , x]W (x)dx = ’
En (f ) = f [x0 , . . . , xn , x]W (x)dx
dx
a a
b
(n+2)
f (ξ(x))

= W (x)dx.
(n + 2)!
a

In deriving the formula above we have used the following identity (see Exercise
4)
d
f [x0 , . . . , xn , x] = f [x0 , . . . , xn , x, x]. (9.22)
dx
Since W (x) > 0 for a < x < b (see [IK66], p. 309), using the mean-value theorem
we obtain
b bx
f (n+2) (ξ) f (n+2) (ξ)
En (f ) = ’ W (x)dx = ’ ωn+1 (t) dt dx (9.23)
(n + 2)! (n + 2)!
a aa

where ξ lies within (a, b). Exchanging the order of integration, letting s = x0 +„ h,
for 0 ¤ „ ¤ n, and recalling that a = x0 , b = xn , yields
b bb

(s ’ x0 ) . . . (s ’ xn )dxds
W (x)dx =
a as

<< . .

. 54
( : 95)



. . >>