Simpson's Rule: the Formula and How it Works

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 x_j = a+jh for j = 0,1,…,n-1,n with h=(b-a)/n ; in particular, x₀ = a and x_n = 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 x_j = a+jh for j = 0,1,…,n-1,n with h=(b-a)/n; in particular, x₀ = a and x_n = b.