Integration of function containing if-test

[Torquil Macdonald Sørensen]

Søren Hauberg wrote:
> tor, 05 02 2009 kl. 18:17 +0100, skrev Torquil Macdonald Sørensen:
>> Hi!
>>
>> I am having some problems with integration of the following function 
>> (this is really a simplification of my function, but the problem I 
>> illustrate is the same).
>>
>> I have a function nu(x) that is defined using an if-test. It has a 
>> special value (0) at x = 0, and is defined in terms of an ordinary 
>> function f(x) everywhere else. Therefore I have:
>>
>> function y = nu(x)
>> 	if (x == 0)
>> 		y = 0;
>> 	else
>> 		y = f(x);
>> 	endif
>> endfunction
>>
>> When I pass it a 0 it returns a zero. But when I pass it a vector [0 1] 
>> it doesn't work, because the if-test is apparently not executed on each 
>> element of the vector individually, but on the vector x as a whole. 
>> Since x = [0 1] is not equal to 0, the first part of the if-test is 
>> never entered, and therefore f(x) is evaluated at both 0 and 1, thereby 
>> returning the wrong result. In my case it returns [ NaN 1], since f(0) = 
>> Nan and f(1) = 1.
> 
> I'm not quite I understand your problem, but:
> 
> 1) Can you just do
> 
>    function y = nu(x)
>      if any (x == 0)
>        y = 0;
>      else
>        y = f(x);
>      endif
>    endfunction
> 
> ? Notice the 'any'.

Thanks Søren, actually this method didn't work for me as is stands, 
because I need it to return a vector even if one of the input values are 
0. This returns "0" of at least one element of x is 0.

> 2) Since you only have a finite number of points where 'x' is zero, then
> I would have thought you could just ignore them when doing integration.
> Can't you just integrate 'f (x)' directly?

Yes, I decided to do this and it works fine. I had looked at the help 
for quadl, not quad, so I didn't think of the possibility to specify a 
point where the function is "singular", as described in "help quad". 
Actually the function I want to integrate is well-defined at x = 0, but 
its defining expression in octave is given  as x*x*g(x) where g(0) = 
NaN. And since 0*NaN is also NaN, I need to treat it as a special case.

Thanks again
Torquil Sørensen