注意
以前のこちらの記事では NumPy を利用した秤量計算について解説しました。
本記事では SymPy を利用した秤量計算について解説します。NumPyの場合と比較してより一般性の高い場合に適用できますので、もし利用する場合には以前のNumPy版よりも本記事で紹介するコードの方を推奨します。
・2023/10/07 追記
最終的なPythonコードにおいて、 表示される質量が正しくない場合があるという重大な欠陥 が見つかったので以下のように修正しました。
変更前
reactant_mass = [solution * weight for (solution, weight) in zip(solutions.values(), reactant_weights)]
変更後
reactant_mass = [solutions[symbol] * weight for (symbol, weight) in zip(symbols, reactant_weights)]
前提
複数種類の試薬（反応物）から生成物を作りたいときに、各試薬をどれくらいの質量で混ぜればいいのかを計算する状況を考えます。例えば、
生成物
| 化学式 | 質量 |
|---|---|
| $ \mathrm{La} \mathrm{O} _{1-x} \mathrm{F} _{x} \mathrm{Bi} _{1-y} \mathrm{Pb} _{y} \mathrm{S} _2 $ | $\mathrm{2 \ g}$ |
試薬
| 化学式 | 分子量 |
|---|---|
| $ \mathrm{La_2 S_3} $ | 373.99 |
| $ \mathrm{Bi_2 O_3} $ | 465.958 |
| $ \mathrm{Bi_2 S_3} $ | 514.14 |
| $ \mathrm{Pb F_2} $ | 245.2 |
| $ \mathrm{Bi F_3} $ | 265.975 |
| $ \mathrm{Bi} $ | 208.98 |
としたときに、
- 生成物に合う各試薬の分子数の比を求め、
- 分子数あたりの質量である分子量と生成物の質量を使って各試薬の質量を求める。
という2stepで計算することになります。2番目は簡単ですが、1番目では連立方程式を解く必要性がでてきます。
1番目はどういうことかというと、
x_1 \mathrm{La_2 S_3} + x_2 \mathrm{Bi_2 O_3} + \cdots + x_6 \mathrm{Bi} \rightarrow \mathrm{La} \mathrm{O} _{1-x} \mathrm{F} _{x} \mathrm{Bi} _{1-y} \mathrm{Pb} _{y} \mathrm{S} _2
が成立するような$x_1, x_2, \cdots, x_6$を求めるということです。
$\mathrm{La}$に注目すると、左辺、右辺それぞれの$\mathrm{La}$の個数は$2 x_1$、$1$であり、これが等しくならなければならないため、$2 x_1 = 1$が成立します。他の試薬についても行うと、
\left\{ \,
\begin{aligned}
2 x_1 &= 1 \\
3 x_2 &= 1-x \\
2 x_4 + 3 x_5 &= x \\
2 x_2 + 2 x_3 + x_5 + x_6 &= 1 - y \\
x_4 &= y \\
3 x_1 + 3 x_3 &= 2
\end{aligned}
\right. \\
という連立方程式ができます。これを頑張って手で解くと、
\begin{align}
x_1 &= \frac{1}{2} \\
x_2 &= \frac{1}{3} - \frac{1}{3}x \\
x_3 &= \frac{1}{6} \\
x_4 &= y \\
x_5 &= \frac{1}{3}x - \frac{2}{3} y \\
x_6 &= \frac{1}{3}x - \frac{1}{3} y\\
\end{align}
となりました。この試料系で組成中の$x, y$を変化させるだけならExcelでどうにかできますが、これを一々色々な試料系で手で解くのはかなり面倒ですし、解く過程やExcelに式を書き起こす際にも間違いが生じ得ます。
SymPyを利用した単純な連立方程式の例
上に示したような連立方程式は SymPy を使って簡単に解くことができます。まずは単純な例として次の連立方程式を考えます。
\left\{ \,
\begin{aligned}
2 x_1 + x_2 &= 11 \\
- x_1 + 3 x_2 &= -23
\end{aligned}
\right. \\
\begin{aligned}
\\
\therefore \ x_1 = 8, \ x_2 = -5
\end{aligned}
これをSymPyを使って解くコードが以下です。
import sympy
# 変数の定義
x_1 = sympy.Symbol('x_1')
x_2 = sympy.Symbol('x_2')
# 式の用意
eq_1 = 2 * x_1 + x_2 - 11
eq_2 = - x_1 + 3 * x_2 + 23
# 式を変数について解く
sols = sympy.solve([eq_1, eq_2], [x_1, x_2])
print(sols)
{x_1: 8, x_2: -5}
このコードの流れは次のとおりです。
-
sympy.Symbol('~~')で変数を定義する。 - 定義した変数を使って連立方程式を作る。
このとき、式の右辺か左辺のどちらかが0になるように移項する必要があります。今の例では、右辺の項を左辺に移項しました。 -
sympy.solve(~~, ~~)で連立方程式を解かせる。
1つ目の引数には連立方程式をリストの形で入力します。2つ目の引数にはどの変数について解くかをリストの形で入力します。今の例では、x_1とx_2の両方について解きたいため、両方を指定しました。
sympy.solve(~~, ~~)で連立方程式を解かせる際、どの変数について解くかを指定できるため、指定していない変数を残したまま解くことができます。例えば、新たな変数$a$を加えた、
\left\{ \,
\begin{aligned}
2 x_1 + x_2 &= 11 + a \\
- x_1 + 3 x_2 &= -23
\end{aligned}
\right. \\
\begin{aligned}
\\
\therefore \ x_1 = 8 + \dfrac{3}{7} a, \ x_2 = -5 + \dfrac{1}{7} a
\end{aligned}
という連立方程式を、sympy.solve(~~, ~~)で$a$を指定せずにSymPyで解くと、
import sympy
# 変数の定義
x_1 = sympy.Symbol('x_1')
x_2 = sympy.Symbol('x_2')
a = sympy.Symbol('a')
eq_1 = 2 * x_1 + x_2 - 11 - a
eq_2 = - x_1 + 3 * x_2 + 23
solutions = sympy.solve([eq_1, eq_2], [x_1, x_2])
print(solutions)
{x_1: 3*a/7 + 8, x_2: a/7 - 5}
と、人間が手で解析的に解くように、解x_1とx_2の結果にaを変数として残したまま解くことができました。NumPyのような連立方程式を数値的にしか解けないライブラリでは解に変数を残すことができず、あらかじめその変数の数値を指定しなければなりません。一方、SymPyでは変数の数値をあらかじめ指定せずに、変数を変数として残したまま解くことができます。
SymPyを利用して試薬の係数を求める
結論から言うと、試薬の係数を求めるコードは以下のようになります。
import sympy
# 反応物（試薬）
reactant = [
{'La': 2, 'S': 3}, # La2S3
{'Bi': 2, 'O': 3}, # Bi2O3
{'Bi': 2, 'S': 3}, # Bi2S3
{'Pb': 1, 'F': 2}, # PbF2
{'Bi': 1, 'F': 3}, # BiF3
{'Bi': 1} # Bi
]
# 変数の定義
x = sympy.Symbol('x')
y = sympy.Symbol('y')
# 生成物
product = {'La': 1, 'O': 1 - x, 'F': x, 'Bi': 1 - y, 'Pb': y, 'S': 2} # LaO(1-x)F(x)Bi(1-y)Pb(y)S2
# --- 計算 ---
# 変数 x_1, x_2, ... を定義してリストへ入れる
symbols = [sympy.Symbol('x_' + str(i)) for i in range(1, 1 + len(reactant))]
# 式（equations）の用意
equations = {k: -v for k, v in product.items()} # -b を代入
for (symbol, react) in zip(symbols, reactant):
for k, v in react.items():
equations[k] += v * symbol # equation に 1 * x1 などを加える
solutions = sympy.solve(list(equations.values()), symbols)
print(solutions)
{x_1: 1/2, x_2: 1/3 - x/3, x_3: 1/6, x_4: y, x_5: x/3 - 2*y/3, x_6: x/3 - y/3}
解説
反応物（試薬）の定義
# 反応物（試薬）
reactant = [
{'La': 2, 'S': 3}, # La2S3
{'Bi': 2, 'O': 3}, # Bi2O3
{'Bi': 2, 'S': 3}, # Bi2S3
{'Pb': 1, 'F': 2}, # PbF2
{'Bi': 1, 'F': 3}, # BiF3
{'Bi': 1} # Bi
]
このコードでの試薬の定義は、試薬それぞれの化学式を辞書型で記し、それらをリストにまとめるという方法を採っています。
辞書型というのは、
a = {'La': 2, 'S': 3}
print(a['La'])
2
のように、定義した文字列に対応する数値などを返してくれる機能を持つ型です。
生成物（試薬）の定義
# 変数の定義
x = sympy.Symbol('x')
y = sympy.Symbol('y')
# 生成物
product = {'La': 1, 'O': 1 - x, 'F': x, 'Bi': 1 - y, 'Pb': y, 'S': 2} # LaO(1-x)F(x)Bi(1-y)Pb(y)S2
生成物の化学式に含まれる変数を先に定義し、その後に生成物を定義します。こちらも化学式を辞書型で記して定義します。NumPyの場合と異なる点は、この時点で$y = 0.08$などと具体的な変数の数値を指定する必要がないということです。
解きたい変数の定義
# 変数 x_1, x_2, ... を定義してリストへ入れる
symbols = [sympy.Symbol('x_' + str(i)) for i in range(1, 1 + len(reactant))]
[x_1, x_2, x_3, x_4, x_5, x_6]
リスト内包表記 を使って、定義した変数をリストに格納しています。len(reactant)でreactantの大きさ（長さ）が返されます。今の場合は6種類の試薬があるので、len(reactant)は6となり、range(1, 7)は1から6までの連番を返します。こうして作られた変数x_1, x_2, … , x_6は試薬の化学式の係数になり、これらが解きたい変数になります。
連立方程式の定義
equations = {k: -v for k, v in product.items()} # -b を代入
{'La': -1,
'O': x - 1,
'F': -x,
'Bi': y - 1,
'Pb': -y,
'S': -2}
辞書内包表記 を使って、連立方程式の右辺だけを辞書型で定義しています。単純に捉えると、productの各要素にマイナスを掛けただけです。元々、方程式は各元素について立てられるものであり、辞書の名前がどの元素についての方程式なのかを表しています。この後、equationsに連立方程式の左辺を加えていきます。
左辺の加算
for (symbol, react) in zip(symbols, reactant):
for k, v in react.items():
equations[k] += v * symbol # equation に 1 * x1 などを加える
{'Bi': 2*x_2 + 2*x_3 + x_5 + x_6 + y - 1,
'F': -x + 2*x_4 + 3*x_5,
'La': 2*x_1 - 1,
'O': x + 3*x_2 - 1,
'Pb': x_4 - y,
'S': 3*x_1 + 3*x_3 - 2}
各試薬reactとその係数symbolについて、対応する元素の方程式にその元素の個数とsymbolの積を加えていきます。for (symbol, react) in zip(symbols, reactant):でsymbolとreactにそれぞれx_1と{'La': 2, 'S': 3}、x_2と{'Bi': 2, 'O': 3}、... 、x_6と{'Bi': 1}というように試薬とその係数が順番に代入されてループしていきます。
そしてfor k, v in react.items():でkとvにそれぞれLaと2、Sと3というようにreact内の元素とそれに対応する数値が順番に代入されてループしていきます。
これらの2重ループについて、equations[k] += v * symbolで、kの方程式にvとsymbolの積を加えます。
これで連立方程式ができあがります。
連立方程式を解く
solutions = sympy.solve(list(equations.values()), symbols)
[2*x_1 - 1,
x + 3*x_2 - 1,
-x + 2*x_4 + 3*x_5,
2*x_2 + 2*x_3 + x_5 + x_6 + y - 1,
x_4 - y,
3*x_1 + 3*x_3 - 2]
[x_1, x_2, x_3, x_4, x_5, x_6]
sympy.solve(~~, ~~)の1つ目の引数では、辞書型の各要素の名前（'Bi'など）が邪魔なので除去します。
質量を求める
連立方程式を解いた後は、各試薬の分子量を使って計算するだけです。具体的には次のような流れで求められます。
- 解いた係数（
x_1, …）のそれぞれに試薬の分子量を掛ける。 - 1.の和で1.を割ることで、和を$\mathrm{1 \ g}$に規格化する。
- 2.に生成物の目標の質量を掛ける。
最終的なPythonコード
上のコードに質量の計算以外にも色々改良して、以下のコードになりました。
import sympy
def main():
# 変数の定義
x = sympy.Symbol('x')
y = sympy.Symbol('y')
# 後で代入する代入値
subs = {x : 0.5, y : 0.08}
# 反応物（試薬） 組成 分子量
reactant = [
[{'La': 2, 'S': 3}, 373.99], # La2S3
[{'Bi': 2, 'O': 3}, 465.958], # Bi2O3
[{'Bi': 2, 'S': 3}, 514.14], # Bi2S3
[{'Pb': 1, 'F': 2}, 245.2], # PbF2
[{'Bi': 1, 'F': 3}, 265.975], # BiF3
[{'Bi': 1}, 208.98] # Bi
]
# 生成物 組成
product = {'La': 1, 'O': 1 - x, 'F': x, 'Bi': 1 - y, 'Pb': y, 'S': 2} # LaO(1-x)F(x)Bi(1-y)Pb(y)S2
product_mass = 2 # 生成物の質量(g)
# 計算
solve(reactant, product, product_mass, subs=subs, decimal_place=4)
# decimal_place 小数第何位まで表示するか
return
def solve(reactant, product, product_mass=1, subs={}, decimal_place=10):
def dict_to_str(d):
return ''.join( [str(k) if v == 1 else str(k) + str(v).replace(' ', '') for k, v in d.items()] ) # 半角スペースを削除
reactant_dicts, reactant_weights = [r[0] for r in reactant], [r[1] for r in reactant] # それぞれ組成・分子量の部分のみ取り出す
symbols = [sympy.Symbol('x_' + str(i)) for i in range(1, 1 + len(reactant_dicts))] # 変数 x_1, x_2, ... を定義してリストへ入れる
# 反応式の表示
str_chem_lefts = [str(symbol) + ' ' + dict_to_str(reactant) for (symbol, reactant) in zip(symbols, reactant_dicts)]
str_chem = ' + '.join(str_chem_lefts) + ' -> ' + dict_to_str(product)
print('# Reaction Formula \n \t', str_chem)
# 式（equation）の用意
equations = {k: -v for k, v in product.items()} # -b を代入
for (symbol, reactant) in zip(symbols, reactant_dicts):
for k, v in reactant.items():
assert k in equations.keys(), f'生成物に含まれない元素（＝{k}）が試薬に含まれています。'
equations[k] += v * symbol # equation に 1 * x1 などを加える
solutions = sympy.solve(list(equations.values()), symbols)
assert len(solutions) > 0, f'解が存在しません'
print('# Solution \n \t', solutions) # 解の表示
solutions = {k: v.subs(subs) if isinstance(v, sympy.Basic) else v for k, v in solutions.items()} # 解に変数を代入
reactant_mass = [solutions[symbol] * weight for (symbol, weight) in zip(symbols, reactant_weights)]
sum_reactant_mass = sum(reactant_mass)
reactant_mass = [v * product_mass / sum_reactant_mass for v in reactant_mass] # 規格化してproduct_massを掛ける
# 質量の表示
print(f'# Mass \n \t product = {dict_to_str(product):<8} {subs}, product_mass = {product_mass} g \n \t \t reactant mass')
for i, (r_mass, r_dict) in enumerate(zip(reactant_mass, reactant_dicts)):
print(f'\t {i + 1:>3} \t {dict_to_str(r_dict):<8} {float(r_mass):>6.{decimal_place}f} g')
return reactant_mass
if __name__ == '__main__':
main()
# Reaction Formula
x_1 La2S3 + x_2 Bi2O3 + x_3 Bi2S3 + x_4 PbF2 + x_5 BiF3 + x_6 Bi -> LaO1-xFxBi1-yPbyS2
# Solution
{x_1: 1/2, x_2: 1/3 - x/3, x_3: 1/6, x_4: y, x_5: x/3 - 2*y/3, x_6: x/3 - y/3}
# Mass
product = LaO1-xFxBi1-yPbyS2 {x: 0.5, y: 0.08}, product_mass = 2 g
reactant mass
1 La2S3 0.8710 g
2 Bi2O3 0.3617 g
3 Bi2S3 0.3992 g
4 PbF2 0.0914 g
5 BiF3 0.1404 g
6 Bi 0.1363 g
このコードの特徴
・reactantとproductの部分を自由にいじることで、試薬と生成物がどんな組成でどんな元素の数であってもこのコード1つで対応することができ、事前に手計算せずに即座に秤量計算ができます。
・生成物の化学式に変数を導入し、その変数を残したまま試薬の係数を解くことができます。これはNumPyの場合では不可能でした。
・生成物の元素数より試薬の数が少なく（正方行列でない）、かつ解が存在するような場合でも正しく解くことができます。これもNumPyの場合では正方行列となる場合しか解けないため不可能でした。