Simpson's rule

In numerical analysis, Simpson's rule (named after Thomas Simpson) is a way to get an approximation of an integral:

\int_{a}^{b} f(x)\, dx.

Basics

Simpson's rule works by approximating

f(x)

by the quadratic polynomial

P(x)

which takes the same values as

f(x)

at a, b, and the midpoint m=(a+b)/2. One can use Lagrange polynomial interpolation to find an expression for this polynomial,

P(x)=f(a)\frac{(x-m)(x-b)}{(a-m)(a-b)}+

f(m)\frac{(x-a)(x-b)}{(m-a)(m-b)}+ f(b)\frac{(x-a)(x-m)}{(b-a)(b-m)} .

Simpson's rule then follows by an easy (albeit tedious) calculation:

\int_{a}^{b} f(x) \, dx\approx \int_{a}^{b} P(x) \, dx =\frac{b-a}{6}\left + 4f\left(\frac{a+b}{2}\right)+f(b)\right .

The error in approximating an integral by Simpson's rule is

-\frac{h^5}{90}f^{(4)}(\xi),

with $h=(b-a)/2$ and $\xi$ some number between $a$ and $b.$

Composite Simpson's rule

We see that Simpson's rule provides an adequate approximation if the interval of integration $b$ is small, which does not happen most of the time. The obvious solution is to split the interval of integration in small subintervals, apply Simpson's rule on each subinterval, and add up the results. In this way one obtains the composite Simpson's rule

\int_a^b f(x) \, dx\approx

\frac{h}{3}\bigg[f(x_0)+2\sum_{j=1}^{n/2-1}f(x_{2j})+ 4\sum_{j=1}^{n/2}f(x_{2j-1})+f(x_n) \bigg],

where $n$ is the number of subintervals in which one splits $b$ with $n$ an even number, $h=(b-a)/n$ is the length of each subinterval, and $x_i=a+ih$ for $i=0, 1, ..., n-1, n$ , in particular, $x_0=a$ and $x_n=b.$ Alternatively, the above can be written as:

\int_a^b f(x) \, dx\approx

\frac{h}{3}\bigg^*.

The maximum error associated with the composite Simpson's rule can be found using the following formula:

-\frac{h^4}{180}(b-a)f^{(4)}(\xi),

where $h$ is the "step length", given by $h=(b-a)/n.$

Python implementation of Simpson's rule

Here is an implementation of Simpson's rule in Python.


def simpson_rule(f, a, b):
  "Approximate the definite integral of f from a to b by Simpson's rule."
  c = (a + b) / 2.0
  h3 = abs(b - a) / 6.0
  return h3 * (f(a) + 4.0*f(c) + f(b))
 Calculates integral of sin(x) from 0 to 1
from math import sin
print simpson_rule(sin, 0, 1)

Integrating sin x from 0 to 1 with this code gives 0.4598622... whereas the true value is 1 − cos 1 = 0.45969769413... .

Matlab implementation of Composite Simpson's rule

%Define the function to integrate (using anonymous functions) f = @(x) x^2; %Set the interval to integrate a = -1; b = 1; %set the number of panels to compute and their length n = 100; h = (b-a)/n; %split up the interval into subintervals x = ^*; %note that matlab matrices are indexed starting at 1 sum = f(x(1)); for i=2:2:n sum = sum + 4*f(x(i)); end; for i=3:2:n-1 sum = sum + 2*f(x(i)); end; %prints out the result f_integrated = (h/3)*(sum + f(x(n+1)))

References

External links

Numerical integration