In [ ]:
from __future__ import division, print_function
from operator import mul
from functools import reduce
from itertools import product
import chempy
from chempy.chemistry import Species, Equilibrium
from chempy.equilibria import EqSystem, NumSysLog, NumSysLin
import numpy as np
import sympy as sp
sp.init_printing()
import matplotlib.pyplot as plt
%matplotlib inline
sp.__version__, chempy.__version__
First, let's define our substances. ChemPy can parse chemical formulae into chempy.Substance instances with information on charge, composition, molar mass (and pretty-printing):
In [ ]:
substance_names = ['H+', 'OH-', 'NH4+', 'NH3', 'H2O']
subst = {n: Species.from_formula(n) for n in substance_names}
assert [subst[n].charge for n in substance_names] == [1, -1, 1, 0, 0], "Charges of substances"
print(u'Composition of %s: %s' % (subst['NH3'].unicode_name, subst['NH3'].composition))
Let's define som e initial concentrations and governing equilibria:
In [ ]:
init_conc = {'H+': 1e-7, 'OH-': 1e-7, 'NH4+': 1e-7, 'NH3': 1.0, 'H2O': 55.5}
x0 = [init_conc[k] for k in substance_names]
H2O_c = init_conc['H2O']
w_autop = Equilibrium({'H2O': 1}, {'H+': 1, 'OH-': 1}, 10**-14/H2O_c)
NH4p_pr = Equilibrium({'NH4+': 1}, {'H+': 1, 'NH3': 1}, 10**-9.26)
equilibria = w_autop, NH4p_pr
[(k, init_conc[k]) for k in substance_names]
In [ ]:
eqsys = EqSystem(equilibria, subst)
eqsys
We can solve the non-linear system of equations by calling the root method (the underlying representation uses pyneqsys:
In [ ]:
x, sol, sane = eqsys.root(init_conc)
x, sol['success'], sane
This system is quite easy to solve for, if we have convergence problems we can try to solve a transformed system. As an example we will use NumSysLog:
In [ ]:
logx, logsol, sane = eqsys.root(init_conc, NumSys=(NumSysLog,))
logx, logsol['success'], sane
In thise case they give essentially the same answer:
In [ ]:
x - logx
We can create symbolic representations of these systems of equations:
In [ ]:
ny = len(substance_names)
y = sp.symarray('y', ny)
i = sp.symarray('i', ny)
K = Kw, Ka = sp.symbols('K_w K_a')
w_autop.param = Kw
NH4p_pr.param = Ka
ss = sp.symarray('s', ny)
ms = sp.symarray('m', ny)
In [ ]:
numsys_log = NumSysLog(eqsys, backend=sp)
f = numsys_log.f(y, list(i)+list(K))
f
In [ ]:
numsys_lin = NumSysLin(eqsys, backend=sp)
numsys_lin.f(y, i)
In [ ]:
A, ks = eqsys.stoichs_constants(False, backend=sp)
[reduce(mul, [b**e for b, e in zip(y, row)]) for row in A]
In [ ]:
from pyneqsys.symbolic import SymbolicSys
subs = list(zip(i, x0)) + [(Kw, 10**-14), (Ka, 10**-9.26)]
numf = [_.subs(subs) for _ in f]
neqs = SymbolicSys(list(y), numf)
neqs.solve([0, 0, 0, 0, 0], solver='scipy')
In [ ]:
j = sp.Matrix(1, len(f), lambda _, q: f[q]).jacobian(y)
init_conc_j = {'H+': 1e-10, 'OH-': 1e-7, 'NH4+': 1e-7, 'NH3': 1.0, 'H2O': 55.5}
xj = eqsys.as_per_substance_array(init_conc_j)
jarr = np.array(j.subs(dict(zip(y, xj))).subs({Kw: 1e-14, Ka: 10**-9.26}).subs(
dict(zip(i, xj))))
jarr = np.asarray(jarr, dtype=np.float64)
np.log10(np.linalg.cond(jarr))
In [ ]:
j.simplify()
j
In [ ]:
eqsys.composition_balance_vectors()
In [ ]:
numsys_rref_log = NumSysLog(eqsys, True, True, backend=sp)
numsys_rref_log.f(y, list(i)+list(K))
In [ ]:
np.set_printoptions(4, linewidth=120)
scaling = 1e8
for rxn in eqsys.rxns:
rxn.param = rxn.param.subs({Kw: 1e-14, Ka: 10**-9.26})
In [ ]:
x, res, sane = eqsys.root(init_conc, rref_equil=True, rref_preserv=True)
x, res['success'], sane
In [ ]:
x, res, sane = eqsys.root(init_conc, x0=eqsys.as_per_substance_array(
{'H+': 1e-11, 'OH-': 1e-3, 'NH4+': 1e-3, 'NH3': 1.0, 'H2O': 55.5}))
res['success'], sane
In [ ]:
x, res, sane = eqsys.root(init_conc, x0=eqsys.as_per_substance_array(
{'H+': 1.7e-11, 'OH-': 3e-2, 'NH4+': 3e-2, 'NH3': 0.97, 'H2O': 55.5}))
res['success'], sane
In [ ]:
init_conc
In [ ]:
nc=60
Hp_0 = np.logspace(-3, 0, nc)
def plot_rref(**kwargs):
fig, axes = plt.subplots(2, 2, figsize=(16, 6), subplot_kw=dict(xscale='log', yscale='log'))
return [eqsys.roots(init_conc, Hp_0, 'H+', plot_kwargs={'ax': axes.flat[i]}, rref_equil=e,
rref_preserv=p, **kwargs) for i, (e, p) in enumerate(product(*[[False, True]]*2))]
In [ ]:
res_lin = plot_rref(method='lm')
[all(_[2]) for _ in res_lin]
In [ ]:
for col_id in range(len(substance_names)):
for i in range(1, 4):
plt.subplot(1, 3, i, xscale='log')
plt.gca().set_yscale('symlog', linthreshy=1e-14)
plt.plot(Hp_0, res_lin[0][0][:, col_id] - res_lin[i][0][:, col_id])
plt.tight_layout()
In [ ]:
eqsys.plot_errors(res_lin[0][0], init_conc, Hp_0, 'H+')
In [ ]:
init_conc, eqsys.ns
In [ ]:
res_log = plot_rref(NumSys=NumSysLog)
In [ ]:
eqsys.plot_errors(res_log[0][0], init_conc, Hp_0, 'H+')
In [ ]:
res_log_lin = plot_rref(NumSys=(NumSysLog, NumSysLin))
eqsys.plot_errors(res_log_lin[0][0], init_conc, Hp_0, 'H+')
In [ ]:
from chempy.equilibria import NumSysSquare
res_log_sq = plot_rref(NumSys=(NumSysLog, NumSysSquare))
eqsys.plot_errors(res_log_sq[0][0], init_conc, Hp_0, 'H+')
In [ ]:
res_sq = plot_rref(NumSys=(NumSysSquare,), method='lm')
eqsys.plot_errors(res_sq[0][0], init_conc, Hp_0, 'H+')
In [ ]:
x, res, sane = eqsys.root(x0, NumSys=NumSysLog, rref_equil=True, rref_preserv=True)
x, res['success'], sane
In [ ]:
x, res, sane = eqsys.root(x, NumSys=NumSysLin, rref_equil=True, rref_preserv=True)
x, res['success'], sane