Switch between two sets of coupled ODE's systems using Event

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."

Answer 1

[Bjorn Gustavsson]

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

[Jose Bolivar]

Hej Björn

Thank you for your suggestion, I have corrected that. 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."

What I'm trying to do with the piece of code above is to "update" the values of t and u. That has worked before when I was working with single (not coupled) ODE's. My problem now is that I don't know how to code it when u is a vector with two columns.

Do you know what can be done there?

[Bjorn Gustavsson]

I think it should be as simple as:

t = [t;at(2:end)];
u = [u;ay(2:end,:)];

HTH

[Jose Bolivar]

Hej again Björn

Thank you once more for your reply.

Sadly it does not solve the problem, the current error says "Error using vertcat. Dimensions of matrices being concatenated are not consistent."

It feels like I should have one code line for u(1) and another one for u(2) since both are changing. Or maybe the error is somewhere else.

[Bjorn Gustavsson]

The problem you have now is that you initially set u to u0, which makes it a 2-by-1 array, so when you try to concatenate it with the output from ode45 which is n x 2, then you get a problem. Try this modification:

u = u0';

HTH

[Jose Bolivar]

Hi again

When I do so I get even more errors. This time I get errors from ode45 which I did not get before.

As I said before, I'm almost sure my problem is how I define "u" in the loop. Matlab also says the problem is there:

Error using vertcat
Dimensions of matrices being concatenated are not consistent.

Error in Carbonpack (line XXX)
u = [u;ay(2:end,:)];

[Bjorn Gustavsson]

This runs through:

%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 1000];
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], u0,opt);
  
  t = [t; at(2:end)];
  u = [u; ay(2:end,:)];
  u0 = u(end,:);
  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,'-')


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(~,u)
c_max = 833./1e6; %mol/L
val = c_max-u(1);
isterm = 1;
dir=0;
end
function [val, isterm, dir]=Event2(~,u)
c_min = 46./1e6; %mol/L
val = c_min-u(1);
isterm = 1;
dir=0;
end

Don't know what else to say

[Jose Bolivar]

Hej Björn!!

Thank you very much, I actually was about to write a new message when I saw you sent this one. I could figure it out and I got the same plot as you did with your code. ODE45 is working now but something is off with the coditions for toggling on and off the ODEs since the plot reaches a value way higher than c_max and reaches steady state instead of activating the other ODE and decreasing the concentration.

I really do not want to disturb you anymore, I think this is a minor stuff that I hope I can figure out soon.

Thank you so much for your help! :D

[Bjorn Gustavsson]

Notice that I've now changed the script, to call ode45 with an initial condition, and to update that one after each run, and fixed the events-functions too.

[Jose Bolivar]

This works!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Taaaaaack!

[Bjorn Gustavsson]

Good,

You're welcome, Ingen orsak, De nada(?)