Mfrhtyj

In some chemical systems containing poorly soluble salts, the concentration of various species can be calculated by solving systems of algebraic equations.

E.g. for a system containing an excess of salt $AC$, if $S_{AC}$ is the number of moles of undissolved salt divided by the volume of the physical system, $A$ is the concentration of dissociated species 'A' and $C$ is the concentration of dissociated species 'C', the following system must be solved:

$left{ matrix {{N_A-A-S_{AC}=0}\{N_C-C-S_{AC}=0}\{A cdot C-K_{AC}=0}} right.$

where $N_A$ and $N_C$ are known real constants expressing the total nominal concentration of species 'A' and 'C', respectively.

The problem is that the 3rd equation is only applicable when there actually is some solid $AC$, i.e. when $S_{AC}>0$, which only happens when $A cdot C-K_{AC} geq 0$. Otherwise, the system reduces to:

$left{ matrix {{N_A-A=0}\{N_C-C=0}\{S_{AC}=0}} right.$

I would like to solve this type of system numerically, using a general method that captures both the case where the solid is present and the case where the solid is not present.

The analytical solution is possible in this simple case:

$left{ matrix {{S_{AC}=max(0,frac {N_A + N_C - sqrt {(N_A-N_C)^2 +4 cdot K_{AC}}} {2})}\{C=N_C-S_{AC}} \{A=N_A-S_{AC}} } right.$

which is consistent, because $N_A + N_C - sqrt {(N_A-N_C)^2 +4 cdot K_{AC}}>0$ when $N_A cdot N_C > K_{AC}$, so in all cases all concentrations are positive.

However, there are often other chemical equilibria to consider (i.e. additional equations, often nonlinear), potentially different salts (i.e. multiple equations subject to conditions), thus the situation can become very complicated (or even impossible) to handle analytically.

I tried various approaches for the general numerical solution, where the last equation is included in all cases but is only 'used' when necessary, e.g.:

$left{ matrix {{N_A-A-S_{AC}=0}\{N_C-C-S_{AC}=0}\{S_{AC} cdot (A cdot C-K_{AC})=0}} right.$

which does not work because the 3rd equation is always satisfied when $S_{AC}=0$, so the solver always concludes that there is no solid.

This other version instead:

$left{ matrix {{N_A-A-S_{AC}=0}\{N_C-C-S_{AC}=0}\{if S_{AC}=0 then 0=0 else A cdot C-K_{AC}=0}} right.$

works for this simple example when the initial value of $S_{AC}$ is $>0$, but often fails for more complicated examples, e.g. with 2 insoluble salts $AC$ and $AK$ sharing a common dissociated species:

$left{ matrix {{N_A-A-S_{AC}-S_{AK}=0}\{N_C-C-S_{AC}=0}\{N_K-K-S_{AK}=0}\{if S_{AC}=0 then 0=0 else A cdot C-K_{AC}=0}\{if S_{AK}=0 then 0=0 else A cdot K-K_{AK}=0}} right.$

Solvers that use multivariate Newton-like algorithms complain about 'singular matrix' for some initial values of the variables, and I don't know how to avoid that.

I tried something like this too:

$left{ matrix {{N_A-A-S_{AC}-S_{AK}=0}\{N_C-C-S_{AC}=0}\{N_K-K-S_{AK}=0}\{max(0,A cdot C-K_{AC})=0}\{max(0,A cdot K-K_{AK})=0}} right.$

which worked even worse :(

I described the problem in programming terms in this post, which has got no replies so far.

Can you suggest a (better) way to handle this kind of systems in a numerical context?

Thanks!