Event in ode45 doesn't work

Hello everyone,

I am trying to solve a system of differential equations using ode45 and I want to make it stop calculating when one variable becomes negative. So I created an event that will stop the solver, but apparently it doesn't work. I tried to reduce the tolerance of the ode45 but again no luck. If someone knows how to fix it please let me know.

function [value, isterminal, direction] = myEvent(t,c)
value      = (c(1) == 0) ;
isterminal = 1;   % Stop the integration
direction  = 0;

i      = 1500;        %Current density A/m2
q_so4  = 0.7;         %Current efficiency of sulfates
q_na   = 0.7;         %Current efficiency of sodium
tw_aem = 1;           %Transport number of water in AEM
tw_cem = 4;           %Transport number of water in CEM
N      = 100;         %Number of cells
A      = 2;           %Area of one membrane m2
V_salt = 5;               %Volume of salt compartment m3
V_acid = 5;               %Volume of acid compartment m3
V_base = 5;               %Volume of base compartment m3

function dxdt = mass_balance(t,x,i,q_so4,q_na,tw_so4,tw_na,A,N,V_salt,V_acid,V_base)

%% Constants
zso4  = -2;               %Charge of sulfate
zna   = 1;                %Charge of sodium
F     = 96485.3;          %Faraday constant sA/mol

%% Partial molar volumes
Vh   = 0;    %Partial molar volume of protons m3/mol
Vna  = 0;    %Partial molar volume of sodium m3/mol
Voh  = 0;    %Partial molar volume of hydroxides m3/mol
Vso4 = 0;    %Partial molar volume of sulfates m3/mol
Vw   = 0;

% Vh   = -5.5*10^(-6);    %Partial molar volume of protons m3/mol
% Vna  = -6.7*10^(-6);    %Partial molar volume of sodium m3/mol
% Voh  = -0.2*10^(-6);    %Partial molar volume of hydroxides m3/mol
% Vso4 = 25*10^(-6);      %Partial molar volume of sulfates m3/mol
% Vw   = 1.81*10^(-5);

%% Fluxes
Nw_na  = tw_na*i/F;               %Water flux in CEM mol/m2/s
Nw_so4 = tw_so4*i/F;              %Water flux in AEM mol/m2/s
Nna    = q_na*i/(F*abs(zna));     %Sodium flux in CEM mol/m2/s
Nso4   = q_so4*i/(F*abs(zso4));   %Sulfate flux in CEM mol/m2/s
Nh     = i/F;                     %Hydrogen cations flux in BIPM mol/m2/s
Noh    = i/F;                     %Hydroxide cations flux in BIPM mol/m2/s
Nw_el  = Nh/2;                    % Water flux due to electrolysis mol/m2/s

cso4_salt = x(1);
cna_salt  = x(2);
cso4_acid = x(3);
cna_base  = x(4);
ch        = x(5);
coh       = x(6);
cw_salt   = x(7);
cw_acid   = x(8);
cw_base   = x(9);
% V_salt    = x(10);
% V_acid    = x(11);
% V_base    = x(12);


dxdt = [ (-N*A*Nso4)/V_salt;               % Salt compartment for sulfate
         (-N*A*Nna)/V_salt;                % Salt compartment for sodium
         (N*A*Nso4)/V_acid;                % Acid compartment for sulfate
         (N*A*Nna)/V_base;                 % Base compartment for sodium
         (N*A*Nh)/V_acid;                  % Acid compartment for hydrogen cations
         (N*A*Noh)/V_base;                 % Base compartment for hydroxides
         (-N*A*(Nw_so4 + Nw_na))/V_salt;   % Salt compartment for water
         (N*A*(Nw_so4 - Nw_el))/V_acid;    % Acid compartment for water
         (N*A*(Nw_na - Nw_el))/V_base];    % Base compartment for water
         

%% Solution of mass balances

hours = 2.2;                        %Operating hours
K     = 25;                         %Number of points
tspan = linspace(0,3600*hours,K);   %Time width
c0 = [2000;                    %Initial conditions
      4000;
      100;
      100;
      200;
      100;
      40000;
      55555;
      55555;

options = odeset('Events', @myEvent,'RelTol',1e-6,'AbsTol',1e-6);
[t,c] = ode45(@(t,c) mass_balance(t,c,i,q_so4,q_na,tw_aem,tw_cem,A,N,V_salt,V_acid,V_base),tspan,c0,options);

[Cris LaPierre]

What isn't working? None of the values are negative over the defined tspan.

Answer 1

[Cris LaPierre]

You've pasted your code all out of order, but once I ordered it correctly, it runs as expected.

EDIT: I did change the first line of the myEvents function to be value = c(1);

The event does not stop the code because none of the values are negative within the defined tspan. However, if I double the tspan, the integration stops around 9,189 hours.

i      = 1500;        %Current density A/m2
q_so4  = 0.7;         %Current efficiency of sulfates
q_na   = 0.7;         %Current efficiency of sodium
tw_aem = 1;           %Transport number of water in AEM
tw_cem = 4;           %Transport number of water in CEM
N      = 100;         %Number of cells
A      = 2;           %Area of one membrane m2
V_salt = 5;               %Volume of salt compartment m3
V_acid = 5;               %Volume of acid compartment m3
V_base = 5;               %Volume of base compartment m3


% Solution of mass balances
hours = 4.2;                        %Operating hours
K     = 25;                         %Number of points
tspan = linspace(0,3600*hours,K);   %Time width
c0 = [2000;                    %Initial conditions
    4000;
    100;
    100;
    200;
    100;
    40000;
    55555;
    55555];

options = odeset('Events', @myEvent,'RelTol',1e-6,'AbsTol',1e-6);
[t,c] = ode45(@(t,c) mass_balance(t,c,i,q_so4,q_na,tw_aem,tw_cem,A,N,V_salt,V_acid,V_base),tspan,c0,options);

plot(t,c)

t(end)

ans = 9.1891e+03

function dxdt = mass_balance(t,x,i,q_so4,q_na,tw_so4,tw_na,A,N,V_salt,V_acid,V_base)

%% Constants
zso4  = -2;               %Charge of sulfate
zna   = 1;                %Charge of sodium
F     = 96485.3;          %Faraday constant sA/mol

%% Partial molar volumes
Vh   = 0;    %Partial molar volume of protons m3/mol
Vna  = 0;    %Partial molar volume of sodium m3/mol
Voh  = 0;    %Partial molar volume of hydroxides m3/mol
Vso4 = 0;    %Partial molar volume of sulfates m3/mol
Vw   = 0;

% Vh   = -5.5*10^(-6);    %Partial molar volume of protons m3/mol
% Vna  = -6.7*10^(-6);    %Partial molar volume of sodium m3/mol
% Voh  = -0.2*10^(-6);    %Partial molar volume of hydroxides m3/mol
% Vso4 = 25*10^(-6);      %Partial molar volume of sulfates m3/mol
% Vw   = 1.81*10^(-5);

%% Fluxes
Nw_na  = tw_na*i/F;               %Water flux in CEM mol/m2/s
Nw_so4 = tw_so4*i/F;              %Water flux in AEM mol/m2/s
Nna    = q_na*i/(F*abs(zna));     %Sodium flux in CEM mol/m2/s
Nso4   = q_so4*i/(F*abs(zso4));   %Sulfate flux in CEM mol/m2/s
Nh     = i/F;                     %Hydrogen cations flux in BIPM mol/m2/s
Noh    = i/F;                     %Hydroxide cations flux in BIPM mol/m2/s
Nw_el  = Nh/2;                    % Water flux due to electrolysis mol/m2/s

cso4_salt = x(1);
cna_salt  = x(2);
cso4_acid = x(3);
cna_base  = x(4);
ch        = x(5);
coh       = x(6);
cw_salt   = x(7);
cw_acid   = x(8);
cw_base   = x(9);
% V_salt    = x(10);
% V_acid    = x(11);
% V_base    = x(12);


dxdt = [ (-N*A*Nso4)/V_salt;               % Salt compartment for sulfate
    (-N*A*Nna)/V_salt;                % Salt compartment for sodium
    (N*A*Nso4)/V_acid;                % Acid compartment for sulfate
    (N*A*Nna)/V_base;                 % Base compartment for sodium
    (N*A*Nh)/V_acid;                  % Acid compartment for hydrogen cations
    (N*A*Noh)/V_base;                 % Base compartment for hydroxides
    (-N*A*(Nw_so4 + Nw_na))/V_salt;   % Salt compartment for water
    (N*A*(Nw_so4 - Nw_el))/V_acid;    % Acid compartment for water
    (N*A*(Nw_na - Nw_el))/V_base];    % Base compartment for water
end

function [value, isterminal, direction] = myEvent(t,c)
value      = c(1); % <--------- #########  Note - this line has been modified #######
isterminal = 1;   % Stop the integration
direction  = 0;
end

Answer 2

[Steven Lord]

You don't want to use == in defining the first output of your event function. That will make the value of your event function 0 (false) except when or if the value passed into the event function is exactly, down to the last bit equal to 0. That could be highly unlikely.

Instead, your event function should return a value that you want to detect when it crosses 0 (as opposed to hitting 0 exactly.) See the example in the "Writing an Event Function" section on this documentation page. Note that the definition of the position output in appleEventsFcn is just y(1). The ODE solver will detect when this value changes sign between evaluations (going from positive to negative or vice versa) and then will work to determine when that function crossed 0. In this example that represents the height of the apple above the ground being 0, or more simply when the apple hits the ground (or maybe Isaac Newton's head.)

Alternately you might want to specify the NonNegative option in your odeset call. See this documentation page for two examples of this option.