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

2 Ansichten (letzte 30 Tage)
Farangis am 23 Jan. 2023
Kommentiert: Davide Masiello am 23 Jan. 2023
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)');
##### 2 Kommentare1 älteren Kommentar anzeigen1 älteren Kommentar ausblenden
Farangis am 23 Jan. 2023
I have a vector of temperatures versus time.

Melden Sie sich an, um zu kommentieren.

### Akzeptierte Antwort

Davide Masiello am 23 Jan. 2023
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
##### 4 Kommentare3 ältere Kommentare anzeigen3 ältere Kommentare ausblenden
Davide Masiello am 23 Jan. 2023
No problem, please consider to click on the Accept button if you think it was useful. Cheers.

Melden Sie sich an, um zu kommentieren.

### Kategorien

Find more on Ordinary Differential Equations in Help Center and File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by