```python
from collections import OrderedDict
import warnings
from IPython.display import display
import sympy as sp
import matplotlib.pyplot as plt
from chempy.chemistry import Equilibrium, ReactionSystem, Substance
from chempy.kinetics.integrated import binary_rev
from chempy.kinetics.ode import get_odesys
%matplotlib inline
sp.init_printing()
```


```python
t, kf, kb, X, Y, Z = sp.symbols('t k_f k_b X Y Z')
binary_rev(t, kf, kb, Y, Z, X, sp)
```


```python
eq = Equilibrium.from_string('Fe+3 + SCN- = FeSCN+2; 10**2.065')
rsys = ReactionSystem(eq.as_reactions(kf=900))
rsys
```


```python
odesys, extra = get_odesys(rsys)
```


```python
def plot_against_analytic(odsy, c0, tend=1, integrator='cvode', atol=1e-12, rtol=1e-10):
    xout, yout, info = odsy.integrate(tend, c0, integrator=integrator, atol=atol, rtol=rtol)
    ref = binary_rev(xout[1:], rsys.rxns[0].param, rsys.rxns[1].param, c0['FeSCN+2'], c0['Fe+3'], c0['SCN-'])
    {k: v for k, v in info.items() if not k.startswith('internal')}
    plt.figure(figsize=(14, 4))
    plt.subplot(1, 2, 1)
    plt.plot(xout[1:], ref, alpha=0.3, linewidth=3)
    _ = odsy.plot_result()
    plt.subplot(1, 2, 2)
    plt.plot(xout[1:], yout[1:, odsy.names.index('FeSCN+2')] - ref)
    ref[:3], ref[-3:]
    return {k: v for k, v in info.items() if not k.startswith('internal')}
```


```python
conc = {'Fe+3': 10e-3, 'SCN-': 2e-3, 'FeSCN+2': 2e-4}
plot_against_analytic(odesys, conc)
```

Below are the derivations for ``binary_rev``:


```python
x, a, b, c, U, V, S = sp.symbols('x a b c U V S')
u, v = 2*a*x + b, sp.sqrt(b**2 - 4*a*c)  # b**2 > 4*a*c
Prim = sp.log((U-V)/(U+V))/V  # primitive recip. 2nd order polynomial
prim = sp.log((u-v)/(u+v))/v  # primitive recip. 2nd order polynomial
prim.diff(x).simplify()
```


```python
y = Y + X - x
z = Z + X - x
```


```python
dxdt = kf*y*z - kb*x
dxdt
```


```python
dxdt.expand()
```


```python
expr2 = dxdt.expand().collect(x)
expr2
```


```python
coeffs = expr2.as_poly(x).coeffs()
coeffs
```


```python
eq = (sp.exp(Prim) - sp.exp(t + S))
eq
```


```python
sp.Eq(Prim*V, (t+S)*V)
```


```python
eq = sp.exp(Prim*V) - sp.exp((t+S)*V)
eq
```


```python
eq2 = eq.subs({U: u})
eq2
```


```python
sol, = sp.solve(eq2, x)
sol
```


```python
s, = sp.solve(sol.subs(t, 0) - X, S)
s
```


```python
sol2 = sol.subs(S, s).simplify()
sol2
```


```python
R, W = sp.symbols('R W')
r = b + 2*X*a
w = sp.exp(-V*t)
sol3 = sol2.subs({w: W, r: R}).simplify().collect(W)
sol3
```


```python
exprs = [
    sp.Eq(x, sol3),
    sp.Eq(R, r),
    sp.Eq(W, w),
    sp.Eq(V, v),
    sp.Eq(a, coeffs[0]),
    sp.Eq(b, coeffs[1]),
    sp.Eq(c, coeffs[2]),
]
for expr in exprs:
    display(expr)

```


```python
master = exprs[0]
for e in exprs[1:]:
    master = master.subs(e.lhs, e.rhs)
master
```


```python
cses, (reform,) = sp.cse(master)
reform
```


```python
cses
```


```python

```


```python
ret = 'return ' + str(reform.rhs).replace('k_f', 'kf').replace('k_b', 'kb')
for key, val in cses:
    print('%s = %s' % (key, str(val).replace('k_f', 'kf').replace('k_b', 'kb').replace('exp', 'be.exp').replace('sqrt', 'be.sqrt')))
print(ret)
```


```python
with warnings.catch_warnings():
    warnings.filterwarnings('error')
    binary_rev(1.0, 1e10, 1e-4, 0, 3e-7, 2e-7)  # check that our formulation avoids overflow
```

Below we will test ``get_native``:


```python
from chempy.kinetics._native import get_native
```


```python
print({k: v.composition for k, v in rsys.substances.items()})
rsys.substances = OrderedDict([(k, Substance.from_formula(v.name)) for k, v in rsys.substances.items()])
print({k: v.composition for k, v in rsys.substances.items()})
```


```python
nsys = get_native(rsys, odesys, 'gsl')
```


```python
plot_against_analytic(nsys, conc, integrator='gsl', atol=1e-16, rtol=1e-16)
```


```python

```

b888721c7fa2cedacdd70e6191114889e824a02f
ipynb
Jupyter Notebook
examples/_integrated.ipynb
bertiewooster/chempy
115adc1d570aa1631baff4374f3128ce23fa7776
BSD-2-Clause
2015-10-30T03:41:05.000Z
2022-03-31T05:01:17.000Z
examples/_integrated.ipynb
bertiewooster/chempy
115adc1d570aa1631baff4374f3128ce23fa7776
BSD-2-Clause
2015-11-03T13:31:23.000Z
2022-03-31T16:46:19.000Z
examples/_integrated.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