``````

In [ ]:

import sympy as sp
sp.init_printing()

``````

Analyitc solution for CSTR of 2 A -> B

``````

In [ ]:

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

``````
``````

In [ ]:

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)
]

``````
``````

In [ ]:

exprs = analytic(*symbs)
exprs

``````
``````

In [ ]:

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

``````
``````

In [ ]:

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

``````
``````

In [ ]:

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

``````
``````

In [ ]:

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

``````
``````

In [ ]:

%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

``````
``````

In [ ]:

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()

``````
``````

In [ ]:

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))
]

``````
``````

In [ ]:

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)
#     ]

``````
``````

In [ ]:

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

``````
``````

In [ ]:

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])

``````
``````

In [ ]:

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')
)

``````
``````

In [ ]:

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

``````
``````

In [ ]:

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

``````
``````

In [ ]:

print(list(rsys.substances))
rsys

``````
``````

In [ ]:

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

``````
``````

In [ ]:

symbs[:-2]

``````
``````

In [ ]:

_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)

``````
``````

In [ ]:

def get_analytic(result):
ref = calc_analytic(
result.xout,
{k: res.named_dep(k)[0] for k in result.odesys.names},
{k: res.named_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

``````
``````

In [ ]:

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

``````
``````

In [ ]:

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

``````
``````

In [ ]:

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

``````
``````

In [ ]:

from chempy.kinetics.integrated import binary_irrev_cstr

``````
``````

In [ ]:

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

``````
``````

In [ ]:

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

``````
``````

In [ ]:

# Pseduo reversible reaction (batch reactor)
#
#     A  +  B   <->   C
# t  a0    b0-x     c0 + x
import sympy
kf, kb, a0, b0, c0, x = sympy.symbols('k_{\\rm\\ f} k_{\\rm\\ b} a0 b0 c0 x', real=True, nonnegative=True)
rf = kf*a0*(b0-x)
rb = kb*(c0+x)
dxdt = (rf - rb).expand().factor(x)
dxdt

``````
``````

In [ ]:

integral = sympy.integrate(1/dxdt, (x, 0, x))
integral

``````
``````

In [ ]:

sol, = sympy.solve(integral - t, x)
sol

``````
``````

In [ ]:

print(sympy.ccode(sol.subs({kf: 'kf', kb: 'kb', a0: 'major', b0: 'minor', c0: 'prod'})))

``````