Numerical Integration

Professor Ph.D. Jorge Bacca

Universidad Industrial de Santander

Colombia

2024

High Dimensional Signal Processing Group

[Link]
jbacquin@[Link]

Professor Jorge Bacca Numerical methods 2024 1 / 43

Outline

1 Introduction to Quadrature

2 Composite Trapezoidal and Simpson’s Rule

Professor Jorge Bacca Numerical methods 2024 2 / 43

Introduction to Quadrature
We now approach the subject of numerical integration. The goal is to
approximate the definite integral of f (x) over the interval [a, b] by
evaluating f (x) at a finite number of sample points.

Definition 1
Suppose that a = x0 < x1 < · · · < xM = b. A formula of the form.
M
X
Q[f ] = wk f (xk ) = w0 f (x0 ) + w1 f (x1 ) + · · · + wM f (xM ) (1)
k=0

with the property that Z b

f (x)dx = Q[f ] + E[f ] (2)
a

is called a numerical integration or quadrature formula. The term E[f ]

 M
is called the truncation error for integration. The values xk k=0 are
 M
called the quadrature nodes, and wk k=0 are called the weights.

Professor Jorge Bacca Numerical methods 2024 3 / 43

Introduction to Quadrature

Depending on the application, the nodes {xk } are chosen in various

ways.
For the trapezoidal rule, Simpson’s rule, and Boole’s rule, the
nodes are chosen to be equally spaced.
For Gauss-Legendre quadrature, the nodes are chosen to be
zeros of certain Legendre polynomials.

Professor Jorge Bacca Numerical methods 2024 4 / 43

Introduction to Quadrature

The derivation of quadrature formulas is sometimes based on polyno-

mial interpolation.
Recall that there exists a unique polynomial PM (x) of degree ≤ M pass-
ing through the M + 1 equally spaced points
{(xk , f (xk ))}M
k=0 . When this polynomial is used to approximate f (x) over
[a, b],and then the integral of f (x) is approximated by the integral of
PM (x), the resulting formula is called a Newton-Cotes quadrature for-
mula (see figure). When the sample points x0 = a and xM = b are used,
it is called a closed Newton-Cotes formula.

Professor Jorge Bacca Numerical methods 2024 5 / 43

Introduction to Quadrature
y y

1.5 1.5
y=f(x) y=f(x)
1.0 1.0

0.5 0.5

x x
0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0

(a) (b)
y y

1.5 1.5
y=f(x) y=f(x)
1.0 1.0

0.5 0.5

x x
0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0

(c) (d)

Figure 1: (a) The trapezoidal rule integrates y = P1 (x) over [x0 , x1 ] = [0.0, 0.5]. (b) Simpson’s
rule integrates y = P2 (x) over [x0 , x1 ] = [0.0, 1.0]. (c) Simpson’s 38 rule integrates y = P3 (x) over
[x0 , x3 ] = [0.0, 1.5]. (d) Boole’s rule integrates y = P4 (x) over [x0 , x4 ] = [0.0, 2.0].
Professor Jorge Bacca Numerical methods 2024 6 / 43
Introduction to Quadrature

Theorem 1 (Closed Newton-Cotes Quadrature Formula).

Assume that xk = x0 + kh are equally spaced nodes and fk = f (xk ). The
first four closed Newton-Cotes quadrature formulas are
Z x1
h
f (x) dx ≈ ( f0 + f1 ) (Trapezoidal rule), (3)
x 2
Z 0x2
h
f (x) dx ≈ ( f0 + 4f1 + f2 ) (Simpson’s rule), (4)
x 3
Z 0x3
3h 3
f (x) dx ≈ ( f0 + 3f1 + 3f2 + f3 ) (Simpson’s rule), (5)
x 8 8
Z 0x4
2h
f (x) dx ≈ (7f0 + 32f1 + 12f2 + 32f3 + 7f4 ) (Boole’s rule) (6)
x0 45

Professor Jorge Bacca Numerical methods 2024 7 / 43

Introduction to Quadrature

Corollary 1 (Newton-Cotes Precision).

Assume that f (x) is sufficiently differentiable; then E[ f ] for
Newton-Cotes quadrature involves an appropriate higher derivative.
The trapezoidal rule has degree of precision n = 1. If f ∈ C2 [a, b], then
x1
h3
Z
h
f (x) dx = (f0 + f1 ) − f (2) (c). (7)
x0 2 12

Simpson’s rule has degree of precision n = 3. If f ∈ C4 [a, b], then

x2
h5
Z
h
f (x) dx = (f0 + 4f1 + f2 ) − f (4) (c). (8)
x0 3 90

Professor Jorge Bacca Numerical methods 2024 8 / 43

Introduction to Quadrature

Corollary 1 (Newton-Cotes Precision).

3
Simpson’s 8 rule has degree of precision n = 3. If f ∈ C4 [a, b], then
x3
3h5 (4)
Z
3h
f (x) dx = (f0 + 3f1 + 3f2 + f3 ) − f (c). (9)
x0 8 80

Boole’s rule has degree of precision n = 5. If f ∈ C6 [a, b], then

x4
8h7 (6)
Z
2h
f (x) dx = (7f0 + 32f1 + 12f2 + 32f3 + 7f4 ) − f (c). (10)
x0 45 945

Professor Jorge Bacca Numerical methods 2024 9 / 43

Introduction to Quadrature
Proof of Theorem 1.
Start with the Lagrange polynomial PM (x) based on x0 , x1 , ..., xM that
can be used to approximate f (x) :
M
X
f (x) ≈ PM (x) = fk LM,k (x), (11)
k=0

where fk = f (xk ) for k = 0, 1, ..., M. An approximation for the integral is

obtained by replacing the integrand f (x) with the polynomial PM (x).
This is the general method for obtaining a Newton-Cotes integration
formula:
Z xM Z xM
f (x) dx ≈ PM (x) dx
x0 x0
M
! M Z
Z xM X X xM 
= fk LM,k (x) dx = fk LM,k (x) dx (12)
x0 k=0 k=0 x0

M Z
X xM  M
X
= LM,k (x) dx fk = wk fk .
k=0 x0 k=0

Professor Jorge Bacca Numerical methods 2024 10 / 43

Introduction to Quadrature

The details for the general computations of the coefficients of wk in (13)

are tedious. We shall give a sample proof of Simpson’s rule, which is
the case M = 2. This case involves the approximating polynomial
(x − x1 )(x − x2 ) (x − x0 )(x − x2 ) (x − x0 )(x − x1 )
P2 (x) = f0 + f1 + f2 . (13)
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )

Since f0 , f1 , and f2 are constants with respect to integration, the

relations in (13) lead to
Z x2 Z x2 Z x2
(x − x1 )(x − x2 ) (x − x0 )(x − x2 )
f (x) dx ≈ f0 dx + f1 dx
x0 x0 (x0 − x1 )(x0 − x2 ) x0 (x1 − x0 )(x1 − x2 )
Z x2
(x − x0 )(x − x1 )
+f2 dx. (14)
x0 (x2 − x0 )(x2 − x1 )

Professor Jorge Bacca Numerical methods 2024 11 / 43

Introduction to Quadrature

We introduce the change of variable x = x0 + ht with dx = h dt to assist

with the evaluation of the integrals in (15). The new limits of integration
are from t = 0 to t = 2. The equal spacing of the nodes xk = x0 + kh
leads to xk − xj = (k − j)h and x − xk = h(t − k), which are used to
simplify (15) and get

Professor Jorge Bacca Numerical methods 2024 12 / 43

Introduction to Quadrature

Z x2 Z 2 Z 2
h(t − 1)h(t − 2) h(t − 0)h(t − 2)
f (x) dx ≈ f0 h dt + f1 h dt (15)
x0 0 (−h)(−2h) 0 (h)(−h)
Z 2
h(t − 0)h(t − 1)
+f2 h dt
0 (2h)(h)
Z 2 Z 2 Z 2
h 2 h 2
= f0 (t − 3t + 2) dt − f1 h (t − 2t)dt + f2 (t2 − t) dt
2 0 0 2 0
 3 2
 t=2  3  t=2
h t 3t t
= f0 − + 2t − f1 h − t2
2 3 2 t=0 3 t=0
 t=2
h t3 t2

+f2 −
2 3 2 t=0
     
h 2 −4 h 2
= f0 − f1 h + f2
2 3 3 2 3
h
= (f0 + 4f1 + f2 )
3
and the proof is complete.
Professor Jorge Bacca Numerical methods 2024 13 / 43
Introduction to Quadrature
Example 1 Consider the function f (x) = 1 + e−x sin(4x), the equally
spaced quadrature nodes x0 = 0.0, x1 = 0.5, x2 = 1.0, x3 = 1.5, and
x4 = 2.0, and the corresponding function values
f0 = 1.00000, f1 = 1.55152, f2 = 0.72159, f3 = 0.93765, and f4 = 1.13390.
Apply the various quadrature formulas (4) through (7).

Professor Jorge Bacca Numerical methods 2024 14 / 43

Introduction to Quadrature
Example 1 Consider the function f (x) = 1 + e−x sin(4x), the equally
spaced quadrature nodes x0 = 0.0, x1 = 0.5, x2 = 1.0, x3 = 1.5, and
x4 = 2.0, and the corresponding function values
f0 = 1.00000, f1 = 1.55152, f2 = 0.72159, f3 = 0.93765, and f4 = 1.13390.
Apply the various quadrature formulas (4) through (7).
The step size is h = 0.5, and the computations are
Z 0.5
0.5
f (x) dx ≈ (1.00000 + 1.55152) = 0.63788
0 2
Z 1.0
0.5
f (x) dx ≈ (1.00000 + 4(1.55152) + 0.72159) = 1.32128
0 3
Z 1.5
3(0.5)
f (x) dx ≈ (1.00000 + 3(1.55152) + 3(0.72159) + 0.93765)
0 8
= 1.64193
Z 2.0
2(0.5)
f (x) dx ≈ (7(1.00000) + 32(1.55152) + 12(0.72159)
0 45
+32(0.93765) + 7(1.13390)) = 2.29444.

Professor Jorge Bacca Numerical methods 2024 14 / 43

Introduction to Quadrature

It is important to realize that the quadrature formulas (4) through (7)

applied in the illustration above give approximations for definite integrals
over different intervals. The graph of the curve y = f (x) and the areas
under the Lagrange polynomials y = P1 (x), y = P2 (x), y = P3 (x), and
y = P4 (x) are shown in Figure 1(a) through (d), respectively.
In Example 1 we applied the quadrature rules with h = 0.5. If the end-
points of the interval [a, b] are held fixed, the step size must be adjusted
for each rule. The step sizes are h = b − a, h = (b − a)/2, h = (b − a)/3,
and h = (b − a)/4 for the trapezoidal rule, Simpson’s rule, Simpson’s
3
8 rule, and Boole’s rule, respectively. The next example illustrates this
point.

Professor Jorge Bacca Numerical methods 2024 15 / 43

Introduction to Quadrature

Example 2 Consider the integration of the function

f (x) = 1 + e−x sin(4x) over the fixed interval [a, b] = [0, 1]. Apply the
various formulas (4) through (7).

Professor Jorge Bacca Numerical methods 2024 16 / 43

Introduction to Quadrature

Example 2 Consider the integration of the function

f (x) = 1 + e−x sin(4x) over the fixed interval [a, b] = [0, 1]. Apply the
various formulas (4) through (7).
For the trapezoidal rule, h = 1 and
Z 1
1
f (x) dx ≈ (f (0) + f (1))
0 2
1
= (1.00000 + 0.72159) = 0.86079.
2

Professor Jorge Bacca Numerical methods 2024 16 / 43

Introduction to Quadrature

Example 2 Consider the integration of the function

f (x) = 1 + e−x sin(4x) over the fixed interval [a, b] = [0, 1]. Apply the
various formulas (4) through (7).
For the trapezoidal rule, h = 1 and
Z 1
1
f (x) dx ≈ (f (0) + f (1))
0 2
1
= (1.00000 + 0.72159) = 0.86079.
2
For Simpson’s rule, h = 1/2, and we get
Z 1
1/2 1
f (x) dx ≈ (f (0) + 4f ( ) + f (1))
0 3 2
1
= (1.00000 + 4(1.55152) + 0.72159) = 1.32128.
6

Professor Jorge Bacca Numerical methods 2024 16 / 43

Introduction to Quadrature

3
For Simpson’s 8 rule, h = 1/3, and we obtain
Z 1
3(1/3) 1 2
f (x) dx ≈ (f (0) + 3f ( ) + 3f ( ) + f (1))
0 8 3 3
1
= (1.00000 + 3(1.69642) + 3(1.23447) + 0.72159) = 1.31440
8

Professor Jorge Bacca Numerical methods 2024 17 / 43

Introduction to Quadrature

3
For Simpson’s 8 rule, h = 1/3, and we obtain
Z 1
3(1/3) 1 2
f (x) dx ≈ (f (0) + 3f ( ) + 3f ( ) + f (1))
0 8 3 3
1
= (1.00000 + 3(1.69642) + 3(1.23447) + 0.72159) = 1.31440
8
For Boole’s rule, h = 1/4, and the result is
Z 1
2(1/4) 1 1 3
f (x) dx ≈ (7f (0) + 32f ( ) + 12f ( ) + 32f ( ) + 7f (1))
0 45 4 2 4
1
= (7(1.00000) + 32(1.65534) + 12(1.55152)
90
+32(1.06666) + 7(0.72159)) = 1.30859.

Professor Jorge Bacca Numerical methods 2024 17 / 43

Introduction to Quadrature

The true value of the definite integral is

Z 1
21e − 4 cos(4) − sin(4)
f (x) dx = = 1.3082506046426...,
0 17e

and the approximation 1.30859 from Boole’s rule is best. The area
under each of the Lagrange polynomials P1 (x), P2 (x), P3 (x), and P4 (x)
is shown in Figure 2(a) through (d), respectively.

Professor Jorge Bacca Numerical methods 2024 18 / 43

Introduction to Quadrature
y y
y=f(x) y=f(x)
1.5 1.5

1.0 1.0

0.5 0.5

x x
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
y y
y=f(x) y=f(x)
1.5 1.5

1.0 1.0

0.5 0.5

x x
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

Figure 2: (a) The trapezoidal rule used over [0, 1] yields the approximation 0.86079. (b)
Simpson’s rule used over [0, 1] yields the approximation 1.32128. (c) Simpson’s 83 rule used over
[0, 1] yields the approximation 1.31440. (d) Boole’s rule used over [0, 1] yields the approximation
1.30859.

Professor Jorge Bacca Numerical methods 2024 19 / 43

Introduction to Quadrature

To make a fair comparison of quadrature methods, we must use the

same number of function evaluations in each method. Our final example
is concerned with comparing integration over a fixed interval [a, b] using
exactly five function evaluations fk = f (xk ),for k = 0, 1, ..., 4 for each
method. when the trapezoidal rule is applied on the four subintervals
[x0 , x1 ], [x1 , x2 ], [x2 , x3 ], and [x3 , x4 ], it is called a composite trapezoidal
rule:
Z x4 Z x1 Z x2 Z x3 Z x4
f (x) dx = f (x) dx + f (x) dx + f (x) dx + f (x) dx
x0 x0 x1 x2 x3
h h h h
≈ (f0 + f1 ) + (f1 + f2 ) + (f2 + f3 ) + (f3 + f4 ) (16)
2 2 2 2
h
= (f0 + 2f1 + 2f2 + 2f3 + f4 ).
2

Professor Jorge Bacca Numerical methods 2024 20 / 43

Introduction to Quadrature

Simpson’s rule can also be used in this manner. When Simpson’s rule
is applied on the two subintervals [x0 , x2 ] and [x2 , x4 ], it is called a
composite Simpson’s rule:
Z x4 Z x2 Z x4
f (x) dx = f (x) dx + f (x) dx
x0 x0 x2
h h
≈ (f0 + 4f1 + f2 ) + (f2 + 4f3 + f4 ) (17)
3 3
h
= (f0 + 4f1 + 2f2 + 4f3 + f4 ).
3
The next example compares the values obtained with (17), (18), and
(7).

Professor Jorge Bacca Numerical methods 2024 21 / 43

Introduction to Quadrature
Example 3 Consider the integration of the function
f (x) = 1 + e−x sin(4x) over [a, b] = [0, 1]. Use exactly five function
evaluations and compare the results from the composite trapezoidal
rule, composite Simpson rule, and Boole’s rule.

Professor Jorge Bacca Numerical methods 2024 22 / 43

Introduction to Quadrature
Example 3 Consider the integration of the function
f (x) = 1 + e−x sin(4x) over [a, b] = [0, 1]. Use exactly five function
evaluations and compare the results from the composite trapezoidal
rule, composite Simpson rule, and Boole’s rule.
The uniform step size is h = 1/4. The composite trapezoidal rule (17)
produces
Z 1
1/4 1 1 3
f (x) dx ≈ (f (0) + 2f ( ) + 2f ( ) + 2f ( ) + f (1))
0 2 4 2 4
1
= (1.00000 + 2(1.65534) + 2(1.55152) + 2(1.06666) + 0.72159)
8
= 1.28358.

Professor Jorge Bacca Numerical methods 2024 22 / 43

Introduction to Quadrature
Example 3 Consider the integration of the function
f (x) = 1 + e−x sin(4x) over [a, b] = [0, 1]. Use exactly five function
evaluations and compare the results from the composite trapezoidal
rule, composite Simpson rule, and Boole’s rule.
The uniform step size is h = 1/4. The composite trapezoidal rule (17)
produces
Z 1
1/4 1 1 3
f (x) dx ≈ (f (0) + 2f ( ) + 2f ( ) + 2f ( ) + f (1))
0 2 4 2 4
1
= (1.00000 + 2(1.65534) + 2(1.55152) + 2(1.06666) + 0.72159)
8
= 1.28358.
Using the composite Simpson’s rule (18), we get
Z 1
1/4 1 1 3
f (x) dx ≈ (f (0) + 4f ( ) + 2f ( ) + 4f ( ) + f (1))
0 3 4 2 4
1
= (1.00000 + 4(1.65534) + 2(1.55152) + 4(1.06666) + 0.72159)
12
= 1.30938.
Professor Jorge Bacca Numerical methods 2024 22 / 43
Introduction to Quadrature

We have already seen the result of Boole’s rule in Example 2:

Z 1
2(1/4) 1 1 3
f (x) dx ≈ (7f (0) + 32f ( ) + 12f ( ) + 32f ( ) + 7f (1))
0 45 4 2 4
= 1.30859.

Professor Jorge Bacca Numerical methods 2024 23 / 43

Introduction to Quadrature

We have already seen the result of Boole’s rule in Example 2:

Z 1
2(1/4) 1 1 3
f (x) dx ≈ (7f (0) + 32f ( ) + 12f ( ) + 32f ( ) + 7f (1))
0 45 4 2 4
= 1.30859.

The true value of the integral is

Z 1
21e − 4 cos(4) − sin(4)
f (x) dx = = 1.3082506046426...,
0 17e

and the approximation 1.30938 from Simpson’s rule is much better than
the value 1.28358 obtained from the trapezoidal rule. Again, the
approximation 1.30859 from Boole’s rule is closest. Graphs for the
areas under the trapezoids and parabolas are shown in Figure 3(a)
and (b), respectively.

Professor Jorge Bacca Numerical methods 2024 23 / 43

Introduction to Quadrature

y y

1.5 y=f(x) 1.5 y=f(x)

1.0 1.0

0.5 0.5

t t
0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00
(a) (b)

Figure 3: (a) The composite trapezoidal rule yields the approximation 1.28358. (b) The
composite Simpson rule yields the approximation 1.30938.

Professor Jorge Bacca Numerical methods 2024 24 / 43

Introduction to Quadrature
3
Example 4 Determine the degree of precision of Simpson’s 8 rule.
It will suffice to apply Simpson’s 38 rule over the interval [0, 3] with the
five test functions f (x) = 1, x, x2 , x3 ,and x4 . For the first four functions,
Simpson’s 38 rule is exact.
Z 3
3
1 dx = 3= (1 + 3(1) + 3(1) + 1)
0 8
Z 3
9 3
x dx = = (0 + 3(1) + 3(2) + 3)
0 2 8
Z 3
3
x2 dx = 9= (0 + 3(1) + 3(4) + 9)
0 8
Z 3
81 3
x3 dx = = (0 + 3(1) + 3(8) + 27).
0 4 8
The function f (x) = x4 is the lowest power of x for which the rule is not
exact. Z 3
243 99 3
x4 dx = ≈ = (0 + 3(1) + 3(16) + 81).
0 5 2 8
Therefore, the degree of precision of Simpson’s 3/8 rule is n = 3.
Professor Jorge Bacca Numerical methods 2024 25 / 43
Composite Trapezoidal and Simpson’s Rule
An intuitive method of finding the area under the curve y = f (x) over
[a, b] is by approximating that area with a series of trapezoids that lie
above the intervals {[xk , xk+1 ]}.
Theorem 2 (Composite trapezoidal rule)
Suppose that the interval [a, b] is subdivided into M intervals [xk , xk+1 ]
of width h = (b − a)/M by using the equally spaced nodes xk = a + kh,
for k = 0, 1, · · · , M. The composite trapezoidal rule for M
subintervals can be expressed in any of three equivalent ways
M
hX
T(f , h) = (f (xk−1 ) + f (xk )) (1a)
2
k=1
h
T(f , h) = (f0 + 2f1 + 2f2 + 2f3 + · · · + 2fM−2 + 2fM−1 + fM ) (1b)
2
M−1
h X
T(f , h) = (f (a) + f (b)) + h (f (xk )) (1c)
2
k=1
Professor Jorge Bacca Numerical methods 2024 26 / 43
Composite Trapezoidal and Simpson’s Rule

This is an approximation to the integral of f (x) over [a, b] and we write

Z b
f (x)dx ≈ T(f , h) (2)
a

Proof. Apply the trapezoidal rule over each subinterval [xk−1 , xk ]. Use
the additive property of integral for subintervals:
Z b M Z xk M
X hX
f (x)dx = f (x)dx ≈ (f (xk−1 ) + f (xk )) (3)
a xk−1 2
k=1 k=1

Since h/2 is a constant, the distributive law of addition can be applied

to obtain (1a). Formula (1b) is the expanded version of (1a). Formula
(1c) shows how to group all the intermediate terms in (1b) that are
multiplied by 2.

Professor Jorge Bacca Numerical methods 2024 27 / 43

Composite Trapezoidal and Simpson’s Rule
√
Example 5 Consider f (x) = 2 + sin(2 x). Use the composite
trapezoidal rule with 11 sample points to compute an approximation to
the integral of f (x) taken over [1, 6]. To generate 11 sample points we
use M = 10 and h = (6 − 1)/10 = 1/2. Using formula (1c), the
computation is

Professor Jorge Bacca Numerical methods 2024 28 / 43

Composite Trapezoidal and Simpson’s Rule

√
Figure 4: Approximating the area under the curve y = 2 + sin(2 x) with the composite
trapezoidal rule.

Professor Jorge Bacca Numerical methods 2024 29 / 43

Composite Trapezoidal and Simpson’s Rule
Theorem 3 (Composite Simpson Rule)
Suppose that [a, b] is subdivided into 2M subintervals [xk , xk+1 ] of equal
width h = (b − a)/2M by using xk = a + kh for k = 0, 1, · · · , 2M. The
composite Simpson rule for 2M subintervals can be expressed in
any of three equivalent ways:
M
hX
S(f , h) = (f (x2k−2 ) + 4f (x2k−1 ) + f (x2k )) (4a)
3
k=1
h
S(f , h) = (f0 + 4f1 + 2f2 + 4f3 + · · · + 2fM−2 + 4fM−1 + fM ) (4b)
3
M−1 M
h 2h X 4h X
S(f , h) = (f (a) + f (b)) + f (x2k ) + f (x2k−1 ) (4c)
3 3 3
k=1 k=1

This is an approximation to the integral of f (x) over [a, b], and we write
Z b
f (x)dx ≈ S(f , h) (5)
a
Professor Jorge Bacca Numerical methods 2024 30 / 43
Composite Trapezoidal and Simpson’s Rule

Proof. Apply Simpson’s rule over each subinterval [xk−1 , xk ]. Use the
additive property of the integral for subintervals:
Z b M Z 2xk M
X hX
f (x)dx = f (x)dx ≈ (f (x2k−2 ) + 4f (x2k−1 ) + f (x2k )) (6)
a x2k−2 3
k=1 k=1

Since h/3 is a constant, the distributive law of addition can be applied

to obtain (4a). Formula (4b) is the expanded version of (4a). Formula
(4c) shows how to group all the intermediate terms in (4b) that are
multiplied by 2 and those that are multiplied by 4.

Professor Jorge Bacca Numerical methods 2024 31 / 43

Composite Trapezoidal and Simpson’s Rule
√
Example 6 Consider f (x) = 2 + sin(2 x). Use the composite Simpson
rule with 11 sample points to compute an approximation to the integral
of f (x) taken over [1, 6]. To generate 11 sample points we must use
M = 5 and h = (6 − 1)/10 = 1/2. Using formula (4c), the computation is

Professor Jorge Bacca Numerical methods 2024 32 / 43

Composite Trapezoidal and Simpson’s Rule

√
Figure 5: Approximating the area under the curve y = 2 + sin(2 x) with the composite Simpson
rule.

Professor Jorge Bacca Numerical methods 2024 33 / 43

Composite Trapezoidal and Simpson’s Rule

Error Analysis

The significance of the next two results is to understand that the er-
ror terms ET (f , h) and ES (f , h) for the composite trapezoidal rule and
composite Simpson rule are of the order O(h2 ) and O(h4 ), respectively.
This shows that the error for Simpson’s rule converges to zero faster
than the error for the trapezoidal rule as the step size h decreases to
zero. In cases where the derivatives of f (x) are known, the formulas

−(b − a)f (2) (c)h2 −(b − a)f (4) (c)h4

ET (f , h) = and ES (f , h) =
12 180
Can be used to eliminate the number of subintervals required to achieve
a specified accuracy.

Professor Jorge Bacca Numerical methods 2024 34 / 43

Composite Trapezoidal and Simpson’s Rule
Corollary 2 (Trapezoidal Rule: Error analysis)
Suppose that [a, b] is subdivided into M subintervals [xk , xk+1 ] of width
h = (b − a)/M. The composite trapezoidal rule
M−1
h X
T(f , h) = (f (a) + f (b)) + h f (xk ) (7)
2
h=1

is an approximation of the integral

Z b
f (x)dx = T(f , h) + ET (f , h). (8)
a

Furthermore, if f ∈ C2 [a, b], there exits a value c with a < c < b so that
the error term ET (f , h) has the form

−(b − a)f (2) (c)h2

ET (f , h) = = O(h2 ). (9)
12
Professor Jorge Bacca Numerical methods 2024 35 / 43
Composite Trapezoidal and Simpson’s Rule
Corollary 3 (Simpson’s Rule: Error Analysis)
Suppose that [a, b] is subdivided into 2M subintervals [xk, xk + 1] of
equal width h = (b − a)/2M. The composite Simpson rule
M−1 M
h 2h X 4h X
S(f , h) = (f (a) + f (b)) + f (x2k ) + f (x2k−1 ) (14)
3 3 3
k=1 k=1

is an approximation to the integral

Z b
f (x)dx = S(f , h) + ES (f , h). (15)
a

Furthermore, if f ∈ C4 [a, b], there exists a value c with a < c < b so that
the error term ES (f , h) has the form

−(b − a)f (4) (c)h4

ES (f , h) = = O(h4 ). (18)
180
Professor Jorge Bacca Numerical methods 2024 36 / 43
Composite Trapezoidal and Simpson’s Rule
√
Example 7 Consider f (x) = 2 + sin(2 x). Investigate the error when
the composite trapezoidal rule is over [1, 6] and the number of
subintervals is 10, 20, 40, 80 and 160.

√
Table 1: Composite trapezoidal rule for f (x) = 2 + sin(2 x) over [1, 6]

Table 1 shows the approximations T(f , h). The antiderivate of f (x) is

√
√ √ sin(2 x)
F(x) = 2x − x cos (2 x) + ,
2
and the true value of definite integral is
Z 6
f (x)dx = F(x)|x=6
x=1 = 8.1834792077.
1

Professor Jorge Bacca Numerical methods 2024 37 / 43

Composite Trapezoidal and Simpson’s Rule

This value was used to compute the values ET (f , h) = 8.1834792077 −

T(f , h) in Table 1. It is important to observe that when h is reduced by
a factor of 1/2 the successive errors ET (f , h) are diminished by approx-
imately 1/4. This confirms that the order is of O(h2 ).

Example 8 Find the number M and the step size h so that the error
ET (f , h) for the Rcomposite trapezoidal rule is less than 5 × 10−9 for the
7
approximation 2 dx/x ≈ T(f , h). The integral is f = 1/x and the first two
derivates are f ′ (x) = −1/x2 and f (2) (x) = 2/x3 . The maximum value of
|f (2) (x)| taken over [2, 7] occurs at the end point x = 2, and those we
have the bound |f (2) (c) ≤ |f (2) (2)| = 1/4, for 2 ≤ c ≤ 7. This is used
with formula (9) to obtain

| − (b − a)f (2) (c)h2 | (7 − 2) 14 h2 5h2

|ET (f , h)| = ≤ = . (17)
12 12 48

Professor Jorge Bacca Numerical methods 2024 38 / 43

Composite Trapezoidal and Simpson’s Rule
The step size h and number M satisfy the relation h = 5/M, and this is
used in (17) to get the relation
125
≤ 5 × 10−9
|ET (f , h)| ≤ (18)
48M 2
Now rewrite (18) so that is easier to solve for M:
25
× 109 ≤ M 2 . (19)
48
Solving (19), we find that 22821.77 ≤ M. Since M must be an inte-
ger, we choose M = 22, 822 and the corresponding step size is h =
5/22, 822 = 0.000219086846. When the composite trapezoidal rule is im-
plemented with this main function evaluations, there is a possibility that
the rounded-off function evaluation will produce a significant amount of
error. When the computation was perform, the result was
 
5
T f, = 1.252762969,
22, 822
Professor Jorge Bacca Numerical methods 2024 39 / 43
Composite Trapezoidal and Simpson’s Rule

R7
which compares favorably with the true value 2 dx/x = ln(x)|x=7 x=2 =
1.252762968. The error is smaller than predicted because the bound
1/4 for |f (2) (c)| was used. Experimentation shows that it takes about
10.001 function evaluations to achieve the desired accuracy of 5 × 10−9 ,
and when the calculation is performed with M = 10, 000, the result is
 
5
T f, = 1.252762973.
10, 000

Professor Jorge Bacca Numerical methods 2024 40 / 43

Composite Trapezoidal and Simpson’s Rule

The composite trapezoidal rule usually requires a large number of func-

tion evaluations to achieve and accurate answer. This is contrasted in
the next example with Simpson’s rule, which will require significantly
fewer evaluations

Example 9 Find the number M and the step size h so that the error
ES (f , h) for theRcomposite Simpson rule is less than 5 × 10−9 for the
7
approximation 2 dx/x ≈ S(f , h). The integral is f (x) = 1/x and f (4) (x) =
24/x5 . The maximum value of |f (4) (x)| taken over [2, 7] occurs at the end
point x = 2, and those we have the bound |f (4) (c)| ≤ |f (4) (2)| = 3/4, for
2 ≤ c ≤ 7. This is used with formula (16) to obtain

| − (b − a)f (4) (c)h4 | (7 − 2) 34 h4 h4

|ES (f , h)| = ≤ = . (20)
180 180 48

Professor Jorge Bacca Numerical methods 2024 41 / 43

Composite Trapezoidal and Simpson’s Rule
The step size h and number M satisfy the relation h = 5/2M, and this is
used in (20) to get the relation
625
|ES (f , h)| ≤ ≤ 5 × 10−9 (21)
768M 4
Now rewrite (21) so that is easier to solve for M:
125
× 109 ≤ M 4 . (22)
768
Solving (22), we find that 112.95 ≤ M. Since M must be an integer,
we choose M = 113, and the corresponding step size is h = 5/226 =
0.02212389381. When the composite Simpson rule was perform, the
result was  
5
S f, = 1.252762969,
226
R7
which agrees with 2 dx/x = ln(x)|72 = 1.252762968.
Professor Jorge Bacca Numerical methods 2024 42 / 43
Composite Trapezoidal and Simpson’s Rule

Experimentation shows that it takes about 129 function evaluations to

achieve the desired accuracy of 510−9 and one the calculation is per-
formed with M = 64 the result is
 
5
S f, = 1.252762973.
128

So we see that the composite Simpson rule using 229 evaluations of

f (x) and the composite trapezoidal rule using 22823 evaluations of f (x)
achieve the same accuracy. In Example 9, Simpson’s rule required
1
about 100 the number of function evaluations.

Professor Jorge Bacca Numerical methods 2024 43 / 43