How can I use ODE45 continuously?

Question

Heejun Lee 2021-7-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/887384-how-can-i-use-ode45-continuously

编辑： Chunru 2021-7-28

采纳的回答： Chunru

在 MATLAB Online 中打开

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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Chunru 2021-7-28

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/887384-how-can-i-use-ode45-continuously#answer_755359

编辑：Chunru 2021-7-28

在 MATLAB Online 中打开

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 个评论
显示 1更早的评论隐藏 1更早的评论

Chunru 2021-7-28

Good point!!! See the update above.

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!!

请先登录，再进行评论。