Solving a system of differential equations using ODE45 by switching between two functions depending upon the solution obtained in last iteration.
7 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Ameer Hamza
le 4 Mai 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.
Réponse acceptée
Ameer Hamza
le 4 Mai 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 commentaires
Plus de réponses (0)
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!