Numerical Integration

Numerical Integration: Numerical integration is the processes of computing the value

of a definite integral from a set of numerical values of the integrand. When applied to the
integration of a function of a single variable, the process is sometimes called mechanical
quadrature; when applied to the computation of a double integral. of a function of two
independent variables it is called mechanical cubature.
The problem of numerical integration, like that of numerical differentiation, is solved by
representing the integrand by an interpolation formula and then integrating this formula
between the desired limits.

General Quadrature Formula for Equidistant Ordinates:

Given a set of data points ( x0 , y0 ), ( x1 , y1 ), ....., ( xn , y n ) of a function y = f (x) , where
f (x) is known explicitly, it is required to compute the value of the definite integrals
b
I =  ydx (1)
a

We desire in this section a general formula for numerical integration using Newton’s
forward difference formula. Let interval [a, b] de derived into n equal subintervals such
that a = x0  x1  x2  ......xn = b. Clearly xn = x0 + nh . Hence the integral becomes
x0 + nh
I =  ydx (2)
x0

Approximate y by Newton’s forward formula, we obtain

  u (u − 1) 2 u (u − 1)(u − 2) 3 u (u − 1)(u − 2)(u − 3) 4 
xn 
y 0 + uy 0 +  y 0 +  y 0 +  y 0 + 
2! 3! 4!
I=   dx
x0 u (u − 1)(u − 2)(u − 3)(u − 4) 5 − − − − −
  y0 +
u (u 1)(u 2)(u 3)(u 4)(u 5)
 y0 + ..
6
 5! 6! 

 (u − 1) 2
2
(u − 3u + 2u ) 3
3 2
(u − 6u + 11u − 6u ) 4
4 3 2

xn  y 0 + uy 0 +  y 0 +  y 0 +  y 0 + 
2 6 24
I=  5  dx
x0  (u − 10u + 35u − 50u + 24u )  y − + − + −  
4 3 2 5 6 5 4 3 2 6
(u 15 u 85u 225u 274u 120 u ) y

0
+ 0

 120 720 
x − x0
Now, u = and x = x0 + hu, dx = hdu , we have
h
 (u − 1) 2
2
(u − 3u + 2u ) 3
3 2
(u − 6u + 11u − 6u ) 4
4 3 2

n 0
y + u y 0 +  y 0 +  y 0 +  y 0 + 
2 6 24
I =  5  hdu
0  (u − 10u + 35u − 50u + 24u )  y − + − + −  
4 3 2 5 6 5 4 3 2 6
(u 15 u 85u 225u 274u 120 u ) y

0
+ 0

 120 720 
 n2  n 3 n 2  2 y 0  n 4 2   y0
3
 n 5 3n 4 11n 3 2   y0
4


ny 0 + y 0 +  −   
+ −n −n 
3
 
+ − + − 3n   + 
xn
 2  3 2  2  4  6  5 2 3  24 
 =h  6
  5  7  6 
 n − 2n 5 + 35n − 50n + 12n 2   y 0 +  n − 15n + 17 n 5 − 225n + 274n − 60n 2   y 0 
x0 4 3 6 4 3

 6 4 3  120  7
  6 4 3  720 
 
(3)
From this general formula (3), we can obtain a variety of quadrature formulas by putting
n = 1, 2, 3,....etc.
(i) Trapezoidal Rule: putting n = 1 in the general formula (3) and neglecting all
differences above the first, we get
    h
y 0  = h  y 0 + ( y1 − y 0 ) =  y 0 + y1 
x1 1 1
 ydx =h  y 0 +
x0  2   2  2
For the next interval [ x1 , x2 ] ,we get in like manner
h
 y1 + y 2 
x2
 ydx =
x1 2
and so on. For the last interval [ xn−1 , xn ] , we get
h
 y n−1 + y n 
xn
 ydx =
xn −1 2
Combining all these expressions, we obtain the rule
h
 y0 + 2( y1 + y 2 + .... + y n−1 ) + y n 
xn
 ydx =
x0 2
This is known as composite trapezoidal rule.
The geometrical significance of this rule is that the curve y = f (x) is replaced by n
straight lines joining the points ( x0 , y0 ) and ( x1 , y1 ); ( x1 , y1 ) and ( x2 , y 2 ),...., ( xn−1 , y n−1 )
and ( xn , y n ). The area bounded by the curve y = f (x) , the ordinates x = x0 and x = xn
and the x -axis is then approximately equivalent to the sum of the areas of the n
trapeziums obtained.
(ii) Simpson’s 1/3 Rule: putting n = 2 in the general formula (3) and neglecting all
differences above the second, we get
x2  8 2 y 0   1 
 ydx =h  0
2 y + 2 y 0 + ( − 2)  = h 2 y 0 + 2( y1 − y 0 ) + ( y 2 − 2 y1 + y 0 )
x0
 3 2   3 

  ydx = =  y 0 + 4 y1 + y 2 
x2 h
x0 3
For the next interval [ x2 , x4 ] ,we get in like manner
h
 y 2 + 4 y3 + y 4 
x4
 ydx =
x2 3
and so on. For the last interval [ xn−2 , xn ] , we get
h
 y n−2 + 4 y n−1 + y n 
xn
 ydx =
xn − 2 3
Combining all these expressions, we obtain the rule
h
 y0 + 4( y1 + y3 + .... + y n−1 ) + 2( y 2 + y 4 + .... + y n−2 ) + y n 
xn
 ydx =
x0 3
This is known as composite Simpson’s 1/3 rule. It should be noted that the rule requires
the division of the whole range into an even number of subintervals of width h .

(iii) Simpson’s 3/8 Rule: putting n = 3 in the general formula (3) and neglecting all
differences above the third, we get
x3
 9 3 2 1 3 
 ydx =h 3 y 0 + y 0 +  y 0 +  y 0 
x0  2 4 8 
 9 3 1 
= h 3 y 0 + ( y1 − y 0 ) + ( y 2 − 2 y1 + y 0 ) + ( y3 − 3 y 2 + 3 y1 + y 0 )
 2 4 8 
3h
 y0 + 3 y1 + 3 y 2 + y3 
x3
 ydx = =
x0 8
For the next interval [ x3 , x6 ] , we get in like manner
3h
 y3 + 3 y 4 + 3 y 4 + y 6 
x6
 ydx =
x3 8
and so on. For the last interval [ xn−3 , xn ] , we get
3h
 yn−3 + 3 yn−2 + 3 yn−1 + yn 
xn
 ydx =
xn − 3 8
Combining all these expressions, we obtain the rule
3h
 y0 + 3( y1 + y 2 + .... + y n−1 ) + 2( y3 + y6 + .... + y n−3 ) + y n 
xn
 ydx =
x0 8
This is known as composite Simpson’s 3/8 rule.
(iv) Boole’s Rule: putting n = 4 in the general formula (3) and neglecting all differences
above the fourth, we get
2h
7 y0 + 32 y1 + 12 y 2 + 32 y3 + 7 y 4 
x4
 ydx =
x0 45
This is known as boole’s rule.
(v) Weddle’s Rule: putting n = 6 in the general formula (3) and neglecting all differences
above the six, we get
3h
 y0 + 5 y1 + y 2 + 6 y3 + y 4 + 5 y5 + y6 
x6
 ydx =
x0 10
For the next interval [ x6 , x12 ] ,we get in like manner
3h
 y6 + 5 y7 + y8 + 6 y9 + y10 + 5 y11 + y12 
x12
 ydx =
x6 10
and so on. For the last interval [ xn−3 , xn ] , we get
 3h
 y n−6 + 5 y n−5 + y n−4 + 6 y n−3 + y n−2 + 5 y n−1 + y n 
xn
 ydx =
xn − 6 10
Combining all these expressions, we obtain the rule
xn 3h  y 0 + 5( y1 + y 7 + .... + y n −5 ) + ( y 2 + y8 + .... + y n −4 ) + 6( y3 + y9 + .... + y n −3 ) 
 ydx =
x0 10 + ( y 4 + y10 + .... + y n −2 ) + 5( y5 + y11 + .... + y n −1 ) + 2( y 6 + y12 + .... + y n −6 ) + y n 

This is known as composite Weddle’s rule.

Error Analysis:
(b − a ) 2
(i) Error for Trapezoidal rule = − h f ( x )
12
(b − a ) 4 iv
(ii) Error for Simpson’s 1/3 rule = − h f (x)
180
(b − a ) 4 iv
(iii) Error for Simpson’s 3/8 rule = − h f (x)
80
2(b − a ) 6 vi
(iv) Error for Boole’s rule = − h f (x)
945
(b − a ) 6 vi
(v) Error for Weddle’s rule = − h f (x)
840
Where f n (x ) is the largest value of the n-th derivative.

1 1
Example1: Evaluate the definite integral I =  dx using a suitable numerical
01 + x
2

integration formula dividing the whole range of integration into six equal parts.
b − a 1− 0 1
Solution: h = = =
n 6 6
The vales of the integrand are tabulated as below:
x x0 = 0 x1 = 1 / 6 x2 = 2 / 6 x3 = 3 / 6 x4 = 4 / 6 x5 = 5 / 6 x6 = 1

1 /(1 + x 2 ) 1 36/37 36/40 36/45 36/52 36/61 1/2

 (i) By Trapezoidal rule:
1 h
 y0 + 2( y1 + y2 + y3 + y4 + y5 ) + y6 
1
 dx =
01 + x
2
2
=
1
1 + 2(36 / 37 + 36 / 40 + 36 / 45 + 36 / 52 + 36 / 52 + 36 / 61 + 1/ 2
26
= 0.78424078
(ii) By Simpson’s 1/3 rule:
1 h
 y 0 + 2( y 2 + y 4 ) + 4( y1 + y3 + y5 ) + y 6 
1
 dx =
01+ x
2
3
=
1
1 + 2(36 / 40 + 36 / 52) + 4(36 / 37 + 36 / 45 + 36 / 61) + 1 / 2
26
= 0.78539801
(iii) By Simpson’s 3/8 rule:
1 3h
 y0 + 2( y3 ) + 3( y1 + y2 + y4 + y5 ) + y6 
1
 dx =
01 + x
2
8
=
3
1 + 2(36 / 45) + 3(36 / 37 + 36 / 40 + 36 / 52 + 36 / 61) + 1/ 2
8 6
= 0.78539592
(iv) By Weddle’s rule:
1 3h
 y0 + 5 y1 + y2 + 6 y3 + y4 + 5 y5 + y6 
1
 dx =
01 + x
2
10
=
1
1 + 5(36 / 37) + 36 / 40 + 6(36 / 45) + 5(36 / 52) + 36 / 61 + 1/ 2
26
= 0.78539962

Exact:
1

1
01 + x
2
dx = 
tan −1
x1
0 = /4

For Simpson’s 1/3 rule:

 / 4 = 0.78539801
or ,  = 3.14159204
Error for Trapezoidal rule = 0.78539816 - 0.78424078 = 0.00115738
Error for Simpson’s 1/3 rule = 0.78539816 - 0.78539801 = 0.00000015

Error for Simpson’s 3/8 rule = 0.78539816 - 0.78539592 = 0.00000224

Error for Weddle’s rule = 0.78539816 - 0.78539962 = 0.00000146

Numerical Double Integration:

x j +1 xi +1
I =   f ( x, y )dxdy
xj xi

Where xi +1 = xi + h and y j +1 = y j + k .

(i) Trapezoidal Rule: n = 1 .

h y j +1
I =  [ f ( xi , y ) + f ( xi +1 , y )]dy
2 yj
hk
= [ f ( xi , y j ) + f ( xi +1 , y j ) + f ( xi , y j +1 ) + f ( xi +1 , y j +1 )] , Where f i , j = f ( xi , x j )
4
hk
= [ f i , j + f i +1, j + f i , j +1 + f i +1, j +1 ]
4
(ii) Trapezoidal Rule: n = 2 .
x j + 2 xi + 2
I=   f ( x, y ) dxdy
xj xi
 h y j +2
I=  [ f ( xi , y ) + 2 f ( xi +1 , y ) + f ( xi + 2 , y )]dy
2 yj
hk
= [ f ( xi , y j ) + 2 f ( xi +1 , y j ) + f ( xi + 2 , y j ) + 2( f ( xi , y j +1 ) + 2 f ( xi +1 , y j +1 ) + f ( xi + 2 , y j +1 ))
4
+ f ( xi , y j + 2 ) + 2 f ( xi +1 , y j + 2 ) + f ( xi + 2 , y j + 2 )]
hk
= [ f i , j + f i + 2, j + f i , j + 2 + f i + 2, j + 2 + 2( f i +1, j + f i , j +1 + f i +1, j + 2 + f i + 2, j +1 ) + 4 f i +1, j +1 ]
4
(iii) Simpson’s 1/3 Rule: n = 2 .
x j + 2 xi + 2
I=   f ( x, y ) dxdy
xj xi

h y j +2
I=  [ f ( xi , y ) + 4 f ( xi +1 , y ) + f ( xi + 2 , y )]dy
3 yj
hk
= [ f ( xi , y j ) + 4 f ( xi +1 , y j ) + f ( xi + 2 , y j ) + 4( f ( xi , y j +1 ) + 4 f ( xi +1 , y j +1 ) + f ( xi + 2 , y j +1 ))
9
+ f ( xi , y j + 2 ) + 4 f ( xi +1 , y j + 2 ) + f ( xi + 2 , y j + 2 )]
hk
= [ f i , j + f i + 2, j + f i , j + 2 + f i + 2, j + 2 + 4( f i +1, j + f i , j +1 + f i +1, j + 2 + f i + 2, j +1 ) + 16 f i +1, j +1 ]
9

Example1:
11 x+ y
I =  e dxdy
00

By Trapezoidal Rule: n = 2 .
Using Trapezoidal and Simpson’s rules with h = k = 0.5 , we have the following table of
x+ y
values of e .
y x 0 0.5 1.0
0 1 1.6487 2.7183
0.5 1.6487 2.7183 4.4817
1 2.7183 4.4817 7.3891
 hk
I= [ f i , j + f i + 2, j + f i , j + 2 + f i + 2, j + 2 + 2( f i +1, j + f i , j +1 + f i +1, j + 2 + f i + 2, j +1 ) + 4 f i +1, j +1 ]
4
0.25
= [1 + 2.7183 + 2.7183 + 7.3891 + 2(1.6487 + 1.6487 + 4.4817 + 4.4817) + 4  2.7183]
4
12.3050
= = 3.0762
4
By Simpson’s 1/3 Rule: n = 2 .
hk
I= [ f i , j + f i + 2, j + f i , j + 2 + f i + 2, j + 2 + 4( f i +1, j + f i , j +1 + f i +1, j + 2 + f i + 2, j +1 ) + 16 f i +1, j +1 ]
9
0.25
= [1 + 2.7183 + 2.7183 + 7.3891 + 4(1.6487 + 1.6487 + 4.4817 + 4.4817) + 16  2.7183]
9
26.59042
= = 2.9545
9
11 x+ y
Exact: I =   e dxdy = 2.9525
00

Error for Trapezoidal rule = 2.9525 - 2.3.0762 = 0.1237

Error for Simpson’s 1/3 rule = 2.9525 - 2.9545 = 0.0020