该方程化简后将是四次方程，使用公式求解并不现实。我们使用牛顿法进行求解。这里不详细讲牛顿法的原理了，只给出公式：
（令方程为f(x)=0）
x=x0-f(x0)/f’(x0)
其中，式子右侧的x0是方程的原近似解，左侧的x是经过此步牛顿法求解后新的方程近似解，相比x0更接近方程的真实解。经过数步迭代后，新近似解x与原近似解x0的差值（对于pH计算而言，可以用相对差值代替绝对差值）小于一个事先指定的，非常接近0的数时，即可以认为已经得到了方程精度极高的数值解。
原方程左边对[H+]求导，可得（为方便起见令x=[H+]）：
f’(x)=1+Kw/(x^2)-(2/(20+V))*((((-Ka1/(x^2)-4*Ka1*Ka2/(x^3))*(1+Ka1/x+Ka1*Ka2/(x^2)))-((Ka1/x+2*Ka1*Ka2/(x^2))*(-Ka1/(x^2)-2*Ka1*Ka2/(x^3))))/((1+Ka1/x+Ka1*Ka2/(x^2))^2))
然后按上面的公式就可以进行计算了。为方便起见，我写了一个计算程序，源代码会在后面贴出。

要想计算更为准确，应当考虑离子强度的影响。本次讨论中，离子的活度系数使用德拜-休克尔公式计算：
γ=10^(-0.512*z^2*I^0.5/(1+0.00328*a*I^0.5))
a为以pm为单位的离子半径。I为离子强度，本例中公式为：
I=0.5*([H+]+[Na+]+[OH-]+[HC2O4-]+4*[C2O42-])
根据二元弱酸分布系数公式，得出：
[HC2O4-]=2/((20+V)*([H+]/Ka1+1+Ka2/[H+]))
[C2O42-]=2/((20+V)*([H+]^2/(Ka1*Ka2)+[H+]/Ka2+1))
浓度平衡常数会受离子强度的影响（即“盐效应”），具体表现为：
Ka1（浓度平衡常数）=Ka1（活度平衡常数）/(γ(HC2O4-)*γ(H+))
Ka2（浓度平衡常数）=Ka2（活度平衡常数）*γ(HC2O4-)/(γ(C2O42-)*γ(H+))
Kw（浓度平衡常数）=Kw（活度平衡常数）/(γ(H+)*γ(OH-))
这就陷入了一个死循环：想要解出各物种的平衡浓度需要知道离子强度值，但是想要知道离子强度值就需要知道各物种的平衡浓度。为此，我们先假定一个离子强度的值（例如I=0mol/L），以此值计算出各物种平衡浓度，再以得出的各物种平衡浓度计算新的离子强度值，之后再计算新一轮的各物种平衡浓度……直到某一轮时，新计算出的离子强度与上一轮的离子强度差值小于某一个接近0的数时，就认为已经计算出了体系真实的离子强度和各物种平衡浓度。
最后，考虑离子强度时pH的定义为：pH=-lg(γ(H+)*[H+])。

代码是python写的。我写的代码比较丑陋，大家凑合一下看。两个收敛限分别设置的10^(-10)和10^(-6)，第一个收敛限为相对偏差。
import math //载入数学模块以计算对数
V = 0.00000000E+00 //加入NaOH溶液的体积，可自行设定
x = 1.00000000E-14 //氢离子浓度的初猜
I = 0.00000000E+00
Ka1 = 5.62341325E-02
Ka2 = 5.12861383E-05
Kw = 1.00000000E-14
yH = 1.00000000E+00
yOH = 1.00000000E+00
yHC2O4 = 1.00000000E+00
yC2O4 = 1.00000000E+00
pH = 7.00000000E+00
ypH = 7.00000000E+00
x1 = 0.00000000E+00
a = 0.00000000E+00
b = 0.00000000E+00
c = 0.00000000E+00
d = 0.00000000E+00
f = 0.00000000E+00
g = 0
h = 0
print("加入滴定剂体积",V)
print("氢离子浓度",x)
print("离子强度",I)
print("第一电离常数",Ka1)
print("第一电离常数的负对数",-(math.log(Ka1))/(math.log(10)))
print("第二电离常数",Ka2)
print("第一电离常数的负对数",-(math.log(Ka2))/(math.log(10)))
print("水的离子积",Kw)
print("水的离子积的负对数",-(math.log(Kw))/(math.log(10)))
print("氢离子活度常数",yH)
print("氢氧根离子活度常数",yOH)
print("草酸氢根离子活度常数",yHC2O4)
print("草酸根离子活度常数",yC2O4)
print("浓度定义下的pH",pH)
print("活度定义下的pH",ypH)
for j in range(64): //离子强度迭代，大循环
for i in range(64): //牛顿法解方程，小循环
print("Iteration: ",i+1)
print("********************")
print("原氢离子浓度：",x)
a = x+0.1*V/(20+V)-Kw/x-(Ka1/x+2*Ka1*Ka2/(x**2))*2/((1+Ka1/x+Ka1*Ka2/(x**2))*(20+V))
b = 1+Kw/(x**2)-(2/(20+V))*((((-Ka1/(x**2)-4*Ka1*Ka2/(x**3))*(1+Ka1/x+Ka1*Ka2/(x**2)))-((Ka1/x+2*Ka1*Ka2/(x**2))*(-Ka1/(x**2)-2*Ka1*Ka2/(x**3))))/((1+Ka1/x+Ka1*Ka2/(x**2))**2))
c = x-a/b
print("迭代后氢离子浓度：",c)
d = ((c-x)**2)**(1/2)
print("绝对误差：",d)
print("相对误差：",d*100/c,"%")
if d/c < 1.00000000E-10:
x = c
break
x = c
if i == 63: //如果出错，则退出程序
print("****ERROR:Converged Failed.")
print("****错误:收敛失败.")
h = 1
if h == 1:
break
exit(0)
print("")
print("!!!!氢离子浓度为：",x)
pH = -(math.log(x))/(math.log(10))
print("浓度定义下的pH",pH)
f = 0.5*(x+0.1*V/(20+V)+Kw/x+2/((20+V)*(x/Ka1+1+Ka2/x))+(4*2)/((20+V)*((x**2)/(Ka1*Ka2)+x/Ka2+1)))
yH = 10**(-0.512*((f)**(1/2)/(1+0.00328*900*(f)**(1/2))))
ypH = -(math.log(x*yH))/(math.log(10))
print("活度定义下的pH",ypH)
print("原离子强度:",I)
print("本轮迭代后离子强度",f)
print("delta I:",((f-I)**2)**(1/2))
if ((f-I)**2)**(1/2) < 1.00000000E-06:
I = f
break
I = f
yOH=10**(-0.512*((I)**(1/2)/(1+0.00328*300*(I)**(1/2))))
yHC2O4=10**(-0.512*((I)**(1/2)/(1+0.00328*500*(I)**(1/2))))
yC2O4=10**(-0.512*4*((I)**(1/2)/(1+0.00328*500*(I)**(1/2))))
Kw=10**(-14)/((10**(-0.512*((I)**(1/2)/(1+0.00328*900*(I)**(1/2)))))*(10**(-0.512*((I)**(1/2)/(1+0.00328*300*(I)**(1/2))))))
Ka1=10**(-1.25)/((10**(-0.512*((I)**(1/2)/(1+0.00328*900*(I)**(1/2)))))*(10**(-0.512*((I)**(1/2)/(1+0.00328*500*(I)**(1/2))))))
Ka2=10**(-4.29)/((10**(-0.512*((I)**(1/2)/(1+0.00328*900*(I)**(1/2)))))*(10**(-0.512*4*((I)**(1/2)/(1+0.00328*500*(I)**(1/2))))))
if j != 63:
print("")
print("已设置好平衡常数，进入下一圈离子强度迭代")
print("")
if j == 63:
print("")
print("离子强度迭代未能在指定圈数内收敛，请谨慎使用输出结果！")
g = 1
print("")
print("体系已达到化学平衡.")
print("")
print("")
print("考虑离子强度的pH值已计算完毕 考虑离子强度的pH值已计算完毕 考虑离子强度的pH值已计算完毕 考虑离子强度的pH值已计算完毕 考虑离子强度的pH值已计算完毕")
print("")
print("")
print("加入滴定剂体积",V)
print("氢离子浓度",x)
print("化学平衡状态下的离子强度",I)
print("经化学平衡状态下离子强度校正后的第一电离常数",Ka1)
print("经化学平衡状态下离子强度校正后的第一电离常数的负对数",-(math.log(Ka1))/(math.log(10)))
print("经化学平衡状态下离子强度校正后的第二电离常数",Ka2)
print("经化学平衡状态下离子强度校正后的第二电离常数的负对数",-(math.log(Ka2))/(math.log(10)))
print("经化学平衡状态下离子强度校正后的水的离子积",Kw)
print("经化学平衡状态下离子强度校正后的水的离子积的负对数",-(math.log(Kw))/(math.log(10)))
print("化学平衡状态下的氢离子活度常数",yH)
print("化学平衡状态下的氢氧根离子活度常数",yOH)
print("化学平衡状态下的草酸氢根离子活度常数",yHC2O4)
print("化学平衡状态下的草酸根离子活度常数",yC2O4)
print("浓度定义下的pH",pH)
print("活度定义下的pH",ypH)
if g == 1:
print("离子强度迭代未能在指定圈数内收敛，请谨慎使用输出结果！")
print("")
print("程序运行完成.")

将计算结果作图，得：（离子强度进行了体积校正，即在算出的离子强度的基础上乘(20+V)/20）


与分析化学课程中的结论“草酸只能一步滴定掉两个氢离子”符合。离子强度的变化规律近似于直线，在化学计量点处斜率发生变化。
如果想要区分开草酸的两个氢离子使得它们能被分步滴定，有一个方法就是电导滴定，即根据滴定过程中溶液的电导变化推知化学计量点的位置。详细的在这里不介绍了。

看不懂帮顶，但你确定吧里会有人看吗

d

d

最严谨的一集