Simpson's rule for integration- finding area and centroid, program does not work.

I have written a program for finding the area and centroid for a mathematical function, unfortunately it does not seem to work, and i am getting various errors.

Here is my code in full: main-

    function main()
    global a b h n temp
    a=(-0.5);
    b=0.5;
    n=100;
    h=(b-a)/n;

[area, temp]= simpson(fun, fun2);
area=2*area;
Xc=(1/area)*temp;
disp(area)
disp(Xc)
end

Simpson's method, the action itself-

function [ area, temp] = simpson( fun, fun2 )
%Simpson's 1/3 integration method and calculating required values
global a b h
[se1, se2]=sumev(h);
[so1, so2]=sumod(h);
area= (h/3)*(fun(a)+fun(b)+(2*se1)+(4*so1));
temp= (h/3)*(fun2(a)+fun2(b)+(2*se2)+(4*so2));

end

summing even intervals-

function [ se1, se2] = sumev( h )
%Looping function for summing the values in even intervals for
%each function
global n
for j=0:2:n-2
    v=(2*j+2)*h;
    se1=se1+fun(v);
    se2=se2+fun2(v);
end

summing odd intervals-

function [ so1, so2 ] = sumod( h )
%Looping function for summing the values in odd intervals for
%each function
global n
for j=0:2:n-1
    v=(2*j+1)*h;
    so1=so1+fun(v);
    so2=so2+fun2(v);

end

The first function-

function [ q ] = fun(v)
%The mathematical function to integrate for the area
q=(2*(v^2)+0.2*v+0.08)*sqrt(1-4*(v^2));

end

Second function, used to find the centroid-

function [ e ] = fun2(v)
%The mathematical function to integrate for the area
e=v*(2*(v^2)+0.2*v+0.08)*sqrt(1-4*(v^2));

end

As far as i see, it's a very simple program, i don't see any missed letters or variable names, yet i'm getting these errors when i try to run main-

Error using fun (line 3) Not enough input arguments.

Error in main (line 8) [area, temp]= simpson(fun, fun2);

Any ideas and help debugging would be much appreciated

Answer 1

[John D'Errico]

I think you do not understand functions.

You have defined fun1 and fun2 as taking input arguments. READ THE ERROR MESSAGE!

Error using fun (line 3) Not enough input arguments.

Error in main (line 8) [area, temp]= simpson(fun, fun2);

Then you call them with no inputs.

[area, temp]= simpson(fun, fun2);

What did you expect to happen? Should MATLAB guess what to pass in? Of course not. You need to tell it.

[Rob Dunham]

You are completely correct that i do not understand functions in matlab, as this is only the second time i'm ever using it, which is why i am here asking for help from people much more experienced. I am fully aware that the syntax in this program is probably incredibly flawed.

after changing to

[area, temp]= simpson(fun, fun2);

i am getting: Undefined function or variable "se1".

Error in sumev (line 7) se1=se1+fun(v);

Error in simpson (line 4) [se1, se2]=sumev(h);

Error in main (line 8) [area, temp]= simpson(fun(h), fun2(h));

se1 is defined as a result from a function in simpson, and as an output argument in sumev, so i don't see where the "undefined" error is coming from.

in sumev line 7- i do not see what is the problem with that statement. the function is well defined and returns a single answer for a single value. v is properly defined in sumev.

and of course main, line 8- it doesn't give specification about the error anymore so i don't know what to look for

[John D'Errico]

(Shiver) Global variables all over the place. These make code pure hell to debug, read, follow, understand.

I see that fun and fun2 are identical in their comments, so if one of them is supposed to do something different you might want to change the comments. It took me awhile before I noticed that extra v* in one of them. And why not use descriptive names? Instead of fun and fun2, is there a good reason why you would not use areafun and centroidfun? Or pick some other descriptive names! You code will become better once it becomes more readable, because then you are less likely to make errors in the writing, and you are more able to debug it.

Similarly, when I see lines like this:

[se1, se2]=sumev(h);
[so1, so2]=sumod(h);

my head hurts. se1, se2? so1, so2? Make your code readable, like a good book. You will not be charged for extra characters.

The error from line 8 tells you that fun and fun2 expect an argument. You cannot call them with no arguments. Maybe you needed to call them as

[area, temp]= simpson(fun(h), fun2(h));

fun expects an argument, as does fun2. Thniking that v is a variable that was defined somewhere else in some other function does not tell matlab what to pass into fun.

For example, consider the (built-in) function sum. It is designed to take an argument. Whatever you pass in, it will return the sum. So if I do this:

H = 1:5;
sum(H)

then MATLAB will return the sum of what is in H. If I now ttry this

sum(1:2:15)

it will return a different sum. A function operates on whatever you pass into it! Your own functions operate exactly the same way.

I don't know your intent in this code without trying to re-write your entire code. But, perhaps I might do it like this:

function [farea,fcent] = mysimpsonsrule(f,a,b,n)
% compute the area and centroid of a function f
% given a lower limit a, upper limit b,
% and total number of points n (must be odd)
if a >= b
  error('must have a < b')
end
if (n < 3) || (rem(n,1) ~= 0) || (rem(n,2) ~= 1)
  error('n must be >= 3, an even integer')
end

% interval between points. (n points, so n-1 intervals).
h = (b - a)/(n-1);

% simpson's rule weights
W = 4 - 2*rem(1:n,2);
W([1 end]) = 1;
W = W*h/3;

% what points will we evaluate the function at?
x = linspace(a,b,n);

% and now evaluate the function at those points
fx = f(x);

% compute function area (many ways to write this.)
farea = sum(W.*fx);

% geometric centroid is just the same op, but with
% x in there and then divide by area.
fcent = sum(W.*x.*fx)/farea;

So if I save that as an m-file, then we can see how it works:

fun = @(v) (2*(v.^2)+0.2*v+0.08).*sqrt(1-4*(v.^2));

[fa,fc] = mysimpsonsrule(fun,-.5,.5,11)
fa =
         0.154568591225133
fc =
        0.0598861070612269

[fa,fc] = mysimpsonsrule(fun,-.5,.5,101)
fa =
         0.160816961670664
fc =
        0.0609457541641495

[fa,fc] = mysimpsonsrule(fun,-.5,.5,1001)
fa =
         0.161000663314883
fc =
        0.0609746753979993

It seems to be converging properly. We can verify the above results by a symbolic computation:

syms v
sfun = (2*(v.^2)+0.2*v+0.08).*sqrt(1-4*(v.^2));

FA = int(sfun,-.5, .5)
FA =
(41*pi)/800

vpa(FA)
ans =
0.16100662349647690347121047339307

FC = int(sfun*v,-.5,.5)/FA
FC =
5/82

vpa(FC)
ans =
0.060975609756097560975609756097561

[Rob Dunham]

Thank you for taking so much time to try and help- a lot of advice in here, and i'll be applying it all from now on. Especially in regards to code readability.

A little update- my program was actually ok, all i had to do was define so/se 1 and 2 as zeroes for an initial value, and leave the parenthesis in "simpson" empty (both when calling the function and in the function line itself)

The reason i am using global variables and splitting every simple calculation into a function of its own is that i have to comply with the instructions, though it makes absolutely no sense, and makes things messier.

Anyways- everything works now, code also changed to be much easier to read, notes left in the code with more explanations. Thanks, i appreciate it.

[John D'Errico]

You can get away without the globals. If you want to pass n around, for example, simply pass it in as an argument to the function that needs it. Or use nested functions. A nested function can see the variables from the workspace of the caller.

But definitely head for readability as much as possible. I think about it in terms of the person who might inherit my code one year, if for example, I get run down by the cross-town bus. (Ok, there are some crazed people online that might be intentionally driving that bus, aiming it at me.) Would my successor want to inherit unreadable code? Even if the code works, if it is impossible to modify in the future because nobody can figure out how it works, then my legacy will be for naught.

Internal comments are a great way to make you code more readable. I like a target of one line of comment for every line of code, at least in complex code. Nobody says you need to hit that target all the time of course. And use descriptive variable names and function names. These make your code self-documenting, and always easier to read and follow.

So always look at your code as if you are reading it for the first time, as if you are the person who will be taking it over. If you do, your code will get better.