Modelling Simpson's Rule for Numerical Integration in Python

Simpson's rule is a very powerful tool to perform numerical integration. To maximise accuracy while requiring fewer divisions, Simpson's rules calculate the integrals using weighting factors. The trapezoidal rule only considers two points, 𝑥_𝑖 and 𝑥_𝑖+1, to estimate a trapezoid's area, but Simpson's rules use more than two points (or many strips) in each round (refer figure shown below). The values of 𝑓(𝑥) at the various points are updated by weighing various criteria in order to decrease the error

Simpson’s 1/3 Rule

This approach calculates the area of two strips simultaneously, so three different values of "x" are taken into account at each turn. To cover the full domain, the number of strips, n, must be exactly even.

For the first two strips, the Simpson's 1/3 rule formula is as follows −

$$\mathrm{A=\frac{1}{3}h(f(x_0)+4f(x_1)+f(x_2))}$$

The total number of strip pairings will therefore be −

$$\mathrm{I=\frac{1}{3}h((x_0)+4f(x_1)+f(x_2)))+\frac{1}{3}h(f(x_2)+4f(x_3)+f(x_4))+...}$$

$$\mathrm{+\frac{1}{3}h(f(x_{n-4})+4f(x_{n-3})+f(x_{n-2}))}$$

$$\mathrm{+\frac{1}{3} h(f(x_{n-2})+4f(x_{n-1})+f(x_n))}$$

The streamlined form can be written as follows for simple programming −

$$\mathrm{I=\frac{1}{3}h{(f(x_0)+f(x_n))+4(f(x_1)+f(x_3)+...+f(x_{n-3})+f(x_{n-1}))+2(f(x_2)+f(x_4)+...+f(x_{n-4})+f(x_{n-2}))}}$$

Finally, putting the terms in summation form will yield −

$$\mathrm{I=\frac{1}{3}h{(f(a)+f(b))+\sum^{n-1}_{i=1,3,5}4f(x_i)+\sum^{n-2}_{2,4,6}2f(x_i)}}$$

Here, it's important to keep in mind that "n" is even.

Let's suppose we want to perform the following integration −

$$\mathrm{\int^{\frac{\pi}{2}}_0x\:sin(x)dx}$$

Now the programming procedure can be written as −

Importing the module −

# Importing module for array
from numpy import * 

Defining the function whose integration has to be performed −

# Defining the function
def f(x):
   return x*sin(x)

Enter left and right limits (𝑎, 𝑏) of the function and the number of trapezoids in which the area is to be divided i.e. 𝑛 (𝒏 should be even).

# Defining function domain (i.e. Left and Right limits)
a=0
b=pi/2
# Defining number of trapezoids
n=4

Evaluate the width of each trapezoid (ℎ = (𝑏 − 𝑎)/𝑛)

# Evaluating the width of trapezoids
h=(b-a)/n

Create an array of 𝑥. As the number of trapezoids is n so, the number of 𝑥's will be 𝑛 + 1, i.e., the "𝑥" array will be having 𝑛 + 1 elements. This is due to the fact that the index of NumPy array starts with 0 not with 1.

# Creating array of x
# (Note: For n trapezoids (no. of x's/nodes = n+1))
x=linspace(a,b,n+1)

First evaluate $\mathrm{\frac{1}{2}(𝑓(𝑎) + 𝑓(𝑏))}$ and assign the result to the variable "I".

# First term of I
I=0.5*(f(a)+f(b))

Run a loop for the evaluation Equation-1 for 𝑥 = 1 to 𝑥 = 𝑛 (as there are 𝑛 + 1 nodes). The important thing to note here is that when the index j in the main loop is even then $\mathrm{2𝑓(𝑥_𝑖)}$ will be evaluated, else $\mathrm{4𝑓(𝑥_𝑖)}$ term will be evaluated.

# Summing over remaining n-2 trapezoids
for j in range(1,n):
   if j%2==0:
      I=I+2*f(x[j])
   else:
      I=I+4*f(x[j])

Finally, outside the loop multiply "I" by h/3 and assign it to "I" and print the result

I=(h/3)*I
print(f'I = {round(I,6)}')

The full Python program will be −

from numpy import *
def f(x):
   return x*sin(x)
# Left and Right limits
a=0
b=pi/2
# Number of trapezoids
n=4
# Width of trapezoid
h=(b-a)/n
# array of x
# For n trapezoids (no. of nodes=n+1)
x=linspace(a,b,n+1)
# First term of I
I=(f(a)+f(b))
# Summing over remaining n-2 trapezoids
for j in range(1,n):
   if j%2==0:
      I=I+2*f(x[j])
   else:
      I=I+4*f(x[j])
I=(h/3)*I
print(f'I = {round(I,6)}')

Output

For n=4, the program returns I = 0.999591,which is very close to the exact result.

Simpson’s 3/8 Rule

The only difference between it and Simpson's 1/3 rule is that three strips are taken into account, which results in four points being included each time. The first three strips' formula is as follows −

$$\mathrm{A=\frac{3}{8}h(f(x_0)+3f(x_1)+3f(x_2)+f(x_3))}$$

The total of all strip pairings will therefore be.

$$\mathrm{I=\frac{3}{8}[(f(x_0)+3f(x_1)+3f(x_2)+f(x_3))+(f(x_3)+3f(x_4)+3f(x_5)+f(x_6))+(f(x_{n-3})+3f(x_{n-2}+3f(x_n-1)+f(x_n))]}$$

Put the terms in summary form at the end to make programming simpler.

$$\mathrm{I=\frac{3}{8}[(f(a)+f(b))+\sum^{n-2}_{i=1,4,7}3[f(x_i)+f(x_{i+1})]+\sum^{n-3}_{i=3,6,9}2f(x_i)]}$$

The quantity, n, of strips must be a multiple of three. We will use two loops, one for the second term and one for the third term, since it is obvious from the formula above that the second term starts with 1 and the third term starts with 3. By three, the loop counter will increase.

Here we will be taking the same problem which we have taken in the previous 1/3 case −

$$\mathrm{\int^\frac{\pi}{2}_0x\:sin(x)dx}$$

The only difference in the algorithm from the 1/3 case is in the way the Equation-2 is solved. If you observe it carefully then you will understand that index of the second term on RHS starts with 1, whereas the index of the third term starting with 3. So, the loop has to be split in two parts as follows −

# Summing over remaining n-2 trapezoids
j=1
while j<n:
   I=I+3*(f(x[j])+f(x[j+1]))
   j=j+3
j=3
while j<n:
   I=I+2*f(x[j])
   j=j+3
I=(3*h/8)*I

The full Python program is as follows −

from numpy import *
def f(x):
   return x*sin(x)
# Left and Right limits
a=0
b=pi/2

# Number of trapezoids
n=6

# Width of trapezoid
h=(b-a)/n

# array of x
# For n trapezoids (no. of nodes=n+1)
x=linspace(a,b,n+1)

# First term of I
I=(f(a)+f(b))

# Summing over remaining n-2 trapezoids
j=1
while j<n:
   I=I+3*(f(x[j])+f(x[j+1]))
   j=j+3
j=3
while j<n:
   I=I+2*f(x[j])
   j=j+3
I=(3*h/8)*I
print(f'I = {round(I,6)}')

Output

For "n=6", the program returns −

I = 0.999819

Conclusion

In this tutorial, we have learned that how to model Simpsons rule in Python to perform numerical integration. The algorithms for both 1/3 and 3/8 rules were elaborated and programmed.