Solve rate equations with different reaction orders using SciPy ode

Question

The image below illustrates the kinetic scheme I am trying to model.

My first example focuses on the conversion of the wood-oil to non-volatiles and volatiles. The rate equations are as follows:

$$ r_w = \frac{dC_w}{dt} = -(K_{nv} + K_v)\,C_w \\ r_{nv} = \frac{dC_{nv}}{dt} = K_{nv}\,C_w \\ r_{v} = \frac{dC_{v}}{dt} = K_{v}\,C_w \\ $$

To solve for the concentrations, I use the odeint solver in SciPy as demonstrated in the Python code below.

dt = 0.0001                           # time step, delta t
tmax = 25                             # max time, s
t = np.linspace(0, tmax, num=tmax/dt) # time vector
nt = len(t)                           # total number of time steps

def rates(c, t):
    """
    w = wood-oil as conc[0]
    nv = non-volatiles as conc[1]
    v = volatiles as conc[2]
    """
    Knv = 0.3
    Kv = 0.8
    rw = -(Knv + Kv) * c[0]
    rnv = Knv * c[0]
    rv = Kv * c[0]
    return [rw, rnv, rv]

cc = sp.odeint(rates, [1, 0, 0], t)

The above code displays the following plot.

In my next approach I attempt to account for the reactions that consume the non-volatiles and volatiles; the rate equations are displayed below.

$$ r_w = \frac{dC_w}{dt} = -(K_{nv} + K_v)\,C_w \\ r_{nv} = \frac{dC_{nv}}{dt} = K_{nv}\,C_w - K_{r1}\,{C_{nv}}^{2.5} - K_{cr}\,{C_{nv}}^{0.9} - K_{c1}\,{C_{nv}}^{1.1} \\ r_{v} = \frac{dC_{v}}{dt} = K_{v}\,C_w + K_{cr}\,{C_{nv}}^{0.9} - K_{d}\,{C_{v}}^{0.9} - K_{a}\,{C_{v}}^{1.4} - K_{g}\,{C_{v}}^{0.8} - K_{r2}\,{C_{v}}^{0.7} $$

And the code I use to solve for the concentrations is shown below.

def dCdt(c, t):
    """
    w = wood-oil as conc[0]
    nv = non-volatiles as conc[1]
    v = volatiles as conc[2]
    """
    Knv = 0.3
    Kv = 0.8
    Kr1 = 9.2e7;    r1 = 2.5
    Kcr = 4.1e-5;   cr = 0.9
    Kc1 = 3.7e5;    c1 = 1.1
    Kd = 8.0e-4;    d = 0.9
    Ka = 6.1e-6;    a = 1.4
    Kg = 1.8e-4;    g = 0.8
    Kr2 = 37.0e5;   r2 = 0.7
    rw = -(Knv + Kv)*c[0]
    rnv = Knv*c[0] - Kr1*c[1]**r1 - Kcr*c[1]**cr - Kc1*c[1]**c1
    rv = Kv*c[0] + Kcr*c[1]**cr - Kd*c[2]**d - Ka*c[2]**a - Kg*c[2]**g - Kr2*c[2]**r2
    return [rw, rnv, rv]

cc2 = sp.odeint(dCdt, [1, 0, 0], t)

Unfortunately, this approach does not work and I get an error from Python about invalid value encountered in double_scalars. It appears that the odeint solver isn't capable of handling the different reaction orders.

So my next approach is to solve the system with the SciPy ode solver. The code for this is shown below:

def dcdt(t, c):
    """
    w = wood-oil as c[0]
    nv = non-volatiles as c[1]
    v = volatiles as c[2]
    """
    Knv = 0.3
    Kv = 0.8
    Kr1 = 9.2e7;    r1 = 2.5
    Kcr = 4.1e-5;   cr = 0.9
    Kc1 = 3.7e5;    c1 = 1.1
    Kd = 8.0e-4;    d = 0.9
    Ka = 6.1e-6;    a = 1.4
    Kg = 1.8e-4;    g = 0.8
    Kr2 = 37.0e5;   r2 = 0.7
    rw = -(Knv + Kv)*c[0]
    rnv = Knv*c[0] - Kr1*c[1]**r1 - Kcr*c[1]**cr - Kc1*c[1]**c1
    rv = Kv*c[0] + Kcr*c[1]**cr - Kd*c[2]**d - Ka*c[2]**a - Kg*c[2]**g - Kr2*c[2]**r2
    return [rw, rnv, rv]

# store concentrations
Coil = np.ones(nt)      # bio-oil concentration
Cnvol = np.zeros(nt)    # non-volatiles concentration
Cvol = np.zeros(nt)     # volatiles concentration

# Setup the ode integrator where 'dopri5' is Runge-Kutta 4th order
r = sp.ode(dcdt).set_integrator('dopri5', nsteps=10000)
r.set_initial_value([1, 0, 0], 0)

# integrate the odes for each time step then store the results
k = 1
while r.successful() and r.t < tmax-dt:
    r.integrate(r.t+dt)
    Coil[k] = r.y[0]
    Cnvol[k] = r.y[1]
    Cvol[k] = r.y[2]
    k += 1

Unfortunately, the ode approach does not work and I receive a warning about the system being stiff. Does anyone have suggestions on how to solve this system of rate equations in Python when the reaction order is not one?

Answer 1

First of all, there seems to be an inconsistency between the equation you wrote for $r_w$ and the code (a minus sign).

Then, there are two issues with your problem:

1. As indicated in the comments, you're system is very stiff, so you should use the vode or lsoda integrator using the 'bdf' option or use the Radau method implemented in the Assimulo package.

2. Second, the model is ill-behaved. If you use the "poor mans debugging method" (i.e. print statements), you will observe that the calculated rates make the concentrations go negative. This will lead to problems when evaluating the rates in the next step since there are powers $< 1$ involved. This results in NaNs appearing and the integrator doesn't know what to do with them.

The solution is to double-check your rate equations, double-check your initial conditions. Maybe "capping" the rate to the current value of the concentration (to allow it to reach 0.0 but not go below which is unphysical).

It will also help the numerical routine if you provide it with the analytical Jacobian (which for this $3x3$ is not too hard.

For your convenience, the code for Assimulo:

def dCdt(t,c):
       blahblah

t0 = 0.0
y0 = np.array([1.0,0.0,0.0])
tmax = 25.0

from assimulo.solvers import Radau5ODE
from assimulo.problem import Explicit_Problem

prb = Explicit_Problem(dCdt,y0)
sim = Radau5ODE(prb)
sim.reset()
t, y = sim.simulate(tmax)