MATLAB: ODE15s Events function not working when switching between different ODE systems

eventsodeode systemode15s

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

Best Answer

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?
Related Question