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;