Re: Integration of function containing if-test

[Ivan Sutoris]

On Thu, Feb 5, 2009 at 6:17 PM, Torquil Macdonald Sørensen
<address@hidden> wrote:
> 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.
>
> So I get problems when integrating this function e.g. with quadl, if the
> list of values at which it is evaluated contains 0.
>
> Is there an elegant way to define this function nu(x) so that is
> actually returns a 0 for each zero element of its input vector argument
> that is given to it by quadl?
>
> Best regards
> Torquil Sørensen
> Norway

Hi

as Soren has written, for integration, value at one point should not
matter. Anyway, about your question: if "f" accepts vector arguments,
you can just compute f(x) for each element and then set y to zero for
those elements where x==0:

    function y = nu(x)
        y = f(x);
        y(x==0) = 0;
    endfunction

Alternatively, you can always vectorize your function by computing
value at each element individually with "for" cycle:

    function y = nu(x)
        y = zeros(size(x));
        for i=1:numel(x)
            if (x(i) == 0)
                y(i) = 0;
            else
                y(i) = f(x(i));
            endif
        endfor
    endfunction

Regards
Ivan Sutoris