How can I use ODE45 continuously?

1 vue (au cours des 30 derniers jours)
Heejun Lee
Heejun Lee le 28 Juil 2021
Modifié(e) : Chunru le 28 Juil 2021
This is the code I initially made.
[T,Y] = ode45(@hw,[0 100],[0.01 10]);
And this is function hw
//
function dy = hw(t,y)
D=0.01;
um=0.2;
Ks=0.05;
X0=0.01;
Sf=10;
dy = zeros(2,1);
dy(1) = -D*y(1) + (um*y(2))./(Ks+y(2))*y(1);
dy(2) = Sf * D - D*y(2) - (um*y(1)*y(2))./(0.4*(Ks+y(2)));
//
For this time I want to use D = 0.01 for t=0~50, and D =0.015 for t=50~100.
And value of X0, Sf for t=50~100 should be the value of X(which is y(1)), S(which is y(2)) at t=50.
How can I make this work?

Réponse acceptée

Chunru
Chunru le 28 Juil 2021
Modifié(e) : Chunru le 28 Juil 2021
Just change D depending on t. Run the solver over the whole time period.
[Update: See @Walter Roberson remarks below. The result looks strange.]
[Update 2: Change the ode options, ie. smaller time step. Three methods gives similar results. It seems that default time step is not good enough for this problem. The branching seems cause no problem as well.]
opts = odeset('MaxStep', 0.001);
[T,Y] = ode45(@hw1,[0 100],[0.01 10], opts);
h(1)=subplot(131); plot(T, Y);
title('D=0.01');
[T,Y] = ode45(@hw,[0 100],[0.01 10], opts);
h(2)=subplot(132); plot(T, Y);
title('D: Branching');
[T1,Y1] = ode45(@hw1,[0 50],[0.01 10], opts);
[T2 Y2] = ode45(@hw2,[50+eps 100],Y1(end,:), opts);
T=[T1;T2]; Y=[Y1; Y2];
h(3)=subplot(133); plot(T, Y)
title('Two ODEs')
%linkaxes(h, 'xy')
function dy = hw(t,y)
%D=0.01;
if t<=50
D = 0.01;
else
D = 0.015;
end
um=0.2;
Ks=0.05;
X0=0.01;
Sf=10;
dy = zeros(2,1);
dy(1) = -D*y(1) + (um*y(2))./(Ks+y(2))*y(1);
dy(2) = Sf * D - D*y(2) - (um*y(1)*y(2))./(0.4*(Ks+y(2)));
end
function dy = hw1(t,y)
D = 0.01;
um=0.2;
Ks=0.05;
X0=0.01;
Sf=10;
dy = zeros(2,1);
dy(1) = -D*y(1) + (um*y(2))./(Ks+y(2))*y(1);
dy(2) = Sf * D - D*y(2) - (um*y(1)*y(2))./(0.4*(Ks+y(2)));
end
function dy = hw2(t,y)
D = 0.015;
um=0.2;
Ks=0.05;
X0=0.01;
Sf=10;
dy = zeros(2,1);
dy(1) = -D*y(1) + (um*y(2))./(Ks+y(2))*y(1);
dy(2) = Sf * D - D*y(2) - (um*y(1)*y(2))./(0.4*(Ks+y(2)));
end
  3 commentaires
Chunru
Chunru le 28 Juil 2021
Good point!!! See the update above.
Heejun Lee
Heejun Lee le 28 Juil 2021
Thanks both of you. The answer and comments are super awesome that I could find out how to do this and what function I have to use and study. Have a nice day!!

Connectez-vous pour commenter.

Plus de réponses (0)

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by