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

Question

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. H₂O = H⁺ + OH^- (concentration of H₂O is not considered)

B. NH₃ + H₂O = NH₄⁺ + OH^- (concentration of H₂O is not considered)

Equilibrium constants:

[H⁺][OH⁻] = 10^-14

[NH₄⁺][OH⁻]/[NH₃] = 1.77×10⁻⁵

So we have 2 reactions (A and B) and 4 substances (H⁺, OH⁻, NH₃ and NH₄⁺), 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:

• H₂O = H⁺ + OH^-
• NH₃ + H₂O = NH₄⁺ + OH^-
• CO₂ + H₂O = H⁺ + HCO₃^-
• HCO₃^- = H⁺ + CO₃^2-

which includes 4 reactions and 7 substances?

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.