MATLAB Answers

How can I use ODE45 continuously?

1 view (last 30 days)
Heejun Lee
Heejun Lee on 28 Jul 2021
Edited: Chunru on 28 Jul 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?

Accepted Answer

Chunru
Chunru on 28 Jul 2021
Edited: Chunru on 28 Jul 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 Comments
Heejun Lee
Heejun Lee on 28 Jul 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!!

Sign in to comment.

More Answers (0)

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by