Simpson's rule is a method for numerical integration. In other words, it's the numerical approximation of definite integrals.

Simpson's rule is as follows:

sim1

In it,

  • f(x) is called the integrand
  • a = lower limit of integration
  • b = upper limit of integration

Simpson's 1/3 Rule

sim01

As shown in the diagram above, the integrand f(x) is approximated by a second order polynomial; the quadratic interpolant being P(x).

The approximation follows,

sim3

Replacing (b-a)/2 as h, we get,

sim4

As you can see, there is a factor of 1/3 in the above expression. That’s why, it is called Simpson’s 1/3 Rule.

If a function is highly oscillatory or lacks derivatives at certain points, then the above rule may fail to produce accurate results.

A common way to handle this is by using the composite Simpson's rule approach. To do this, break up [a,b] into small subintervals, then apply Simpson's rule to each subinterval. Then, sum the results of each calculation to produce an approximation over the entire integral.

If the interval [a,b] is split up into n subintervals, and n is an even number, the composite Simpson's rule is calculated with the following formula:

sim7

where xj = a+jh for j = 0,1,…,n-1,n with h=(b-a)/n ; in particular, x0 = a and xn = b.

Example in C++:

To approximate the value of the integral given below where n = 8:

sim9
#include<iostream>
#include<cmath>
using namespace std;

float f(float x)
{
	return x*sin(x);	//Define the function f(x)
}

float simpson(float a, float b, int n)
{
	float h, x[n+1], sum = 0;
	int j;
	h = (b-a)/n;
	
	x[0] = a;
	
	for(j=1; j<=n; j++)
	{
		x[j] = a + h*j;
	}
	
	for(j=1; j<=n/2; j++)
	{
		sum += f(x[2*j - 2]) + 4*f(x[2*j - 1]) + f(x[2*j]);
	}
	
	return sum*h/3;
}

int main()
{
	float a,b,n;
	a = 1;		//Enter lower limit a
	b = 4;		//Enter upper limit b
	n = 8;		//Enter step-length n
	if (n%2 == 0)
		cout<<simpson(a,b,n)<<endl;
	else
		cout<<"n should be an even number";
	return 0;
}

Simpson's 3/8 Rule

Simpson's 3/8 rule is similar to Simpson's 1/3 rule, the only difference being that, for the 3/8 rule, the interpolant is a cubic polynomial. Though the 3/8 rule uses one more function value, it is about twice as accurate as the 1/3 rule.

Simpson’s 3/8 rule states :

sim6

Replacing (b-a)/3 as h, we get,

sim5

Simpson’s 3/8 rule for n intervals (n should be a multiple of 3):

sim8

where xj = a+jh for j = 0,1,…,n-1,n with h=(b-a)/n; in particular, x0 = a and xn = b.