Why Matlab integral function does not work properly?

8 views (last 30 days)
Rob
Rob on 23 Aug 2013
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!!!!

Answers (1)

Mike Hosea
Mike Hosea on 24 Aug 2013
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

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!