```python
import sympy as sm
import matplotlib.pyplot as plt
import numpy as np
from chempy import ReactionSystem
from chempy.units import to_unitless, SI_base_registry as si, default_units as u, default_constants as const
from chempy.kinetics.ode import get_odesys
from chempy.kinetics.rates import RampedTemp
sm.init_printing()
%matplotlib inline
```


```python
t, t0, A, B, C1 = sm.symbols('t t0 A B C1')
```


```python
y = C1/sm.E**((A*(((t + t0)*(-B + t + t0))/sm.E**(B/(t + t0)) - B**2*sm.Ei(-(B/(t + t0)))))/2)
```


```python
y
```


```python
(y.diff(t)/y).simplify().expand().simplify().factor().powsimp(force=True)
```


```python
y.subs(t, 0)
```


```python
yunit0 = y.subs(C1, C1/y.subs(t, 0)).simplify()
yunit0
```


```python
from scipy.special import expi
f = sm.lambdify([t, t0, A, B], yunit0, modules=['numpy', {'Ei': expi}])
```


```python
R = 8.314472
T_K = 290
kB = 1.3806504e-23
h = 6.62606896e-34
dH = 80e3
dS = 10
rsys1 = ReactionSystem.from_string("""
NOBr -> NO + Br; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol)
""".format(dH=dH, dS=dS))
kref = 20836643994.118652*T_K*np.exp(-(dH - T_K*dS)/(R*T_K))
kref
```


```python
_A = kB/h*np.exp(dS/R)
_B = dH/R
```


```python
f(np.array([0, 1, 5, 20]), 290, _A, _B)
```


```python
NOBr0_M = 0.7
init_cond = dict(
    NOBr=NOBr0_M*u.M,
    NO=0*u.M,
    Br=0*u.M
)
t = 20*u.second
```


```python
def integrate_and_plot(rsys):
    odes, extra = get_odesys(rsys, unit_registry=si, constants=const, substitutions={
        'temperature': RampedTemp([T_K*u.K, 1*u.K/u.s])})
    fig, all_axes = plt.subplots(2, 3, figsize=(14, 6))
    for axes, odesys in zip(all_axes, [odes, odes.as_autonomous()]):
        res = odesys.integrate(t, init_cond, integrator='cvode')
        t_sec = to_unitless(res.xout, u.second)
        NOBr_ref = NOBr0_M*f(t_sec, T_K, _A, _B)
        cmp = to_unitless(res.yout, u.M)
        ref = np.empty_like(cmp)
        ref[:, odesys.names.index('NOBr')] = NOBr_ref
        ref[:, odesys.names.index('Br')] = NOBr0_M - NOBr_ref
        ref[:, odesys.names.index('NO')] = NOBr0_M - NOBr_ref
        axes[0].plot(t_sec, cmp)
        axes[1].plot(t_sec, cmp - ref)
        res.plot_invariant_violations(ax=axes[2])
        assert np.allclose(cmp, ref)
        print({k: v for k, v in res.info.items() if not k.startswith('internal')})    
```


```python
integrate_and_plot(rsys1)
```


```python
rsys2 = ReactionSystem.from_string("""
NOBr -> NO + Br; MassAction(EyringHS([{dH}*J/mol, {dS}*J/K/mol]))
""".format(dH=dH, dS=dS))
```


```python
integrate_and_plot(rsys2)
```

b5de45dcf79aeabca76a8fc15118ae282e8d8649
ipynb
Jupyter Notebook
examples/_eyring_kinetics_ramped_temperature.ipynb
bertiewooster/chempy
115adc1d570aa1631baff4374f3128ce23fa7776
BSD-2-Clause
2015-10-30T03:41:05.000Z
2022-03-31T05:01:17.000Z
examples/_eyring_kinetics_ramped_temperature.ipynb
bertiewooster/chempy
115adc1d570aa1631baff4374f3128ce23fa7776
BSD-2-Clause
2015-11-03T13:31:23.000Z
2022-03-31T16:46:19.000Z
examples/_eyring_kinetics_ramped_temperature.ipynb
bertiewooster/chempy
115adc1d570aa1631baff4374f3128ce23fa7776
BSD-2-Clause
2016-06-06T19:55:48.000Z
2022-03-19T23:39:13.000Z
Qwen/Qwen-72B
1. YES
2. YES