# Coupled ODEs by ode45 and could not get the figures

1 visualizzazione (ultimi 30 giorni)
Hao Li il 21 Lug 2024
Modificato: David Goodmanson il 21 Lug 2024
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')
##### 4 CommentiMostra 2 commenti meno recentiNascondi 2 commenti meno recenti
William Rose il 21 Lug 2024
Very nice non-obvious observation by @Sam Chak.
David Goodmanson il 21 Lug 2024
Modificato: David Goodmanson il 21 Lug 2024
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.

Accedi per commentare.

### Risposte (1)

Torsten il 21 Lug 2024
% 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
##### 0 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

### Categorie

Scopri di più su Time Series Objects in Help Center e File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by