Numerical Integration

There are two main reasons for us to do numerical integration

(i) Analytical integration may be impossible or infeasible

(ii) We may wish to integrate tabulated data rather than known functions.

Here we are interested in finding the definite integration of f(x) between two limits say a and
b. further to note the definite integration provides the area under the curve y=f(x) bounded
between the x-axis and the ordinates at x=a and x=b.

Newton’s Cotes Formula

Computation of integral when f(x) is given in tabular form is called numerical integration
𝑏
where geometrically ∫𝑎 𝑓(𝑥) 𝑑𝑥 represents the area under the curve 𝑦 = 𝑓(𝑥), x-axis and the
ordinates at x=a and x=b.

Let the interval [𝑎, 𝑏] = [𝑥0 , 𝑥𝑛 ] be subdivided into n equal sub-intervals of width h.
𝑥−𝑥0
Then 𝑥 = 𝑥0 + 𝑛ℎ, given 𝑛 = 𝑎𝑛𝑑 𝑑𝑥 = ℎ𝑑𝑛.
ℎ

𝑏 𝑛
∫ 𝑓(𝑥) 𝑑𝑥 = ℎ ∫ 𝑓(𝑥0 + 𝑛ℎ) 𝑑𝑛 (𝐴)
𝑎 0

We know that the Newton’s forward formula is given by

𝑛(𝑛 − 1) 2 𝑛(𝑛 − 1)(𝑛 − 2) 3

𝑓(𝑥) = 𝑦0 + 𝑛∆𝑦0 + ∆ 𝑦0 + ∆ 𝑦0
2! 3!
𝑛(𝑛 − 1)(𝑛 − 2)(𝑛 − 3) 4
+ ∆ 𝑦0 . . .
4!
Putting in (A), we have
𝑏 𝑛
∫ 𝑓(𝑥) 𝑑𝑥 = ℎ ∫ 𝑓(𝑥0 + 𝑛ℎ) 𝑑𝑛
𝑎 0
𝑛
𝑛(𝑛 − 1) 2 𝑛(𝑛 − 1)(𝑛 − 2) 3
= ℎ ∫ [𝑦0 + 𝑛∆𝑦0 + ∆ 𝑦0 + ∆ 𝑦0 + . . . ] 𝑑𝑛
0 2! 3!

Integrating and simplifying, we get

 𝑛
∫ 𝑓(𝑥0 + 𝑛ℎ) 𝑑𝑛
0
𝑛 𝑛(2𝑛 − 3) 2 𝑛(𝑛 − 2)2 3
= 𝑛ℎ[𝑦0 + ∆𝑦0 + ∆ 𝑦0 + ∆ 𝑦0
2 12 24
+ . . . (𝐵)

This formula is known as Newton’s Cotes Quadrature Formula.

Particular cases of Newton’s Cotes formula n=1, n=2, n=3, n=4, and n=6 are called Trapezoidal
Rule, Simpson’s 1/3 Rule, Simpson’s 3/8 Rule, Boole’s Rule and Weddle’s Rule respectively

Trapezoidal Rule

Putting n=1 in equation (B) and taking curve through (𝑥0 , 𝑦0 ) and (𝑥1 , 𝑦1 ) as a straight line i.e. a
polynomial of first degree so that the differences of order higher than first become zero, we get
𝑥1 =𝑥0 +ℎ
1 ℎ
∫ 𝑓(𝑥) 𝑑𝑥 = ℎ [𝑦0 + ∆𝑦0 ] = [𝑦0 + 𝑦1 ]
𝑥0 2 2

Similarly,
𝑥2
ℎ
∫ 𝑓(𝑥) 𝑑𝑥 = [𝑦 + 𝑦2 ]
𝑥1 2 1
𝑥3
ℎ
∫ 𝑓(𝑥) 𝑑𝑥 = [𝑦 + 𝑦3 ],
𝑥2 2 2

.
𝑥𝑛
ℎ
∫ 𝑓(𝑥) 𝑑𝑥 = [𝑦 + 𝑦𝑛 ]
𝑥𝑛−1 2 𝑛−1

Now adding above integrals, we have

𝑏 𝑥𝑛
ℎ
𝐼 = ∫ 𝑓(𝑥) 𝑑𝑥 = ∫ 𝑓(𝑥) 𝑑𝑥 = [𝑦 + 2(𝑦1 + 𝑦2 + . . . 𝑦𝑛−1 ) + 𝑦𝑛 ]
𝑎 𝑥0 2 0

This is known as general Trapezoidal Rule or Trapezoidal Rule

 1 1
Example 01 Evaluate ∫0 𝑑𝑥 by analytical method and trapezoidal rule taking n=10 and
1+𝑥 2
discuss the results.

Solution
1
1
∫ 𝑑𝑥 = |tan−1 𝑥|10
0 1 + 𝑥2

= tan−1 1 − tan−1 0
𝜋
= − 0 = 0.7854
4
here

𝑏−𝑎 1−0
ℎ= = = 0.1
𝑛 10
X 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Y Yo=1 (Y1)0.99 0.96 0.92 0.86 0.80 0.74 0.67 0.61 0.55 0.50

1
0.1
𝐼 = ∫ 𝑓(𝑥) 𝑑𝑥 = [𝑦 + 2(𝑦1 + 𝑦2 + 𝑦3 + 𝑦4 + 𝑦5 . . . 𝑦9 ) + 𝑦10 ]
0 2 0

= 0.05[1 + 2(0.99 + 0.96 + 0.92 + 0.86 + 0.80 + 0.74 + 0.67 + 0.61 + 0.55) + 0.5]

= 0.785

Simpson’s 1/3 Rule

Substituting n=2 in equation B and taking the curve through (𝑥0 , 𝑦0 ), (𝑥1 , 𝑦1 ) and (𝑥2 , 𝑦2 ) as a
polynomial of degree two (parabola) so that the differences of order higher than second vanish,
we get
𝑥2 =𝑥0 +2ℎ
1 ℎ
∫ 𝑓(𝑥) 𝑑𝑥 = 2ℎ [𝑦0 + ∆𝑦0 + ∆2 𝑦0 ] = [𝑦0 + 4𝑦1 + 𝑦2 ]
𝑥0 6 3

Note:
∆𝑦0 = 𝑦1 − 𝑦0 , ∆2 𝑦0 = 𝑦0 − 2𝑦1 + 𝑦2

Similarly
𝑥4 =𝑥0 +4ℎ
ℎ
∫ 𝑓(𝑥) 𝑑𝑥 = [𝑦 + 4𝑦3 + 𝑦4 ]
𝑥2 3 2
 𝑥6
ℎ
∫ 𝑓(𝑥) 𝑑𝑥 = [𝑦 + 4𝑦5 + 𝑦6 ]
𝑥4 3 4

.
𝑥𝑛
ℎ
∫ 𝑓(𝑥) 𝑑𝑥 = [𝑦 + 4𝑦𝑛−1 + 𝑦𝑛 ]
𝑥𝑛−2 3 𝑛−2

Summing all these results, we obtain

𝑏 𝑥𝑛
𝐼 = ∫ 𝑓(𝑥) 𝑑𝑥 = ∫ 𝑓(𝑥) 𝑑𝑥
𝑎 𝑥0
ℎ
= [𝑦 + 4(𝑦1 + 𝑦3 + . . . 𝑦𝑛−1 ) + 2(𝑦2 + 𝑦4 + . . . 𝑦𝑛−2 ) + 𝑦𝑛 ]
3 0
This is known as general formula of Simpson’s 1/3 Rule. ( n must be an even)
1 1
Evaluate ∫0 𝑑𝑥 by analytical method and simpson’s 1/3 rule taking n=10 and discuss the
1+𝑥 2
results.

Solution
1
1
∫ 𝑑𝑥 = |tan−1 𝑥|10
0 1 + 𝑥2

= tan−1 1 − tan−1 0
𝜋
= − 0 = 0.7854
4
here

𝑏−𝑎 1−0
ℎ= = = 0.1
𝑛 10
X 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Y yo1 (Y1)0.99 0.96 0.92 0.86 0.80 0.74 0.67 0.61 0.55 0.50

1
0.1
𝐼 = ∫ 𝑓(𝑥) 𝑑𝑥 = [𝑦 + 4(𝑦1 + 𝑦3 + 𝑦5 + 𝑦7 + 𝑦9 ) + 2(𝑦2 + 𝑦4 + 𝑦6 + 𝑦8 ) + 𝑦10 ]
0 3 0
 = 0.03[1 + 4(0.99 + 0.92 + 0.80 + 0.67 + 0.55) + 2(0.96 + 0.86 + 0.74 + 0.61) + 0.5]

= 0.7068

Simpson’s 3/8 Rule

Putting n=3 in Newton’s Cote’s formula derived above and taking the curve through (𝑥0 , 𝑦0 ),
(𝑥1 , 𝑦1 ), (𝑥2 , 𝑦2 ) and (𝑥3 , 𝑦3 ) as a polynomial of degree two so that the differences of order
higher than three vanish, we get
𝑥3 =𝑥0 +3ℎ
3 3 1 3ℎ
∫ 𝑓(𝑥) 𝑑𝑥 = 3ℎ [𝑦0 + ∆𝑦0 + ∆2 𝑦0 + ∆3 𝑦0 ] = [𝑦 + 3𝑦1 + 3𝑦2 + 𝑦3 ]
𝑥0 2 2 8 8 0

Similarly
𝑥6
3ℎ
∫ 𝑓(𝑥) 𝑑𝑥 = [𝑦 + 3𝑦4 + 3𝑦5 + 𝑦6 ]
𝑥3 8 3
𝑥9
3ℎ
∫ 𝑓(𝑥) 𝑑𝑥 = [𝑦 + 3𝑦7 + 3𝑦8 + 𝑦9 ]
𝑥6 8 6

.
𝑥𝑛
3
∫ 𝑓(𝑥) 𝑑𝑥 = [𝑦 + 3𝑦𝑛−2 + 3𝑦𝑛−1 + 𝑦𝑛 ]
𝑥𝑛−3 8 𝑛−3

Summing all these results, we obtain

𝑏 𝑥𝑛
𝐼 = ∫ 𝑓(𝑥) 𝑑𝑥 = ∫ 𝑓(𝑥) 𝑑𝑥
𝑎 𝑥0
3ℎ
= [𝑦 + 3(𝑦1 + 𝑦2 + . . . 𝑦𝑛−1 ) + 2(𝑦3 + 𝑦6 + . . . 𝑦𝑛−3 ) + 𝑦𝑛 ]
8 0
This is known as Simpson’s 3/8 Rule. ( n must be multiple of 3)

The following table shows current flow in a circuit at different time interval

T 0 3 6 9 12 15 18 21 24 27
I 0 2 5 10 19 24 39 50 66 86
Find the total charge Q in the circuit when t=27 seconds.