How to Incorporate a Piecewise Function into an ODE System?

The function p(t) in my system of ordinary differential equations is a piecewise function. I have tried two different approaches for incorporating the piecewise function in Matlab, however, none have worked.

For my first approach, I defined p(t) in the ODE file (a nested function):

    function dy  =  ODE_v1(t,y)

        dy(1,1) =  j(t)*y(1)-k(t)*y(1)-0.25*p(t)*y(1);
        dy(2,1) =  j(t)*y(1)-k(t)*y(2)-0.50*p(t)*y(2);
        dy(3,1) =  j(t)*y(2)-k(t)*y(3)-0.75*p(t)*y(3);

        % nested functions

             function y  =  j(t) 
                 y = 0.31*sin(t*2*pi/365);
             end

             function y  =  k(t)
                 y = 0.91*sin(t*4*pi/365);
             end

             function y = p(t)
                 y = ...
                     (0*t             ) .* (t >= 0 & t < 500) + ...
                     (0*t+0.50        ) .* (t == 500) + ...
                     (0*t             ) .* (t > 500 & t < 750) + ...
                     (0*t+0.50        ) .* (t == 750) + ...
                     (0*t             ) .* (t > 750);
             end
    end

For my second approach, I defined p(t) in the ODE solver.

    tf = 1000;
    dt = 0.01;
    t = 0:dt:tf;

    y0 = [300;200;100];

    p = @(t) q(t);         % p(t) is defined here, where q(t) is defined in a separate file q.m

    callODE = @(t,y) ODE_v1(t,y,p);    % In the ODE file: function dy = ODE_v1(t,y,p)
    [t,y] = ode45(callODE,t,y0);

The file q.m:

    function y = q(t)
    y = ...
       (0*t             ) .* (t >= 0 & t < 500) + ...
       (0*t+0.50        ) .* (t == 500) + ...
       (0*t             ) .* (t > 500 & t < 750) + ...
       (0*t+0.50        ) .* (t == 750) + ...
       (0*t             ) .* (t > 750);
    end

For both approaches, Matlab considered p(t) as p(t) = 0. Are there other approaches for incorporating a piecewise function? Thank you!

[Torsten]

Numerical solvers can't handle impulse responses.

Best wishes

Torsten.

Answer 1

[NING LIU]

Is your problem solved? I have the same kind of problem...

[Walter Roberson]

None of the ode*() numeric functions can handle general piecewise functions or impulse functions .

The situations that they can handle are:

1. The condition resolves exactly the same way for every invocation during a single call to the ode*() function because the tspan is manipulated to ensure that the conditions are the same for the entire stretch; or
2. The condition resolves exactly the same way for every invocation during a single call to the ode*() function because event functions are used to stop execution at boundaries; or
3. The conditions are carefully matched so that the first and second derivatives of the code are continuous at each condition boundary. For example, using interp1() with the default 'linear' is not good enough, but using interp1() with 'spline' would be acceptable mathematically

[NING LIU]

Thanks for your comments. I tried to create the piecewise function with Heaviside, and then integrate it into ode(). Fortunately it works!

[Walter Roberson]

Heaviside() is primarily for symbolic work; you would use dsolve() for that, not one of the ode() functions. However, dsolve() is not very likely to find closed form solutions for piecewise functions.

If you defined a function with heaviside or with using conditional computation along the line of sin(x).*(t>2) and you used one of the ode*() functions, then whatever result you got out was very likely wrong. The functions such as ode45() do not always notice that the computation is meaningless, but the mathematical operations they carry out are only valid under the conditions I described above. You might not have received an error message, but you are unlikely to have received the correct answer.