How o define two different data set of (t , T) for an ode solver that it's variables are (t , C)?

How o define two different data set of (t , T) for an ode solver that it's variables are (t , C)?

I have an ode function that gives t (time) and C (Consentration) and there are k values that changes versus T (Temperature) and also temperature changes versus t. I have two different data set of (t , T) and I have to solve the ode equations to obtain C for both (t , T) data set in a one code.

Finally, I need to optimize k and Ea values for the two data set (which will be different from each other) and "alpha" (which should be the same for each data set).

So, I think I do not know how to define these two (t , T) data set for ode solver. I would be thankful if anyone could help me.

Here is my ode:

function dcdt = Cell_deg_60C(t,c,k0_vec,Ea_vec,alpha)

R = 8.314; % J/K.mol

%% For T=60 C

% I need to define T data in a vector like it T = [320.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15]

% and For T=80 C ---> T= [321.15 334.15 348.15 353.15 353.15 353.15 353.15 353.15 353.15 353.15]

%% here I just defined one data set of T, because I have no idea how to define the second one

if t <= 720

T = 297.15 + ((320.15-297.15)*(t/540)); % K , s

else

T = 333.15;

end

%% ------------------------------------------------------------------------------------------------

k01 = k0_vec(1);

k02 = k0_vec(2);

k03 = k0_vec(3);

Ea1 = Ea_vec(1);

Ea2 = Ea_vec(2);

Ea3 = Ea_vec(3);

k1 = (k01.*exp(-Ea1./(R.*T)));

k2 = (k02.*exp(-Ea2./(R.*T)));

k3 = (k03.*exp(-Ea3./(R.*T)));

%% D=1 G=2 M=3 O=4

dcdt = zeros(4,1);

dcdt(1) = - k1.*c(1) - k3.*c(1);

dcdt(2) = k1.*c(1);

dcdt(3) = k1.*c(1) + 2*k3.*c(1) - k2.*c(3);

dcdt(4) = alpha*k2*c(3);

dcdt = [dcdt(1);dcdt(2);dcdt(3);dcdt(4)];

% here is odesolver mfile

c0 = zeros(4,1);

c0(1) = 0.14; % mol/l

c0(2) = 0.00045;

c0(3) = 0.01529;

c0(4) = 0.00101;

k0_vec = [0.63, 0.03, 0.951];

Ea_vec = [17000, 10000, 100000];

alpha = 1.44;

tspan = 60*[

9

12

27

42

57

72

87

102

117

132

]; % second

[t,c] = ode45(@(t,c) Cell_deg_60C(t,c,k0_vec,Ea_vec,alpha),tspan,c0);

plot (t/60,c);

legend('Disac','GISA','Monosac','Otheracids')

xlabel('time (min)');

ylabel('C at T=60C (mol/L)');

[Davide Masiello]

Hi @Farangis, could you please clarify what the dependence of T on time is?

Do you have a vector of temperatures versus time?

Or do you have a T(t) analytical function?

[Farangis]

Hi, @Davide Masiello.

I have a vector of temperatures versus time.

Answer 1

[Davide Masiello]

I think your code is fine, below you can see a slightly different approach that uses interpolation to define the T(t) dependence.

Regarding optimization, you'd have to provide experimental concentration profiles to carry it out.

% Parameters and ICs
c0      = [0.14 0.00045 0.01529 0.00101];                                   % Initial concentrations, mol/L
k0      = [0.63, 0.03, 0.951];                                              % Pre-exponential factors
Ea      = [17000, 10000, 100000];                                           % Activation energies
alpha   = 1.44;                                                             % Time vector, s

% Experimental temperatures and interpolation
time    = 60*[9 12 27 42 57 72 87 102 117 132];
TEMP1   = [320.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15];
TEMP2   = [321.15 334.15 348.15 353.15 353.15 353.15 353.15 353.15 353.15 353.15];
T1      = @(t)interp1(time,TEMP1,t);
T2      = @(t)interp1(time,TEMP2,t);

[t,c] = ode45(@(t,c) Cell_deg_60C(t,c,k0,Ea,alpha,T1),[time(1),time(end)],c0);
figure(1)
plot (t/60,c);
title('Kinetics with first temperature evolution')
xlabel('time (min)');
ylabel('C at T=60C (mol/L)');
legend('Disac','GISA','Monosac','Otheracids','Location','Best')

[t,c] = ode45(@(t,c) Cell_deg_60C(t,c,k0,Ea,alpha,T2),[time(1),time(end)],c0);
figure(2)
plot (t/60,c);
title('Kinetics with second temperature evolution')
xlabel('time (min)');
ylabel('C at T=80C (mol/L)');
legend('Disac','GISA','Monosac','Otheracids','Location','Best')

function dcdt = Cell_deg_60C(t,c,k0,Ea,alpha,T1)
R = 8.314;                                                                  % Universal gas constantJ/K.mol
k = (k0.*exp(-Ea/(R*T1(t))));                                               % Rate constants

dcdt(1,1) = -k(1)*c(1)-k(3)*c(1);
dcdt(2,1) = k(1)*c(1);
dcdt(3,1) = k(1)*c(1)+2*k(3)*c(1)-k(2)*c(3);
dcdt(4,1) = alpha*k(2)*c(3);
end

[Farangis]

thank you very much for your help.

I want just use T vector instead of T interpolation equation. Can I just define temperatures like this TEMP1 = [320.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15 333.15];

TEMP2 = [321.15 334.15 348.15 353.15 353.15 353.15 353.15 353.15 353.15 353.15]; and ignore the interpolation equation?

[Davide Masiello]

I don't think so, because the ODE solver has an adaptive time step. This means that it will integrate the equations at time steps different from those where your temperatures are measured, i.e. you need to interpolate the temperature values to be able to predict the temperature at any time is used by ode45.

Alternatively, you could use an analytical function T(t), if you know what that is.

[Farangis]

Thank you very much for the answer.

All the best!

[Davide Masiello]

No problem, please consider to click on the Accept button if you think it was useful. Cheers.