# How to solve chemical equilibrium equations $\bf{N}^T\bf{X}+\bf{C}=\bf{Y}$ and $\bf{N}\ln\bf{Y}=\ln\bf{K}$

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:

$$\mathbf{C}=\begin{pmatrix}c_\text{H+}\\c_\text{OH-}\\c_\text{NH3}\\c_\text{NH4+}\end{pmatrix}=\begin{pmatrix}0\\0\\0.1\\0\end{pmatrix}$$

$$\mathbf{N}=\begin{pmatrix}\nu_\text{A,H+}&\nu_\text{A,OH-}&0&0\\0&\nu_\text{B,OH-}&\nu_\text{B,NH3}&\nu_\text{B,NH4}\end{pmatrix}=\begin{pmatrix}1&1&0&0\\0&1&-1&1\end{pmatrix}$$

$$\mathbf{K}=\begin{pmatrix}k_A\\k_B\end{pmatrix}=\begin{pmatrix}10^{-14}\\1.77\times10^{-5}\end{pmatrix}$$

The solutions are:

$$\mathbf{X}=\begin{pmatrix}x_A\\x_B\end{pmatrix}=\begin{pmatrix}7.57\times10^{-12}\\1.32\times10^{-3}\end{pmatrix}$$

$$\mathbf{Y}=\begin{pmatrix}y_\text{H+}\\y_\text{OH-}\\y_\text{NH3}\\y_\text{NH4+}\end{pmatrix}=\begin{pmatrix}7.57\times10^{-12}\\1.32\times10^{-3}\\0.0987\\1.32\times10^{-3}\end{pmatrix}$$

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?

• I take it $\ln\mathbf{Y}$ is $\ln$ applied to each entry of the vector $\mathbf{Y}$? Jul 6, 2021 at 14:01
• Yes. Are there any simple notations for $(\ln y_1, \ln y_2, ..., \ln y_n)^T$? Jul 6, 2021 at 14:06
• 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
– Sal
Jul 6, 2021 at 15:58
• No. The second equation actually includes m equations and n variables, and generally n > m. Jul 6, 2021 at 16:11

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.

• 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. Jul 6, 2021 at 16:16
• Post your values to verify. Jul 6, 2021 at 16:18
• What's the result? I'm not familiar with Mathematica and I use SciPy's fsolve and encountered some problems: stackoverflow.com/questions/68302523 Jul 8, 2021 at 13:37
• 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]. Jul 8, 2021 at 13:47
• 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. Jul 8, 2021 at 14:27