Simpson's Rule: The Formula And How It Works - FreeCodeCamp

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.