MATLAB: Switch between two sets of coupled ODE’s systems using Event

coupled odeeventodeode systemode45odeset

Hi everyone!
I have two sets of coupled ODE that I want to switch on and off depending on the reached value when solving.
My system looks like this:
%diff equations
%u(1)=concentration 

%u(2)=density

pak1=@(t,u) [Ka.*(c_sat_179-u(1))-(q.*u(2)); (u(2).*mu_max.*u(1))./(K_s+u(1))];
pak2=@(t,u) [-(q.*u(2)); (u(2).*mu_max.*u(1))./(K_s+u(1))];
I want to start by solving pak1, the initial concentration is c_min. When a certain calue c_max is achieved, I want to turn off pak1 and start solving pak2, when c_min is achieved again I want to turn off pak2 and switch on pak 1 and so on. I could manage to swith between single (not coupled) ODE's but I'm struggling now, I attach the entire code.
clear all
close all
clc
%Values-------------------------------------
p_tot=150.*10000; %Pa

p_par=p_tot.*0.179; %Pa
H=2938.4; %Pa m^3 mol^-1
c_sat_100=(p_tot./H)./1000; %mol/L





c_sat_179=(p_par./H)./1000; %mol/L
c_max=833./1e6; %mol/L
Ka=1.2;
tspan = [0 100];
tstart = tspan(1);
t=tstart;
tend = tspan(end);
c_min = 46./1e6; %mol/L
u0=[c_min;1/1000]; %initial conditions for ODE's
u=u0;
%X=g/m^3
q=0.9./(1000*3600); %mol / (g s)
K_s=(2.2e-5*1000)/44; %mol/m^3
mu_max=1.3/3600; %1/s
%diff equations----------------------------
%u(1)=concentration 
%u(2)=density
pak1=@(t,u) [Ka.*(c_sat_179-u(1))-(q.*u(2)); (u(2).*mu_max.*u(1))./(K_s+u(1))];
pak2=@(t,u) [-(q.*u(2)); (u(2).*mu_max.*u(1))./(K_s+u(1))];
fcn=pak1;
opt=odeset('AbsTol',1e-10,'RelTol',1e-12,'Events',@Event1);
% 
ateout = [];
ayeout = [];
aieout = [];
while (t(end) < tend)
[at, ay, ate, aye, aie] = ode45(fcn, [t(end), tend], u,opt);
t = cat(1, t, at(2:end));
u(1) = cat(1, u(1), ay(2:end,1));
u(2) = cat(1, u(2), ay(2:end,2));
  
ateout = [ateout; ate];
ayeout = [ayeout; aye];
aieout = [aieout; aie];
if u(end,1) >= c_max
fcn=pak2;
opt=odeset('AbsTol',1e-10,'RelTol',1e-12,'Events',@Event2);
elseif u(end,1) <= c_min 
fcn=pak1;
opt=odeset('AbsTol',1e-10,'RelTol',1e-12,'Events',@Event1);
end
end
plot(t,u(1).*1e6,'-o')
set(gca, 'Nextplot','add','box','on','FontSize',14)
xlabel('Time (s)','Interpreter','LaTeX'); 
ylabel('CO$_2$ concentration ($\mu$mol L$^{-1}$)','Interpreter','LaTeX');
%Event functions-----------------
function [val, isterm, dir]=Event1(t,u)
c_max=833./1e6; %mol/L
val=[c_max-u(end,1)];
isterm=[1];
dir=0;
end
function [val, isterm, dir]=Event2(t,u)
c_min = 46./1e6; %mol/L
val=[c_min-u(end,1)];
isterm=[1];
dir=0;
end
My problem is apparently here
t = cat(1, t, at(2:end));
u(:,1) = cat(1, u(end,1), ay(2:end, :));
u(:,2) = cat(1, u(end,2), ay(2:end, :));
There I get following error:
"Error using cat. Dimensions of matrices being concatenated are not consistent."

Best Answer

The problem is you event-functions. You should (most likely) make the check for the first component of u, since the event-function tests on the "current" value of u, not the last element of the entire solution. Something like this:
%Event functions-----------------
function [val, isterm, dir]=Event1(t,u)
c_max=833./1e6; %mol/L

val=[c_max-u(1)];
isterm=[1];
dir=0;
end
function [val, isterm, dir]=Event2(t,u)
c_min = 46./1e6; %mol/L
val=[c_min-u(1)];
isterm=[1];
dir=0;
end
HTH
Related Question