As we saw in the previous notebook, the rate of chemical reactions can be described by ordinary differential equations.
In this notebook we will look at a classic example in chemical kinetics: Robertson's example. It is a 3-species reaction system:
$$ A \overset{k_1}{\rightarrow} B \\ B + C \overset{k_2}{\rightarrow} A + C \\ 2 B \overset{k_3}{\rightarrow} B + C $$where A, B and C represent three different chemical species (e.g. reactive molecules dissolved in water). The system is interesting from a numerical point of view because the stepping needs to be performed using an implicit method (requiring the Jacobian) when solving for large times.
The rate of each process follows the law of mass action, i.e. the rate is proportional to the concentration of each reacting species (to the power of their multiplicity). The proportionality constant is known as the rate constant of the reaction ($k_1,\ k_2\ \&\ k_3$ in our case). If we denote the rate of each reaction:
$$ r_1 = k_1[A] \\ r_2 = k_2[B][C] \\ r_3 = k_3[B]^2 $$$[A],\ [B],\ [C]$ denotes the concentration of respective species. We can now formulate a system of ordinary differential equations describing how the concentrations evolve over time:
$$ \frac{d[A]}{dt} = r_2 - r_1 \\ \frac{d[B]}{dt} = r_1 - r_2 - r_3 \\ \frac{d[C]}{dt} = r_3 $$We will now express these differential equations (and their Jacobian) symbolically using SymPy. We saw in the previous notebook how we can use matrices with integer coefficients to describe the system of ODEs. In practice, those matrices, will consist of mostly zeros. So dictionaries are a better suited data structure (using dictionaries is one way to represent sparse matrices).
Let us therefore write a function, which creactes SymPy expressions from dictionaries mapping strings (names) to integers (coefficients). Here is one possible representation of Robertson's example:
In [ ]:
reactions = [
# (coeff, r_stoich, net_stoich)
('k1', {'A': 1}, {'B': 1, 'A': -1}),
('k2', {'B': 1, 'C': 1}, {'A': 1, 'B': -1}),
('k3', {'B': 2}, {'B': -1, 'C': 1})
]
names = 'A B C'.split()
In [ ]:
%load_ext scipy2017codegen.exercise
Use either the %exercise or %load magic to get the exercise / solution respectively:
In [ ]:
%exercise exercise_symbolic.py
To complete the above exercise you may want to remember what the law of mass action looks like from the previous notebook:
$$ \frac{dc_i}{dt} = \sum_j S_{ij} r_j \\ r_j = k_j\prod_l c_l^{R_{jl}} $$where $k_j$, S and R corresponds to coeff, net_stoich and r_stoich respectively.
In [ ]:
sym.init_printing()
ydot, y, k = mk_exprs_symbs(reactions, names)
ydot
$\dot{\mathbf{y}}$ now represent our ODE system, where $\mathbf{y}$ is our state vector (concentrations). We will need a callback to evaluate $\dot{\mathbf{y}}$ when we integrate this ODE system numerically (using scipy.integrate.odeint). As we have seen SymPy can provide us with this callback:
In [ ]:
t = sym.symbols('t') # not used in this case.
f = sym.lambdify((y, t) + k, ydot)
In [ ]:
import numpy as np
from scipy.integrate import odeint
In [ ]:
tout = np.logspace(-6, 6)
k_vals = (0.04, 1e4, 3e7) # from the literature
y0 = [1, 0, 0]
yout, info = odeint(f, y0, tout, k_vals, full_output=True)
In [ ]:
import matplotlib.pyplot as plt
%matplotlib inline
In [ ]:
plt.loglog(tout, yout)
plt.legend(names)
print("The Jacobian was evaluated %d times." % info['nje'][-1])
If we look closer at the info-dictionary we will see that odeint (or rather LSODA which is the unerlying package) switched method from an explicit Adams method to an implicit Backward Differentiation Formula (BDF). It is common for chemical kinetics problems that the problem becomes stiff.
By default, the solver will approximate the elements in the Jacobian matrix by taking finite differences of $\mathbf{f}$. This is often works quite satisfactorily, but for larger systems it sometimes fails. A more robust (and faster) approach is to provide a callback which evaluates an analytic Jacobian. Using SymPy we can do this quite effortlessly:
In [ ]:
%exercise exercise_lambdify_jac.py
In [ ]:
yout, info = odeint(f, y0, tout, k_vals, full_output=True, Dfun=J_cb)
In [ ]:
plt.loglog(tout, yout)
plt.legend(names)
print("The Jacobian was evaluated %d times." % info['nje'][-1])
We see that the solver needed to evaluate the Jacobian fewer times (due to it being essentially exact this time around). For larger systems the impact of an analytic Jacobian is often even greater (being the difference between a failed and successful integration).
Benchmarking with and without the analytic Jacobian callback:
In [ ]:
%timeit odeint(f, y0, tout, k_vals)
In [ ]:
%timeit odeint(f, y0, tout, k_vals, Dfun=J_cb)