Why Matlab integral function does not work properly?

Hello,

I have defined function myfun as: function f=myfun(x) if x>1 f=1./(x.*x); else f=4./sqrt(x); end end

This function is continuous and has continuous first derivative. When I run

integral(@myfun,0,1)+integral(@myfun,1,inf).

ans =

    9.0000

As expected. However, when I enter >> integral(@myfun,0,inf)

ans =

   72.0123.

Which is surprising!!!!

Answer 1

[Mike Hosea]

The problem is that "if" is not vectorized the way you would expect. I think

if x > 0

means

if all(x(:) > 0)

So this was probably going to go into the else clause a lot, even when it shouldn't. Here are 4 ways of doing it right.

function doit
integral(@myfun0,0,inf)
integral(@myfun1,0,inf)
integral(@myfun2,0,inf)
integral(@myfun3,0,inf)

function f = myfun_scalar(x)
% Only pass scalar x to this function.
if x > 1 
    f=1/(x*x); 
else
    f=4/sqrt(x);
end

function f = myfun0(x)
% Vectorize with ARRAYFUN.
f = arrayfun(@myfun_scalar,x);

function f = myfun1(x)
% Use logical indexing.
f = zeros(size(x));
f(x > 1) = 1./(x(x > 1).^2);
f(~(x > 1)) = 4./sqrt(x(~(x > 1)));

function f = myfun2(x)
% Use FIND.
f = zeros(size(x));
idx = find(x > 1);
f(idx) = 1./(x(idx).*x(idx));
idx = find(~(x > 1));
f(idx) = 4./sqrt(x(idx));

function f = myfun3(x)
% Use a loop.
f = zeros(size(x));
for k = 1:numel(x)
    if x(k) > 1
        f(k) = 1/(x(k)*x(k));
    else
        f(k) = 4/sqrt(x(k));
    end
end