```python
import sympy as sp
sp.init_printing()
```

Analyitc solution for CSTR of 2 A -> B


```python
symbs = t, f, fcA, fcB, IA, IB, k, c1, c2 = sp.symbols('t f phi_A phi_B I_A I_B k c1 c2', real=True, positive=True)
symbs
```


```python
def analytic(t, f, fcA, fcB, IA, IB, k, c1, c2):
    u = sp.sqrt(f*(f + 4*fcA*k))
    v = u*sp.tanh(u/2*(t - c1))
    return [
        (-f + v)/(2*k),
        sp.exp(-f*t)*c2 + (f + 2*(fcA + fcB)*k - v)/(2*k)
    ]
```


```python
exprs = analytic(*symbs)
exprs
```


```python
exprs0 = [expr.subs(t, 0) for expr in exprs]
exprs0
```


```python
sol = sp.solve([expr - c0 for expr, c0 in zip(exprs0, [IA, IB])], [c1, c2])
sol
```


```python
exprs2 = [expr.subs(dict(zip([c1, c2], sol[1]))) for expr in exprs]
exprs2
```


```python
[expr.subs(t, 0).simplify() for expr in exprs2]
```


```python
%matplotlib inline
import numpy as np
from chempy import ReactionSystem
from chempy.kinetics.ode import get_odesys
rsys = ReactionSystem.from_string("OH + OH -> H2O2; 'k'")
cstr, extra = get_odesys(rsys, include_params=False, cstr=True)
fr, fc = extra['cstr_fr_fc']
print(cstr.names, cstr.param_names)
cstr.exprs
```


```python
args = tout, c0, params = np.linspace(0, .17), {'OH': 2, 'H2O2': 3}, {'k': 5, fc['OH']: 42, fc['H2O2']: 11, fr: 13}
res = cstr.integrate(*args)
res.plot()
```


```python
def analytic_alt(t, f, fcA, fcB, IA, IB, k):
    u = sp.sqrt(f*(f + 4*fcA*k))
    q = sp.atanh(-u*(sp.sqrt(f) + 2*k*IA)/(f*(f + 4*fcA*k)))
    v = u*sp.tanh(u/2*t - q)
    w = sp.exp(-f*t)/2/k
    y = 2*k*(fcA + fcB)
    return [
        (-f + v)/(2*k),
        w*(sp.exp(f*t)*(f + y - v) - y + 2*k*(IA + IB))
    ]
```


```python
def analytic_alt2(t, f, fcA, fcB, IA, IB, k, n):
    one_point_five = sp.S(3)/2
    a0, b0, x, y = fcA, fcB, IA, IB
    Sqrt, Tanh, ArcTanh, E = sp.sqrt, sp.tanh, sp.atanh, sp.E
    return [
        (-f + Sqrt(f)*Sqrt(f + 8*a0*k)*
     Tanh((Sqrt(f)*Sqrt(f + 8*a0*k)*t - 2*ArcTanh((-(f**one_point_five*Sqrt(f + 8*a0*k)) - 4*Sqrt(f)*k*Sqrt(f + 8*a0*k)*x)/(f**2 + 8*a0*f*k)))/2)
    )/(4*k),
    
    (-8*b0*k + 8*b0*E**(f*t)*k + E**(f*t)*f*n - 4*a0*k*n + 4*a0*E**(f*t)*k*n + 4*k*n*x + 8*k*y - 
    E**(f*t)*Sqrt(f)*Sqrt(f + 8*a0*k)*n*Tanh((Sqrt(f)*Sqrt(f + 8*a0*k)*
         (t - (2*ArcTanh((-(f**one_point_five*Sqrt(f + 8*a0*k)) - 4*Sqrt(f)*k*Sqrt(f + 8*a0*k)*x)/(f**2 + 8*a0*f*k)))/(Sqrt(f)*Sqrt(f + 8*a0*k))))
        /2))/(8*E**(f*t)*k)
    ]
#     return [
#     (-f + Sqrt(f)*Sqrt(f + 8*a0*k)*
#    Tanh((Sqrt(f)*Sqrt(f + 8*a0*k)*t - 2*ArcTanh((-(f**one_point_five*Sqrt(f + 8*a0*k)) - 4*Sqrt(f)*k*Sqrt(f + 8*a0*k)*x)/(f**2 + 8*a0*f*k)))/2)
#   )/(4*k),
# (E**(f*t)*f - 4*a0*k - 8*b0*k + 4*a0*E**(f*t)*k + 8*b0*E**(f*t)*k + 4*k*x + 8*k*y - 
#   E**(f*t)*Sqrt(f)*Sqrt(f + 8*a0*k)*Tanh((Sqrt(f)*Sqrt(f + 8*a0*k)*
#        (t - (2*ArcTanh((-(f**one_point_five*Sqrt(f + 8*a0*k)) - 4*Sqrt(f)*k*Sqrt(f + 8*a0*k)*x)/(f**2 + 8*a0*f*k)))/(Sqrt(f)*Sqrt(f + 8*a0*k))))
#       /2))/(8*E**(f*t)*k)
#     ]
        
```


```python
n = sp.Symbol('n')
exprs_alt = analytic_alt2(*symbs[:-2], n)
exprs_alt
```


```python
cses, expr_cse = sp.cse([expr.subs({fcA: sp.Symbol('fr'), fcB: sp.Symbol('fp'), f: sp.Symbol('fv'),
                                     IA: sp.Symbol('r'), IB: sp.Symbol('p')}) for expr in exprs_alt])
```


```python
print(
    '\n'.join(['%s = %s' % (lhs, rhs) for lhs, rhs in cses] + ['return (\n    %s\n)' % str(expr_cse)[1:-1]]).replace(
        '3/2', 'three/2').replace(
        'exp', 'be.exp').replace(
        'sqrt', 'be.sqrt').replace(
        'atanh', 'ATANH').replace(
        'tanh', 'be.tanh\n    ').replace(
        'ATANH', 'atanh')
)
```


```python
[expr.subs(t, 0).simplify() for expr in exprs_alt]
```


```python
[expr.diff(t).subs(t, 0).simplify() for expr in exprs_alt]
```


```python
print(list(rsys.substances))
rsys
```


```python
print(cstr.names, cstr.param_names)
cstr.exprs
```


```python
symbs[:-2]
```


```python
_cb = sp.lambdify(symbs[:-2] + (n,), exprs_alt)
def calc_analytic(xout, y0, p):
    return _cb(xout, p[fr], p[fc['OH']], p[fc['H2O2']], y0['OH'], y0['H2O2'], p['k'], 1)
```


```python
def get_analytic(result):
    ref = calc_analytic(
        result.xout,
        {k: res.get_dep(k)[0] for k in result.odesys.names},
        {k: res.get_param(k) for k in result.odesys.param_names})
    return np.array([ref[{'OH': 0, 'H2O2': 1}[k]] for k in result.odesys.names]).T
```


```python
yref = get_analytic(res)
print(yref.shape)
```


```python
res.plot()
res.plot(y=yref)
```


```python
res.plot(y=res.yout - yref)
```


```python
from chempy.kinetics.integrated import binary_irrev_cstr
```


```python
def get_analytic2(result):
    ref = binary_irrev_cstr(result.xout, result.get_param('k'), result.get_dep('OH')[0],
                           result.get_dep('H2O2')[0], result.get_param(fc['OH']), result.get_param(fc['H2O2']),
                           result.get_param(fr))
    return np.array([ref[{'OH': 0, 'H2O2': 1}[k]] for k in result.odesys.names]).T
```


```python
res.plot(y=res.yout - get_analytic2(res))
```

68efe3a9973f0570150c175ba3cef510f8ef23a9
ipynb
Jupyter Notebook
chempy/kinetics/tests/_derive_analytic_cstr_bireac.ipynb
Narsil/chempy
ac7217f45a8cfe3b11ca771f78f0a04c07708818
BSD-2-Clause
chempy/kinetics/tests/_derive_analytic_cstr_bireac.ipynb
Narsil/chempy
ac7217f45a8cfe3b11ca771f78f0a04c07708818
BSD-2-Clause
chempy/kinetics/tests/_derive_analytic_cstr_bireac.ipynb
Narsil/chempy
ac7217f45a8cfe3b11ca771f78f0a04c07708818
BSD-2-Clause
2022-03-21T09:01:48.000Z
2022-03-21T09:01:48.000Z
Qwen/Qwen-72B
1. YES
2. YES