How can I use ODE45 continuously?

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?

Answer 1

[Chunru]

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

[Walter Roberson]

The mathematics behind ode45() say that you will get the wrong answer if have if with different branches taken during any one call to ode45(), except in the case that the second derivative of all equations used in the code are continuous. You are effectively creating an impulse around D, and the first derivative of an impulse is not continuous, so the answer you get from the above will be wrong.

You should instead break up the system into two different ode45() calls, one for [0 50] and the other for [50 100]. The output of the first call should become the input to the second call, but the second call should be arranged to use the different D value.

[Chunru]

Good point!!! See the update above.

[Heejun Lee]

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