Simpson's Rule

Simpson's rule is a repetitive process for estimating the definite integral of a function (the area underneath the curve). Most calculators use a modification of this rule. As mentioned earlier, it is repetitive -- and it must be repeated for an even number of iterations. It is generally more accurate than the othe repetitive approximation of the definite integral -- the Trapezoidal Rule.

What is this wonderful rule? Here it is, using Integ for the integral sign (the formatting might be slighly messed up, due to the difficulty of putting this in pASCII].):

Integ(a, b, f(x) dx) ~=
b-a [f(x0)+4f(x1)+2(fx2)+4(fx3)...+4f(xn-a)+f(xn)]
3n

The idea behind all this is that we use portions of the graph in the form y=cx^2+dx+e for constants c d and e.

The error from Simposon's rule can be found with the following formula:

M(b-a)^5
180n^4

M here is a positive real number such that |f⁽⁴⁾(x)| <= M for every x in [a,b].

Simpson's rule is still very useful for finding definite integrals using calculators. Many programs that will compute a definite integral using Simpson's rule are available from ticalc.org or hpcalc.org.

It is interesting to note that Simpson's Rule is the best polynomial approximation method for definite integrals in the general case. It might be intuitive to assume that higher-degree polynomial approximations to an integral would yield stricter error bounds, but this is not the case. I can personally vouch for this fact, as I was forced in high school calculus through the painful, ink-intensive process of deriving a general cubic approximation formula and its (more lax) error bound solely to stress the point that it just doesn't get any better than Simpson's rule.

---

To clarify: the polynomial approximations in question are nth-degree fits to (n+1)-point intervals, not general polynomials. The idea is that you take n+1 points evenly spaced out on the curve and fit a degree n polynomial to exactly those points. The whole point of this construction is that it is entirely algebraic. While a high-order Taylor approximation may often yeild better results, calculus is required to construct the approximating polynomial in the first place.

Numerical method for integrating an arbitrary function real->real. This is like the well-known trapezium rule, except that the function is approximated in each subinterval as a quadratic, rather than a straight line. This makes it better. Symbolically expressed as:

_a∫^b f(x) dx ≅ (h/3)_i=0∑^m-1(f(a+2ih)+4f(a+(2i+1)h)+f(a+(2i+2)h))

Where m is the number of subintervals that the interval b-a is divided into. h is half the size of those sub-intervals, ie 2h=(b-a)/m

The constant 1/3 multiplier falls out of the integration of a quadratic. Don't question it (or consult a derivation of this rule from first principles).

The expression that is being summed is the integral of the quadratic used to approximate the function in the subinterval i (Read that again if necessary. A picture would help. I recommend that you consult http://metric.ma.ic.ac.uk/integration/techniques/definite/numerical-methods/simpsons-rule/, the page that helped me to understand this concept. In addition, it has a quick derivation of this sucker).

Maths buffs will note that this has nothing whatsoever to do with Pascal's Triangle, despite having a symmetric coefficient pattern.

In case you are curious, the polynomial is chosen to agree with the original function at the start of the interval, the middle of the interval, and the end of the interval.

Shonky matlab code follows. It has been tested with GNU octave 2.1.36 and matlab 6. If you do insist on running it (In violation of the licence: don't), you will need to place it in a file, then invoke the function.

---

% This may or may not calculate simpson's rule.
% You are not licensed to use this code for anything at all
% If you do use it, I disclaim all liability. I do not warrant that this code does anything at all
% If it does something when you run it, it's your own fault for submitting this to whatever 
% interpreter you used.
function retval = simpson(func, m, a, b)
    bmina = (b-a);
    h = bmina/(2*m);
    hover3 = bmina/(6*m);
    retval = 0;
    for i = 0:(m-1)
        atwoih = a+2*i*h;
        % These are calculated by addition to avoid loss of precision
        atwoih1 = atwoih+h;
        atwoih2 = atwoih1+h;
        retval = retval + feval(func, atwoih) + 4*feval(func, atwoih1) + feval(func, atwoih2);
    end
    %This line was missing in previous versions. It is rather crucial.
    retval = retval*hover3;

---