Coupled ODEs by ode45 and could not get the figures

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')

[Umar]

Hi Hao,

You did a good job because solving ode is not easy job in Matlab. It took me a while to fix the bug in your code. I had to analyze your code with out plot and trying to display results after executing function in matlab, it lead me to error, which was declaration of function % Solve ODEs [T, Y] = ode45(@(t, Y) complexSystem(t, Y), tspan, Y0); The syntax supposed to be [t,y] = ode45(odefun,tspan,y0) For more information, please refer to https://www.mathworks.com/help/matlab/ref/ode45.html So, after updating this part of code, and using snippet code for plot, which still had error, plot(T, Y) should be replaced with plot(t,Y). Please try to execute code using modified code below and this should resolve your problem now. Also, for now use

% Initial conditions Y0 = [298; 0.75; 0.5; 0.5; 0.5; 0.5]; % Time span tspan = [0 10];

If you use 10000 for tspan, it will take longer for Matlab to execute and will never plot the function. Just another tip to keep in mind.

% Define the ODE function function dYdt = complexSystem(t, Y) % Constants and parameters Aa = 1e14; Ea = 2.24e-19; Ac = 4e13; Ec = 2.03e-19; As = 1e15; Es = 2.24e-19; Aec = 1e12; Eec = 1.4e-19; Kb = 1.38e-23; z0 = 0.033; 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);

    % Extract variables from Y
    T = Y(1); Xa = Y(2); Xc = Y(3); Xs = Y(4); Z = Y(5); Soc = Y(6);

    % ODEs
    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)); 

    % Heat generation rates
    Qa = -ma * ha * dXadt;
    Qc = -mc * hca * dXcdt;
    Qs = -ma * hs * dXsdt;
    Qec = C * dSocdt;
    Q = Qa + Qc + Qs + Qec;

    % Temperature change
    dTdt = Q;

    % Output the derivatives
    dYdt = [dTdt; dXadt; dXcdt; dXsdt; dZdt; dSocdt];
end

% Initial conditions Y0 = [298; 0.75; 0.5; 0.5; 0.5; 0.5]; % Time span tspan = [0 10];

% Solve the ODEs [t, Y] = ode45(@complexSystem, tspan, Y0);

% Plot the results

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')

Please see attached snippet code.

[Sam Chak]

Hi @Hao Li,

It is important to note that the value of the parameter 'Kb' is exceedingly small. When dividing by a very small value, you must be aware of the potential numerical implications and the resulting behavior of the system.

Given the extremely small magnitude of 'Kb,' ranging around , would you expect the initial response of the (rate of change of temperature) to be as high as the order of ? Please provide your assessment and analysis of the expected system response in this scenario.

[William Rose]

Very nice non-obvious observation by @Sam Chak.

[David Goodmanson]

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.

Answer 1

[Torsten]

% 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