3.8 KiB

Raw Blame History

title
Simpson's Rule

Simpson's Rule

In numerical analysis, Simpson's rule is a method for numerical integration (numerical approximation of definite integrals).

Simpson's rule approximates the integration of the form,

where,

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. One common way of handling this problem is by breaking up the interval [a,b] into a number of small subintervals. Simpson's rule is then applied to each subinterval, with the results being summed to produce an approximation for the integral over the entire interval. This sort of approach is termed the composite Simpson's rule.

Suppose that the interval [a,b] is split up into n subintervals, with n being an even number. Then, the composite Simpson's rule is given by,

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:

Approximate the value of the integral given below, taking n = 8.

Implementation of Simpson's 1/3 Rule in C++ is as follows :

#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 is that, here, the interpolant is a cubic polynomial. The 3/8 rule is about twice as accurate as the 1/3 rule, but it uses one more function value.

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) :