User defined function with vector inputs returns incorrect values?

Hello

I wrote a user-defined function m-file to use the Trapezoidal Rule in order to estimate the integral of a function over a certain interval [a b]. The inputs will include the function (as an anonymous function handle), the interval, and the number of segments to be divided. Output is the integral estimate. Here is the m-file so far:

function [I] = trapezoid(fun,Int,n)
% trapezoid: Implements the trapezoidal rule on a function over a
% specified interval for n number of segments. n can be a scalar or a
% vector.

%input:
%fun = function to be integrated
%Int = [a b] = integration limits
%n = number of segments
%output:
%I = integral estimate

a = Int(:,1);
b = Int(:,2);

I = zeros(1,length(n)); %preallocate solution vector

if length(n) > 1
    for i = 1:length(n)
        I(i) = trapezoid(fun,Int,n(i));
    end
end

x = a;
h = (b - a)./n;
s = fun(a);
en = n-1;
for j = 1:en
    x = x + h;
    s = s + 2.*fun(x);
end
s = s + fun(b);
I = s.*(h./2);

So far, when using the example anonymous function from the command window

fun = @(x) 4*x.^3 - 42*x.^2 + 300

over the interval [0 10], when a scalar value is used for n, the function m-file works and gives verifyable answers. But when the number of segments input is a vector, such as n=[2 3 4 5 10 15 20 30 40 50], the first value of the solution vector is correct, but the others are incorrect (the solution vector that outputs from n=[2 3 4 5 10 15 20 30 40 50] is

 [-250.0000  604.9383  750.0000  728.0000  462.0000  321.6790  245.0000  165.1605  124.3594   99.6704]

, when the correct values should be

[-250, -666.6667, -812.5000, -880, -970, -986.6667, -992.5000, -996.6667, -998.1250, -998.8000]

for the trapezoidal rule).

How can I get my program m-file to accept the vector for the number of segments input AND return the correct results? I have tried vectorizing & adding "." before multiplication/division operators. How can I get the correct output values? Thank you

Answer 1

[Stephen23]

if length(n) > 1
    ...
    return % add this!
end

The reason is simple: you detect when n is non-scalar, and inside that loop you call the function for each n value... but then once you have finished this you simply throw away all of that work by continuing on with the rest of the function (which therefore redefines I, making your loop and self-calling function totally superfluous!). Instead, once you have finished with the loop, you simply need to exit the function, as I showed above using return.

I also recommend replacing length with numel.

Note that a conceptually simpler version of your function would simply use a local or a nested function, e.g.:

function I = trapezoid(fun,a,b,n)
I = zeros(size(n));
for k = 1:numel(n)
    I(k) = subfun(fun,a,b,n(k));
end
end
function I = subfun(fun,a,b,n)
x = a;
h = (b - a)/n;
s = fun(a);
for j = 1:n-1
    x = x + h;
    s = s + 2*fun(x);
end
s = s + fun(b);
I = s*(h/2);
end