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:

In it,

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

Simpson's 1/3 Rule

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,

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

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:

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:

``````#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 :

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

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

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