Using IF condition with ODE

Question

Faidzul Hakim 2020-3-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/513901-using-if-condition-with-ode

关闭： MATLAB Answer Bot 2021-8-20

I would need some help to check on the proper utilization of conditions inside ODE equations.

Below can be found the code. Graphes were managed to be plotted but they were not as I expected them to be.

----------------------------------------------------------------------------------------------------------------------------------------------------------------

%Constants
Iapp=40*60;
K(1)=0.001;                
K(2)=10e-4;
K(3)=1;
K(4)=5;
K(5)=0.8;
% Initial Conditions
y0 = [3.949270942 7e-4 0.000281838];
% Time inverval
t_exp = [0 30 60 120 180 240 300];
% Solver;
options = odeset('AbsTol',1e-6);
[t,y] = ode15s(@manu,t_exp,y0,options,Iapp,K);
% Graphics:
             figure(1)
            %Figure 1-Ca
subplot(311)
plot(t,y(:,1),'r')
legend('Ca2+ theo')
title('Ca2+')
xlabel('treatment time / min')
ylabel('Concentration / M')
hold on
grid
            %Figure 2-CO3
subplot(312)
plot(t,y(:,2),'r')
legend('CO32- theo')
title('CO32-')
xlabel('treatment time / min')
ylabel('Concentration / M')
hold on
grid
            %Figure 3-OH
subplot(313)
plot(t,y(:,3),'r')
legend('OH- theo')
title('OH-')
xlabel('treatment time / min')
ylabel('Concentration / M')
hold on
grid
function dydt = manu(t,y,Iapp,K)
if y(2)>=K(2) 
Ca = -K(1)*y(1);
CO3 = -K(4)*y(2);
OH = 0;
elseif y(2)<=K(2)
Ca = 0;
CO3 = K(3);
OH = Iapp/96500*K(5);
end
dydt =  [Ca; CO3; OH];
end

9 个评论
显示 7更早的评论隐藏 7更早的评论

Torsten 2020-3-30

In this case Walter's suggestion to use Matlab's event option is the way to go.

Or - if the equations remain that simple - you can integrate using pencil and paper.

Faidzul Hakim 2020-3-30

Thanks Torsten for your feedback.

The real equations are a little bit complicated than that. I've picked simple cases to check where the issue came from.

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

Answer 1

Walter Roberson 2020-3-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/513901-using-if-condition-with-ode#answer_422799

You need to use Events in the options to terminate the ode integration, and then restart from that time. None of the ode* functions can deal with discontinuities in any of the derivatives (and not in two more derivatives beyond either.) The ode* functions do not always notice the discontinuities, but their mathematical models rely upon there not being any discontinuities, so even if they do not complain then the results are invalid.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Faidzul Hakim 2020-3-30

Thanks Walter.

Faidzul Hakim 2020-3-31

编辑：Faidzul Hakim 2020-3-31

在 MATLAB Online 中打开

Hi again Walter,

I've tried using the Events options, I seem to be stucked where I would like to redefine the initial values for all my three ODEs. Could you please show me how I can extract my value of dy(1)/dt, dy(2)/dt and dy(3)/dt when my condition of dy(2)/dt achieves 'zero'?

The illustration is given in photo below and the code can be found further down:

K(1)=0.001;                
K(2)=8.5079e-4;
K(3)=1;
K(4)=5;
K(5)=0.02;
% Initial Conditions
y0 = [3.949270942 7e-4 0.000281838];
% Experimental values
t_exp = [0 30 60 120 180 240 300];
% Solver;
options = odeset('AbsTol',1e-6,'Events',@criticEvents);
[t,y,te,ye,ie] = ode15s(@manu1,t_exp,y0,options,K);
if (y(:,2)-K(2))<0
    [t,y,te,ye,ie] = ode15s(@manu1,t_exp,y0,options,K);
    
elseif (y(:,2)-K(2))>=0
    y0=[ye(:,1) ye(:,2) ye(:,3)];
    [t,y,te,ye,ie] = ode15s(@manu2,t_exp,y0,options,K);
end
% Graphics:
             figure(1)
            %Figure 1-M
subplot(311)
plot(t,y(:,1),'r')
legend('M theo')
title('M')
xlabel('time')
ylabel('unit')
hold on
grid
            %Figure 2-N
subplot(312)
plot(t,y(:,2),'r')
legend('N theo')
title('N')
xlabel('time')
ylabel('unit')
hold on
grid
            %Figure 3-P
subplot(313)
plot(t,y(:,3),'r')
legend('P theo')
title('P')
xlabel('time')
ylabel('unit')
hold on
grid
function dydt = manu2(t,y,K)
M = -K(1)*y(1);
N = -K(4)*y(2);
P = 0;
dydt =  [M; N; P];
end
function dydt = manu1(t,y,K)
M = 0;
N = K(3);
P = K(5);
dydt =  [M; N; P];
end
function [critic,isterminal,direction] = criticEvents(t,y)
critic=(y(:,2)-K(2));
isterminal=1;
direction=+1;
end

Answer 2

darova 2020-4-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/513901-using-if-condition-with-ode#answer_423308

在 MATLAB Online 中打开

Maybe you don't need event function for this case. I just add persistent variable to your ode function

function main
clear functions
% your main code
end
function dydt = manu(t,y,Iapp,K)
persistent flag
if isempty(flag)
    flag = false;
end
if y(2)>=K(2) || flag
    % variables
    flag = true;
elseif y(2)<=K(2)
    Ca = 0;
    CO3 = K(3);
    OH = Iapp/96500*K(5);
end
dydt =  [Ca; CO3; OH];
end