¯

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

2π

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