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
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

Sign in to comment.

More Answers (0)

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!