TRIMESTER 2 TMA1301 COMPUTATIONAL METHODS CHAPTER 3

CHAPTER 3 NUMERICAL INTEGRATION

In this Chapter, you will learn:

• the definite integral,

• three methods to do the numerical integration when the integral is complicated.
• error analysis

1. DEFINITE INTEGRAL
b
The exact value of a definite integral ∫a
f ( x)dx can be computed only when the function
f ( x ) is integrable in finite terms.

Example 1: Find the following definite integral.

2
2 x3
∫ x dx =
2
0 3 0
3
2 03
= −
3 3
8
=
3

Example 2: Find the following definite integral.

2
∫0
3x 2 + e x dx =

The definite integral of a function has an interpretation as the area under a curve.
Whenever the function f ( x ) cannot be exactly integrated in finite terms or the evaluation
of its integral is complicated, we can use the numerical integration. The numerical
integrations that we are going to discuss in this chapter are Trapezoidal rule, Simpson’s
rule and Romberg algorithm.

KCY/NOORSHAHIDA/THL Page | 1
TRIMESTER 2 TMA1301 COMPUTATIONAL METHODS CHAPTER 3

2. TRAPEZOIDAL RULE

A1 A2 An

a = x1 x2 x3 xn-1 xn xn+1= b

Let CD be the curve and its equation is y = f ( x ) . To evaluate

b
∫
a
f ( x)dx , we subdivide

the interval from a to b into n equal parts, so that each part is equal to h =
(b − a ) .
n

The area of the first trapezium, A1 =

h
[ f (x1 ) + f (x 2 )]
2

The area of the second trapezium, A2 =

h
[ f (x2 ) + f (x3 )]
2

And the area of nth trapezium, An =

h
[ f (xn ) + f (xn+1 )]
2
Total area,
AT = ∫ f ( x )dx ≈
h
{ f (x1 ) + f (x n+1 ) + 2[ f (x 2 ) + f (x3 ) + ... + f (x n )]} where x n = a + (n − 1)h, n = 1,2,..
b

a 2
The above formula is known as the Composite Trapezoidal rule for numerical integration.
Example 3: Dividing the range into 4 equal parts, find the approximate value of
1 1
∫0 1 + x 2 dx by using Composite Trapezoidal rule.
1− 0
n = 4; f ( x ) =
1
; h= = 0.25
1+ x 2
4

AT =
0.25
{ f (0) + f (1) + 2[ f (0.25) + f (0.5) + f (0.75)]}
2
=
0.25
(6.262353)
2
= 0.782794

KCY/NOORSHAHIDA/THL Page | 2
TRIMESTER 2 TMA1301 COMPUTATIONAL METHODS CHAPTER 3

Example 4: Compute an approximate value for integral given in Example 2 by using the
Composite Trapezoidal rule with 5 points. Then compare with the actual value of the
integral and find its absolute error.

2−0
Integral in Example 2: ∫ 3 x 2 + e x dx ; 5 points ⇒ n = 4 ; f ( x ) = 3 x 2 + e x ; h =
2
= 0.5
0 4

AT =
0.5
{ f (0) + f (2) + 2[ f (0.5) + f (1.0) + f (1.5)]}
2
=
0.5
(59.08644)
2
= 14.77161
2
Actual value, I = ∫ 3 x 2 + e x dx = __________________
0

Absolute Error = Actual value − Approximated value

2.1 Error Analysis

In this section, error incurred in using the Trapezoidal rule to estimate an integral is
analyzed.

Theorem on Precision of Trapezoidal Rule

If f ' ' exists and is continuous on the interval [a, b] and if the Composite Trapezoidal
rule T with uniform spacing h is used to estimate the integral I = ∫ f ( x )dx , then for
b

some ζ in (a, b ) ,

I −T = −
1
(b − a )h 2 f ' ' (ζ ) = ο (h 2 )
12

KCY/NOORSHAHIDA/THL Page | 3
TRIMESTER 2 TMA1301 COMPUTATIONAL METHODS CHAPTER 3

Proof:
Firstly, prove the theorem when a = 0, b = 1 and h = 1. In this case, we have to show that

∫ f (x )dx − 2 [ f (0) + f (1)] = − 12 f ' ' (ζ )

1 1 1
(1)
0

Let p be the polynomial of degree 1 that interpolates f at 0 and 1. Then p is given by

p( x ) = f (0) + [ f (1) − f (0 )]x

∫ f (x )dx = f (0) + 2 [ f (1) − f (0)] = 2 [ f (0) + f (1)]

1 1 1
Hence,
0

By the following theorem,

Theorem on Interpolation Errors I

If p is the polynomial of degree at most n that interpolates f at the n + 1 distinct nodes
x0 , x1 ,..., xn belonging to an interval [a, b ] and if f (n+1) is continuous, then for each x in
[a, b], there is a ξ in (a, b ) for which
n
f (x ) − p(x ) = f (n+1) (ξ )∏ ( x − xi )
1
(n + 1)! i =0

we have (n = 1, x0 = 0, x1 = 1)

f (x ) − p(x ) = f ' ' [ξ ( x )]x ( x − 1)

1
(2)
2

where ξ ( x ) depends on x in (0,1) . From Equation (2), it follows that

f ( x )dx − ∫ p( x )dx = f ' ' [ξ ( x )]x ( x − 1)dx

1 1 1 1
∫0 0 2 ∫0

That f ' ' [ξ ( x )] is continuous can be proven by solving Equation (2) for f ' ' [ξ ( x )] and
verifying the continuity. Notice that x ( x − 1) does not change sign in the interval [0,1] .
Hence, by the Mean-Value Theorem for Integrals (refer to page 5), there is a point x = s
for which

f ' ' [ξ ( x )]x( x − 1)dx = f ' ' [ξ (s )]∫ x( x − 1)dx = − f ' ' (ξ )
1 1 1
∫
0 0 6
By putting all these equations together, Equation (1) is obtained. From Equation (1), by
making a change of variable, the Basic Trapezoidal rule with its error term is obtained.

KCY/NOORSHAHIDA/THL Page | 4
TRIMESTER 2 TMA1301 COMPUTATIONAL METHODS CHAPTER 3

b−a
∫ f (x )dx = 2 [ f (a ) + f (b )] − 12 (b − a ) f ' ' (ξ )
1 1 3
(3)
0

b−a
∫ f (x )dx − 2 [ f (a ) + f (b)] = − 12 (b − a ) f ' ' (ξ )
1 1 3
0

This is the Trapezoidal rule and error term for the interval [a, b ] with only one subinterval,
which is the entire interval. Thus, the error term is ο h 3 , where h = b − a . Here, ξ is in ( )
(a, b ) .
Now, let the interval [a, b] be divided into n equal subintervals by points x0 , x1 ,..., xn with
spacing h. Applying Equation (3) to subinterval [xi , xi +1 ] , we have

f ( x )dx =
h
[ f (xi ) + f (xi+1 )] − 1 h 3 f ' ' (ξ i )
xi +1
∫
xi 2 12

where xi < ξ i < xi +1 . We use this result over the interval [a, b ] , obtaining the Composite
Trapezoidal rule
n −1
h n−1 h 3 n−1
f ( x )dx = ∑ ∫ f ( x )dx = ∑ i [ ( ) ( )] ∑ f ' ' (ξ i )
b xi +1
∫a
i =0
xi 2 i =0
f x + f x i +1 −
12 i =0
(4)

The final term in Equation (4) is the error term and it can be simplified in the following
b−a
way: Since h = , the error term for the Composite Trapezoidal rule is
n

h 3 n−1 b − a 2  1 n−1  b−a 2

− ∑ f ' ' (ξ i ) = − h  ∑ f ' ' (ξ i ) = − h f ' ' (ζ )
12 i =0 12  n i =0  12

 1 n −1 
The average  ∑ f ' ' (ξ i ) lies between the least and the greatest values of f ' ' on the
 n i =0 
interval (a, b ) . Hence, by the Intermediate-Value Theorem, it is f ' ' (ζ ) for some point
ζ in (a, b ) . This completes the proof of the error formula.

Mean-Value Theorem for Integrals Intermediate-Value Theorem

Let f be continuous on [a, b] and assume If the function g is continuous on an

that g is Riemann-integrable on [a, b] . If interval [a, b] , then for each c between g (a )
g ( x ) > 0 on [a, b] , then there exists a point and g (b ) , there is a point ξ in [a, b] for
which g (ξ ) = c .
=
ξ such that a ≤ξ ≤b and
f ( x )g ( x )dx = f (ξ )∫ g ( x )dx
b b
∫a a

KCY/NOORSHAHIDA/THL Page | 5
TRIMESTER 2 TMA1301 COMPUTATIONAL METHODS CHAPTER 3

Error formula for Trapezoidal rule:

1. Basic Trapezoidal rule: −

1
(b − a )3 f ' ' (ξ ) for a < ξ < b
12

b−a 2
2. Composite Trapezoidal rule: − h f ' ' (ξ ) for a < ξ < b
12
1
− x2
Example 5: If the Composite Trapezoidal rule is to be used to compute ∫e
0
dx with an
1
error of at most × 10 −4 , how many points should be used? Reminder:
2
n subintervals = (n+1) points
b−a 2
The error formula is − h f ' ' (ζ ) .
12
2 2
( )
f ( x ) = e − x , f ' ( x ) = −2 xe − x , f ' ' ( x ) = 4 x 2 − 2 e − x
2

Thus, f ' ' ( x ) ≤ 2 on the interval [0,1] .

b−a 2 b−a 2
h f ' ' (ζ ) = h f ' ' (ζ ) = h 2 (2) = h 2
1 1
−
12 12 12 6

1 1 1
To have an error of at most × 10 −4 , we require h 2 ≤ × 10 −4 or h ≤ 3 × 10 −2 .
2 6 2
h = 1 / n : so n ≥ 58 . Hence, 59 or more points will certainly produce the desired accuracy.

1 sin x
Example 6: How many subintervals are needed to approximate ∫0 x
dx with error not to
1
exceed × 10 −5 using the Composite Trapezoidal rule? Here, the integrand,
2
f ( x ) = x −1 sin x , is defined to be 1 when x is 0.

b−a 2
The error formula is − h f ' ' (ζ ) .
12

Now, establish a bound on f ' ' ( x ) for x in the range [0,1] . Taking derivatives in the usual
way is not satisfactory because each term contains x with a negative power, and it is
difficult to find an upper bound on f ' ' ( x ) . However using Taylor series, we have

KCY/NOORSHAHIDA/THL Page | 6
TRIMESTER 2 TMA1301 COMPUTATIONAL METHODS CHAPTER 3

x 2 x 4 x6 x8
f (x ) = 1 − + − + − ...
3! 5! 7! 9!
2 x 4 x 3 6 x 5 8x 7
f ' (x ) = − + − + − ...
3! 5! 7! 9!
2 3× 4x 2 5 × 6x 4 7 × 8x6
f ' ' (x ) = − + − + − ...
3! 5! 7! 9!

Thus, on the interval [0,1] , f ' ' ( x ) cannot exceed

1
because
2

2 3× 4 5 × 6 7 × 8 1 1 1 1 1 1  1
+ + + + ... < + +  + + + ...  <
3! 5! 7! 9! 3 10 24  2 4 8  2

Therefore, the error term (b − a )h 2 f ' ' (ζ ) / 12 cannot exceed h 2 / 24 .

1 2 1  1 2
− (b − a )h 2 f ' ' (ζ ) / 12 = (b − a )h 2 f ' ' (ζ ) / 12 = h  = h
12  2  24

1 1 2 1
To have an error of at most ×10 −5 , we require h ≤ × 10 −5 or h ≤ 1.2 × 10 −2 .
2 24 2

( )
n ≥ 1 / 1.2 10 2 = 91.29 . This analysis shows that 92 subintervals are needed.

3. SIMPSON’S RULE

A1 A2 An

a = x1 x2 x3 x2n-1 x2n x2n+1= b

Let CD be the curve and its equation is y = f ( x ) . To evaluate

b
∫ a
f ( x)dx , we subdivide

the interval from a to b into 2n equal parts, so that each part is equal to h =
(b − a ) .
2n

By Simpson’s parabolic rule,

KCY/NOORSHAHIDA/THL Page | 7
TRIMESTER 2 TMA1301 COMPUTATIONAL METHODS CHAPTER 3

the first area, A1 = h[ f ( x1 ) + 4 f ( x 2 ) + f ( x3 )]

1
3

the second area, A2 = h[ f ( x3 ) + 4 f ( x 4 ) + f ( x5 )]

1
3

and the nth area, An = h[ f ( x 2 n −1 ) + 4 f ( x 2 n ) + f ( x 2 n +1 )]

1
3

Total area,
b
h  f ( x1 ) + f ( x 2 n +1 ) + 4[ f ( x 2 ) + f ( x 4 ) + ... + f ( x 2 n )]
f ( x )dx ≈   where x n = a + (n − 1)h, n = 1,2,..
AS =
∫a 3 + 2[ f ( x3 ) + f ( x5 ) + ... + f ( x 2 n −1 )] 

Example 7: Dividing the range into 4 equal parts, find the approximate value of
1 1
∫0 1 + x 2 dx by using Simpson’s rule.
11
AS =  { f (0) + f (1) + 4[ f (0.25) + f (0.75)] + 2 f (0.5)}
3 4

= (9.42470588)
1
12
= 0.785392

Example 8: Compute an approximate value for integral given in Example 2 by using the
Simpson’s rule with 5 points. Then compare with the actual value of the integral and find
its absolute error.

2−0
+ e x dx ; 5 points ⇒ n = 4 ; f ( x ) = x 2 ; h =
2
∫ 3x = 0.5
2
Integral in Example 2:
0 4

AS =
1
(0.5){ f (0) + f (2) + 4[ f (0.5) + f (1.5)] + 2[ f (1.0)]}
3
=
0.5
(86.34726)
3
= 14.39121
2
Actual value, I = ∫ 3 x 2 + e x dx = __________________
0

Absolute Error = Actual value − Approximated value

KCY/NOORSHAHIDA/THL Page | 8
TRIMESTER 2 TMA1301 COMPUTATIONAL METHODS CHAPTER 3

3.1 Error Analysis

A numerical integration rule over two equal subintervals with partition points
a = x1 , a + h and a + 2h = b = x2 n+1 is the widely used Basic Simpson’s rule:

f ( x )dx ≈ [ f (a ) + +4 f (a + h ) + f (a + 2h )]
a+2h h
∫
a 3
(5)

The error term in Simpson’s rule can be established by using the Taylor series:

f (a + h ) = f + hf '+
1 2 1 1
h f ' '+ h 3 f ' ' '+ h 4 f (4 ) + ...
2! 3! 4!

Now replacing h by 2h:

24
f (a + 2h ) = f + 2hf '+2h 2 f ' '+ h 3 f ' ' '+ h 4 f (4 ) + ...
4
3 4!

where the functions f , f ' , f ' ' ,... on the right-hand side are evaluated at a.

Using these two series, we obtain

f (a ) + 4 f (a + h ) + f (a + 2h ) = 6 f + 6hf '+4h 2 f ' '+2h 3 f ' ' '+

20 4 (4 )
h f + ...
4!

and thereby, we have

h
[ f (a ) + 4 f (a + h ) + f (a + 2h )] = 2hf + 2h 2 f '+ 4 h 3 f ' '+ 2 h 4 f ' ' '+ 20 h 5 f (4 ) + ... (6)
3 3 3 3 ⋅ 4!

Hence, we have a series for the right-hand side of Equation (5). Now let’s find one for the
left-hand side. The Taylor series for F (a + 2h ) is

25
F (a + 2h ) = F (a ) + 2hF ' (a ) + 2h 2 F ' ' (a ) + h 3 F ' ' ' (a ) + h 4 F (4 ) (a ) + F (5 ) (a ) + ...
4 2
3 3 5!

Let F ( x ) = ∫ f (t )dt . By the Fundamental Theorem of Calculus, F ' = f . We observe that

x

F (a ) = 0 and F (a + 2h ) is the integral on the left-hand side of Equation (5). Since

F ' ' = f ' , F ' ' ' = f ' ' , and so on, we have

a+2 h 2 5 5 (4 )
f ( x )dx = 2hf + 2h 2 f '+ h 3 f ' '+ h 4 f ' ' '+
4 2
∫
a 3 3 5 ⋅ 4!
h f + ... (7)

Subtracting Equation (6) from Equation (7), we obtain

KCY/NOORSHAHIDA/THL Page | 9
TRIMESTER 2 TMA1301 COMPUTATIONAL METHODS CHAPTER 3

5
f ( x )dx = [ f (a ) + 4 f (a + h ) + f (a + 2h )] − h f (4 ) − ...
a+2 h h
∫a 3 90

We can rewrite the Basic Simpson’s rule over the interval [a, b ] as

(b − a )  a+b 
∫ f (x )dx ≈ 6  f (a ) + 4 f  2  + f (b )
b

h f (ξ ) for some ξ in (a, b ).

1 5 (4 )
with error term −
90

Suppose that the interval [a, b ] is subdivided into an even number of subintervals, say n,
each of width h = (b − a ) / n . Then the partition points are xi = a + ih for 0 < i < n , where
= =

n is divisible by 2. Now from basic calculus, we have

n/2
f ( x )dx = ∑ ∫ f ( x )dx
b a + 2 ih
∫a
i =1
a + 2 (i −1)h

Using the Basic Simpson’s rule, we have, for the right-hand side,
n/2
≈∑
h
{ f (a + 2(i − 1)h ) + 4 f (a + (2i − 1)h ) + f (a + 2ih )}
i =1 3
( n / 2 )−1 (n / 2 ) ( n / 2 )−1
h 
=  f (a ) + ∑ f (a + 2ih ) + 4 ∑ f (a + (2i − 1)h ) + ∑ f (a + 2ih ) + f (b )
3 i =1 i =1 i =1 

Thus, we obtain the Composite Simpson’s rule over n (even) subintervals

(n / 2 ) (n / 2 ) / 2
h 
f ( x )dx ≈ [ ( ) ( )] ∑ ( ( ) ) ∑ f (a + 2ih ) where h = (b − a ) / n .
b
∫a

3
f a + f b + 4
i =1
f a + 2i − 1 h + 2
i =1 

The error term is −

1
(b − a )h 4 f (4 ) (ξ ) for some ξ in (a, b ).
180

Error formula for Simpson’s rule:

h f (ξ ) for some ξ in (a, b )

1 5 (4 )
1. Basic Simpson’s rule: −
90

2. Composite Simpson’s rule: −

1
(b − a )h 4 f (4 ) (ξ ) for some ξ in (a, b )
180

KCY/NOORSHAHIDA/THL Page | 10
TRIMESTER 2 TMA1301 COMPUTATIONAL METHODS CHAPTER 3

2
Example 9: Apply the Basic Simpson’s rule to approximate ∫ ln xdx , and find an upper
1

bound for the error in your approximations.

∫ ln xdx ≈ 3 [ f (a ) + 4 f (a + h) + f (a + 2h)]
2 h
1

=
0.5
[ f (1) + 4 f (3 / 2) + f (2)]
3
0.5  3 
=  ln 1 + 4 ln + ln 2 ≈ 0.3858
3  2 

h f (ξ ) where 1 < ξ < 2 . Since

1 5 (4 )
The error for Basic Simpson’s rule is −
90

f (4 ) (x ) = − , f (4 ) ( x ) = − 4 = 4 . Thus, f (4 ) ( x ) ≤ 6 on the interval [1,2] .

6 6 6
4
x x x

The error is at most − h f (ξ ) ≤

1 5 (4 ) (0.5) (6) = 1 ≈ 0.0021
5

90 90 480

4. ROMBERG ALGORITHM

The Romberg algorithm produces a triangular array of numbers, all of which are
b
numerical estimates of the definite integral ∫
a
f ( x)dx . The array is denoted here by the
notation
m
n 0 1 2 3 …

0 R(0,0 )
1 R(1,0 ) R(1,1)
2 R(2,0) R(2,1) R(2,2)
3 R(3,0) R(3,1) R(3,2 ) R(3,3)
M M M M M O
R(n,0 ) R(n,1) R(n,2 ) R(n,3) L R(n, m )

The first column of this table contains estimates of the integral obtained by the recursive
trapezoid formula with decreasing values of the step size. Explicitly, R(n,0 ) is the result

KCY/NOORSHAHIDA/THL Page | 11
TRIMESTER 2 TMA1301 COMPUTATIONAL METHODS CHAPTER 3

of applying the trapezoid rule with 2 equal subintervals. The first of them, R(0,0 ) , is
n

obtained with just one trapezoid:

R (0,0 ) =
1
(b − a )[ f (a ) + f (b )]
2

Note that R(n,0 ) is obtained easily from R(n − 1,0 ) , that is

n −1
R (n,0) =
1 2
R (n − 1,0) + h ∑ f [a + (2k − 1)h] ; h =
(b − a ) , n ≥ 1
2 k =1 2n

The second and successive columns in the Romberg array are generated by the
extrapolation formula

R(n, m ) = R(n, m − 1) +
1
[R(n, m − 1) − R(n − 1, m − 1)]
4 −1
m

with n ≥ 1 , m ≥ 1 and m ≤ n .

Example 10: If R (8,2 ) = 8 and R (7,2 ) = 1 , what is R (8,3) ?

R(8,3) = R(8,2 ) +
1
[R(8,2) − R(7,2)] = 73
63 9

1 1
Example 11: Compute ∫ 1+ x
0 2
dx by the Romberg algorithm using 4 equal parts and

show your calculations correct to SIX decimal places. Then, compare with the values
obtained from the methods of Composite Trapezoidal rule (Example 3) and Simpson’s
rule (Example 7).

n 0 1 2 3

0 0.750000

1 0.775000 0.783333

2 0.782794

KCY/NOORSHAHIDA/THL Page | 12