酸碱滴定、络合滴定、氧化还原滴定和沉淀滴定是分析化学中化学分析的核心内容^[1]。虽然4种滴定反应的本质不同，且可能伴随着各种副反应，滴定过程较为复杂，但是利用林邦的副反应思想均可以将这4种滴定反应统一为M + L = N^[1-4]。在此基础上，可以将这4种滴定曲线统一为两种数学模型。模型的统一可以使学生更好地理解滴定分析的原理和本质。同时，目前计算机技术的快速发展为解析这些复杂的分析化学模型提供了可能。尤其是MATLAB软件能够提供高效的数值计算功能，能使用户从繁杂的数学运算分析中解脱出来；特别是其友好的用户界面及接近数学表达式的自然化语言，使用户易于学习和掌握；而且具有完备的图形处理功能，实现计算结果和编程的可视化。它应该是化学工作者解析复杂的滴定分析模型的强有力的工具。因此，本文在统一4种滴定分析数学模型的基础上，利用MATLAB软件解析滴定曲线方程模型，并绘制滴定曲线。

以标准NaOH溶液滴定HCl及Cl^-滴定Ag⁺为例，滴定反应分别为：

酸碱滴定：H⁺ + OH- = H₂O

沉淀滴定：Ag⁺ + Cl^- = AgCl

酸碱滴定及沉淀滴定反应可统一为：

M + L = N

其中M分别代表H⁺和Ag⁺，L分别代表OH^-和Cl^-，N分别代表H₂O和AgCl。在酸碱滴定和沉淀滴定反应中N为纯水和沉淀，[N] = 1。因此，考虑到体积变化，两种滴定曲线方程可以统一为^{[1, 5]}：

$\begin{array}{l}{K_{\rm{t}}}{\left[{\rm{M}} \right]^2} + {K_{\rm{t}}}c_{\rm{M}}^0\frac{{\alpha - 1}}{{\alpha + 1}}\left[{\rm{M}} \right] - 1{\rm{ = }}0\\\left[{\rm{M}} \right]{\rm{ = }}\frac{{ - {K_{\rm{t}}}c_{\rm{M}}^0\frac{{\alpha - 1}}{{\alpha + 1}} + \sqrt {{{\left( {{K_{\rm{t}}}c_{\rm{M}}^0\frac{{\alpha - 1}}{{\alpha + 1}}} \right)}^2} + 4{K_{\rm{t}}}} }}{{2{K_{\rm{t}}}}}\end{array}$

滴定常数K_t表示滴定反应进行的程度，在酸碱滴定中为K_t = 1/K_w = 10¹⁴，而在沉淀滴定中为K_t = 1/K_sp。滴定分数a = c_M/c_L。

当不考虑副反应时，利用MATLAB可以非常容易地解析方程并绘制滴定曲线，程序如下：

a = 0:0.001:2;

Kt = 10^14;

C0 = 0.1

d = (a - 1)./(a + 1);

b = Kt*C0*d;

e = b.^2;

g = sqrt(e + 4*Kt);

M = (-b + g)./(2*Kt);

pH = -log10(M)

plot(a,pM,'.')

将1/K_w = 10¹⁴，1/K_sp = 10^9.75、10^12.31、10^16.08分别赋值于K_t，即可绘制出0.1 mol·L^-1 NaOH滴定0.1 mol·L^-1HCl的酸碱滴定曲线，以及0.1 mol·L^-1 Cl^-、Br^-、I^-分别滴定0.1 mol·L^-1 Ag⁺的沉淀滴定曲线(图1)。

在统一模型中强酸碱滴定与沉淀滴定唯一的差异是K_t，而将1/K_w及不同的1/K_sp赋值于K_t，即可绘制与文献^{[1, 5, 6]}理论计算结果相一致的滴定曲线。

当强碱滴定弱酸时，由林邦副反应思想，强碱滴定弱酸的反应式可以写为：

$\begin{array}{*{20}{c}}{{\rm{O}}{{\rm{H}}^ - } + {{\rm{H}}^{\rm{ - }}}{\rm{ = }}{{\rm{H}}_2}{\rm{O}}}\\ + \\{{{\rm{A}}^ - }}\\\parallel \\{{\rm{HA}}}\end{array}$

则条件滴定常数为：

${K'_{\rm{t}}} = \frac{1}{{\left[ {{{\rm{H}}^{{{\rm{ + }}^\prime }}}} \right]\left[ {{\rm{O}}{{\rm{H}}^{\rm{ - }}}} \right]}}{\rm{ = }}\frac{{{K_{\rm{t}}}}}{{{\alpha _{{\rm{H}}\left( {\rm{A}} \right)}}}}$

考虑体积变化：

${c_{{\rm{HA}}}}{\rm{ = }}\frac{{c_{{\rm{HA}}}^0}}{{1 + \alpha }}$

其中c_HA为滴定过程中弱酸浓度，[H^+′]为未参与中和反应的质子的总浓度。式中α_H(A)为质子的副反应系数，c_HA⁰为弱酸的原始浓度。

根据林邦副反应思想，将K_t′代入到强碱滴定强酸的滴定曲线方程^{[3, 4]}，得强碱滴定弱酸的滴定曲线方程为：

${K'_{\rm{t}}}{\left[ {{{\rm{H}}^{ + '}}} \right]^2} + {K_{\rm{t}}}c_{{\rm{HA}}}^{\rm{0}}\frac{{\alpha - 1}}{{\alpha + 1}}\left[ {{{\rm{H}}^{ + '}}} \right] - 1{\rm{ = }}0$

其中[H^+′] = α_H(A)[H⁺]，${K_{\rm{t}}}'{\rm{ = }}\frac{{{K_{\rm{t}}}}}{{{\alpha _{{\rm{H}}\left( {\rm{A}} \right)}}}}$；且

$\begin{array}{l}{\alpha _{{\rm{H}}\left( {\rm{A}} \right)}} = 1 + \frac{{\left[ {{{\rm{A}}^ - }} \right]}}{{{K_{\rm{a}}}}} = \frac{{{K_{\rm{a}}} + \left[ {{{\rm{H}}^ + }} \right] + c_{{\rm{HA}}}^0}}{{{K_{\rm{a}}} + \left[ {{{\rm{H}}^ + }} \right]}}\\\left[ {{{\rm{H}}^{ + '}}} \right] = \frac{{ - {{K'}_{\rm{t}}}c_{\rm{M}}^0\frac{{\alpha - 1}}{{\alpha + 1}} + \sqrt {{{\left( {{{K'}_{\rm{t}}}c_{\rm{M}}^0\frac{{\alpha - 1}}{{\alpha + 1}}} \right)}^2} + 4{{K'}_{\rm{t}}}} }}{{2{{K'}_{\rm{t}}}}};{\rm{pH}} = - \lg \frac{{\left[ {{{\rm{H}}^{ + '}}} \right]}}{{{\alpha _{{\rm{H}}\left( {\rm{A}} \right)}}}}\end{array}$

由于副反应系数α_H(A)随滴定过程而变化，因此解析该滴定曲线方程较为复杂。本文利用一个简单的迭代：

${\alpha _{{\rm{H}}\left( {\rm{A}} \right)}}\left( i \right) = \frac{{{K_{\rm{a}}} + \left[{{{\rm{H}}^ + }} \right]\left( {i - 1} \right) + c_{{\rm{HA}}}^0}}{{{K_{\rm{a}}} + \left[{{{\rm{H}}^ + }} \right]\left( {i - 1} \right)}}$

很容易解决了α_H(A)随滴定过程变化的问题，准确地绘制了强碱滴定弱酸的滴定曲线。

利用MATLAB解析该方程并绘制滴定曲线程序如下：

a = 0:0.001:2;

Kt = 10^14;

Ka = 1.8*10^(-5);

b = a - 1;

c = a + 1;

d = b./c;

for i = 2:1:2001

   H(1) = 10^(-2.87);

   alfa(i) = (Ka + H(i - 1) + 0.1/c(i))/(Ka + H(i - 1));

   Ktsh(i) = Kt/alfa(i);

   e1(i) = Ktsh(i)*0.1*d(i);

   f1(i) = e1(i)^2;

   g1(i) = sqrt(f1(i) + 4*Ktsh(i));

   h(i) = (-e1(i) + g1(i))/(2*Ktsh(i));

   H(i) = h(i)/alfa(i);

end

pH = -log10(H)

plot(a,pH,'.')

由图2(A)可以看出，其滴定突跃范围与文献^[1]中的结果相一致，表明该模型能够准确绘制出强碱滴定弱酸时的滴定曲线，且无需解复杂的一元三次方程^[5]。另外，将10^-3、10^-5、10^-7、10^-9、10^-11分别赋值于K_a，可以非常容易地绘制出强碱滴定不同弱酸的滴定曲线(图2(B))。

对于氧化还原反应，采用林邦副反应思想，把一个半反应看作是另一个半反应的副反应，便可和络合滴定反应统一为M + L = N。其中L可为Y或e^-。而滴定曲线方程统一为^{[1, 5]}：

$\begin{array}{l}{K_{\rm{t}}}{\left[{\rm{M}} \right]^2} + \left( {{K_{\rm{t}}}c_{\rm{M}}^0\left( {\alpha - 1} \right)/\left( {\alpha + 1} \right) + 1} \right)\left[{\rm{M}} \right] - c_{\rm{M}}^0{\rm{ = }}0\\\left[{\rm{M}} \right]{\rm{ = }}\frac{{ - \left( {{K_{\rm{t}}}c_{\rm{M}}^0\left( {\alpha - 1} \right)/\left( {\alpha + 1} \right) + 1} \right) + \sqrt {{{\left( {{K_{\rm{t}}}c_{\rm{M}}^0\left( {\alpha + 1} \right) + 1} \right)}^2} + 4{K_{\rm{t}}}c_{\rm{M}}^0/\left( {\alpha + 1} \right)} }}{{2{K_{\rm{t}}}}}\end{array}$

以0.01 mol·L^-1 Y滴定0.01 mol·L^-1的Ca²⁺为例，利用MATLAB解析该方程并绘制滴定曲线，程序如下：

a = 0:0.001:2;

Kt = 10^10.69;

b = a - 1;

c = a + 1;

d = b./c;

c0 = 0.01;

e1 = Kt*c0*d + 1;

f1 = e1.^2;

g1 = sqrt(f1 + 4*Kt*c0./c);

M = ( - e1 + g1)/(2*Kt);

pM = - log10(M);

plot(a,pM,'.')

由图3可以读出EDTA滴定Ca²⁺的突跃范围为(5.299,7.692)，与文献^[7]结果相一致。因此利用统一的滴定模型能够准确绘制出络合滴定曲线。将不同K_MY′赋值给K_t，赋予c_M⁰不同的初始浓度，则分别得到络合滴定曲线(图4)。

氧化还原滴定时，以Ce⁴⁺滴定Fe²⁺为例：

$\begin{array}{*{20}{c}}{{\rm{C}}{{\rm{e}}^{4 + }} + {{\rm{e}}^ - } = {\rm{C}}{{\rm{e}}^{3 + }}}\\ + \\{{\rm{F}}{{\rm{e}}^{3 + }}}\\\parallel \\{{\rm{F}}{{\rm{e}}^{2 + }}}\end{array}$

主反应常数为：

${K_{\rm{t}}} = {10^{\frac{{E_{{\rm{C}}{{\rm{e}}^{4 + }}/{\rm{C}}{{\rm{e}}^{3 + }}}^\Theta }}{{0.059}}}}$

不考虑副反应时，滴定曲线方程为：

${K_{\rm{t}}}{\left[{{\rm{C}}{{\rm{e}}^{4 + }}} \right]^2} + \left( {{K_{\rm{t}}}{c_{{\rm{C}}{{\rm{e}}^{4 + }}}}\left( {\alpha - 1} \right)/\left( {\alpha + 1} \right) + 1} \right)\left[{{\rm{C}}{{\rm{e}}^{4 + }}} \right] - {c_{{\rm{C}}{{\rm{e}}^{4 + }}}}{\rm{ = }}0$

或

${K_{\rm{t}}}{\left[{{\rm{F}}{{\rm{e}}^{2 + }}} \right]^2} + \left( {{K_{\rm{t}}}{c_{{\rm{F}}{{\rm{e}}^{2 + }}}}\left( {\alpha - 1} \right)/\left( {\alpha + 1} \right) + 1} \right)\left[{{\rm{F}}{{\rm{e}}^{2 + }}} \right] - {c_{{\rm{F}}{{\rm{e}}^{2 + }}}}{\rm{ = }}0$

考虑副反应时，有：

$\begin{array}{l}{\alpha _{\rm{e}}} = \frac{{\left[ {{{\rm{e}}^{ - '}}} \right]}}{{\left[ {{{\rm{e}}^ - }} \right]}} = \frac{{\left[ {{{\rm{e}}^ - }} \right] + \left[ {{\rm{F}}{{\rm{e}}^{2 + }}} \right]}}{{\left[ {{{\rm{e}}^ - }} \right]}} = 1 + \frac{{\left[ {{\rm{F}}{{\rm{e}}^{2 + }}} \right]}}{{\left[ {{{\rm{e}}^ - }} \right]}}\\{K_{{\rm{Fe}}}} = \frac{{\left[ {{\rm{F}}{{\rm{e}}^{2 + }}} \right]}}{{\left[ {{\rm{F}}{{\rm{e}}^{3 + }}} \right]\left[ {{{\rm{e}}^ - }} \right]}}\\\frac{{\left[ {{\rm{F}}{{\rm{e}}^{2 + }}} \right]}}{{\left[ {{{\rm{e}}^ - }} \right]}} = {K_{{\rm{Fe}}}}\left[ {{\rm{F}}{{\rm{e}}^{3 + }}} \right]\\{\alpha _{\rm{c}}} = 1 + {K_{{\rm{Fe}}}}\left[ {{\rm{F}}{{\rm{e}}^{3 + }}} \right]\\{{K'}_{\rm{t}}} = \frac{{{K_{\rm{t}}}}}{{{\alpha _{\rm{e}}}}} = \frac{{{K_{\rm{t}}}}}{{1 + {K_{{\rm{Fe}}}}\left[ {{\rm{F}}{{\rm{e}}^{3 + }}} \right]}}\\{{K'}_{\rm{t}}} = {\left[ {{\rm{C}}{{\rm{e}}^{4 + }}} \right]^2} + {{K'}_{\rm{t}}}{c_0}\left( {\alpha - 1} \right)/\left( {\alpha + 1} \right) + 1\left[ {{\rm{C}}{{\rm{e}}^{4 + }}} \right] - {c_0}/\left( {\alpha + 1} \right) = 0\end{array}$$\begin{array}{l}{\alpha _{\rm{e}}} = \frac{{\left[{{{\rm{e}}^{ - '}}} \right]}}{{\left[{{{\rm{e}}^ - }} \right]}} = \frac{{\left[{{{\rm{e}}^ - }} \right] + \left[{{\rm{F}}{{\rm{e}}^{2 + }}} \right]}}{{\left[{{{\rm{e}}^ - }} \right]}} = 1 + \frac{{\left[{{\rm{F}}{{\rm{e}}^{2 + }}} \right]}}{{\left[{{{\rm{e}}^ - }} \right]}}\\{K_{{\rm{Fe}}}} = \frac{{\left[{{\rm{F}}{{\rm{e}}^{2 + }}} \right]}}{{\left[{{\rm{F}}{{\rm{e}}^{3 + }}} \right]\left[{{{\rm{e}}^ - }} \right]}}\\\frac{{\left[{{\rm{F}}{{\rm{e}}^{2 + }}} \right]}}{{\left[{{{\rm{e}}^ - }} \right]}} = {K_{{\rm{Fe}}}}\left[{{\rm{F}}{{\rm{e}}^{3 + }}} \right]\\{\alpha _{\rm{c}}} = 1 + {K_{{\rm{Fe}}}}\left[{{\rm{F}}{{\rm{e}}^{3 + }}} \right]\\{K_t}' = \frac{{{K_{\rm{t}}}}}{{{\alpha _{\rm{e}}}}} = \frac{{{K_{\rm{t}}}}}{{1 + {K_{{\rm{Fe}}}}\left[ {{\rm{F}}{{\rm{e}}^{3 + }}} \right]}}\\{K_t}' = {\left[ {{\rm{C}}{{\rm{e}}^{4 + }}} \right]^2} + {\rm{ }}\left( {{K_t}'{c_0}\left( {\alpha - 1} \right)/\left( {\alpha + 1} \right) + 1} \right)\left[ {{\rm{C}}{{\rm{e}}^{4 + }}} \right] - {c_0}/\left( {\alpha + 1} \right) = 0\end{array}$

或

$\begin{array}{l}{{K'}_{\rm{t}}}{\left[ {{\rm{F}}{{\rm{e}}^{2 + }}} \right]^2} + {{K'}_{\rm{t}}}{c_0}\left( {\alpha - 1} \right)/\left( {\alpha + 1} \right) + 1\left[ {{\rm{F}}{{\rm{e}}^{2 + }}} \right] - {c_0}/\left( {\alpha + 1} \right) = 0\\\left[ {{\rm{C}}{{\rm{e}}^{4 + }}} \right]{\rm{或}}\left[ {{\rm{F}}{{\rm{e}}^{2 + }}} \right]{\rm{ = }}\frac{{ - {{K'}_{\rm{t}}}{c_0}\left( {\alpha - 1} \right)/\left( {\alpha + 1} \right) + 1 + \sqrt {{{K'}_{\rm{t}}}{c_0}\left( {\alpha - 1} \right)/\left( {\alpha + 1} \right) + {1^2} + 4{K_{\rm{t}}}{c_0}/\left( {\alpha + 1} \right)} }}{{2{{K'}_{\rm{t}}}}}\end{array}$

由于副反应系数α_e随滴定过程[Fe³⁺]的变化而变化，因此解析该滴定曲线方程较为复杂。本文利用一个简单的迭代α_e(i) = 1 + K_Fe[Fe³⁺](i)，计算出变化的α_e。利用MATLAB则能够准确地解析滴定曲线方程，并绘制了滴定曲线，程序如下：

a = 0:0.001:2;

kt = 10^(1.44/0.059);

kfe = 10^(0.68/0.059);

b = a - 1;

c = a + 1;

d = b./c;

Fe2(1) = 0.1;

for i = 2:1:2001

   alfa(i) = 1 + kfe*(0.1/c(i) - Fe2(i - 1));

   ktsh(i) = kt/alfa(i);

   e1(i) = ktsh(i)*0.1*d(i) + 1;

   f1(i) = e1(i)^2;

   g1(i) = sqrt(f1(i) + 4*ktsh(i)*0.1/c(i));

   Fe2(i) = (-e1(i) + g1(i))/(2*ktsh(i));

   Fe3(i) = 0.1/c(i) - Fe2(i);

   Ce4(i) = Fe2(i);

   Ce3(i) = Fe3(i);

   E(i) = 1.44 + 0.059*log10(Ce4(i)/Ce3(i));

   E1(i) = 0.68 + 0.059*log10(Fe3(i)/Fe2(i));

end

plot(a,E,'.',a,E1,'.')

所绘制的Ce⁴⁺与Fe²⁺相互滴定的滴定曲线如图5所示，可以看出由统一模型所绘制的曲线与文献^{[1, 5]}的理论计算结果相一致。

酸碱滴定与沉淀滴定曲线的统一模型为${K'_{\rm{t}}}{\left[ {\rm{M}} \right]^2} + {K'_{\rm{t}}}c_{\rm{M}}^0\frac{{\alpha - 1}}{{\alpha + 1}}\left[ {\rm{M}} \right] - 1{\rm{ = }}0$，络合滴定与氧化还原滴定曲线的统一模型为{K'_{\rm{t}}}{\left[ {\rm{M}} \right]^2} + {\rm{ }}\left( {{{K'}_{\rm{t}}}c_{\rm{M}}^0\left( {\alpha - 1} \right)/\left( {\alpha + 1} \right) + 1} \right)\left[ {\rm{M}} \right] - c_{\rm{M}}^0{\rm{ = }}0。通过简单迭代法计算滴定过程中变化的副反应系数，利用MATLAB解析滴定曲线统一模型，能够精确地绘制出4种滴定曲线，且方法简单易操作。