Composite Simpsons 1/3 matlab code

[haz]

hey guys,
So I'm working on a MATLAB function that uses simpsons 1/3 rule to find an integral.
This is what I have done so far, but I'm not 100% confident.. I seem to get double when I use greater segments.
If anyone would be able to have a look at give me some tips that would be very much appreciated.

function [ I ] = simpsons3( func,a,b,n )
%Finds estimate integral from a to b of function using simpson 1/3 rule

% % INPUTS
% % func = function
% % a = lower limit
% % b = upper limit
% % n=number of segments used for integration
% % OUTPUTS
% % I = integral estimate

h = (b-a)/(n-1);
sumE = 0;
for i = a+2*h:h:b-2*h
sumE = sumE + func(i);
end
sumO = 0;
for i = a+h:h:b-h
sumO = sumO + func(i);
end
I = (h/3)*(func(a)+4*sumE+2*sumO+func(b));

 

[DrClaude]

So I'm working on a MATLAB function that uses simpsons 1/3 rule to find an integral.
This is what I have done so far, but I'm not 100% confident.. I seem to get double when I use greater segments.

Always good to test with a function for which you know the analytical result and compare.

function [ I ] = simpsons3( func,a,b,n )
%Finds estimate integral from a to b of function using simpson 1/3 rule

% % INPUTS
% % func = function
% % a = lower limit
% % b = upper limit
% % n=number of segments used for integration
% % OUTPUTS
% % I = integral estimate

h = (b-a)/(n-1);

Minor point, but if n is the number of intervals, then

h = (b-a)/n;

sumE = 0;
for i = a+2*h:h:b-2*h
sumE = sumE + func(i);
end
sumO = 0;
for i = a+h:h:b-h
sumO = sumO + func(i);
end

First, it is a bad idea to use real values for a for loop counter. Use integers and calculate the independent variable from h at each step.

Second, you are doubly counting the interior points. Each sum should add points that are 2h apart, not h.

Third, such an implementation of the Simpson rule only works if the number of intervals is even (or n odd according to your definition of h), otherwise it will return a result with a greater error.

 

[haz]

Thankyou very much for your reply.
I actually was able to make these adjustments before reading your reply but you have given me confidence so thankyou!
I am however confused as to why it is a bad idea to set up my for loop counters as I have done. Would you be able to tell be me why it is bad and give me examples as too when this would fail.

 

[DrClaude]

Because of rounding errors, you can end up with unexpected behaviour at the termination of a for loop. Say you have a for loop that counts to 15, but when you add the step sizes h, you end up with the last values being 15.00000001. This being > 15, the loop will exit without the last point being counted.

 

[alphy]

%simpson 1/3 rule formula
disp(I);

% % The simpson 1/3 rule requires an even number of segments, so I added an
% % if statement to check if n is even. If it's not, I use the trapezoidal
% % rule for the last segment.
if mod(n,2) == 0
I = I + (h/2)*(func(b-h)+func(b));
disp(I);
else
I = I + (h/3)*(func(b-h)+4*func(b-h)+func(b));
disp(I);
end

Dear colleague,

Your code seems to be on the right track! However, I have a few suggestions that might help improve your function.

Firstly, I noticed that you are using a for loop to calculate the sum of even and odd segments separately. This can be simplified by using vectorization. You can create a vector for the even and odd segments, and then use the sum function to calculate their sums. This will not only make your code more efficient, but also easier to read and debug.

Secondly, I would recommend adding some error handling to your code. For example, checking if the input values for a, b, and n are valid (e.g. n must be even, a and b must be within the function's domain), and displaying an error message if they are not. This will make your function more robust and user-friendly.

Lastly, I would suggest using more descriptive variable names. This will make your code easier to understand and maintain in the long run.

Overall, it looks like you have a good understanding of the Simpson's 1/3 rule and how to implement it in MATLAB. Keep up the good work and don't hesitate to reach out for further assistance. Best of luck with your project!