How can I use ODE45 continuously?

1 次查看(过去 30 天)
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?

采纳的回答

Chunru
Chunru 2021-7-28
编辑:Chunru 2021-7-28
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 个评论
Chunru
Chunru 2021-7-28
Good point!!! See the update above.
Heejun Lee
Heejun Lee 2021-7-28
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!!

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by