Adding delay to a system of differential equations using the Heaviside function

I have a system of differential equations which discribes a certain scenario of drug administration. I want to model the this drug administration using the heaviside function as follows:

I have a seperate function "first_term.m" to create the first part of the equation and another function "second_term.m" to create the second part . And then there is another function "add_RHS.m" to combine both these terms and pass it to ModelRHS(t,x,param).

Here's my add_RHS.m function that defines model equations.

if some condition > 0
    dxdt(j) = dxdt(j) + first_term + heaviside(t-1);
end

if some condition < 0
    dxdt(j) = dxdt(j) + second_term + heaviside(t-1);
end

I have added added an event function to define the event of adding drugs and added the Heaviside term to the above function "add_RHS.m". But my output doen't look like what I wanted it to be. Can someone please help me with this?

editparams; %file that contains parameters
Tend = 100;
Nt = 100;
% Define RHS functions
RHS = @(t,x)RHS(t,x,param); 
%Execution-----------------------------------------------------------------
x0 = [0.004, 0.05, 0.1]; %Initial condition
t = linspace(0,Tend,Nt); %TSPAN
% Options with event function
options = odeset('Events', @heavi);
[t, A] = ode45(RHS, t, x0,options);

% Event function
    function [value, isterminal, direction] = heavi(t, x)
        value = heaviside(t-1);
        isterminal = 1;
        direction = 0;
    end
  

[Torsten]

If you know in advance the time when the last term will be added to your model equations, why don't you simply integrate up to t=tau without this term and afterwards restart the integration with this term included ?

[UserCJ]

Thanks Torsten for your suggession. Can you please show that in a sample code for me?

Answer 1

[Torsten]

  tau = 1.0;
  
  fun1 = @(t,y) y(1);
  tstart = 0.0;
  tend = tau;
  tspan = [tstart,tend];
  y0 = 1;
  [T1,Y1] = ode45(fun1,tspan,y0);  % Integate from 0 up to tau
  
  fun2 = @(t,y) -y(1);
  tstart = tau;
  tend = 2.0;
  tspan = [tstart,tend];
  y0 = Y1(end,1);
  [T2,Y2] = ode45(fun2,tspan,y0);  % Integrate from tau to Tend
  
  T = [T1(1:end-1);T2];
  Y = [Y1(1:end-1,1);Y2];
  plot(T,Y)                        % Plot result from 0 to Tend

[UserCJ]

This makes sense to me! Thanks Torsten! But there are situations where we do not know Tau and we need to model it as a general system as I mentioned in my question. If that's the case, how do I change the code for that?

[UserCJ]

@Torsten Also, do I need to change my equations in fun2 if I want to run it from tau to Tend? I can see you have used -y1 for fun2. I have a system of three equations and I can't think of anyway of changing the system . In my case RHS comes with a function handle.

[Torsten]

Copy RHS to RHS1, RHS2. Implement the equations you use from 0 to tau in RHS1, implement the equations you use from tau to tend in RHS2. Call the integrator at the first place (from 0 to tau) with RHS1, at the second place (from tau to tend) with RHS2.

If tau is unknown, then the "value" in the event function will depend on a solution variable or a combination of solution variables, I guess. Prepare the code as written above. After the event took place, control is given back to the calling program. Then call ode45 anew with RHS2 and starting values at the time of the event (that are returned by ode45).

[UserCJ]

@Torsten Thanks for your comment. Can you please post a sample code for the latter part " After the event took place, control is given back to the calling program. Then call ode45 anew with RHS2 and starting values at the time of the event (that are returned by ode45)". Thanks!

[Torsten]

  function main

    tend = 2.0;
  
    fun1 = @(t,y) y(1);
    tstart = 0.0;
    tend = tend;
    tspan = [tstart,tend];
    y0 = 1;
    options = odeset('Events',@event);
    [T1,Y1,te,ye] = ode45(fun1,tspan,y0,options);  % Integate from 0 up to the point where y(1) = e
  
    fun2 = @(t,y) -y(1);
    tstart = te;
    tend = tend;
    tspan = [tstart,tend];
    y0 = ye(1);
    [T2,Y2] = ode45(fun2,tspan,y0);               % Integrate from the point where y(1) = e to Tend
  
    T = [T1(1:end-1);T2]
    Y = [Y1(1:end-1,1);Y2]
    plot(T,Y)                                     % Plot result from 0 to Tend
  
  end

  function [value, isterminal, direction] = event(t, y)
    value = exp(1) - y(1);
    isterminal = 1;
    direction = 0;
  end

[UserCJ]

Thanks Torsten. I copied RHS to RHS1 and RHS2 and modified your code accordingly. I can see the first part with RHS1 is running and then it returns the following error

Index exceeds the number of array elements. Index must not exceed 1

[Torsten]

Since I don't see your code, I can't help.

But I think the reason behind it is somehow that your system has more than one solution variable y(1).

Of course, this error can be repaired.

[UserCJ]

yes it has three solution variables. Oh thanks! I fixed it! Thanks so much for your help!

[UserCJ]

Hi Torsten, how can I assign this Heaviside function(and the event function) to only one of my variables, say x2? I have a function handle that defines the RHS of the system of equations. I want to assign this delay only to the equation dx2/dt.

[Torsten]

Interrupt the integration with RHS1 as usual at the time when the ODE function for x2 changes.

In RHS2, leave the equations for x1 and x3 as they are and only change the RHS for x3.

Or did I misunderstand your question ?

[UserCJ]

No you have got my point correctly.

"Interrupt the integration with RHS1 as usual at the time when the ODE function for x2 changes."

Does that mean changing the tend in both integrations like below?

    tend = 2.0;
  
    fun1 = @(t,y) y(1);
    tstart = 0.0;
    tend1 = some_tend;
    tspan = [tstart,tend1];
    y0 = 1;
    options = odeset('Events',@event);
    [T1,Y1,te,ye] = ode45(fun1,tspan,y0,options);  % Integate from 0 up to the point where y(1) = e
  
    fun2 = @(t,y) -y(1);
    tstart = te;
    tend2 = tend;
    tspan = [tstart,tend2];
    y0 = ye(1);
    [T2,Y2] = ode45(fun2,tspan,y0);               % Integrate from the point where y(1) = e to Tend
  
    T = [T1(1:end-1);T2]
    Y = [Y1(1:end-1,1);Y2]
    plot(T,Y)  

[Torsten]

tend = 2.0;

RHS1 = @(t,y) [y(1);y(2)]
tstart = 0.0;
tend = tend;
tspan = [tstart,tend];
y0 = [1;1];
options = odeset('Events',@event);
[T1,Y1,te,ye] = ode45(RHS1,tspan,y0,options);  % Integate from 0 up to the point where y(2) = e

RHS2 = @(t,y) [y(1);-y(2)];
tstart = te;
tend = tend;
tspan = [tstart,tend];
y0 = ye;
[T2,Y2] = ode45(RHS2,tspan,y0);               % Integrate from the point where y(2) = e to Tend with a modified equation for y(2)

T = [T1(1:end-1,:);T2]
Y = [Y1(1:end-1,:);Y2]
plot(T,Y)  

end

function [value, isterminal, direction] = event(t, y)
  value = exp(1) - y(2);   % interrupt integration when y(2) reaches e, no matter what y(1) is
  isterminal = 1;
  direction = 0;
end

[UserCJ]

Thanks Torsten! Will try it and check.

[UserCJ]

Hi Torsten, I'm just trying to understand the relationship between y(2) in RHS1 and y(2) in RHS2. Did you randomly put it like -y(2) in RHS2 and in the event function or is that the way that we want to assign our funtion when we define a delay?

Also, I'm creating my model as a combination of a few function files as mentioned in the original equation (first_term.m, second_term.m and add_rhs.m). So if I need to change only one equation to set for the delay, how can I do it? Previously, I directly changed add_rhs.m with a heaviside term(see the original question above), but now I don't need to add the Heaviside function to all three equations. I change only y2 (or x2).

Thanks so much for your help!

[Torsten]

Did you randomly put it like -y(2) in RHS2 and in the event function or is that the way that we want to assign our funtion when we define a delay?

I assumed that the equation for solution component y(2) should change when y(2) is equal to e.

So I set the event to trigger this point (value = y(2) - exp(1)).

Then I assumed that at this point when y(2) reaches e, the equation for y(2) should change to dy(2)/dt = - y(2). That's what I implemented in RHS2.

The equation for y(1) remains the same in both phases of integration.

Also, I'm creating my model as a combination of a few function files as mentioned in the original equation (first_term.m, second_term.m and add_rhs.m). So if I need to change only one equation to set for the delay, how can I do it? Previously, I directly changed add_rhs.m with a heaviside term(see the original question above), but now I don't need to add the Heaviside function to all three equations. I change only y2 (or x2).

I don't know how you arranged your code. Make sure that the correct equations are solved in the different phases of integration. The easiest way to do so is to have two function files RHS1.m and RHS2.m for the two phases of integration that supply the valid equations for these phases. In these two function files, you may call other functions that are equal for both phases so that RHS1 and RHS2 only work as some kind of collector routines.