ode45 taking too long to solve set of ODEs.

The title sums it up well. I'm using the ode45 function in order to solve a set of ODEs. However, the script fails and I get an error, as it ends up stopping itself before using up all of the PCs memory. I have used a similar script successfully to solve a simpler set. So I was wondering if there is a way to maybe make the code do less iterations so it can actually solve the set? My code is:

function dYfuncvecdW = PBR_ODE(W,Yfuncvec)
F_EtOH = Yfuncvec(1);
F_DEE = Yfuncvec(2);
F_Ete = Yfuncvec(3);
F_AA = Yfuncvec(4);
F_But = Yfuncvec(5);
F_BuOH = Yfuncvec(6);
F_Bte = Yfuncvec(7);
F_H2O = Yfuncvec(8);
F_H2 = Yfuncvec(9);
y = Yfuncvec(10);
%% Total Flowrate
Ft = F_EtOH + F_DEE + F_Ete + F_AA + F_But + F_BuOH + F_Bte + F_H2;
%% Defined Values
Pto = 4000; %kPa
alpha = 0.001; %kg^-1
Fto = 1290; %kmol/hr
%% Stoichiometry
P_EtOH = Pto * (F_EtOH/Ft) * y;
P_DEE = Pto * (F_DEE/Ft) * y;
P_Ete = Pto * (F_Ete/Ft) * y;
P_AA = Pto * (F_AA/Ft) * y;
P_But = Pto * (F_But/Ft) * y;
P_BuOH = Pto * (F_BuOH/Ft) * y;
P_Bte = Pto * (F_Bte/Ft) * y;
P_H2O = Pto * (F_H2O/Ft) * y;
P_H2 = Pto * (F_H2/Ft) * y;
%% Rate Constants
k1_EtOH = 0.005;
k2_DEE = 0.1455;
k3_EtOH = 0.0427;
k4_AA = 45294;
k5_EtOH = 0.0488;
k6_BuOH = 0.0288;
%% Equilibrium Constants
K_EQ1 = 0.998;
K_EQ2 = 6.71;
K_EQ3 = 9.03; 
K_EQ4 = 0.148;
K_EQ5 = 20860;
K_EQ6 = 39502;
%% Net Rates
rEtOH = - k1_EtOH*(P_EtOH^2 - (P_DEE*P_H2O)/K_EQ1) - k3_EtOH*(P_EtOH - (P_AA*P_H2)/K_EQ3) - k5_EtOH*(P_EtOH^2 - (P_BuOH*P_H2O)/K_EQ5) + k2_DEE*(P_DEE - (P_EtOH*P_Ete)/K_EQ2);
rDEE = - k2_DEE*(P_DEE - (P_EtOH*P_Ete)/K_EQ2) + (k2_DEE/2)*(P_EtOH^2 - (P_DEE*P_H2O)/K_EQ1);
rEte = + k2_DEE*(P_DEE - (P_EtOH*P_Ete)/K_EQ2);  
rAA = - k4_AA*((P_AA^2 * P_H2) - (P_But*P_H2O)/K_EQ4) + k3_EtOH*(P_EtOH - (P_AA*P_H2)/K_EQ3);
rBut =  (k4_AA/2)*((P_AA^2 * P_H2) + (P_But*P_H2O)/K_EQ4); 
rBuOH = - k6_BuOH*(P_BuOH - (P_Bte*P_H2O)/K_EQ6) + (k5_EtOH/2)*(P_EtOH^2 - (P_BuOH*P_H2O)/K_EQ5);
rBte =  k6_BuOH*(P_BuOH - (P_Bte*P_H2O)/K_EQ6);
rH2O =  (k1_EtOH/2)*(P_EtOH^2 - (P_DEE*P_H2O)/K_EQ1) + (k4_AA/2)*((P_AA^2 * P_H2) - (P_But*P_H2O)/K_EQ4) + (k5_EtOH/2)*(P_EtOH^2 - (P_BuOH*P_H2O)/K_EQ5) + k6_BuOH*(P_BuOH - (P_Bte*P_H2O)/K_EQ6);
rH2 = k3_EtOH*(P_EtOH - (P_AA*P_H2)/K_EQ3);
%% Selectivity
if (W > 0.0001)
    S_BuOH = F_BuOH / (F_DEE + F_AA);
else
    S_BuOH = 0;
end
%% Differential Equations
dF_EtOH_dW = rEtOH;
dF_DEE_dW = rDEE;
dF_Ete_dW = rEte;
dF_AA_dW = rAA;
dF_But_dW = rBut;
dF_BuOH_dW = rBuOH;
dF_Bte_dW = rBte;
dF_H2O_dW = rH2O;
dF_H2_dW = rH2;
dydW = 0 - (alpha / 2 / y * Ft / Fto);
dYfuncvecdW = [dF_EtOH_dW; dF_DEE_dW; dF_Ete_dW; dF_AA_dW; dF_But_dW; dF_BuOH_dW; dF_Bte_dW; dF_H2O_dW; dF_H2_dW; dydW]

;

Answer 1

[Torsten]

Try to use ode15s instead of ode45.

And prescribe the vector of output times. Otherwise ode45 will save all of its successful integration steps which can make your PC collapse.

Answer 2

[Sam Chak]

Hi @H,

A simple test shows that the code works with ode15s() solver, and it produces results.

Take note: If the initial values are zeros(1, 10), then it produces NaN. So, you may need to check if the equations are correct such that one or more states must NOT approach zero.

tspan         = [0 1];
Yfuncvec0     = 0.1*ones(1, 10);    % initial values
[W, Yfuncvec] = ode15s(@PBR_ODE, tspan, Yfuncvec0);

plot(W, Yfuncvec(:,5), 'linewidth', 1.5), grid on, % legend('show')
xlabel('W'), ylabel('State #5')

function dYfuncvecdW = PBR_ODE(W,Yfuncvec)
    F_EtOH = Yfuncvec(1);
    F_DEE = Yfuncvec(2);
    F_Ete = Yfuncvec(3);
    F_AA = Yfuncvec(4);
    F_But = Yfuncvec(5);
    F_BuOH = Yfuncvec(6);
    F_Bte = Yfuncvec(7);
    F_H2O = Yfuncvec(8);
    F_H2 = Yfuncvec(9);
    y = Yfuncvec(10);
    %% Total Flowrate
    Ft = F_EtOH + F_DEE + F_Ete + F_AA + F_But + F_BuOH + F_Bte + F_H2;
    %% Defined Values
    Pto = 4000; %kPa
    alpha = 0.001; %kg^-1
    Fto = 1290; %kmol/hr
    %% Stoichiometry
    P_EtOH = Pto * (F_EtOH/Ft) * y;
    P_DEE = Pto * (F_DEE/Ft) * y;
    P_Ete = Pto * (F_Ete/Ft) * y;
    P_AA = Pto * (F_AA/Ft) * y;
    P_But = Pto * (F_But/Ft) * y;
    P_BuOH = Pto * (F_BuOH/Ft) * y;
    P_Bte = Pto * (F_Bte/Ft) * y;
    P_H2O = Pto * (F_H2O/Ft) * y;
    P_H2 = Pto * (F_H2/Ft) * y;
    %% Rate Constants
    k1_EtOH = 0.005;
    k2_DEE = 0.1455;
    k3_EtOH = 0.0427;
    k4_AA = 45294;
    k5_EtOH = 0.0488;
    k6_BuOH = 0.0288;
    %% Equilibrium Constants
    K_EQ1 = 0.998;
    K_EQ2 = 6.71;
    K_EQ3 = 9.03; 
    K_EQ4 = 0.148;
    K_EQ5 = 20860;
    K_EQ6 = 39502;
    %% Net Rates
    rEtOH = - k1_EtOH*(P_EtOH^2 - (P_DEE*P_H2O)/K_EQ1) - k3_EtOH*(P_EtOH - (P_AA*P_H2)/K_EQ3) - k5_EtOH*(P_EtOH^2 - (P_BuOH*P_H2O)/K_EQ5) + k2_DEE*(P_DEE - (P_EtOH*P_Ete)/K_EQ2);
    rDEE = - k2_DEE*(P_DEE - (P_EtOH*P_Ete)/K_EQ2) + (k2_DEE/2)*(P_EtOH^2 - (P_DEE*P_H2O)/K_EQ1);
    rEte = + k2_DEE*(P_DEE - (P_EtOH*P_Ete)/K_EQ2);  
    rAA = - k4_AA*((P_AA^2 * P_H2) - (P_But*P_H2O)/K_EQ4) + k3_EtOH*(P_EtOH - (P_AA*P_H2)/K_EQ3);
    rBut =  (k4_AA/2)*((P_AA^2 * P_H2) + (P_But*P_H2O)/K_EQ4); 
    rBuOH = - k6_BuOH*(P_BuOH - (P_Bte*P_H2O)/K_EQ6) + (k5_EtOH/2)*(P_EtOH^2 - (P_BuOH*P_H2O)/K_EQ5);
    rBte =  k6_BuOH*(P_BuOH - (P_Bte*P_H2O)/K_EQ6);
    rH2O =  (k1_EtOH/2)*(P_EtOH^2 - (P_DEE*P_H2O)/K_EQ1) + (k4_AA/2)*((P_AA^2 * P_H2) - (P_But*P_H2O)/K_EQ4) + (k5_EtOH/2)*(P_EtOH^2 - (P_BuOH*P_H2O)/K_EQ5) + k6_BuOH*(P_BuOH - (P_Bte*P_H2O)/K_EQ6);
    rH2 = k3_EtOH*(P_EtOH - (P_AA*P_H2)/K_EQ3);
    %% Selectivity
    if (W > 0.0001)
        S_BuOH = F_BuOH / (F_DEE + F_AA);
    else
        S_BuOH = 0;
    end
    %% Differential Equations
    dF_EtOH_dW = rEtOH;
    dF_DEE_dW = rDEE;
    dF_Ete_dW = rEte;
    dF_AA_dW = rAA;
    dF_But_dW = rBut;
    dF_BuOH_dW = rBuOH;
    dF_Bte_dW = rBte;
    dF_H2O_dW = rH2O;
    dF_H2_dW = rH2;
    dydW = 0 - (alpha / 2 / y * Ft / Fto);
    dYfuncvecdW = [dF_EtOH_dW; dF_DEE_dW; dF_Ete_dW; dF_AA_dW; dF_But_dW; dF_BuOH_dW; dF_Bte_dW; dF_H2O_dW; dF_H2_dW; dydW];
end