不同条件下溶液中的物质浓度？化学反应平衡？MATLAB求解非线性方程组

一种常见的化学反应平衡关系是离子在溶液中的水解平衡，由于存在平衡常数，所以在给定条件下（温度，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 求解非线性方程组 以化学反应平衡为例
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}$] mol/L和${CO_{3}}^{2-}$浓度范围在[$10^{-9}$, $2\times 10^{-9}$] mol/L时对应的$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 中国

 

特别感谢：

杨翠婷