利用简单的化学及编程知识解决一类水溶液中各粒子浓度的计算

(相关资料图)

于是得出

那么所有存在形式的浓度均可以用 \(c(\ce{PO_4^{3-}})\) 和 \(c(\ce{H^+})\)表示。如果我们用想用 \(c(\ce{H3PO4})\) 来表示，则可以得到 \(c(\ce{PO_4^{3-}})=K_{a_1}K_{a_2}K_{a_3}\times {c(\ce{H^+})^{-3}}\times c(\ce{H3PO4})\)。

不妨设往溶液中加入的磷酸的浓度为 \(c_{总}\)，就用 \(c(\ce{PO4^{3-}})\) 作为未知数，利用物料守恒可以得到方程

显然我们如果知道了 \(c(\ce{H^+})\)，和 \(c_{总}\)，可以轻易地计算出 \(c(\ce{PO4^{3-}})\)，进而计算出其他离子的浓度。

那么我们就可以将四种含 \(\ce{P}\) 粒子的浓度随 \(p\ce H\) 变化的图像画出来，这就是在题目中常见的图。

如果溶液中还有其他的酸或碱，其各个离子的浓度也可以根据这种方法计算出来。

考虑电荷守恒

放在实际的溶液中，其电荷一定是守恒的。我们上方只得到了各粒子浓度关于氢离子浓度的函数，正负电荷并不一定守恒。于是我们现在需要计算出一个合适的氢离子浓度，使得溶液中电荷守恒，那么这种情况就是满足的。由于各粒子浓度均非负数，那么我们应该可以证明氢离子浓度的解是唯一的（我没有证明，感觉可能是唯一的）。

程序设计部分

如何计算出氢离子浓度

我们根据电荷守恒，可以得到一个一元高次方程。如果我们尝试整理系数并利用求根公式求解，那么在程序上可能会比较难以实现。

不妨尝试试出氢离子浓度，计算出各阴阳离子的浓度之后，直接比较等式左右电荷绝对值的差，如果这个差的绝对值小于一个数，我们可以认为这个氢离子浓度是正确的。

但是如果按照浓度逐一尝试，每次为尝试的氢离子浓度增加一定的浓度，那么如果我们要保证精度，就需要尝试非常多次。

可以感性理解出正负电荷量差的绝对值随氢离子浓度的变化是一个连贯的曲线，我们的目标是找到这个曲线上的最小值点，这个点函数值为 \(0\)。这个函数上面可能有无数极值点（还是我感觉的，也有可能是单调的，这个时候可以直接使用二分法解决，显然但粒子种类非常多的时候已经很难用手证明是否在整个定义域上都是单调的了），而我们需要找到那个最小值点。这个时候可以考虑使用模拟退火算法（退火是不是本身就是化学名词）。

模拟退火就是比较常规的写法，记得比较上一个位置的答案而不是目前已经得到的最优的答案（很多模拟退火的竞赛题都是直接比较最优的答案，但是这里不行）。

一个精度问题

还是以磷酸举例，如果我们以 \(c(\ce{PO4^{3-}})\) 为中间变量计算，如果 \(c(\ce{H^+})\) 的浓度非常大，那么得到的 \(c(\ce{PO4^{3-}})\) 会非常小，由于默认的小数存储结构精度有限，一些有效数字可能会被直接省略掉，造成较大的精度丢失。那么我们可以在计算之前先找出在一定氢离子浓度下占据主导地位的粒子，然后用直接用它作为中间变量计算。（也许这点误差真的能够忽略不计吧，但是还是在这里写出来吧）

另一个大问题

磷酸 \(K_{a_1}=6.9\times 10^{-3},K_{a_2}=6.2\times 10^{-8},K_{a_3}=4.8\times 10^{-13}\)，乘起来是一个非常小的数字，超过了 C 语言所能记录的浮点数精度。这意味着在计算过程中极易出现爆精度的问题，尽管可以采用先乘后除的办法，并且这些极小的数所代表的粒子浓度本身就可以忽略不计，但是还是会在一些玄学的地方带来误差，这样在比较两个浓度都非常小的粒子的时候就不能保证准确性了（我瞎猜说）。

解决的办法就是写自己写记录小数的数据结构（即高精度），但是这会浪费非常多的时间，并转化为一个彻头彻尾的代码问题。

实际上，程序中采用的课本上的平衡常数本身就可能存在着误差，所以在这里我们并不考虑解决这个问题。

一些目前版本的bug

• 对于在实际情况中浓度相同的粒子可能计算出的结果会有细微差异。
• 在我的电脑上运行程序的时候载入会有几秒钟的时间（这确实有点不太正常），可能是代码输入部分的实现问题。我不是专业的程序员，完全不会写 IO，请见谅。

一些省流

\(0.1\ce{mol\cdot L^{-1}}\ce{H3PO4}\) 溶液的计算结果：

pH: 1.6375c(H_3PO_4):    0.076956557415177c(H^+):        0.023043504584859c(H_2PO_4^-):  0.023043380585157c(HPO_4^2-):   0.000000061999666c(OH^-):       0.000000000000434c(PO_4^3-):    0.000000000000000

\(0.1\ce{mol\cdot L^{-1}}\ce{CH3COONa}\) 溶液的计算结果：

pH: 8.8785c(Na^+):       0.100000000000000c(CH_3COO^-):  0.099992441663165c(OH^-):       0.000007559670777c(CH_3COOH):   0.000007558336835c(H^+):        0.000000001322809

\(0.1\ce{mol\cdot L^{-1}}\ce{NH3\cdot H2O}\) 溶液的计算结果：

pH: 11.1247c(NH_3-H_2O):  0.098667174249750c(NH_4^+):     0.001332825750250c(OH^-):       0.001332514123592c(H^+):        0.000000000007505

代码

建立结构体记录每种弱酸或碱。结构体内有根据氢离子浓度计算出对应酸根离子们的浓度的函数。

#include #include #include using namespace std;typedef long double ld;typedef pairttfa;//这个 pair 记录答案中的离子名称及相应的浓度#define mp make_pairconst ld K_w=1e-14;mt19937 gen(chrono::system_clock::now().time_since_epoch().count());uniform_real_distributionrd(0,1);struct basecid{int n;ld c;//n 元酸或碱，以及浓度 cvectorname;//名字，按照电荷数的绝对值从小到大排vectorK;//电离平衡常数bool is_acid;//是酸(1)，还是碱(0)，如果是碱直接用将氢离子浓度转化成氢氧根离子浓度用同样的方法计算。inline ld getk(int loc,ld c_H){//具体原理已经讲过，系数*该粒子浓度=总离子浓度，其中系数和氢离子及选择的粒子有关，这里求出这个系数ld val=1,bas=1;for(int i=loc-1;i>=0;--i){bas=bas*c_H/K[i];val+=bas;}bas=1;for(int i=loc+1;i<=n;++i){bas=bas*K[i-1]/c_H;val+=bas;}return val;}inline ld e_it(ld c_H,ld c_all){//计算总电子数，c_all 如果为大于 0 则忽略目前结构体内的浓度，而使用这个浓度if(c_all<0)c_all=c;if(!is_acid)c_H=K_w/c_H;int key=0;ld minv=getk(0,c_H);for(int i=1;i<=n;++i){ld now=getk(i,c_H);if(now=0;--i){bas=bas*c_H/K[i];}bas=mid;for(int i=key+1;i<=n;++i){bas=bas*K[i-1]/c_H;val+=i*bas;}return is_acid?-val:val;}inline vector lis_it(ld c_H,ld c_all){//单独又写了一个函数，将对应氢离子浓度下的各离子的浓度得出if(c_all<0)c_all=c;if(!is_acid)c_H=K_w/c_H;int key=0;ld minv=getk(0,c_H);for(int i=1;i<=n;++i){ld now=getk(i,c_H);if(nowlis(n+1);lis[key]={name[key],mid};for(int i=key-1;i>=0;--i){bas=bas*c_H/K[i];lis[i]={name[i],bas};}bas=mid;for(int i=key+1;i<=n;++i){bas=bas*K[i-1]/c_H;lis[i]={name[i],bas};}return lis;}}book[100];struct normal{//在溶液中浓度不变的粒子string name;int e;ld c;inline ld e_it(ld c_all){if(c_all<0)c_all=c;return e*c_all;}};vectorother;vectormix;inline ld e_forall(double c_H){//计算出该氢离子浓度下的电荷数ld val=0;for(auto x:mix)val+=x.e_it(c_H,-1);for(auto x:other)val+=x.e_it(-1);return val+c_H-K_w/c_H;}inline ld SA(){//模拟退火，要记录 lastans。没怎么调参ld loc=1e-7,ans=abs(e_forall(loc)),lasans=ans;for(ld T=1;T>=1e-14;T*=0.995){ld nloc=max(loc+(rd(gen)*2-1.0)*T,(ld)1e-16);ld nval=abs(e_forall(nloc)),der=nval-lasans;//printf("%Lf %Lf %Lf ",nloc,nval,der);if(der<0)lasans=ans=nval,loc=nloc;else if(exp(-der/T)>rd(gen))lasans=nval,loc=nloc/*,printf("accept")*/;//puts("");}//printf("%.15Lf %.15Lf\n",loc,ans);return loc;}inline void initset(){//这个写法可能有些丑陋，可能降低了效率book[1].n=1,book[1].is_acid=1;book[1].K.push_back(1.75e-5);book[1].name.push_back("CH_3COOH");book[1].name.push_back("CH_3COO^-");book[2].n=1,book[2].is_acid=0;book[2].K.push_back(1.8e-5);book[2].name.push_back("NH_3-H_2O");book[2].name.push_back("NH_4^+");book[3].n=2,book[3].is_acid=1;book[3].K.push_back(4.5e-7);book[3].K.push_back(4.7e-11);book[3].name.push_back("H_2CO_3");book[3].name.push_back("HCO_3^-");book[3].name.push_back("CO_3^2-");book[4].n=3,book[4].is_acid=1;book[4].K.push_back(6.9e-3);book[4].K.push_back(6.2e-8);book[4].K.push_back(4.8e-13);book[4].name.push_back("H_3PO_4");book[4].name.push_back("H_2PO_4^-");book[4].name.push_back("HPO_4^2-");book[4].name.push_back("PO_4^3-");book[5].n=2,book[5].is_acid=1;book[5].K.push_back(1.4e-2);book[5].K.push_back(6.0e-8);book[5].name.push_back("H_2SO_3");book[5].name.push_back("HSO_3^-");book[5].name.push_back("SO_3^2-");book[6].n=1,book[6].is_acid=1;book[6].K.push_back(4.0e-8);book[6].name.push_back("HClO");book[6].name.push_back("ClO-");}inline bool cmp(ttfa x,ttfa y){//这个写法比较丑陋，但是我并不会比较高级的写法return x.second==y.second?x.firsty.second;}void forehead(){cout<<"English is used to avoid compatibility issues! (1.0 by BSE)"<>opt;if(opt=="auto"){int num;ld c;cin>>num>>c;basecid ac=book[num];ac.c=c;mix.push_back(ac);}else if(opt=="fixed"){normal fx;cin>>fx.name>>fx.e>>fx.c;other.push_back(fx);}else if(opt=="new"){basecid ac;cin>>ac.n>>ac.is_acid;for(int i=0;i>K;ac.K.push_back(K);}for(int i=0;i<=ac.n;++i){string tmp;cin>>tmp;ac.name.push_back(tmp);}cin>>ac.c;mix.push_back(ac);}else if(opt=="end")break;else cout<<"Non-conforming input!"<each;each.push_back(mp("H^+",c_H));each.push_back(mp("OH^-",K_w/c_H));for(auto ac:mix){vectortmp=ac.lis_it(c_H,-1);for(auto ppair:tmp){each.push_back(mp(ppair.first,ppair.second));}}for(auto li:other){//这一步虽然可以在输入的时候完成，但是还是在这里写比较清晰一些each.push_back(mp(li.name,li.c));}sort(each.begin(),each.end(),cmp);cout<<"pH: "<

后记

可能会改进一下算法，提升一下精度。最重要的是改良一下输入。

另外如果想要输入其他弱酸或弱碱，请按照下面的形式，看代码应该能知道怎么输入，下例子为 $ 0.1 \ce{mol\cdot L^{-1}}\ce{H2C2O4}$ 溶液。

new 2 15.6e-2 1.5e-4H_2C_2O_4 HC_2O_4^- C_2O_4^2-0.1

输出为

pH: 1.2836c(H^+):        0.052049217002032c(HC_2O_4^-):  0.051750936225208c(H_2C_2O_4):  0.048099923386503c(C_2O_4^2-):  0.000149140388288c(OH^-):       0.000000000000192

先写这么多吧，还有好多想写的没写。如果错误和漏洞请多多反馈，谢谢观看。\(\rm by \; BigSmall\_En\)