Issue with simpsons rule

I have made the code:

function I = simprule(f, a, b, n)
%SIMPRULE Simpsons rule integration.
% I = SIMPRULE(F, A, B, N) returns Simpsons rule approximation 
% for the integral of f(x) from x=A to x=B, using N subintervals,
% where F is a function handle.

h = (b-a) / n;
x = a:h:b; % an array of length n+1
S = 0;
L = 0; 
for l = 1:2:n %generates the odd number array
S = S + 4*f(x(l));
end
for j = 2:2:n % generates the even number array
L = L + 2*f(x(j));
end 
I = (h/3)*(f(a)+ S + L + f(x(n)));

However this doesn't seem to work when estimating f. Any thoughts on what may be the issue?

[Josh McCrow]

Nevermind, the issue was with l = 1:2:n and j = 2:2:n...

Answer 1

[John D'Errico]

What do you mean about "estimating f"? Simpson's rule has absolutely no purpose in estimating a function, merely an integral.

What does not work? Surely you can tell us why you think it is inadequate? In fact, that code SHOULD come close to working to compute an approximate integral of a data series. It has at least two significant problems.

First, this is a dangerous way to write things if h is a floating point number, since h will not be exactly the value you think it is, with error in the least significant bits.

h = (b-a) / n;
x = a:h:b;

Far better is to use linspace:

x = linspace(a,b,n+1);

so n+1 points, uniformly spanning the interval [a,b], with ABSOLUTE ASSURANCE you get the desired result. (I think I recall they may have fixed colon for this problem some years ago, or maybe I am wrong.)

A much bigger problem is what you do at the first and last point.

Look at your loop:

for l = 1:2:n %generates the odd number array
S = S + 4*f(x(l));
end

What does this do when l=1? By the way, using l as a variable is a REALLY BAD IDEA. Note how easily it is to misread l as 1. Your code will become a death trap of debugging one day. If you insist on the use of L, use an upper case L.

Anyway, when l==1, you get 4*f(a). Is that how Simpson's rule starts? Worse, then you add in f(a) at the end. You do the same thing at the end. So the approximate integral as you wrote it will be

(5*f(a) + 2*f(a+h) + 4*f(a+2*h) + ... + 5*f(b))*h/3

I'm pretty darn sure that is the wrong formula. Yes, I know its been a long time since I took a numerical animalysis class, but I am pretty sure. In fact, the weights for a paneled Simpson's rule would be:

(f(a) + 4*f(a+h) + 2*f(a+2*h) + ... + 4*f(b-h) + f(b))*h/3

They should follow the pattern: [1 4 2 4 2 4 ... 2 4 1]

Finally, you want to test that n (the number of sub-intervals) must be even, as if it is odd, then that last panel ends up in a ditch. Here, by the way, is how I might write Simpson's rule.

if rem(n,2)
  error('The Simpson''s sky is falling. n was odd.')
end
W = mod(0:n,2)*2 + 2;
W([1,end]) = 1;
x = linspace(a,b,n);
h = (b-a)/n;
I = h/3*sum(W.*f(x));