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

2 views (last 30 days)
Farangis on 23 Jan 2023
Commented: Davide Masiello on 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 CommentsShow NoneHide None
Davide Masiello on 23 Jan 2023
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 on 23 Jan 2023
I have a vector of temperatures versus time.

Sign in to comment.

### Accepted Answer

Davide Masiello on 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 CommentsShow 2 older commentsHide 2 older comments
Farangis on 23 Jan 2023
Thank you very much for the answer.
All the best!
Davide Masiello on 23 Jan 2023
No problem, please consider to click on the Accept button if you think it was useful. Cheers.

Sign in to comment.

### Categories

Find more on Assembly 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!