Solving a system of differential equations using ODE45 by switching between two functions depending upon the solution obtained in last iteration.

Problem: I have a system of two coupled differential equations but i want to put an conditional statement so that depending upon that condition it evaluates the system with different parameter values.

My attempt: I couldn't do the problem and hence I don't have a MWE but here is what i tried. I defined two separate functions with the desired parameters and somehow I wanted to check for the result i obtained in the last iteration and then using an if block to satisfy my condition i wanted to let the solver know which function to evaluate.

I am aware about these event functions and examples like ballode, but i don't know how to formulate it in that way.

timerange= 0:0.5:200;
IC= [0.1,0.1];%initial conditions
threshold=0.5;
% [t,y] =ode45(@(t,y) fn_1(t,y),timerange, IC);

%Not a MWE but more of how i want the system to behave.
%I want that whenever the solution of first DE exceed the threshold,
%It switches to the fn_2

if y(:,1)>=threshold
    [t,y] =ode45(@(t,y) fn_2(t,y),timerange, IC);
else
    [t,y] =ode45(@(t,y) fn_1(t,y),timerange, IC);
end

%---
    
plot(t,y(:,1),'r');
hold on
plot(t,y(:,2),'b');
xlabel('n')
ylabel('I')
grid on

function rk1 =fn_1(t,y)

r      = 0.5; K = 0.5;
eps    = 0.001;  rho= 0.02;
alpha1 = 2.15; c1= 0.025;
gammaI = 0.02; I0= 0.01;

n= y(1);
I= y(2);
    
rk1(1)= r*n*(1- n/K)-eps*n*I;
rk1(2) = I0 + (rho*I*n)/(alpha1+n) -c1*I*n - gammaI*I;
rk1=rk1(:);
end

function rk1 =fn_2(t,y)

r      = 0.5;   K= 0.5;
eps    = 0.001; rho= 0.02;
alpha1 = 2.15;  c1 = 0.025;
I0     = 0.5;     gammaI = 0.2;

n= y(1);
I= y(2);
    
rk1(1)= r*n*(1- n/K)-eps*n*I;
rk1(2) = I0 + (rho*I*n)/(alpha1+n) -c1*I*n - gammaI*I;
rk1=rk1(:);
end

Any help with the code or a hint in the right direction would be very helpful.

Thank you.

[Ameer Hamza]

Can you write your equations in mathematical form and attach it as an image or use latex to write the equations.

[Vira Roy]

Sure Ameer and thanks for your response. It is a simple equation and can be anything because i really want to understand the concept of doing so,

So n is here something like a tumor which grows according to logistic model and it's growth is hampered by the immune. If you will plot the first function fn_1 you will see tumor growth is fast while parameters in fn_2 are such that immune growth is faster, I want to switch between function from fn_1 to fn_2 whenever my tumor grows beyond the threshold. So basically i want the other parameters to take over and control it. And when it reaches a lower threshold say threshold = 0.1 it again switches to fn_1.

condition being that tumor fraction which starts from 0.1 or anything for that matter.

Thanks.

[Ameer Hamza]

So If I understand correctly, you have two sets of parameters.

First

r      = 0.5; K = 0.5;
eps    = 0.001;  rho= 0.02;
alpha1 = 2.15; c1= 0.025;
gammaI = 0.02; I0= 0.01;

Second

r      = 0.5;   K= 0.5;
eps    = 0.001; rho= 0.02;
alpha1 = 2.15;  c1 = 0.025;
I0     = 0.5;     gammaI = 0.2;

and you want to switch between the two based on the value of 'n'. Right? If n > 0.1 the switch to the second set of parameters and if n < 0.1 and switch back to first.

[Vira Roy]

Yes Ameer, thats right. although i want If n > 0.5 then switch to the second set of parameters and if n < 0.1 and switch back to first.

[Ameer Hamza]

And what about the values between 0.1 and 0.5? Should it maintain the last state? Also, does the condition depends on the value of I?

[Vira Roy]

Yes Ameer thats where i got stumped, I need to check for the initial conditions too before switching. So i need to update the initial conditions before i switch. So basically it should keep doing what it was doing but with the new parameters now. No the value of I is not important, only the value of n is important here and it should check for that.

[Ameer Hamza]

Check my code in the answer. It seems that ever with the second set of parameters. The value of 'n' still grows. It never decays.

Answer 1

[Ameer Hamza]

According to the problem you defined in your comment, this is the correct way to implement it using ode45

timerange = 0:0.5:200;
IC= [0.1,0.1];%initial conditions
threshold=0.5;

[t,y] =ode45(@(t,y) fn_1(t,y), timerange, IC);

plot(t,y(:,1),'r');
hold on
plot(t,y(:,2),'b');
xlabel('n')
ylabel('I')
grid on

function rk1 =fn_1(t, y)
n = y(1);
I = y(2);
if n < 0.1
    r      = 0.5; K = 0.5;
    eps    = 0.001;  rho= 0.02;
    alpha1 = 2.15; c1= 0.025;
    gammaI = 0.02; I0= 0.01;
else
    r      = 0.5;   K= 0.5;
    eps    = 0.001; rho= 0.02;
    alpha1 = 2.15;  c1 = 0.025;
    I0     = 0.5;     gammaI = 0.2;
end

rk1(1) = r*n*(1- n/K)-eps*n*I;
rk1(2) = I0 + (rho*I*n)/(alpha1+n) -c1*I*n - gammaI*I;
rk1 = rk1(:);
end

However, even the 2nd set of parameters does not decrease the value of n.

[Vira Roy]

Ameer, I am not worried about the result per se. What i want to undersatnd is the concept and i may be naive for asking this. This if block you have added in the function. Does it check for the value of n everytime like it should do or it checks once and just keep running anyways. Usually in event functions, we check for the certain event and then make the switch or something. Sorry I am not too well versed with the matlab and hence my question.

Is there anyway to check if the switching happens . I played with values by changing it largely and i see the effect. A little confirmation of the concept would be helpful.

[Ameer Hamza]

It checks the value of 'n' at each iteration and chooses the parameters accordingly. Also, you don't need to worry about using new initial conditions at switching because it is continually solving the differential equation and use the value of n and I from the previous time-step.

[Vira Roy]

Well, thanks Ameer. Turns out the simple solution is usually correct and I was making it complicated uselessly. A lot of thanks Ameer. Do you mind if i play with the code a bit and then accept your answer.

[Ameer Hamza]

Sure. No problem!