]>

## Simpson's Method for numerical integration

Numerical integration uses the concept that a definite integral of a function between two limits is simply the “area under the curve”. The only complication is estimating the integrand function. Of course if you knew the function that was the integral of the integrand function you would not have to do a numerical integral; you could do it analytically. Simpson’s Method breaks the interval between the limits into a number of subintervals that you specify and then approximates the function in each subinterval as a quadratic polynomial that passes through the two end points and the mid-point of the subinterval. You know how to integrate a quadratic polynomial, so then you just add up the integrals of all the quadratic polynomials that approximate the integrands in each subinterval to get the total integral. For smooth well behaved integrands, if you choose enough narrow subintervals, then this is a good approximation to the area under the curve. This is shown graphically in the figure below for the function

$f\left(x\right)=-{\left(x-1\right)}^{3}+2{\left(x-1\right)}^{2}+2$ between the limits $x=0$ and $x=3$ .

If you fit the three points in each subinterval to a quadratic polynomial then integrate each of these and then add them together you get the formula

$\underset{a}{\overset{b}{\int }}f\left(x\right)dx\approx \frac{h}{3}\left[f\left({x}_{0}\right)+4f\left({x}_{1}\right)+2f\left({x}_{2}\right)+4f\left({x}_{3}\right)+\cdots 2f\left({x}_{N-2}\right)+4f\left({x}_{N-1}\right)+f\left({x}_{N}\right)\right]$

where you divide the interval $\left[a,b\right]$ into $N$ subintervals and $N$ is even. You can see the pattern of coefficients 1, 4, 2, 4, 2, 4, 2, …, 2, 4, 1 and ${x}_{0}=a$ and ${x}_{N}=b$ . The subinterval width is $h=\frac{b-a}{N}$ . The beauty of this formula is that the approximation error associated with using finite subinterval width $h$ varies as ${h}^{4}$ Thus if you halve the value of ${h}^{}$ (i.e. choose twice as many subintervals) then the error is reduced by ${2}^{4}$ or 16 over the error in the original estimate. Simpson’s Method is called a fourth order accurate numerical algorithm because of this property.

Simpson’s Method is an example of the type of method used for numerically estimating definite integrals by the int(f,x=a..b,numeric) command in Maple and quad function in Matlab. While these do not use Simpson’s Method specifically, they use more sophisticated methods that are similar in principle.