在分析化学实验设计中，需对物质选取、用量、溶液pH调节等做定量分析，常用的近似公式、经验判别准则具有一定的局限性；而在精确求解过程中，常会遇到高次方程，往往只能得到部分点值关系。在分析化学中引入其他分析工具，如MATLAB已有较多实例。相较于MATLAB，Python作为一种编程语言，能自行决定更多的算法细节、可以实现更多功能，同时Python开源、免费，有更大的优势。

分析化学中涉及的等量关系常是幂次关系x^αy^β，如平衡常数(K_a、K_b、K_sp、K_稳)、电荷守恒(CBE)、物料守恒(MBE)，几乎不涉及其他复杂函数关系。对于某一特定化学反应体系，借助平衡常数，通过消元可以得到一个二元等量关系，将其转化为单变元函数后或直接利用隐式关系，可通过Python绘制其关系曲线图，从而能更加直观地分析各变量之间关系。

由化学反应体系列出各物质的MBE和体系的CBE (理论上，CBE可由MBE推导而来，但CBE更加清晰直观)；利用平衡常数，消元可得所研究的两变量之间的等量关系，将其转化为Python程序语言，便可绘制关系曲线、求值。

四大滴定的滴定曲线绘制、终点误差计算等在中国科学技术大学教材《分析化学》^[1]已有详细介绍。下面以磷酸钙中Bi³⁺、Ca²⁺混合离子滴定的指示剂选取、滴定突跃、理论终点误差为例^[2]，简要分析如何利用Python处理分析化学中的相关问题。测定流程如图1所示。

沉淀完磷酸根离子后，溶液c_Bi³⁺ = 0.01900 mol·L⁻¹，c_Ca²⁺ = 0.01000 mol·L⁻¹，pH = 2.0，此时EDTA的酸效应系数为${\left. {{\alpha _{\rm{Y}}}} \right|_{{\rm{pH = }}2.0}} = {10^{13.51}}$，二甲酚橙理论变色点XO：${\left. {{\rm{pB}}{{\rm{i}}_{\rm{t}}}} \right|_{{\rm{pH = }}2.0}} = 5.40$，颜色由紫红变为亮黄。Ca²⁺、Bi³⁺与EDTA结合的稳定常数为K_CaY = 4.90 × 10¹⁰、K_BiY = 8.71 × 10²⁷。

取20.00 mL待测液，加入60.00 mL蒸馏水，以0.02000 mol·L⁻¹的EDTA滴定。

滴定体系的MBE:

(1) ${\rm{Ca}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left[ {{\rm{CaY}}} \right] + \left[ {{\rm{Ca}}} \right] = {c_{{\rm{Ca}}}} = 0.20/(V + {\rm{ }}80)$

(2) ${\rm{Bi}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left[ {{\rm{BiY}}} \right] + \left[ {{\rm{Bi}}} \right] = {c_{{\rm{Bi}}}} = 0.38/(V + {\rm{ }}80)$

(3) ${\rm{EDTA}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left[ {{\rm{CaY}}} \right] + \left[ {{\rm{BiY}}} \right] + [{\rm{Y''}}] = {c_{{\rm{EDTA}}}} = 0.02V/(V + 80)$

设$[{\rm{Y}}'']$ = [Y]α_Y(H)、${K_1}''$ = K_CaY/α_Y(H)、${K_2}''$ = K_BiY/α_Y(H)

(4)$由可得：[{\rm{Ca}}] = \frac{1}{{1 + [{\rm{Y}}'']{K_1}''}}\frac{{0.2}}{{V + 80}}$

(5)$由可得：[{\rm{Bi}}] = \frac{1}{{1 + [{\rm{Y}}'']{K_2}''}}\frac{{0.38}}{{V + 80}}$

将(4)、(5)代入(3)中可得：$\left( {\frac{{0.2[{\rm{Y}}'']{K_1}''}}{{1 + [{\rm{Y}}'']{K_1}''}} + \frac{{0.38[{\rm{Y}}'']{K_2}''}}{{1 + [{\rm{Y}}'']{K_2}''}}} \right)\frac{1}{{V + 80}} + [{\rm{Y}}''] = \frac{{0.02V}}{{V + 80}}$

(6) $解出V：V = \left( {\frac{{0.2[{\rm{Y}}'']{K_1}''}}{{1 + [{\rm{Y}}'']{K_1}''}} + \frac{{0.38[{\rm{Y}}'']{K_2}''}}{{1 + [{\rm{Y}}'']{K_2}''}} + 80[{\rm{Y}}'']} \right)\frac{1}{{0.02 - [{\rm{Y}}'']}}$

将上述等量关系转化为Python程序：

1. import matplotlib.pyplot as plt

2. import numpy as np

3. K1 = 4.9^*10^**10/(10^**13.51)     #K1''

4. K2 = 8.71^*10^**27/(10^**13.51)     #K2''

5. Vmin = 0     #加入的EDTA体积范围

6. Vmax = 40

7. py = np.linspace(-20, -2, 100000)      # py = -lg(y)，y or y_为[Y'']，作为自变量

8. y_ = []

9. V = []

10. pBi = []

11.

12. for i in py:      #计算V, pBi

13.   y_.append(10^**i)

14.   y = y_[-1]

15.   V.append((0.2^*K1^*y/(1+K1^*y)+0.38^*K2^*y/(1+K2^*y)+80^*y)/(0.02-y))  #即(6)式

16.    pBi.append(-np.log10(0.38/(V[-1]+80)/(1+K2^*y)))     # pBi = -lg(Bi³⁺)，即(5)式

17.

18. def pfind(x, x0):     #查找求值

19.  for i in range(0, len(x)):

20.    if x[i] > = x0:

21.      return i

22.  print("No result!")

23.

24. print("pBi: ", pBi[pfind(V, 19.00)])      #理论终点

25. print("Et: ", (V[pfind(pBi, 5.40)]/19-1)^*100)     #以二甲酚橙为指示剂的理论终点误差

26. print("pBiJump: ", pBi[pfind(V, 18.98)], pBi[pfind(V, 19.02)])

27.          #滴定突跃:19.00mL(1±0.1%)对应的pBi

28.

29. plt.plot(V, pBi, 'black')     #以V为横坐标，pBi为纵坐标，绘制pBi-V关系曲线

30. plt.xlim(Vmin, Vmax)

31. plt.xlabel("＄V＄/mL")

32. plt.ylabel(“＄pBi＄”)

33. plt.show()

计算结果：

滴定突跃(即V_sp·(1 ± 0.1%)范围内的pBi)：$ {p\mathrm{B}\mathrm{i}}_{\mathrm{J}\mathrm{u}\mathrm{m}\mathrm{p}} $= [5.42, 11.43]，在此范围内的指示剂均可。

理论终点(即V_sp = 19.00 mL时的pBi)：pBi = 8.42，二甲酚橙理论变色点XO：${\left. {p{\rm{B}}{{\rm{i}}_{\rm{t}}}} \right|_{{\rm{pH = }}2.0}} = 5.40$。

二甲酚橙变色点在突跃下端附近，并不是最理想的指示剂，若以二甲酚橙为指示剂的理论终点误差：(V_ep和V_sp分别表示理论终点和化学计量点时加入滴定剂的体积)

$Et = \frac{{{V_{{\rm{ep}}}} - {V_{{\rm{sp}}}}}}{{{V_{{\rm{sp}}}}}} \times 100\% = - 0.1\% $，滴定曲线如图2所示。

由滴定曲线图2可知，EDTA滴定混合离子，铋离子有明显突跃，终点误差较小，能准确滴定Bi³⁺。

由以上分析可知，运用Python分析，只需要将等量关系转化为Python语言(大部分是直接输入公式)，并不需要复杂程序结构，所用到的函数较少而且比较固定(利用numpy、matplotlib库)；结果直观、便于分析，同时能精确求解滴定突跃、理论变色点、终点误差。而对于传统近似算法，需要考虑公式适用条件，有较多繁杂数学计算；当体系改变时需要重新计算，过程较繁琐。若通过传统方法得到滴定曲线，需要计算不同EDTA体积时的pBi，过程繁杂。上述程序的重复利用性较好：对于EDTA滴定的大部分问题均有效，当体系的参数改变时，如溶液pH、金属离子浓度、金属离子种类、指示剂种类，只需修改部分常数值，重新运行程序即可，而不用从头计算。

在铬废液回收处理实验设计中^[3]，首先需要考虑Cr(Ⅵ)的还原，下面就还原剂选择及pH的控制做简要分析：

Cr(Ⅵ)在不同pH下的存在形式不同，氧化还原电对不同，主要存在以下平衡：

H₂CrO₄的电离：K_a1 = 0.18，K_a2 = 3.2 × 10⁻⁷

(7) $2\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{2-}+2{\mathrm{H}}^{+}\rightleftharpoons {\mathrm{C}\mathrm{r}}_{2}{\mathrm{O}}_{7}^{2-}+{\mathrm{H}}_{2}\mathrm{O} \ \ \ \ \ {K}_{3} = {10}^{14.69} $

(8) $体系\text{MBE}： {c}_{\mathrm{C}\mathrm{r}}=\left[{\mathrm{H}}_{2}\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}\right]+\left[\mathrm{H}\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{-}\right]+\left[\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{2-}\right]+2\left[{\mathrm{C}\mathrm{r}}_{2}{\mathrm{O}}_{7}^{2-}\right] $

Cr(Ⅵ)各形态浓度：

(9) ${K_{{\rm{a2}}}} = \frac{{[{\rm{CrO}}_4^{2 - }][{{\rm{H}}^{\rm{ + }}}]}}{{{\rm{[HCrO}}_4^ - ]}}，记{K_1} = \frac{1}{{{K_{{\rm{a2}}}}}} = {10^{6.49}} \Rightarrow {\rm{[HCrO}}_4^ - ] = {K_1}[{\rm{CrO}}_4^{2 - }][{{\rm{H}}^{\rm{ + }}}]$

(10) ${K_{{\rm{a1}}}}{K_{{\rm{a2}}}} = \frac{{[{\rm{CrO}}_4^{2 - }]{{[{{\rm{H}}^{\rm{ + }}}]}^2}}}{{{\rm{[}}{{\rm{H}}_2}{\rm{CrO}}_4^{}]}}，记{K_1} = \frac{1}{{{K_{{\rm{a1}}}}{K_{{\rm{a2}}}}}} = {10^{7.24}} \Rightarrow {\rm{[}}{{\rm{H}}_2}{\rm{CrO}}_4^{}] = {K_2}[{\rm{CrO}}_4^{2 - }]{[{{\rm{H}}^{\rm{ + }}}]^2}$

(11) $\left[{\mathrm{C}\mathrm{r}}_{2}{\mathrm{O}}_{7}^{2-}\right]={K}_{3}{\left[\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{2-}\right]}^{2}{\left[{\mathrm{H}}^{+}\right]}^{2} $

将(9)、(10)、(11)代入式(8)得：

(12) ${c}_{\mathrm{C}\mathrm{r}}={K}_{2}\left[\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{2-}\right]{\left[{\mathrm{H}}^{+}\right]}^{2}+{K}_{1}\left[\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{2-}\right]\left[{\mathrm{H}}^{+}\right]+\left[\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{2-}\right]+2{K}_{3}{\left[\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{2-}\right]}^{2}{\left[{\mathrm{H}}^{+}\right]}^{2} \\记a = 2{K}_{3}{\left[{\mathrm{H}}^{+}\right]}^{2} ，b = {K}_{2}{\left[{\mathrm{H}}^{+}\right]}^{2}+{K}_{1}\left[{\mathrm{H}}^{+}\right]+1 \\则，{c}_{\mathrm{C}\mathrm{r}}=a{\left[\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{2-}\right]}^{2}+b\left[\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{2-}\right] $

${\rm{CrO}}_4^{2 - }分布分数：δ({\rm{CrO}}_4^{2 - }) = \frac{{[{\rm{CrO}}_4^{2 - }]}}{{{c_{{\rm{Cr}}}}}} = \frac{1}{{a[{\rm{CrO}}_4^{2 - }] + b}}$

此处与常规酸碱分布分数不同，酸碱电离一般只涉及H⁺、OH⁻数目变化，如磷酸各形态都只含1个P，其分布分数与本身浓度无关；由于式(7)，Cr(Ⅵ)各形态之间有Cr数目的变化，会出现关于Cr高于一次的等量关系，分布分数与Cr浓度直接相关。

设c_Cr = 1.000 mol·L⁻¹，由式(12)可求得[$ \mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{2-} $] (此时也等于其分布分数)与[$ {\mathrm{H}}^{+} $]关系，

(13)$[{\rm{CrO}}_4^{2 - }] = \frac{{ - b + \sqrt {{b^2} + 4a} }}{{2a}}$

再由式(9)、(10)、(11)可求得$ \left[{\mathrm{H}}_{2}\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}\right]\mathrm{、}\left[\mathrm{H}\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{-}\right]\mathrm{、}\left[{\mathrm{C}\mathrm{r}}_{2}{\mathrm{O}}_{7}^{2-}\right] $与[$ {\mathrm{H}}^{+} $]关系。

绘制各物质分布分数与pH的关系曲线如图3所示(程序见补充材料1.1)：

图3清晰地展示了不同pH下Cr(Ⅵ)的存在形式(此处pH = 1.0–6.0时，Cr(Ⅵ)的主要存在形式不止一种)，再根据不同形态下氧化还原电对的电极电势便可判断与各种还原剂反应的能力以及最适pH范围，从而为还原剂选择提供指导。

由以上分析可发现，利用平衡常数和体系MBE可以得到各物质形态分布分数与[H⁺]的关系，通过Python绘制其关系曲线，从而对各物质形态的分布有更加直观的了解，为后续分析提供指导。在传统方法中，难以对全貌有很清晰的刻画，若要得到分布分数与pH的关系曲线，需要进行大量的单点计算，过程十分繁琐。

前文指出：由式(7)的特殊性，Cr(Ⅵ)分布分数与Cr浓度直接相关，记Cr浓度为c (0.1000–6.000 mol·L⁻¹)，

由式(12)：

(14)$[{\rm{CrO}}_4^{2 - }] = \frac{{ - b + \sqrt {{b^2} + 4ac} }}{{2a}} \Rightarrow \delta ({\rm{CrO}}_4^{2 - }) = \frac{{ - b + \sqrt {{b^2} + 4ac} }}{{2ac}}$

再由式(9)、(10)、(11)可求得$ \left[{\mathrm{H}}_{2}\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}\right]、\left[\mathrm{H}\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{-}\right]、\left[{\mathrm{C}\mathrm{r}}_{2}{\mathrm{O}}_{7}^{2-}\right] $与[$ {\mathrm{H}}^{+} $]关系。

绘制Cr(Ⅵ)分布分数与pH、c关系曲面如图4所示(程序见补充材料1.2)。

由图4，不同浓度下，Cr(Ⅵ)分布分数与pH关系曲线极为相似；同一pH下，Cr(Ⅵ)分布分数随Cr浓度变化较大(即图中条纹)。当浓度较小时，Cr(Ⅵ)不同形态的分布分数变化较大。取pH = 7.0时的Cr(Ⅵ)分布分数随Cr浓度关系曲线如图5所示(程序见补充材料1.3)。

从图5中可以看出，当pH = 7.0时，$ {\mathrm{C}\mathrm{r}}_{2}{\mathrm{O}}_{7}^{2-} $分布分数随Cr浓度增大而增大，在较高浓度下为主要存在形式，在极小浓度下约为0；$ \mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{2-}\mathrm{、}\mathrm{H}\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{-} $分布分数随Cr浓度减小而增大，在较小浓度下为主要存在形式。

由图3、4中可以发现，$ {\mathrm{C}\mathrm{r}}_{2}{\mathrm{O}}_{7}^{2-}、\mathrm{H}\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{-} $分布分数较相似。由$ {\mathrm{H}}_{2}\mathrm{C}\mathrm{r}{\mathrm{O}}_{4} $的二级电离和式(7)：

$2\mathrm{H}\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{-}\rightleftharpoons \left(2{\mathrm{H}}^{+}+2\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{2-}\right)\rightleftharpoons {\mathrm{C}\mathrm{r}}_{2}{\mathrm{O}}_{7}^{2-}+{\mathrm{H}}_{2}\mathrm{O}$

$K = \frac{{[{\rm{C}}{{\rm{r}}_2}{\rm{O}}_7^{2 - }]}}{{{{{\rm{[HCrO}}_4^ - ]}^2}}} = {K_{{\rm{a2}}}}{K_{\rm{3}}} = 50.15$

$ {\mathrm{C}\mathrm{r}}_{2}{\mathrm{O}}_{7}^{2-}、\mathrm{H}\mathrm{C}\mathrm{r}{\mathrm{O}}_{4}^{-} $间的平衡关系与pH无关，相互关联，所以有图3、4中的相似性，浓度较大时$ {\mathrm{C}\mathrm{r}}_{2}{\mathrm{O}}_{7}^{2-} $分布分数更大，浓度极小时关系反转，在图4、5中有所体现。

利用Python将变量之间的关系可视化，结果清晰直观，可以从全局把握变量间的关系，不局限于一点、一线，这是传统分析方法难以企及的。

相对于传统方法而言，利用Python分析不需要记忆大量的近似公式及其适用条件，只需要对不同体系的基本原理有所了解，通过MBE、CBE、平衡常数建立等量关系，再将其转化为研究变量间的函数关系，便可通过Python求解和绘制相关曲线。利用Python分析相关问题能有效减少繁杂计算，同时能对问题有更加深入、全面的认识，结果更加清晰直观，从而为分析化学实验设计提供更好的决策依据。

补充材料：可通过链接http://www.dxhx.pku.edu.cn免费下载。