Info

This question is closed. Reopen it to edit or answer.

Using IF condition with ODE

1 view (last 30 days)
Faidzul Hakim
Faidzul Hakim on 30 Mar 2020
Closed: MATLAB Answer Bot on 20 Aug 2021
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 Comments
Torsten
Torsten on 30 Mar 2020
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
Faidzul Hakim on 30 Mar 2020
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.

Answers (2)

Walter Roberson
Walter Roberson on 30 Mar 2020
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 Comments
Faidzul Hakim
Faidzul Hakim on 30 Mar 2020
Thanks Walter.
Faidzul Hakim
Faidzul Hakim on 31 Mar 2020
Edited: Faidzul Hakim on 31 Mar 2020
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

darova
darova on 1 Apr 2020
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

Community Treasure Hunt

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

Start Hunting!