# Simpson’s 1/3rd Rule – Algorithm, Implementation in C With Solved Examples

Numerical Methods & Algorithms / Monday, October 22nd, 2018

# Simpson’s 1/3 rule

In Simpson’s 1/3 rule the interval [a, b] is divided into two equal sub-intervals by the points x0, x1, x2, where h = (b- a)/2 and x2 = x1 +h.

In the previous article we generate Trapezoidal Rule from the general integration formula based on Newton’s forward interpolation formula

$\int_{a}^{b}{f\left( x \right)dx}=nh\left[ {{y}_{0}}+\frac{n}{2}\Delta {{y}_{0}}+\frac{2{{n}^{2}}-3n}{12}{{\Delta }^{2}}{{y}_{0}}+\frac{{{n}^{3}}-4{{n}^{2}}+4n}{24}{{\Delta }^{3}}{{y}_{0}}+…. \right]……..(1)$

Simpson’s 1/3 rule is obtained by putting n = 2 in (1). In this case, the third and higher order differences do not exist.

$\int_{a}^{b}{f\left( x \right)dx}=\int_{{{x}_{0}}}^{{{x}_{n}}}{f\left( x \right)dx}=2h\left[ {{y}_{0}}+\Delta {{y}_{0}}+\frac{1}{6}{{\Delta }^{2}}{{y}_{0}} \right]$

$\Rightarrow \int_{a}^{b}{f\left( x \right)dx}=2h\left[ {{y}_{0}}+\left( {{y}_{1}}-{{y}_{0}} \right)+\frac{1}{6}\left( {{y}_{2}}-2{{y}_{1}}+{{y}_{0}} \right) \right]$

$\therefore \int_{a}^{b}{f\left( x \right)dx}=\frac{h}{3}\left[ {{y}_{0}}+4{{y}_{1}}+{{y}_{2}} \right]$

The above rule is known as Simpson’s 1/3 rule or simply Simpson’s rule.

## Composite Simpson’s 1/3 rule

Let the interval [a, b] be divided into n (an even number) equal sub-intervals by the points x0, x1, x2, …, xn, where xi = x0 + ih, i = 1, 2, …, n. Then

$\int_{a}^{b}{f\left( x \right)dx}=\int_{{{x}_{0}}}^{{{x}_{2}}}{f\left( x \right)dx}+\int_{{{x}_{2}}}^{{{x}_{4}}}{f\left( x \right)dx}+….+\int_{{{x}_{n-2}}}^{{{x}_{n}}}{f\left( x \right)dx}$

$=\frac{h}{3}\left[ {{y}_{0}}+4{{y}_{1}}+{{y}_{2}} \right]+\frac{h}{3}\left[ {{y}_{2}}+4{{y}_{3}}+{{y}_{4}} \right]+….+\frac{h}{3}\left[ {{y}_{n-2}}+4{{y}_{n-1}}+{{y}_{n}} \right]$

$=\frac{h}{3}\left[ {{y}_{0}}+4\left( {{y}_{1}}+{{y}_{3}}+…+{{y}_{n-1}} \right)+2\left( {{y}_{2}}+{{y}_{4}}+…+{{y}_{n-2}} \right)+{{y}_{n}} \right]$

This formula is known as Simpson’s 1/3 composite rule for numerical integration.

## Error in Simpson’s 1/3 rule

The error committed in the formula is given by

${{E}_{1}}=\int_{{{x}_{0}}}^{{{x}_{2}}}{f\left( x \right)dx}-\frac{h}{3}\left[ {{y}_{0}}+4{{y}_{1}}+{{y}_{2}} \right]\simeq -\frac{1}{90}{{h}^{5}}{{f}^{iv}}\left( {{x}_{0}} \right)$

which is the error in the interval [x0, x2].

The total error committed in composite Simpson’s 1/3 rd rule is given by

$E\simeq -\frac{1}{90}{{h}^{5}}\frac{n}{2}{{f}^{iv}}\left( \xi \right)=-\frac{b-a}{180}{{h}^{4}}{{f}^{iv}}\left( \xi \right),{{x}_{0}}<\xi <{{x}_{n}},\left( \because b=a+nh \right)$

where ξ is the point having largest value of the fourth order derivatives, and since the number of sub-intervals is n/2.

## Algorithm of Simpson’s 1/3 rule

Step 1. Start;
Step 2. Input function f(x);
Step 3. Read a,b,n; // the lower and upper limits and number of sub-intervals
Step 4. Compute h=(b-a)/n;
Step 5. Sum = [f(a)-f(a+nh)];
Step 6. for i=1 to n-1 step 2 do
Compute sum = sum + 4*f(a+ih)+2*f(a+(i+1)h);
endfor;
Step 7. Compute result = sum * h/3;
Step 8. Print result;
Step 9. Stop;

## Simpson’s 1/3 rule implementation in C

/* Program Simpson’s 1/3 rule
Program to find the value of integration of a function
f(x) using Simpson’s 1/3 rule. Here we assume that
f(x) = x^3. */

#include<stdio.h>

void main()
{
float f(float);
float a,b,h,sum;
int i,n;
printf("Enter the values of a,b: ");
scanf("%f%f",&a,&b);
printf("Enter the value of n: ");
scanf("%d",&n);
if(n%2!=0)
{
printf("\nNumber of subdivision should be even");
exit(0);
}
h=(b-a)/n;
sum = f(a)-f(a+n*h);
for(i=1;i<n;i++)
sum += 4*f(a+ih)+2*f(a+(i+1)*h);
sum *= h/3;
printf("\nValue of the integration is %f",sum);
}

float f(float x)
{
return(x*x*x);
}

Output

Enter the values of a, b: 0 1

Enter the value of n: 100

Value of the integration is 3.750000

 Example 01

Find the value of

$\int_{1}^{5}{{{\log }_{10}}xdx}$

Taking 8 sub-intervals, correct up to four decimal places by Simpson’s 1/3rd rule.

Solution:

$Let~~f\left( x \right)={{\log }_{10}}x$

$Here~~{{x}_{0}}=1,~~{{x}_{n}}=5,~~n=8$

$\therefore h=\frac{{{x}_{n}}-{{x}_{0}}}{n}=\frac{5-1}{8}=0.5$

The tabulated values of f(x) for different values of x are given below,

 x 1 1.5 2 2.5 3 3.5 4 4.5 5 f(x) 0 0.17609 0.30103 0.39794 0.47712 0.54407 0.60206 0.65321 0.69897

By Simpson’s 1/3rd rule, we have,

$\int_{a}^{b}{f\left( x \right)dx}=\frac{h}{3}\left[ {{y}_{0}}+4\left( {{y}_{1}}+{{y}_{3}}+{{y}_{5}}+{{y}_{n7}} \right)+2\left( {{y}_{2}}+{{y}_{4}}+{{y}_{6}} \right)+{{y}_{8}} \right]$

$\therefore \int_{1}^{5}{{{\log }_{10}}xdx}=\frac{0.5}{3}[0+4\left( 0.17609+0.39794+0.54407+0.65321 \right)$

$+2\left( 0.30103+0.47712+0.60206 \right)+0.69897]$

$\therefore \int_{1}^{5}{{{\log }_{10}}xdx}=1.75744\simeq 1.7574$

Correct upto four decimal places.

 Example 02

Find the value of

$\int_{1.2}^{1.6}{\left( x+\frac{1}{x} \right)dx}$

Taking 4 sub-intervals, correct up to four significant figures.

Solution:

$Let~~f\left( x \right)=\left( x+\frac{1}{x} \right)$

$Here~~{{x}_{0}}=1.2,~~{{x}_{n}}=1.6,~~n=4$

$\therefore h=\frac{1.6-1.2}{4}=0.1$

The tabulated values of f(x) for different values of x are given below:

 x 1.2 1.3 1.4 1.5 1.6 f(x) 2.03333 2.06923 2.11429 2.16667 2.225

By Simpson’s 1/3rd rule, we have,

$\int_{a}^{b}{f\left( x \right)dx}=\frac{h}{3}\left[ {{y}_{0}}+4\left( {{y}_{1}}+{{y}_{3}} \right)+2\left( {{y}_{2}} \right)+{{y}_{4}} \right]$

$\int_{1.2}^{1.6}{\left( x+\frac{1}{x} \right)dx}=\frac{0.1}{3}\left[ 2.033333+4\left( 2.069231+2.166667 \right)+2\left( 2.114284 \right)+2.225 \right]$

$\therefore \int_{1.2}^{1.6}{\left( x+\frac{1}{x} \right)dx}=0.847683\simeq 0.8477$