MATLAB: Can someone tell me why the error is so large for the composite simpsons rule

composite simpson's approximationcomposite simpson's ruleerror bounds

My assignment was to code a composite simpsons rule where the exact value of the integral of xsin(x) from 0 to 1 is 0.301168678 but I keep getting 0.862149054988026 as the approximation and that is not within my error bound. Am I entering the composite simpsons rule wrong???

if true
  % EV = sin(1)-cos(1); %exact value
a=0; %starting point
b=1; %endpoint
f=@(x) x.*sin(x);
%d=f(a);
for j=0:6
  n=10^(j);
  h=(b-a)/n;
  Err_CS= 1/(36*(n)^4); %error bound
  CS=(1/3)*(f(a) + 2.*sum(f((a + 2*h):2*h:(b - 2*h)))+4.*sum(f((a + h):2*h:(b -h)) +f(b)))*h; %composite simpsons rule
fprintf('%8.2e  %14.8e  %14.8e  %14.8e\n', n, CS, abs(CS-EV), Err_CS) 
end
end

Best Answer

Hi Briyahna,

Yes the Simpson's rule expression is wrong, but only in the typo sense of having a misplaced parenthesis. Instead of

CS=(1/3)*(f(a) + 2.*sum(f((a + 2*h):2*h:(b - 2*h)))+4.*sum(f((a + h):2*h:(b -h))  +f(b)))*h;
it should be
CS=(1/3)*(f(a) + 2.*sum(f((a + 2*h):2*h:(b - 2*h)))+4.*sum(f((a + h):2*h:(b -h))) +f(b ))*h;

Even better would be to have (h/3) in front as per the usual Simpson expression rather than the h factor hiding all the way at the end.

Related Solutions

MATLAB: Issue with simpsons rule

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));

MATLAB: Simpson’s Rule

I think you made an error on the line

   for j = 1:n-1

It should be

   for j = 1:k-1

where k = n/2. As it stands now, the x2 reaches a value of a+(2*n-1)*h which is far beyond the range from 0 to 1.

Best Answer

Related Solutions

MATLAB: Issue with simpsons rule

MATLAB: Simpson’s Rule

Related Question