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

3 views (last 30 days)
UserCJ on 13 Mar 2022
Commented: Torsten on 17 Mar 2022
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
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
UserCJ on 14 Mar 2022
Thanks Torsten for your suggession. Can you please show that in a sample code for me?

Torsten on 14 Mar 2022
Edited: Torsten on 14 Mar 2022
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
Torsten on 17 Mar 2022
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.