一種常見的化學反應平衡關係是離子在溶液中的水解平衡，由於存在平衡常數，所以在給定條件下（溫度，pH值等）能夠求出溶液中鹽離子和弱酸根離子的濃度。但是，化學平衡方程本身是非線性的，而且在一般情況下，溶液中存在眾多的平衡體系。使得該問題的求解難度很高。例如：

已知碳酸（$H_{2}CO_{3}$）溶液中存在如下水解平衡關係，計算$H^{+}$濃度為$10^{-5}$mol/L，${CO_{3}}^{2-}$濃度為$10^{-9}$mol/L時，溶液中$H_{2}CO_{3}$的濃度。

 $$H^{+}+OH^{-}\rightleftharpoons H_{2}O  \qquad    K_{w}=10^{-14}$$

 $${HCO_{3}}^{-}+H^{+}\rightleftharpoons H_{2}CO_{3}  \qquad    K_{a1}=4.2\times 10^{-7}$$

 $${CO_{3}}^{2-}+H^{+}\rightleftharpoons {HCO_{3}}^{-}   \qquad   K_{a2}=5.61\times 10^{-11}$$

分析：

該問題共涉及6種物質（$H^{+}$，$OH^{-}$，${CO_{3}}^{2-}$，${HCO_{3}}^{-}$，$H_{2}CO_{3}$，$ H_{2}O$），其中五種物質存在濃度（$H^{+}$，$OH^{-}$，${CO_{3}}^{2-}$，${HCO_{3}}^{-}$，$H_{2}CO_{3}$），除去已知濃度的$H^{+}$和${CO_{3}}^{2-}$之外，還有3種物質的濃度需要求解（$OH^{-}$，${HCO_{3}}^{-}$，$H_{2}CO_{3}$）。顯然，該溶液體系共有三種化學平衡關係，可列出三個化學平衡方程，理論上可以求解出未知的三種物質的濃度。化學平衡方程如下：

$$K_{w}=\left [ H^{+} \right ]\cdot \left [ OH^{-} \right ]$$

$$K_{a1}=\frac{\left [ H^{+} \right ]\cdot \left [ {HCO_{3}}^{-} \right ]}{\left [ H_{2}CO_{3} \right ]}$$

$$K_{a2}=\frac{\left [ H^{+} \right ]\cdot \left [ {CO_{3}}^{2-} \right ]}{\left [ {HCO_{3}}^{-} \right ]}$$

 將已知量帶入上述方程後，可以使用MATLAB中的solve方法求解該非線性方程組，程式碼及說明如下：

clear,clc
%% matlab 求解非線性方程組 以化學反應平衡為例
num=20;
CO3=10^(-9);   % 已知的CO3^2-濃度
H=10^(-5);     % 已知的H+濃度
kw=10^(-14);
ka1=4.2*10^(-7);
ka2=5.61*10^(-11);
% H 和 CO3為已知量，在方程組中視為常數
syms  OH HCO3 H2CO3 exp1  % 需要求解的其他3種物質（OH，HCO3，H2CO3）和需要求解的方程組（exp1）設定為符號變數
exp1=[kw==H*OH,ka1==(H*HCO3)/H2CO3,ka2==(H*CO3)/HCO3];  % 需要求解的方程組
% 非線性方程組求解方法 Y = solve(eqns,vars)，eqns為需要求解的方程組，vars需要求解的變數
% 注意：該問題中共存在3個方程和3個未知數，若eqns中列舉了全部三個方程，則vars中必須全部列舉所有未知變數才能求解
% 例如，在該問題中，保持exp1不變，vars=[OH HCO3]，則會導致無解。
% solution為結構體變數，需要求解的變數儲存在其中，預設儲存為解析解（分數表示）。
% 分別為：solution.OH，solution.HCO3，solution.H2CO3
solution = solve(exp1, [OH HCO3 H2CO3]);
% vpa方法可以將解析解（分數表示）轉化為數值解（小數表示）
% double方法可以將符號變數轉化為數值變數
H2CO3_result=double(vpa(solution.H2CO3)); 

該程式碼的計算結果：$H_{2}CO_{3}$的濃度為0.0042 mol/L

特別的，如果想要研究$H_{2}CO_{3}$與$H^{+}$和${CO_{3}}^{2-}$濃度的關係，也可以用類似的方法實現。

例如，可以計算出$H^{+}$濃度範圍在[$10^{-5}$ $2\times 10^{-5}$]和${CO_{3}}^{2-}$濃度範圍在[$10^{-9}$ $2\times 10^{-9}$]時對應的$H_{2}CO_{3}$濃度，以$H^{+}$濃度為x軸，${CO_{3}}^{2-}$濃度為y軸，$H_{2}CO_{3}$濃度為z軸繪製三維影象進行展示。程式碼如下：

clear,clc
%% matlab 求解非線性方程組 以化學平衡繪圖為例
num=20;    % 每個維度的取樣數量
H=linspace(1*10^(-5),2*10^(-5),num);   % CO3^2-濃度範圍
CO3=linspace(1*10^(-9),2*10^(-9),num); % H+濃度範圍
H2CO3=zeros(num,num);
for indexCO3=1:num
    for indexH=1:num
        [H2CO3(indexCO3,indexH)]=ChemicalEquilibrium(CO3(indexCO3),H(indexH));
    end
end
[X,Y] = meshgrid(CO3,H); %將其x，y軸網格化
Fig = mesh(X,Y,H2CO3); %繪製三維曲面圖
xlabel('H^+濃度 mol/L','fontsize',10,'fontname','宋體');
ylabel('{CO_3}^{2-}濃度 mol/L','fontsize',10,'fontname','宋體');
zlabel('H_2CO_3濃度 mol/L','fontsize',10,'fontname','宋體');

function [H2CO3_result]=ChemicalEquilibrium(CO3,H)
kw=10^(-14);
ka1=4.2*10^(-7);
ka2=5.61*10^(-11);
syms  OH H2CO3 HCO3 exp1
exp1=[kw==H*OH,ka2==(H*CO3)/HCO3,ka1==(H*HCO3)/H2CO3,];
solution = solve(exp1, [OH HCO3 H2CO3]);
H2CO3_result=double(vpa(solution.H2CO3));
end

執行結果如圖所示：

參考資料：

Equations and systems solver - MATLAB solve - MathWorks 中國

特別感謝：

楊翠婷