Solving a system of differential equations using ODE45 by switching between two functions depending upon the solution obtained in last iteration.
7 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
Ameer Hamza
2020 年 5 月 4 日
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.
採用された回答
Ameer Hamza
2020 年 5 月 4 日
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 件のコメント
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!