In [1]:
from __future__ import absolute_import, division, print_function
from collections import defaultdict
from ipywidgets import interact
import matplotlib.pyplot as plt
from chempy import Reaction, Substance, ReactionSystem
from chempy.kinetics.ode import get_odesys
from chempy.kinetics.analysis import plot_reaction_contributions
from chempy.printing.tables import UnimolecularTable, BimolecularTable
from chempy.util.graph import rsys2graph
import sympy
sympy.init_printing()
%matplotlib inline
In [2]:
A, B, C, D = map(Substance, 'ABCD')
One = sympy.S.One
reactions = r0, r1, r2 = [
    Reaction({'A'}, {'B'}, 4*One/100, name='R1: A cons.'),
    Reaction({'B', 'C'}, {'A', 'C'}, 10**(4*One), name='R2: A reform.'),
    Reaction({'B': 2}, {'B', 'C'}, 3*10**(7*One), name='R3: C form.')
]
rsys = ReactionSystem(reactions, (A, B, C))
rsys
Out[2]:
AB 0.04 R1: A cons.
B + CA + C 10000 R2: A reform.
2 BB + C 3⋅107 R3: C form.
In [3]:
uni, not_uni = UnimolecularTable.from_ReactionSystem(rsys)
bi, not_bi = BimolecularTable.from_ReactionSystem(rsys)
assert not (not_uni & not_bi), "There are only uni- & bi-molecular reactions in this set"
In [4]:
uni
Out[4]:
In [5]:
bi
Out[5]:
In [6]:
rsys2graph(rsys, 'robertson.png', save='.')
from IPython.display import Image; Image('robertson.png')
Out[6]:
In [7]:
odesys, extra = get_odesys(rsys, include_params=True)
odesys.exprs
Out[7]:
$\displaystyle \left( - \frac{y_{0}}{25} + 10000 y_{1} y_{2}, \ \frac{y_{0}}{25} - 30000000 y_{1}^{2} - 10000 y_{1} y_{2}, \ 30000000 y_{1}^{2}\right)$
In [8]:
odesys.get_jac()
Out[8]:
$\displaystyle \left[\begin{matrix}- \frac{1}{25} & 10000 y_{2} & 10000 y_{1}\\\frac{1}{25} & - 60000000 y_{1} - 10000 y_{2} & - 10000 y_{1}\\0 & 60000000 y_{1} & 0\end{matrix}\right]$
In [9]:
c0 = defaultdict(float, {'A': 1})
result = odesys.integrate(1e10, c0, integrator='cvode', nsteps=2000)
{k: v for k, v in result.info.items() if not k.startswith('internal')}
Out[9]:
{'njvev': 0,
 'nfev': 1561,
 'n_dls_rhs_evals': 0,
 'n_dls_jac_evals': 19,
 'n_root_evals': 0,
 'n_nonlin_solv_iters': 1560,
 'n_steps': 950,
 'njev': 19,
 'n_nonlin_solv_conv_fails': 0,
 'n_rhs_evals': 1561,
 'n_lin_solv_setups': 212,
 'n_err_test_fails': 92,
 'time_prec': 0.0,
 'time_jtimes': 0.0,
 'time_wall': 0.252828045,
 'time_rhs': 0.23890527399999986,
 'time_quads': 0.0,
 'time_roots': 0.0,
 'time_cpu': 0.252918,
 'time_jac': 0.002769452,
 'steps': array([], dtype=float64),
 'fpes': array([], dtype=int64),
 'orders': array([], dtype=int64),
 'success': True,
 'atol': [1e-08],
 'rtol': 1e-08,
 'mode': 'adaptive'}
In [10]:
extra['rate_exprs_cb'](result.xout, result.yout)
Out[10]:
array([[4.00000000e-02, 0.00000000e+00, 0.00000000e+00],
       [3.99999881e-02, 5.91526936e-14, 2.66417874e-06],
       [3.99999451e-02, 2.10004515e-11, 5.62962875e-05],
       ...,
       [9.19944050e-09, 9.19944016e-09, 2.53889215e-17],
       [8.48996297e-09, 8.48996273e-09, 2.16238493e-17],
       [8.26464726e-09, 8.26464722e-09, 2.04913266e-17]])
In [11]:
result.plot(xscale='log', yscale='log')
Out[11]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff433ad5ac8>
In [12]:
fig, axes = plt.subplots(2, 2, figsize=(14, 6))
plot_reaction_contributions(result, rsys, extra['rate_exprs_cb'], 'AB', axes=axes[0, :])
plot_reaction_contributions(result, rsys, extra['rate_exprs_cb'], 'AB', axes=axes[1, :],
                            relative=True, yscale='linear')
plt.tight_layout()
/drone/src/github.com/bjodah/chempy/chempy/kinetics/analysis.py:86: RuntimeWarning: divide by zero encountered in true_divide
  factor = 1/xyp[1][:, rsys.as_substance_index(sk)] if relative else 1
/drone/src/github.com/bjodah/chempy/chempy/kinetics/analysis.py:90: RuntimeWarning: invalid value encountered in multiply
  ax.plot(varied, factor*y,

We could also have parsed the reactions from a string:

In [13]:
str_massaction = """
A -> B; 'k1'
B + C -> A + C; 'k2'
2 B -> B + C; 'k3'
"""
In [14]:
rsys3 = ReactionSystem.from_string(str_massaction, substance_factory=lambda formula: Substance(formula))
In [15]:
rsys3.substance_names()
Out[15]:
('A', 'B', 'C')
In [16]:
odesys3, extra3 = get_odesys(rsys3, include_params=False, lower_bounds=[0, 0, 0])
extra3['param_keys'], extra3['unique']
Out[16]:
([], OrderedDict([('k1', None), ('k2', None), ('k3', None)]))
In [17]:
odesys3.exprs, odesys3.params, odesys3.names, odesys3.param_names
Out[17]:
((-p_0*y_0 + p_1*y_1*y_2, p_0*y_0 - p_1*y_1*y_2 - p_2*y_1**2, p_2*y_1**2),
 (p_0, p_1, p_2),
 ('A', 'B', 'C'),
 ('k1', 'k2', 'k3'))
In [18]:
def integrate_and_plot(A0=1.0, B0=0.0, C0=0.0, lg_k1=-2, lg_k2=4, lg_k3=7, lg_tend=9):
    plt.figure(figsize=(14, 4))
    tout, yout, info = odesys3.integrate(
        10**lg_tend, {'A': A0, 'B': B0, 'C': C0},
        {'k1': 10**lg_k1, 'k2': 10**lg_k2, 'k3': 10**lg_k3},
        integrator='cvode', nsteps=3000)
    plt.subplot(1, 2, 1)
    odesys3.plot_result(xscale='log', yscale='log')
    plt.legend(loc='best')
    plt.subplot(1, 2, 2)
    plt.plot(tout[tout<.05], yout[tout<.05, odesys3.names.index('B')])
    _ = plt.legend('best')
interact(integrate_and_plot) #, **kw)
Out[18]:
<function __main__.integrate_and_plot(A0=1.0, B0=0.0, C0=0.0, lg_k1=-2, lg_k2=4, lg_k3=7, lg_tend=9)>
In [19]:
# We could also have used SymPy to construct symbolic rates:
import sympy
rsys_sym = ReactionSystem.from_string("""
A -> B; sp.Symbol('k1')
B + C -> A + C; sp.Symbol('k2')
2 B -> B + C; sp.Symbol('k3')
""", rxn_parse_kwargs=dict(globals_={'sp': sympy}), substance_factory=lambda formula: Substance(formula))
odesys_sym, _ = get_odesys(rsys_sym, params=True)
for attr in 'exprs params names param_names'.split():
    print(getattr(odesys_sym, attr))
(-k1*y_0 + k2*y_1*y_2, k1*y_0 - k2*y_1*y_2 - k3*y_1**2, k3*y_1**2)
(k1, k2, k3)
('A', 'B', 'C')
()

For larger systems it is easy to loose track of what substances are actually playing a part, here the html tables can help (note the yellow background color):

In [20]:
rsys.substances['D'] = D
uni, not_uni = UnimolecularTable.from_ReactionSystem(rsys)
uni
Out[20]:
AR1: A cons.
B
C
D
In [21]:
bi, not_bi = BimolecularTable.from_ReactionSystem(rsys)
bi
Out[21]:
ABCD
A
B-R3: C form.R2: A reform.
C--
D---
In [ ]: