Composite Simpson's Rule

Question

I have this code for the Composite Simpson's Rule. However, I have been fiddling with it for quite a while and I can't seem to get it to work.

How can I fix this algorithm?

function out = Sc2(func,a,b,N)
% Sc(func,a,b,N)
% This function calculates the integral of func on the interval [a,b]
% using the Composite Simpson's rule with N subintervals.
x=linspace(a,b,N+1);
% Partition [a,b] into N subintervals
fx=func(x);
h=(b-a)/(2*N);
%define for odd and even sums
sum_even = 0;
for i = 1:N-1
   x(i) = a + (2*i-2)*h;
   sum_even = sum_even + func(x(i));
end

sum_odd = 0;

for i = 1:N+1
   x(i) = a + (2*i-1)*h;
   sum_odd = sum_odd + func(x(i));
end
% Define the length of a subinterval
out=(h/3)*(fx(1)+ 2*sum_even + 4*sum_odd +fx(end));
% Apply the composite Simpsons rule
end

Answer 1

Well for one thing, your h definition is wrong. h stands for the step size of each interval you want to estimate. You are unnecessarily dividing by 2. Remove that 2 in your h definition. You also are evaluating your function at the values of n not x. You should probably remove this statement because you end up not using this in the end.

Also, you are summing from 1 to N+1 or from 1 to N-1 for either the odd or even values, This is incorrect. Remember, you are choosing every other value in an odd interval, or even interval, so this should really be looping from 1 to N/2 - 1. To escape figuring out what to multiply i with, just skip this and make your loop go in steps of 2. That's besides the point though.

I would recommend that you don't loop over and add up the values for the odd and even intervals that way. You can easily do that by specifying the odd or even values of x and just applying a sum. I would use the colon operator and specify a step size of 2 to exactly determine which values of x for odd or even you want to apply to the overall sum.

You also are declaring x to be your n-point interval, yet you are overwriting those values in your loops. You actually don't need that x declaration in your code in that case.

As such, here's a modified version of your function with the optimizations I have in mind:

function out = Sc2(func, a, b, N)

h = (b – a) / N; %// Width of each interval
odd = 1 : 2 : n-1; %// Define odd interval 
xodd = a + h*odd; %// Create odd x values 
even = 2 : 2 : n-2; %// Create even interval
xeven = a + h*even; % Create even x values 

%// Return area
out = (h/3)*(func(a) + 4*sum(func(xodd)) + 2*sum(func(xeven))+ func(b));

---

However, if you want to get your code working, you simply have to change your for loop iteration limits as well as your value of h. You also have to remove some lines of code, and change some variable names. Therefore:

function out = Sc2(func,a,b,N)
% Sc(func,a,b,N)
% This function calculates the integral of func on the interval [a,b]
% using the Composite Simpson's rule with N subintervals.

%// Define width of each segment
h = (b - a) / N; %// Change

%//define for odd and even sums
sum_even = 0;
for i = 2 : 2 : N-2 %// Change 
   x = a + i*h; %// Change
   sum_even = sum_even + func(x);
end

sum_odd = 0;

for i = 1 : 2 : N-1 %// Change
   x = a + i*h %// Change
   sum_odd = sum_odd + func(x);
end

%// Output area
out = (h / 3)*(func(a) + 2*sum_even + 4*sum_odd + func(b)); %// Change

end