# Conditional non-linear systems of equations
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.


```python
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)
$$


```python
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](https://en.wikipedia.org/wiki/Solubility_equilibrium) 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:


```python
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:


```python
get_f(c, init_concs + eq_constants, False)
```


```python
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):


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


```python
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:


```python
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):


```python
%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

0f3b66dbbe1657af8ebc1064412f708cd83841ba
ipynb
Jupyter Notebook
examples/conditional.ipynb
bjodah/pyneqsys
1e8b63ee6e820d33a87b95c057655df03e1e4fc6
BSD-2-Clause
2015-11-12T09:18:36.000Z
2022-03-16T23:25:03.000Z
examples/conditional.ipynb
bjodah/pyneqsys
1e8b63ee6e820d33a87b95c057655df03e1e4fc6
BSD-2-Clause
2016-01-01T22:37:28.000Z
2018-08-15T11:17:48.000Z
examples/conditional.ipynb
bjodah/pyneqsys
1e8b63ee6e820d33a87b95c057655df03e1e4fc6
BSD-2-Clause
2018-01-17T02:53:21.000Z
2021-02-02T22:52:56.000Z
Qwen/Qwen-72B
1. YES
2. YES