MATLAB: Simpson’s Rule

homeworksimpsons rule

I coded Simpson's Rule, but I'm not sure if it's right.
f = @(x) exp(-x.^2);
true = integral(f, 0, 1);
%simpson's rule
n = 128;
k= n/2;
a = 0; b = 1;
h = (b-a)/n;
x = a + h;
sum1 = (h/3)*(f(b)+ f(a));
sum1 =  sum1 + (4*h/3)*f(x);
for j = 1:n-1
    x1 = a+(2*j)*h;
    x2 = a + (2*j+1)*h;
    sum1 = sum1 + (2*h/3)*f(x1)+(4*h/3)*f(x2);
end
error1 = abs(true-sum1);
fprintf(' Simpsons %d : %0.8f \n', n, sum1);
fprintf('Simpsons Error: %0.8f \n', error1);

Best Answer

  • 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.
    • Related Question