DAE ode15s coupled initial conditions / boundary condition

Hi guys,

I try to solve differential-algebraic equations - which is basicly describing the dissociation of some Acid in water with another Base.

The Systems looks like this:

%Initial conditions 
initital_c_C0 = 0;
initital_c_C1 = 0;
initital_c_C2 = 0;
initital_c_H = 0;
initital_c_Na = 0;
initital_c_OH = 0;
initital_c_NaOH = 0;
initital_F_IA = 0.5;
initital_F_NaOH = 1;

%% Definition of differential-algebraic equation system using the so called mass matrix M(t,y)*y' = f(t,y)
N = zeros(9);

for i = 1 : (7)
    N(i,i) = 1;
end
M=sparse(N);

%% Call of the ode-solver
% feeding
startvector = [initital_c_C0; initital_c_C1; initital_c_C2; initital_c_H; initital_c_Na; initital_c_OH; initital_c_NaOH; initital_F_NaOH; initital_F_IA];
disp('Start ODE15S');

t_end = 7;
time_interval = [0:1:t_end];
options = odeset('Mass', M, 'RelTol', 1e-5, 'Abstol', 1e-5, 'MStateDependence', 'weak', 'MassSingular','yes');

And my DAE

%% Startvector
c_C0 = startvector(1);
c_C1 = startvector(2);
c_C2 = startvector(3);
c_H  = startvector(4);
c_OH = startvector(5);
c_Na = startvector(6);
c_NaOH = startvector(7);
F_Na = startvector(8);
F_IA = startvector(9);

%% DAE 
% reaction rates
f(1,:) = f_1^2 * k_1_R * c_C1 * c_H - f_1 * k_1_H * c_C0;
f(2,:) = f_2 * f_1 * k_2_R * c_C2 * c_H - f_1 * k_2_H * c_C1 - f_1^2 * k_1_R * c_C1 * c_H + f_1 * k_1_H * c_C0;
f(3,:) = - f_2 * f_1 * k_2_R * c_C2 * c_H - f_1 * k_2_H * c_C1;
f(4,:) = - f_1^2 * k_1_R * c_C1 * c_H + f_1 * k_1_H * c_C0 - f_2 * f_1 * k_2_R * c_C2 * c_H + f_1 * k_2_H * c_C1 + k_w_H - f_1^2 * k_w_R * c_H * c_OH;
f(5,:) = k_w_H - f_1^2 * k_w_R * c_H * c_OH;
f(6,:) = k_Na_H * c_NaOH - k_Na_R * f_1^2 * c_Na  * c_OH;
f(7,:) = - k_Na_H * c_NaOH + k_Na_R * f_1^2 * c_Na * c_OH;

% Feed mass balance
f(8,:) = c_Na + c_NaOH - F_Na;
f(9,:) = c_C0 + c_C1 + c_C2 - F_IA ;

So what I basicly want to implement is that:

% Feed mass balance
f(8,:) = c_Na + c_NaOH - F_Na;
f(9,:) = c_C0 + c_C1 + c_C2 - F_IA ;

is true for for each time or in other words - that the sum of the components c_C0, c_C1 and c_C2 is always F_IA (constant). Same for c_NA and c_NaOH.

How could I implement this logic in matlab?

For the sake of clarity I did not show the code where I wrote down the constants, K, k, f and so on.

Thanks!

[darova]

If you have some equations

And then you want to add some condition:

You overdefine system. Or you have some uknown constants?

[Marko]

Yes you are right.

So what I did was to use the constraints in a algebraic equation and eliminate one differential equation. This way I still have a determined system.

Anyway, the result is still wrong.

Heres my code:

%Initial conditions 
initital_c_C0 = 0;
initital_c_C1 = 0;
initital_c_C2 = 0;
initital_c_H = 0;
initital_c_Na = 0;
initital_c_OH = 0;
initital_c_NaOH = 0;
%% Definition of differential-algebraic equation system using the so called mass matrix M(t,y)*y' = f(t,y)
N = zeros(9);
for i = 1 : (7)
    N(i,i) = 1;
end
M=sparse(N);
%% Call of the ode-solver
% feeding
startvector = [initital_c_C0; initital_c_C1; initital_c_C2; initital_c_H; initital_c_Na; initital_c_OH; initital_c_NaOH; initital_F_NaOH; initital_F_IA];
disp('Start ODE15S');
t_end = 7;
time_interval = [0:1:t_end];
options = odeset('Mass', M, 'RelTol', 1e-5, 'Abstol', 1e-5, 'MStateDependence', 'weak', 'MassSingular','yes');

function f = equations_solver_mass_balance_pH(t, startvector)

%% Paramater
% activity coefficient - Index is z (charge number) 
f_1 = 1;
f_2 = 1;

% forward reaction rate constants set to arbitrary number
k_1_H  = 1;
k_2_H  = 1;
k_w_H  = 1;
k_Na_H = 1;

% Dissociation constants
K1 = 10^(-3.85);
K2 = 10^(-5.5);
Kw = 10^(-14);
K_Na = 10^(13);

% calculated backward reaction rate
k_1_R  = k_1_H/K1;
k_2_R  = k_2_H/K2;
k_w_R  = k_w_H/Kw;
k_Na_R = k_Na_H/K_Na;

%Feed concentration  
F_IA = 0.5;
F_Na = 1;


%% Startvector
c_C0 = startvector(1);
c_C2 = startvector(2);
c_H  = startvector(3);
c_OH = startvector(4);
c_Na = startvector(5);
c_C1 = startvector(6);
c_NaOH = startvector(7);

%% DAE 
% reaction rates
f(1,:) = f_1^2 * k_1_R * c_C1 * c_H - f_1 * k_1_H * c_C0;
%f(2,:) = f_2 * f_1 * k_2_R * c_C2 * c_H - f_1 * k_2_H * c_C1 - f_1^2 * k_1_R * c_C1 * c_H + f_1 * k_1_H * c_C0;
f(2,:) = - f_2 * f_1 * k_2_R * c_C2 * c_H - f_1 * k_2_H * c_C1;
f(3,:) = - f_1^2 * k_1_R * c_C1 * c_H + f_1 * k_1_H * c_C0 - f_2 * f_1 * k_2_R * c_C2 * c_H + f_1 * k_2_H * c_C1 + k_w_H - f_1^2 * k_w_R * c_H * c_OH;
f(4,:) = k_w_H - f_1^2 * k_w_R * c_H * c_OH;
f(5,:) = k_Na_H * c_NaOH - k_Na_R * f_1^2 * c_Na  * c_OH;

% Feed mass balance
f(6,:) = c_C0 + c_C1 + c_C2 - F_IA ;
f(7,:) = c_Na + c_NaOH - F_Na;
end

[darova]

I think you are overcompicating the problem

%% Startvector
% ...
c_C1 = F_FIA - c_C2 - c_C0;
c_NaOH = F_Na - c_Na;

%% DAE