Here are chemical equilibrium equations:

$\begin{cases} \mathbf{N}^\text{T}\mathbf{X}+\mathbf{C}=\mathbf{Y}\\ \mathbf{N}\ln\mathbf{Y}=\ln\mathbf{K} \end{cases}$

Here $\mathbf{C}=\begin{pmatrix}c_1\\c_2\\\vdots\\c_n\end{pmatrix}$ are initial concentrations (mol/L) of n types of substances;

$\mathbf{N}=\begin{pmatrix} \nu_{11}&\nu_{12}&\cdots&\nu_{1n}\\ \nu_{21}&\nu_{22}&\cdots&\nu_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ \nu_{m1}&\nu_{m2}&\cdots&\nu_{mn}\\ \end{pmatrix}$ are coefficients of m chemical reactions (negative for reagents, positive for products, and zero if the substance is not related to the reaction);

$\mathbf{X}=\begin{pmatrix}x_1\\x_2\\\vdots\\x_m\end{pmatrix}$ are amounts of reactions (concentration changes for each reaction);

$\mathbf{Y}=\begin{pmatrix}y_1\\y_2\\\vdots\\y_n\end{pmatrix}$ are final (equilibrium state) concentrations;

$\mathbf{K}=\begin{pmatrix}k_1\\k_2\\\vdots\\k_m\end{pmatrix}$ are equilibrium constants for each reaction.

I need to solve X (m variables) and Y (n variables) by given C, N and K, either analytic or numeric solutions.

For Example

0.1 mol/L ammonia solution, we have 2 equilibriums:

A. H2O = H+ + OH- (concentration of H2O is not considered)

B. NH3 + H2O = NH4+ + OH- (concentration of H2O is not considered)

Equilibrium constants:

[H+][OH] = 10-14

[NH4+][OH]/[NH3] = 1.77×10−5

So we have 2 reactions (A and B) and 4 substances (H+, OH, NH3 and NH4+), and known values are:




The solutions are:



Here I used the proximation of $y_\text{OH-}=y_\text{NH4+}$, so just solved a quadratic equation.

But how can I solve a more complicated system like:

  • H2O = H+ + OH-
  • NH3 + H2O = NH4+ + OH-
  • CO2 + H2O = H+ + HCO3-
  • HCO3- = H+ + CO32-

which includes 4 reactions and 7 substances?

  • 1
    $\begingroup$ I take it $\ln\mathbf{Y}$ is $\ln$ applied to each entry of the vector $\mathbf{Y}$? $\endgroup$ Jul 6, 2021 at 14:01
  • 1
    $\begingroup$ Yes. Are there any simple notations for $(\ln y_1, \ln y_2, ..., \ln y_n)^T$? $\endgroup$
    – auntyellow
    Jul 6, 2021 at 14:06
  • $\begingroup$ It seems your second equation determines $Y$, so if $N$ has a left inverse we have $\ln Y =N^{-1} \ln K$. With $Y$ known the first equation determines $X$ similarly $\endgroup$
    – Sal
    Jul 6, 2021 at 15:58
  • $\begingroup$ No. The second equation actually includes m equations and n variables, and generally n > m. $\endgroup$
    – auntyellow
    Jul 6, 2021 at 16:11

1 Answer 1


This problem can be handled with nonlinear programming, through a minimization process.

Follows a self explanatory MATHEMATICA script showing the procedure. You can try it with more complicated arrangements.

c = {0, 0, 0.1, 0};
n = {{1, 1, 0, 0}, {0, 1, -1, 1}};
k = {10^-14, 1.77 10^-5};
X = {x1, x2};
Y = {y1, y2, y3, y4};
dif1 = Transpose[n].X - c - Y;
dif2 = n.Log[Y] - Log[k];
obj = Norm[dif1] + Norm[dif2];
vars = Join[X, Y];
restrs = Table[2/10 > vars[[i]] > 10^-14, {i, 1, Length[vars]}];
sol = NMinimize[{obj, restrs}, vars, Method -> "DifferentialEvolution"]
obj /.sol[[2]]

It gives a solution better than yours.

  • $\begingroup$ In my second example NH4HCO3/(NH4)2CO3, I don't know how many solutions there are, and I don't know if MATHIMATICA can find all solutions. $\endgroup$
    – auntyellow
    Jul 6, 2021 at 16:16
  • $\begingroup$ Post your values to verify. $\endgroup$
    – Cesareo
    Jul 6, 2021 at 16:18
  • $\begingroup$ What's the result? I'm not familiar with Mathematica and I use SciPy's fsolve and encountered some problems: stackoverflow.com/questions/68302523 $\endgroup$
    – auntyellow
    Jul 8, 2021 at 13:37
  • $\begingroup$ I ran your codes in Mathematica Cloud and got the result {0.167346,{x1=1.E-14,x2=5.21993E-8,y1=6.92626E-9,,y2=1.44213E-6,y3=0.00711273,y4=0.0846503}}, which is not correct. The accurate result is X= [7.56662656e-12 1.32159290e-03] Y= [7.56662656e-12 1.32159291e-03 9.86784071e-02 1.32159290e-03]. $\endgroup$
    – auntyellow
    Jul 8, 2021 at 13:47
  • $\begingroup$ If you compute the values for dif1and dif2with the values which allegedly is the "answer", the residuals are $\{-7.56663, -1.32127, -0.2, 0., 34.5395, 6.90851\}$ and the MATHEMATICA script furnishes the residuals $\{-6.92625*10^-9, -1.38993*10^-6, -0.107113, -0.0846502, -0.00114382, -0.0308008\}$ so the solutions given by the script are much better than the so called "answer". Regards. $\endgroup$
    – Cesareo
    Jul 8, 2021 at 14:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.