Using IF condition with ODE

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

[darova]

How should check this? Where are original equations?

[Faidzul Hakim]

I assigned just simple equations to my ODEs in order to check the proper-using of IF conditions. I can list down the equations below:

Condition 1 (when CO3 value exceeds CO3 limit):

d[Ca]/dt = -k1*[Ca]

d[CO3]/dt = -k2*[CO3]

d[OH]/dt = 0

Condition 2 (when CO3 value does not exceed CO3 limit):

d[Ca]/dt = 0

d[CO3]/dt = k3

d[OH]/dt = k4

k1 to k4 (all K(n) in code) are arbitrary constants.

[darova]

[darova]

[Faidzul Hakim]

You can omit Ca_exp, because they are the experimental data that I would like to compare with.

Iapp/96500*K(5) is constant. That's why the equation in the comment, I call it k4.

[Torsten]

As written, your integration will never be successful. Once you reach condition 1 from condition 2, y(2) starts to decrease again and condition 2 again becomes active. So you will switch from one condition to the other over and over again.

Maybe you meant that once condition 1 is met, you work with the equations for this condition for the rest of the time ?

[Faidzul Hakim]

Yes Torsten. That is what I intented to do. Any suggestion for improvement?

[Torsten]

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]

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]

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.

[Faidzul Hakim]

Hi Walter. That means, the program can read the condition that I assigned for 2nd element of ODEs using IF but since there is a discontinuity, it is not showing the right result?.

I have never used Events in the options previously. If you wouldn't mind, could you please develop more regarding this matter? Is there an easy example you can provide to clarify?

Only if you dont mind. Thank you btw for your feedback.

[Walter Roberson]

I suggest you read through ballode

https://www.mathworks.com/help/matlab/math/ode-event-location.html#bu7wjcg

[Faidzul Hakim]

Thanks Walter.

[Faidzul Hakim]

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]

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

result

[Faidzul Hakim]

I tried it. It doesn't work.