# Coupled ODEs by ode45 and could not get the figures

57 views (last 30 days)
Hao Li on 21 Jul 2024 at 5:34
Edited: David Goodmanson on 21 Jul 2024 at 19:43
Hi everyone!
I've been learning MATLAB for about six months, so I apologize if I ask any basic questions. In part of my project, I have a system of six coupled ODEs (all functions of t). You can refer to the attached picture for more details.
To obtain the T profile in terms of t, I decided to define the function similarly to how many others do:
function dYdt = complexSystem(t, Y)
Aa = 1e14; % Pre-exponential factor for xa
Ea = 2.24e-19; % Activation energy for xa
Ac = 4e13; % Pre-exponential factor for xc
Ec = 2.03e-19; % Activation energy for xc
As = 1e15; % Pre-exponential factor for xs
Es = 2.24e-19; % Activation energy for xs
Aec = 1e12; % Pre-exponential factor for xec
Eec = 1.4e-19; % Activation energy for xec
Kb = 1.38e-23; % Boltzmann constant
z0 = 0.033; % Initial value for z
ma = 0.13;
ha = 1.714e6;
mc = 0.29;
hca = 3.14e5;
hs = 2.57e5;
C = -3600 * 1.1 * 3.0 * (1-0.12-0.51);
T = Y(1);
Xa = Y(2);
Xc = Y(3);
Xs = Y(4);
Z = Y(5);
Soc = Y(6);
dXadt = -Aa * Xa * exp(-Ea/(Kb*T));
dXcdt = Ac * Xc * exp(1 - Xc) * exp(-Ec / (Kb * T));
dXsdt = -As * Xs * exp((-Z / z0) * (-Es / (Kb * T)));
dZdt = Aa * Xs * exp((-Z / z0) * (-Es / (Kb * T)));
dSocdt = -Aec * (1 - Xc) * Xa * exp(-Eec / (Kb * T));
Qa = -ma * ha * dXadt;
Qc = -mc * hca * dXcdt;
Qs = -ma * hs * dXsdt;
Qec = C * dSocdt;
Q = Qa + Qc + Qs + Qec;
dTdt = Q;
dYdt = [dTdt; dXadt; dXcdt; dXsdt; dZdt; dSocdt];
end
Then, I defined the initial conditions and solved the ODEs using ode45.
% Initial conditions
Y0 = [298; 0.75; 0.04; 0.015; 0.0331; 0]; % [T0; Xa0; Xc0; Xs0; Z0; Soc0]
% Time span
tspan = [0 100000];
% Solve ODEs
[T, Y] = ode45(@(t, Y) complexSystem(t, Y), tspan, Y0);
Lastly, I attempted to plot the results. However, the code took forever to run, and I didn't get any figures from it.
Could my initial conditions be incorrect? Or are my ODEs not converging?
Thank you so much for your time and help!
figure;
plot(T, Y)
legend('T', 'X_a', 'X_c', 'X_s', 'Z', 'S_{oc}')
title('Solution of Corrected Complex System')
xlabel('Time t')
ylabel('State Variables')
William Rose on 21 Jul 2024 at 18:07
Very nice non-obvious observation by @Sam Chak.
David Goodmanson on 21 Jul 2024 at 18:07
Edited: David Goodmanson on 21 Jul 2024 at 19:43
Hello HL,
First of all it would be a good thing to rescale some of the variables, most importantly the time variable, because the relevant times here are extremely small. (That means that the specified time range of 0 to 10000 is too long by 18 orders of magnitude or so). You can bring in planck's constant and do things more formally, but the leading factor in first five equations range from 1e12 to 1e15. If you use a new variable tau with e.g.
t = 1e-14 tau, d/dt = 1e14 d/dtau
and plug that in, then the first differential eqn becomes
dXadtau = -(1e-14*Aa) * Xa * exp(-Ea/(Kb*T));
where the prefactor now = 1 and the ones in the other fou eqns. range from .01 to 10.
There appear to only be four independent variables. dXsdt and dZdt are proportional, so
Aa*dXsdt + As*dZ/dt = 0 --> Aa*Xs + As*Z = constant
Is that the intent? Also the Q equation implies
dQdt = - ma*ha*dXadt - mc*hca*dXcdt - ma*hs*dXsdt + C*dSocdt
so again there is a conserved quantity
-Q - ma*ha*Xa - mc*hca*Xc - ma*hs*Xs + C*Soc = constant
It's somewhat less of an issue but rescaling to a new temperature U, where T = 1e3*U or something on that order, would also help.

Torsten on 21 Jul 2024 at 9:54
% Initial conditions
Y0 = [298; 0.75; 0.04; 0.015; 0.0331; 0]; % [T0; Xa0; Xc0; Xs0; Z0; Soc0]
% Time span
tspan = [0 3000];
% Solve ODEs
[T, Y] = ode15s(@(t, Y) complexSystem(t, Y), tspan, Y0);
figure;
plot(T, Y)
legend('T', 'X_a', 'X_c', 'X_s', 'Z', 'S_{oc}')
title('Solution of Corrected Complex System')
xlabel('Time t')
ylabel('State Variables')
function dYdt = complexSystem(t, Y)
Aa = 1e14; % Pre-exponential factor for xa
Ea = 2.24e-19; % Activation energy for xa
Ac = 4e13; % Pre-exponential factor for xc
Ec = 2.03e-19; % Activation energy for xc
As = 1e15; % Pre-exponential factor for xs
Es = 2.24e-19; % Activation energy for xs
Aec = 1e12; % Pre-exponential factor for xec
Eec = 1.4e-19; % Activation energy for xec
Kb = 1.38e-23; % Boltzmann constant
z0 = 0.033; % Initial value for z
ma = 0.13;
ha = 1.714e6;
mc = 0.29;
hca = 3.14e5;
hs = 2.57e5;
C = -3600 * 1.1 * 3.0 * (1-0.12-0.51);
T = Y(1);
Xa = Y(2);
Xc = Y(3);
Xs = Y(4);
Z = Y(5);
Soc = Y(6);
dXadt = -Aa * Xa * exp(-Ea/(Kb*T));
dXcdt = Ac * Xc * exp(1 - Xc) * exp(-Ec / (Kb * T));
dXsdt = -As * Xs * exp((-Z / z0) * (-Es / (Kb * T)));
dZdt = Aa * Xs * exp((-Z / z0) * (-Es / (Kb * T)));
dSocdt = -Aec * (1 - Xc) * Xa * exp(-Eec / (Kb * T));
Qa = -ma * ha * dXadt;
Qc = -mc * hca * dXcdt;
Qs = -ma * hs * dXsdt;
Qec = C * dSocdt;
Q = Qa + Qc + Qs + Qec;
dTdt = Q;
dYdt = [dTdt; dXadt; dXcdt; dXsdt; dZdt; dSocdt];
end