ODE15s Events function not working when switching between different ODE systems

Hello everyone

This is the third time I ask for help about this project I'm working on and I want to thank you all that have helped me in the past (especially Björn if you are reading this).

I have a pretty complex code that works mostly fine, I have a minor problem that I do not understand. I'm using event functions to switch on and off different sets of ODE systems, this works fine but it doesn't work in the end, where I want the function to stop when the CO2 concentration is zero, instead, matlab continues solving and it gives me negative CO2 concentrations that make no sense.

Possible causes:

• I'm currently solving with ODE15s and it takes a while to solve (~1 min), maybe the problem is the type of ODE solver I'm using but I have tried with (I think) all of them and the solution never converges.
• Something is off with the event function that should stop the calculation at c_CO2=0
• Something is off with my ODE system (which I doubt because the rest before the end looks fine)

I attach the code below, I have commented where the problems could be (Björn if you are reading this, I think this is the last time I'll ask for help since it is the last last step):

clear all
close all
clc


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





tspan = [0 5e4];
tstart = tspan(1);
t=tstart;
tend = tspan(end);


c_0 = (150./1e6).*1000; %mol/m^3 Initial concentration CO2
H2CO3_0=0; %mol/m^3 Initial concentration H2CO3
HCO3_0=6437e-6*1000; %mol/m^3 Initial concentration HCO3-
H_0=0.00000001*1000; %mol/m^3 Initial concentration H+

pH_75=0.00000003162*1000; %mol/m^3 H+ concentration that corresponds to pH=7.5
pH_76=0.00000002512*1000; %mol/m^3 H+ concentration that corresponds to pH=7.6
X0=1*1000; %g/m^3
u0=[c_0;X0;H2CO3_0;HCO3_0;H_0]; %initial conditions for ODE's


%X=g/m^3

Ka=1.2e-3; %1/s
q=(0.9./(1000.*3600)); %mol / (g s)
K_s=3000.*((2.2e-5.*1000)./44); %mol/m^3
mu_max=3.6e-04; %1/s


u=u0';

%diff equations----------------------------
%u(1)=concentration CO2
%u(2)=density
%u(3) concentration H2CO3
%u(4) concentration HCO3-
%u(5) concentration H+ / pH

k_a=18; %s^-1 Backward constant
k_b=0.04; %s^-1 Forward constant
k_21=1e7; %s^-1 Forward constant
k_12=5e10./1000; %m^3/(mol*s) Backward constant

%Ka.*(c_sat_179-u(1))= CO2 dissolution
%pak1: CO2 dissolution: on

pak1=@(t,u) [Ka.*(c_sat_179-u(1))-(q.*u(2))+(k_a.*u(3))-(k_b.*u(1)); %dCO2/dt
    (u(2).*mu_max.*u(1))./(K_s+u(1)); %dX/dt
    ((k_12.*u(4).*u(5))+k_b.*u(1)-(k_21 + k_a).*u(3)); %dH2CO3/dt
    -((k_12.*u(4).*u(5))-k_21.*u(3)); %dHCO3/dt
    -((k_12.*u(4).*u(5))-k_21.*u(3))]; %dH/dt];

%pak2: CO2 dissolution: off
pak2=@(t,u) [-(q.*u(2))+(k_a.*u(3))-(k_b.*u(1)); %dCO2/dt
    (u(2).*mu_max.*u(1))./(K_s+u(1)); %dX/dt
    ((k_12.*u(4).*u(5))+k_b.*u(1)-(k_21 + k_a).*u(3)); %dH2CO3/dt
    -((k_12.*u(4).*u(5))-k_21.*u(3)); %dHCO3/dt
    -((k_12.*u(4).*u(5))-k_21.*u(3))]; %dH/dt];

fcn=pak1;
opt=odeset('AbsTol',1e-10,'RelTol',1e-12,'Events',@Event1);


while (t(end) < tend)
  
  [at, ay] = ode15s(fcn, [t(end), tend], u0,opt); %ODE15s maybe an issue?
  
  t = [t; at(2:end)];
  u = [u; ay(2:end,:)];
  u0 = u(end,:);
  
  if u(end,5) >= pH_75 %change pak below pH = 7.5
    fcn = pak2;
    opt = odeset('AbsTol',1e-10,'RelTol',1e-12,'Events',@Event2);
  elseif u(end,5) <= pH_76 %change pak above pH = 7.6
    fcn=pak1;
    opt=odeset('AbsTol',1e-10,'RelTol',1e-12,'Events',@Event1); 
  elseif u(end,2) <= (Ka.*(c_sat_179-u(1))+k_a.*u(3)-k_b.*u(1))./q; %Possible problem: This should let c_CO2 reach zero, c_CO2 starts to decrease (very good) but the solution doesn't stop at c_CO2=0 but it becomes negative..
    fcn = pak2;
    opt = odeset('AbsTol',1e-10,'RelTol',1e-12,'Events',@Event3);
  end
end

subplot(1,2,1)
plot(t./3600,u(:,1).*(1e6./1000),'-')
set(gca, 'Nextplot','add','box','on','FontSize',14)
xlabel('Time (hr)','Interpreter','LaTeX'); 
ylabel('CO$_2$ concentration ($\mu$mol L$^{-1}$)','Interpreter','LaTeX');
%title('K_a = 2K_a')


subplot(1,2,2)
plot(t./3600,u(:,2)./1000,'-')
set(gca, 'Nextplot','add','box','on','FontSize',14)
xlabel('Time (hr)','Interpreter','LaTeX'); 
ylabel('Algae density (g L$^{-1}$)','Interpreter','LaTeX');

figure

subplot(3,2,1)
plot(t./3600,u(:,3).*(1e6./1000),'-')
hold on
set(gca, 'Nextplot','add','box','on','FontSize',14)
xlabel('Time (hr)','Interpreter','LaTeX'); 
ylabel('H$_2$CO$_3$ concentration ($\mu$mol L$^{-1}$)','Interpreter','LaTeX');

subplot(3,2,2)
plot(t./3600,u(:,4).*(1e6./1000),'-')
hold on
set(gca, 'Nextplot','add','box','on','FontSize',14)
xlabel('Time (hr)','Interpreter','LaTeX'); 
ylabel('HCO$_3^{-1}$ concentration ($\mu$mol L$^{-1}$)','Interpreter','LaTeX');

subplot(3,2,3)
plot(t./3600,-log10(u(:,5)./1000),'-')
hold on
set(gca, 'Nextplot','add','box','on','FontSize',14)
xlabel('Time (hr)','Interpreter','LaTeX'); 
ylabel('pH','Interpreter','LaTeX');

%Event functions-----------------

function [val, isterm, dir]=Event1(~,u)
pH_75=0.00000003162*1000; %mol/m^3
val = pH_75-u(5);
isterm = 1;
dir=0;
end
function [val, isterm, dir]=Event2(~,u)
pH_76=0.00000002512*1000; %mol/m^3
val = pH_76-u(5);
isterm = 1;
dir=0;
end
function [val, isterm, dir]=Event3(~,u) %-----------------------Possible problem u(1)=c_CO2
val = 0-u(1); %Should stop if c_CO2=0 and this value should not be negative
isterm = 1;
dir=0;
end

Answer 1

[Bjorn Gustavsson]

It might be that it is only in event3 that you check for the CO2 concentration. You could also check if you change to the exact even-function from ballode:

function [value,isterminal,direction] = events(t,y)
% Locate the time when height passes through zero in a decreasing direction
% and stop integration.
value = y(1);     % detect height = 0
isterminal = 1;   % stop the integration
direction = -1;   % negative direction

That should be pretty identical to your CO2 concentration going negative?

[Jose Bolivar]

Hello again kind sir

I have tried this, unfortunately it does not work..

I've bee thinking, even if the events function does not work, the CO2 concentration shouldn't be negative.. never. Maybe there is some error in the code, which I don't really think so since everything else looks alright to me... My second (and most probably according to me) guess is that the solver ODE15s goes kaputt towards the end and that is somehow messing with the events function

[Jose Bolivar]

I have also tried to stop the function using event3 based on pH [u(5)] (which is not what I want but just to try) like the other events functions and that did not work either.

[Bjorn Gustavsson]

Check the documentation of the odeNN-functions. Some of them allows for constraining the solution to be positive. That should correspond to setting a field "nonnegative" or something reasonably similar.

I imagine you could combine the different stopping-events, making sure that an event is set as soon as any of the events occur.

[Jose Bolivar]

Good suggestion! I just read this thread: https://se.mathworks.com/matlabcentral/answers/409126-how-to-make-ode15s-not-to-cross-zero

Since I have a mass matrix I'm supposed to use events, NonNegative won't work in my case (I tried also). But that is not bad news, I guess then that the problem is with Event3, it is weird though that Event1 and Event2 work perfectly. maybe something is off in the while loop?

In this part the CO2 will start to decrease:

  elseif u(end,2) <= (Ka.*(c_sat_179-u(1))+k_a.*u(3)-k_b.*u(1))./q; % Let c_CO2 reach zero
    fcn = pak2;
    opt = odeset('AbsTol',1e-10,'RelTol',1e-12,'Events',@Event3);
    

So pak2 is on (good) and Event3 should also be on.. Another possibility is that Event3 is not defined properly:

function [val, isterm, dir]=Event3(~,u)
val = 0-u(1); %u(1) is the CO2 concentration, it should stop when zero is reached
isterm = 1;
dir=0;
end

[Jose Bolivar]

Hi again!

I could fix the negative concentration by changing some parameters, event is still not working but It's not that important as long as the concentration doesn't go below zero. Thank you again!