はじめに
1価の酸性溶液と1価の塩基性溶液を混ぜ合わせると水素イオンと水酸化物イオンが結合して水分子になる中和という反応が生じる。これによって混合溶液全体のpHは変化する。今回は、1価の酸性溶液に塩基性溶液を混合していった場合においてpHがどのように変化するのかをPythonを用いて数値解析を行う。具体的には、まず化学平衡式を立てることにより中和反応を理論モデル化し、水素イオンに対する方程式を立てる。つぎに、その方程式をニュートンラフソン法もしくは二分法という数値解析手法を用いることで解の近似値を求める。そしてそれをpHの定義($pH=-log_{10}[H^+]$)に沿ってプロットすることで以下のような中和滴定曲線を描写する。
中和滴定曲線のモデル方程式
以下のサイトのモデルを参考にした。
一価の弱酸$HA$に強塩基$BOH$を滴下していくという条件で考える。
電荷保存の法則より以下の等式が常に成立する。
[H^+]+[B^+]=[OH^-]+[A^-]
この式を$x=[H^+]$とした方程式に変形することを目的とする。
まず、一価の酸$HA$の化学平衡について考える。
\rm HA\rightleftharpoons H^{+}+A^{-}
このように、100%電離するとは限らないことに注意する。
これにより、化学平衡式は以下のようになる。
K_a=\rm \dfrac {\left[ H^{+}\right] \left[ A^{-}\right] }{\left[ AH\right] }
ただし、初期の酸の濃度$[AH]$は$C_A$とすると反応前後で原子の総量は変化しない（物質収支）より、
C_A=\rm \left[ A^{-}\right] +\left[ AH\right]
となる。したがって、
[A^{-}]=\frac{K_aC_A}{K_a+[H^{+}]}
1価の塩基$BOH$についても同様に考えると以下の式が成立する。
C_B=\rm \left[ B^{+}\right] +\left[ BOH\right]
ただし、$BOH$は強塩基なので近似的に以下の式で表すことができる。
\rm \left[ B^{+}\right]=C_B
ところで、水の電離の式は以下のように表すことができる。
K_w=\rm \left[ H^{+}\right] \left[ OH^{-}\right]
ゆえに、
\left[ OH^{-}\right]=\frac{K_w}{\left[ H^{+}\right]}
以上の考察により、$[B^+],[A^-],[OH^-]$の式を以下の式に代入する。
[H^+]+[B^+]=[OH^-]+[A^-]
したがって、
[H^+]+C_B=\frac{K_w}{\left[ H^{+}\right]}+\frac{K_aC_A}{K_a+[H^{+}]}
これにより、$[H^+]$について以下の3次方程式が成立する。
[H^+]^3+(K_a+C_B)[H^+]^2+(C_B K_a-K_a C_A-K_w)[H^+]-K_w K_a=0
ただし、混合前の酸性、塩基性溶液の初期濃度を$C_{A0},C_{B0}$とする。
また、酸性溶液の体積を$V_A$とし、滴定した塩基性溶液の体積を$V_B$とし、以下の希釈の前提条件で考える。
C_A=C_{A0}\frac{V_A}{V_A+V_B}
C_B=C_{B0}\frac{V_B}{V_A+V_B}
ところで、3次方程式を解の公式を用いて代数的にに解き厳密解を求めるのはとても難しいので、近似解を求めることを考える。
ニュートンラフソン法による方程式の数値解析
アルゴリズム
そこで、上記の$[H^+]$に関する3次方程式の近似解を数値解析(ニュートンラフソン法)によって求めることができるかについて考察する。
詳しい説明は以下の記事を参考にされたい。
本記事では概要のみ述べるものとする。
まず、実数$x$に対する以下の方程式が与えられているとする。
f(x)=0
この方程式を『常識的』に満たしそうな解を経験則的に初期解$x_1$として配置する。
今回の場合は、$x=[H^+]$なので、$1.0\times10^{-14}<x<1.0\times10^{-1}$程度の範囲で初期解を配置するとする。
そこで、以下の漸化式を新たに考える。
x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}
ただし、$f'(x)$は$|h|<<1$として以下のように表すことができるものとする。（微分の定義）
f'(x)=\frac{f(x+h)-f(x)}{(x+h)-x}=\frac{f(x+h)-f(x)}{h}
この漸化式を回して行き、$|f(x_n)|$が一定値以下になったとき、計算を終了するものとする。
プログラム例
例として$x^2-x=0$の近似解を求めるプログラムを作成すると、以下のようになる。
import numpy as np
import matplotlib.pyplot as plt
import math
import japanize_matplotlib
h=0.001
#このプログラムはx^2-x=0を解くプログラムである。
# （解の個数は分からない）
# x^2-x=0の解はx=0,「1」である。
def equation(x):
return x**2-x
#初期値を1.5とする
x=1.5
while abs(equation(x))>1.0*10**-5:
diff_equation = (equation(x+h)-equation(x))/h
x=x-equation(x)/diff_equation
print(x)
#1.000000195408536が出力される
プログラム
さて、これをもとに、中和滴定のプログラムを作成すると以下のようになる。
ただし、化学平衡に関する定数を以下のようにした。
K_w=[H^+][OH^{-}]=1.0\times 10^{-14}[{mol}^2/L^2],K_a=1.0\times 10^{-5}[mol/L]
プログラムコード
import numpy as np
import matplotlib.pyplot as plt
import math
import japanize_matplotlib
h=1.0*10**-10
n=100
K_w=10**-14
K_a=1.75*10**-5
C_A0=0.1
C_B0=0.2
V_A=10
def equation(x,V_A,V_B1):
C_A=(C_A0)*V_A/(V_A+V_B1)
C_B=(C_B0)*V_B1/(V_A+V_B1)
return x**3+(K_a+C_B)*x**2+(K_a*C_B-K_a*C_A-K_w)*x-K_w*K_a
V_B=np.linspace(0.1,10,n)
pH_array=[]
for i in range(n):
x=1.0
V_B1=V_B[i]
while abs(equation(x,V_A,V_B1))>1.0*10**-20:
diff_equation = (equation(x+h,V_A,V_B1)-equation(x,V_A,V_B1))/h
x=x-equation(x,V_A,V_B1)/diff_equation
pH=-math.log10(x)
pH_array.append(pH)
plt.xlabel("塩基溶液の体積")
plt.ylabel("pH")
plt.plot(V_B,pH_array)
plt.savefig("中和滴定曲線.png")
plt.show()
これを実行すると以下のようなグラフが出力される。
このように、中和点近くになるとpHが急激に変化する現象（pHジャンプ）もしっかりと再現できていることが分かる。
二分法によるプログラム
上記のプログラムは、電荷保存の式をこねくり回して出来た式を3次方程式に直して解析したものである。
ところで、3次方程式に直す前の以下の式を上手く利用することは出来ないだろうか？
[H^+]+C_B=\frac{K_w}{\left[ H^{+}\right]}+\frac{K_aC_A}{K_a+[H^{+}]}
したがって、
[H^+]+C_B-(\frac{K_w}{\left[ H^{+}\right]}+\frac{K_aC_A}{K_a+[H^{+}]})=0
これを満たす解を二分法を用いて見つける。
今回は、解の範囲が$0<x<1.0$程度になると予測が付けられるので、最大値を1として最小値を0とした。
なぜ、ニュートンラフソン法ではなく二分法を用いたかは、単純に$x>0$という条件を常に成立させておかないと、対数を含むpH計算でエラーになってしまうからである。これは、ニュートンラフソン法のデメリットである。つまり、ある範囲のみ探索するということがニュートンラフソン法では非常に難しい。一方で二分法なら、初期値は2つ必要だが、そこで区切られた領域のみの探索が可能である。
2分法について
以下、2分法について概略を述べる。
区間$a<x<b$において、単調増加もしくは単調減少を示す関数$f(x)$が存在するとする。
f(x_{ex})=0
となる、$x_{ex}$が区間$a<x<b$に存在するとする。
このとき、
f(x_{min})f(x_{max})<0
を満たすとき、
x_{min}<x_{ex}<x_{max}
となる。したがって、
x_{max}=\frac{x_{min}+x_{max}}{2}
と新たにおいた$x_{max}$において、
f(x_{min})f(x_{max})<0
が成立すれば、
x_{min}<x_{ex}<x_{max}
を満たすことになる。
この操作を繰り返していき、誤差関数が一定値以下になったら試行を終了する。これによって方程式の近似値を求めることができる。
プログラム
import numpy as np
import matplotlib.pyplot as plt
import math
import japanize_matplotlib
h=1.0*10**-15
n=100
K_w=10**-14
K_a=1.75*10**-5
C_A0=0.1
C_B0=0.2
V_A=10
def equation(x,V_A,V_B1):
C_A=(C_A0)*V_A/(V_A+V_B1)
C_B=(C_B0)*V_B1/(V_A+V_B1)
return (x+C_B-(K_w/x+K_a*C_A/(K_a+x)))
V_B=np.linspace(0.1,10,n)
pH_array=[]
for i in range(n):
#初期値
## 最大値
x_p=1.0
## 最小値
x_q=0
#塩基溶液の体積
V_B1=V_B[i]
#二分法
x=(x_p+x_q)/2
while abs(equation(x,V_A,V_B1))>10**-5:
x=(x_p+x_q)/2
if equation(x,V_A,V_B1)*equation(x_p,V_A,V_B1)<0:
x_q=x
else:
x_p=x
pH=-math.log10(x)
pH_array.append(pH)
plt.xlabel("塩基溶液の体積")
plt.ylabel("pH")
plt.plot(V_B,pH_array)
plt.savefig("中和滴定曲線.png")
plt.show()
これを実行すると以下のようなグラフが出力される。
まとめ
今回は、1価の酸性、塩基性溶液の中和滴定について、化学平衡モデルを考えて、イオン濃度に対する方程式を立式した。そして、その方程式は3次方程式になった。これの厳密解を求めるのは難しいため、ニュートンラフソン法を用いて近似解を算出し、その結果をプロットすることで、中和滴定曲線を描写することができた。その曲線は俊敏な挙動を示すpHジャンプも再現できていたので、かなり精度が高いといえる。ただし、今後の課題として水素イオン濃度の初期解を経験則的に与えているので、それによっては精度が悪化してしまう可能性がある。また、ニュートンラフソン法では、正の値のみの探索を行うことが非常に難しいため、今回の場合は、範囲指定ができる2分法が適していると考えられる。