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

8 views (last 30 days)
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.
  7 Comments
Vira Roy
Vira Roy on 4 May 2020
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
Ameer Hamza on 4 May 2020
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.

Sign in to comment.

Accepted Answer

Ameer Hamza
Ameer Hamza on 4 May 2020
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.
  4 Comments
Vira Roy
Vira Roy on 4 May 2020
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.

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!