NUMERICAL INTEGRATION

INTRODUCTION
Numerical integration is used to obtain approximate answers for definite integrals that cannot be
solved analytically.
Numerical integration is a process of finding the numerical value of a definite integral

z bg
b
I = f x dx ,
a

when a function y = f (x) is not known explicitly. But we give only a set of values of the function
y = f (x) corresponding to the same values of x.
To evaluate the integral, we fit up a suitable interpolation polynomial to the given set of values
of f (x) and then integrate it within the desired limits. Here we integrate an approximate interpolation
formula instead of f (x). When this technique is applied on a function of single variable, the process
is called Quadrature.

z bg
b
Suppose we are required to evaluate the definite integral I = f x dx , firstly we have to
a

approximate f (x) by a polynomial φ( x ) of suitable degree.

Then we integrate f(x) within limits [a, b],

z bg zbg
b b
i.e., f x dx ≈ φ x dx ,
a a
the difference

LM f x dx − φ x dx ,OP ,
z
MN b g bg P z
b b

a Q
a

is called the Error of approximation.

GENERAL QUADRATURE FORMULA FOR EQUIDISTANT ORDINATES
Consider an integral

z bg
b
I = f x dx (1)
a

Let f (x) take the values f (x0) = y0, f (x0 + h) = y1, …, f (x + nh) = yn, when x = x0, x = x0
+ h, …, x = x0 + nh respectively.
To evaluate I, we replace f (x) by a suitable interpolation formula. Let the interval [a, b] be
divided into n subintervals with the division points a = x0 < x0 + h < … < x0 + xh = b where the
h is the width of each subinterval. Approximating f (x) by Newton’s forward interpolation formula
we can write the integral (1) as

z bg z FGH b g∆ y I
x0 + nh x0 + nh
u u−1
I =
x0
f x dx =
x0
y0 + u∆y0 +
2!
2
0 JK
+ ... dx ,

(2)
x − x0
since u= ,
h
i.e., x = x0 + uh

⇒ dx = hdu

and x = x0
⇒ u=0

x = x0 + nh
⇒ u = n.
Expression (2) can be written as

z
F y + u∆y + u − u ∆ y + u − 3u + 2u ∆ y I dx +
n

GH JK
2 3 2
I =h 0 0
2
0
3
0
0
2 6

z
F u − 6u + 11u + 6u ∆ y + ...I dx
n
h G+ JK
4 3 2

H
4
0
0
24

L n ∆y + F n − n I ∆ y +
∴ I = h Mny + GH 3 2 JK 2
2 3 2 3

MN 2
0
0 0
 Fn I∆y Fn I OP
GH 4 JK 6
4 3
∆4 y0
GH 5 JK
5
− n3 + n2 0
+ 3 4 11n 3
− n + − 3n 2 + ... .
2 3 24 PQ (3)

The equation (3) is called General Gauss Legendre Quadrature formula, for equidistant ordinates
from which we can generate any Numerical integration formula by assigning suitable positive integral
value to n. Now we deduce four quadrature formulae, namely
(a) Trapezoidal rule (b) Simpson’s one-third rule (c) Simpson’s three-eighths rule and (d)
Weddle’s rule from the general quadrature formula (3).
TRAPEZOIDAL RULE
Substituting n = 1 in the relation (3) and neglecting all differences greater than the first we get

z LM OP
x0 + h

I1 = bg
f x dx = h y0 +
N
1
2
∆y 0
Q
x0

=
h
2
b
2 y0 + y1 − y0 =
h
2
y0 + y1 , g b g
for the first subinterval [x0, x0 + h],
similarly, we get

z bg
x0 + 2 h

l2 = f x dx =
h
2
b
y1 + y2 , g
x0 + h

z bg
x0 + 3h

l3 = f x dx =
h
2
b
y2 + y3 , g
x0 + 2 h

L
L
L

z bg
x0 + nh

In = f x dx =
h
d
yn − 1 + yn , i
b g
x0 + n − 1 h
2

for the other integrals.

Adding I1, I2, ..., In
we get I 1 + I 2 + ... + I n

z z z z
x0 + h x0 + 2 h x0 + 3h x0 + nh

= f ( x ) dx + f ( x ) dx + f ( x ) dx + ... + f ( x ) dx
x0 x0 + h x0 + 2 h x0 + ( n −1) h
 h h h h
= y0 + y1 + y1 + y2 + y2 + y3 + ... + yn − 1 + yn ,
2 2 2 2

z bg
x0 + nh

⇒ f x dx =
h
2
b g d
y0 + yn + 2 y1 + y2 + ... + yn − 1 , i
x0

z bg b g d i
b
h
I = f x dx = y0 + yn + 2 y1 + y2 + ... + yn − 1 . (4)
a
2

The formula (4) is called Trapezoidal rule for numerical integration. The error committed in this
formula is given by

bg b g f ′′bξg ,
3
h3 − b−a
E≈− f ′′ ξ =
12 12n 2
where a = x0 < ξ < xn = b .
Note: Trapezoidal rule can be applied to any number of subintervals odd or even

SIMPSON’S ONE-THIRD RULE

Substituting n = 2 in the General quadrature formula given by (3) and neglecting the third and other
higher order differences
we get

z LM FG 8 − 2IJ ∆ y OP
x0 + 2 h

I1 = bg
f x dx = h 2 y0 + 2 ∆y0 +
N H3 K Q
2
0
x0

LM b g 13 b y gOPQ
N
= h 2 y0 + 2 y1 − y0 + 2 − 2 y1 + y0

h
= y0 + 4 y1 + y2
3

h
∴ I1 = y0 + 4 y1 + y2 ,
3

z
x0 + 2 h

Similarly I2 = bg
f x dx =
h
3
y 2 + 4 y3 + y4 ,
x0 + 2 h

L
L
L
 z
x0 + nh

In/2 = bg
f x dx =
h
3
yn − 2 + 4 yn − 1 + yn .
x0 + b n − 2 g h

Adding I 1 , I 2 , I n / 2 we get

z bg z bg z bg
x0 + 2 h x0 + 2 h x0 + nh

I 1 + I 2 + ... + I n / 2 = f x dx + f x dx + ... + f x dx
x0 x0 + 2 h x0 + b n − 2 g h

h h h
= y0 + 4 y1 + y2 + y2 + 4 y3 + y4 + ... + yn − 2 + 4 yn − 1 + y n ,
3 3 3

h
= ( y0 + yn ) + 4 ( y1 + y3 + y5 + ... + yn −1 ) + 2 ( y2 + y4 + ... + yn − 2 )
3
h
= y0 + yn + 4 × (sum of odd ordinates) + 2 × (sum of even ordinates)
3

=
h
3
b g b
y0 + yn + 4 × sum of the odd ordinates + g
h
3
= b
2 × sum of the even ordinates g
The above rule is known as Simpson’s one-third rule. The error committed in Simpson’s one-
third rule is given by

bg b g bg
5
− nh5 iv b−a
E≈ f ξ =− f iv ξ
180 2880n 4
where a = x0 < ξ < xn = b (for n subintervals of lengths h).
Note:
1. The above formula may written as

z
x0 + nh

I= f ( x ) dx
x0

2. Simpson’s one-third rule can be applied only when the given interval [a, b] is subdivided into even number
of subintervals each of width h and within any two consecutive subintervals the interpolating polynomial
φ(x) is of degree 2.

SIMPSON’S THREE-EIGHTH’S RULE

We assume that within any three consecutive subintervals of width h, the interpolating polynomial
bg
φ x approximating f (x) is of degree 3. Hence substituting n = 3, i.e., the General quadrature formula
and neglecting all the differences above ∆3 , we get
 z LM FG IJ FG 81 − 27 + 9IJ ∆ y OP
x0 + 3h

bg
3
9 9 2
I1 =
x0
f x dx = h 3 y0 +
N 2
∆y0 + 9 −
2 H
∆ y0 +
K H4 K 6 Q
0

LM b g 94 b y g 3
b gOPQ
N
= h 3 y0 + 9 y1 − y0 + 2 − 2 y1 + y0 +
8
y3 − 3 y2 + 3 y1 + y0

3h
= y0 + 3 y1 + 3 y2 + y3 ,
8

z
x0 + 6 h

Similarly I2 = bg
f x dx =
3h
8
y3 + 3 y4 + 3 y5 + y6 ,
x0 + 3h

z
x0 + 9 h

I3 = bg
f x dx =
3h
8
y6 + 3 y7 + 3 y8 + y9 ,
x0 + 6h

L
L
L

z
x0 + nh

In/3 = bg
f x dx =
3h
yn − 3 + 3 y n − 2 + 3 yn − 1 + yn .
b g
x0 + n − 3 h
8

Adding I1 , I 2 , ..., I n / 3 we get

z bg z bg z bg
x 0 + 3h x0 + 6 h x0 + nh

I 1 + I 2 + ... + I n / 3 = f x dx + f x dx + ... + f x dx ,
x0 x 0 + 3h b g
x0 + n − 3 h

3h 3h
⇒ I = y0 + 3 y1 + 3 y2 + y3 + y3 + 3 y4 + 3 y5 + y6 +
8 8

3h
... + y n − 3 + 3 yn − 2 + 3 yn − 1 + y n ,
8

∴ I =
3h
8
b g d
y0 + yn + 3 y1 + y2 + y4 + y5 + ... + yn − 1 + i
d
2 y3 + y6 + ... + yn − 3 . i
Note:
1. Simpson’s three-eighths rule can be applied when the range [a, b] is divided into a number of subintervals,
which must be a multiple of 3.
 =
3h
10
b g b
y0 + yn + y2 + y4 + y8 + y10 + y14 + y16 + ... +

i d
yn − 4 + yn − 2 + 5 y1 + y5 + y7 + y11 + ... + yn − 5 + yn − 1 + i
d i d
6 y3 + y9 + y15 + ... + yn − 3 + 2 y6 + y12 + ... + yn − 6 . i
Note:
1. Weddle’s rule requires at least seven consecutive equispaced ordinates with in the given interval (a, b).
2. It is more accurate than the Trapezoidal and Simpson’s rules.
3. If f (x) is a polynomial of degree 5 or lower, Weddle’s rule gives an exact result.

z
x
x
Example 9.1 Calculate the value dx correct up to three significant figures taking six intervals by Trapezoidal
1+ x
0
rule.
Solution Here we have

bg
f x =
x
1+ x
,

a = 0, b = 1 and n = 6,

b − a 1− 0 1
∴ h= = = .
n 6 6
x 0 1/6 2/6 3/6 4/6 5/6 6/6 = 1
y = f (x) 0.00000 0.14286 0.25000 0.33333 0.40000 0.45454 0.50000
y0 y1 y2 y3 y4 y5 y6

The Trapezoidal rule can be written as

I =
h
2
b g b
y0 + y6 + 2 y1 + y2 + y3 + y4 + y5 g
=
1
12
(0. 00000 + 0.50000 + 2 014286
. b
+ 0.25000 + 0.3333 + 0.40000 + 0.45454 g
= 0.30512 .
∴ I = 0.0305 , correct to three significant figures.

z
1
dx
Example 9.2 Find the value of , taking 5 subinterval by Trapezoidal rule, correct to five significant figures.
1 + x2
0
Also compare it with its exact value.
Solution Here

bg
f x =
1
1 + x2
,
 a = 0, b = 1 and n = 5,
1− 0 1
∴ h= = + 0.2 .
5 5
x 0.0 0.2 0.4 0.6 0.8 1
y = f (x) 1.000000 0.961538 0.832069 0.735294 0.609756 0.500000
y0 y1 y2 y3 y4 y5

Using trapezoidal rule we get

z b g b g
1
dx h
I = = y0 + y5 + 2 y1 + y 2 + y3 + y4
0
1+ x 2
2

=
0.2
2
b g b
1.000000 + 0.500000 + 2 0.961538 + 0.862069 + 0.735294 + 0.609756 g
= 0.783714 ,
∴ I = 0.78373 , correct to five significant figures.
The exact value

z
1
1 1
= dx = tan −1 x
0
1 + x 2 0

π
= tan −1 1 − tan −1 0 = = 0.7853981
4

z
1
1
dx = 0.78540 ,
0
1 + x2

correct to five significant figures.

z
5

Example 9.3 Find the value of log10 xdx, taking 8 subintervals correct to four decimal places by Trapezoidal rule.
1
Solution Here

bg
f x = log10 x ,
a = 1, b = 5 and n = 8,

b −a 5−1
∴ h= = = 0.5.
n 8
x 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
f (x) 0.00000 0.17609 0.30103 0.39794 0.47712 0.54407 0.60206 0.65321 0.69897
y0 y1 y2 y3 y4 y5 y6 y7 y8
 Using Trapezoidal rule we can write

I =
h
2
b g b
y0 + y8 + 2 y1 + y2 + y3 + y 4 + y5 + y6 + y7 g
=
0.5
2
b g b
0.00000 + 0.69897 + 2 0.17609 + 0.30103 + 0.39794 + g
0.5
2
b
2 0.47712 + 0.54407 + 0.60206 + 0.65321 g
= 17505025
.

z
5
∴ I = log10 xdx = 1.75050 .
1

z
0.6

Example 9.4 Find the value e x dx, taking n = 6, correct to five significant figures by Simpson’s one-third rule.
0

Solution We have

bg
f x = ex ,
a = 0, b = 0.6, n = 6.

b − a 0.6 − 0
∴ h= = = 0.1.
n 6
x 0.0 0.1 0.2 0.3 0.4 0.5 0.6
y = f (x) 1.0000 1.10517 1.22140 1.34986 1.49182 1.64872 1.82212
y0 y1 y2 y3 y4 y5 y6

The Simpson’s rule is

I =
h
3
b g b g b
y 0 + y6 + 4 y1 + y3 + y5 + 2 y2 + y 4 g
=
0.1
3
b
1.00000 + 182212
. . g b
+ 4 110517 + 134986
. + 1.64872 + 2 122140
. + 149182
. g b g
=
01
3
.
b g b
2.82212 + 4 4.10375 + 2 2.71322 g b g
= 0.8221186 ≈ 0.82212
∴ I = 0 82212
 x 1.50 1.65 1.80 1.95 2.10 2.25 2.40 2.55 2.70 2.85 3.00

A
y= 1.025 1.081 1.132 1.182 1.249 1.308 1.375 1.438 1.498 1.571 1.632
x

Using Simpson’s rule, we get

z b g b g b g
3.0
A h
dx = y0 + y10 + 4 y1 + y 3 + y5 + y 7 + y9 + 2 y2 + y4 + y6 + y8
1.5
x 3

=
0.15
3
1025
. b
+ 1.632 + g b
4 1081
. + 1182
. + 1308
. + 1438
. + 1571
. + g
0.15
3
2 1132
. b
+ 1249
. + 1.375 + 1.498 g
= 19743
.

z
3
A
∴ dx = 1.9743.
1.5
x

b g z
3
A
Using the formula 0.018 T = dx ,
1.5
x

we get b0.018g T = 19743

.

⇒ T=
1.9743
0.018
= 110 sec approximately b g
∴ T = 110 sec.

z
1
1
Example 9.7 Evaluate dx, by taking seven ordinates.
1 + x2
0

Solution We have
n + 1= 7 ⇒ n = 6
The points of division are

1 2 3 4 5
0, , , , , , 1.
6 6 6 6 6
x 0 1/6 2/6 3/6 4/6 5/6 1

1
y= 1.0000000 0.9729730 0.9000000 0.8000000 0.6923077 0.5901639 0.5000000
1 + x2

1
Here h = , the Simpson’s three-eighths rule is
6