Sometimes when performing modelling work in physical sciences we use different sets of equations to describe our system depending on conditions. Sometimes it is not known beforehand which of those formulations that will be applicable (only after having solved the system of equations can we reject or accept the answer). pyneqsys provides facilities to handle precisely this situation.
In [ ]:
from __future__ import (absolute_import, division, print_function)
from functools import reduce
from operator import mul
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from pyneqsys.symbolic import SymbolicSys, linear_exprs
sp.init_printing()
Let's consider precipitation/dissolution of NaCl: $$ \rm NaCl(s) \rightleftharpoons Na^+(aq) + Cl^-(aq) $$
In [ ]:
init_concs = iNa_p, iCl_m, iNaCl = [sp.Symbol('i_'+str(i), real=True, negative=False) for i in range(3)]
c = Na_p, Cl_m, NaCl = [sp.Symbol('c_'+str(i), real=True, negative=False) for i in range(3)]
prod = lambda x: reduce(mul, x)
texnames = [r'\mathrm{%s}' % k for k in 'Na^+ Cl^- NaCl'.split()]
if the solution is saturated, then the solubility product will be constant:
$$ K_{\rm sp} = \mathrm{[Na^+][Cl^-]} $$in addition to this (conditial realtion) we can write equations for the preservation of atoms and charge:
In [ ]:
stoichs = [[1, 1, -1]]
Na = [1, 0, 1]
Cl = [0, 1, 1]
charge = [1, -1, 0]
preserv = [Na, Cl, charge]
eq_constants = [Ksp] = [sp.Symbol('K_{sp}', real=True, positive=True)]
def get_f(x, params, saturated):
init_concs = params[:3] if saturated else params[:2]
eq_constants = params[3:]
le = linear_exprs(preserv, x, linear_exprs(preserv, init_concs), rref=True)
return le + ([Na_p*Cl_m - Ksp] if saturated else [NaCl])
Our two sets of reactions are then:
In [ ]:
get_f(c, init_concs + eq_constants, False)
In [ ]:
f_true = get_f(c, init_concs + eq_constants, True)
f_false = get_f(c, init_concs + eq_constants, False)
f_true, f_false
We have one condition (a boolean describing whether the solution is saturated or not). We provide two conditionals, one for going from non-saturated to saturated (forward) and one going from saturated to non-saturated (backward):
In [ ]:
from pyneqsys.core import ConditionalNeqSys
cneqsys = ConditionalNeqSys(
[
(lambda x, p: (x[0] + x[2]) * (x[1] + x[2]) > p[3], # forward condition
lambda x, p: x[2] >= 0) # backward condition
],
lambda conds: SymbolicSys(
c, f_true if conds[0] else f_false, init_concs+eq_constants
),
latex_names=['[%s]' % n for n in texnames], latex_param_names=['[%s]_0' % n for n in texnames]
)
In [ ]:
c0, K = [0.5, 0.5, 0], [1] # Ksp for NaCl(aq) isn't 1 in reality, but used here for illustration
params = c0 + K
Solving for inital concentrations below the solubility product:
In [ ]:
cneqsys.solve([0.5, 0.5, 0], params)
no surprises there (it is of course trivial).
In order to illustrate its usefulness, let us consider addition of a more soluable sodium salt (e.g. NaOH) to a chloride rich solution (e.g. HCl):
In [ ]:
%matplotlib inline
ax_out = plt.subplot(1, 2, 1)
ax_err = plt.subplot(1, 2, 2)
xres, sols = cneqsys.solve_and_plot_series(
c0, params, np.linspace(0, 3), 0, 'kinsol',
{'ax': ax_out}, {'ax': ax_err}, fnormtol=1e-14)
_ = ax_out.legend()
note the (expected) discontinuity at $\mathrm{[Na^+]_0 = 2}$ at which point the solution became saturated