🧪

【化学でPython】SymPy：化学反応速度式などの解析解をさくっと導出する

に公開
Python
sympy
chemistry
化学でpython
tech

はじめに

この「化学でPython」シリーズでは、化学の分野で有用な Python ライブラリを紹介しています。

今回紹介するのは、SymPy です。

SymPy とは？

SymPy は、Python で 記号計算（Symbolic Computation） を行うためのライブラリです。

物理化学のレポートや論文作成で、「複雑な微分方程式を解きたい」「反応速度式の一般解を導出したい」といった場面に遭遇したことはありませんか？

手計算だと計算ミスが怖い場面でも、SymPy を「高度な計算用紙」として使うことで、正確かつ高速に数式処理を行うことができます。

具体的には、NumPy などが数値を扱うのに対し、SymPy は x や y といった文字式を文字のまま扱って計算します。

  • 公式サイト: Link
  • GitHub: Link (⭐️ 14k+)

インストール

pip で簡単にインストールできます。Google Colab にはデフォルトでインストールされています。

terminal
pip install sympy

基本的な使い方

まずは、SymPy の基本的な操作感を見てみましょう。

変数を定義し、微分や因数分解を行ってみます。

from sympy import symbols, expand, factor, diff, integrate, sin, exp, init_printing

# 数式をきれいに表示するための設定（Jupyter/Colab用）
init_printing()

# 1. 記号（シンボル）の定義
x, y = symbols('x y')

# 2. 式の展開と因数分解
expr = (x + y)**3
print(f"元の式: {expr}")

expanded = expand(expr)
print(f"展開: {expanded}")

factored = factor(expanded)
print(f"因数分解: {factored}")

# 3. 微分
# f = x^2 * sin(x) を x で微分
f = x**2 * sin(x)
df = diff(f, x)
print(f"微分結果: {df}")

# 4. 積分
# f を x で不定積分
int_f = integrate(f, x)
print(f"積分結果: {int_f}")
実行結果
元の式: (x + y)**3
展開: x**3 + 3*x**2*y + 3*x*y**2 + y**3
因数分解: (x + y)**3
微分結果: x**2*cos(x) + 2*x*sin(x)
積分結果: -x**2*cos(x) + 2*x*sin(x) + 2*cos(x)

実践例: 連鎖反応の速度論的解析

ここでは具体的な化学のタスクとして、「連鎖反応 A \xrightarrow{k_1} B \xrightarrow{k_2} C の濃度変化」 を解析してみましょう。

この反応系は、最初 A だけが存在し、それが B に変わり、さらに B が分解して C になるというモデルです。

中間体 B の濃度がどのように変化するか、また B の濃度が最大になるのはいつか、を手計算で求めるのは少々骨が折れますが、SymPy なら、あっという間です。

1. 微分方程式を立てて解く

まず、反応速度式（微分方程式）を定義し、一般解を求めます。

  • \frac{d[A]}{dt} = -k_1 [A]
  • \frac{d[B]}{dt} = k_1 [A] - k_2 [B]
  • \frac{d[C]}{dt} = k_2 [B]
import sympy as sp
import matplotlib.pyplot as plt
import numpy as np

# 1. シンボルの定義
t = sp.symbols('t', positive=True) # 時間 t > 0
k1, k2 = sp.symbols('k1 k2', positive=True) # 速度定数 k > 0
A0 = sp.symbols('A0', positive=True) # 初期濃度

# 時間の関数として濃度を定義
A = sp.Function('A')(t)
B = sp.Function('B')(t)
C = sp.Function('C')(t)

# 2. 微分方程式の定義
# d[A]/dt = -k1 * [A]
eq_A = sp.Eq(A.diff(t), -k1 * A)

# d[B]/dt = k1 * [A] - k2 * [B]
eq_B = sp.Eq(B.diff(t), k1 * A - k2 * B)

# d[C]/dt = k2 * [B] (今回は物質収支から [C] = A0 - [A] - [B] としてもよいが、練習のため方程式として解く)
eq_C = sp.Eq(C.diff(t), k2 * B)

# 3. 微分方程式を解く (dsolve)
# まず A を解く
sol_A = sp.dsolve(eq_A, A, ics={A.subs(t, 0): A0}) # icsは初期条件 (t=0で A=A0)
print(f"[A](t) = {sol_A.rhs}")

# 求まった A の式を B の方程式に代入して解く
eq_B_subs = eq_B.subs(A, sol_A.rhs)
sol_B = sp.dsolve(eq_B_subs, B, ics={B.subs(t, 0): 0}) # 初期条件: t=0で B=0
print(f"[B](t) = {sol_B.rhs}")

# 求まった B の式を C の方程式に代入して解く
eq_C_subs = eq_C.subs(B, sol_B.rhs)
sol_C = sp.dsolve(eq_C_subs, C, ics={C.subs(t, 0): 0}) # 初期条件: t=0で C=0
print(f"[C](t) = {sol_C.rhs}")
実行結果
[A](t) = A0*exp(-k1*t)
[B](t) = A0*k1*exp(-k2*t)/(k1 - k2) - A0*k1*exp(-k1*t)/(k1 - k2)
[C](t) = A0*k1*k2*(-exp(-k2*t)/k2 + exp(-k1*t)/k1)/(k1 - k2) + A0

これで各成分の濃度変化の厳密解が得られました。[C](t) の式は少し複雑に見えますが、SymPy が自動で整理してくれます。

2. 中間体Bが最大になる時間を求める

次に、中間体 B の濃度が最大になる時間 t_{max} を求めます。

これは \frac{d[B]}{dt} = 0 となる t を求めれば良いはず。

# [B] の式を取り出す
B_expr = sol_B.rhs

# t で微分
dB_dt = sp.diff(B_expr, t)

# 微分係数が 0 になる t を解く
t_max_sols = sp.solve(dB_dt, t)
print(f"t_max の候補: {t_max_sols}")

# 実数解としてきれいな形を取り出す（通常は1つ解が得られます）
t_max = t_max_sols[0]
print(f"t_max = {t_max}")
実行結果
t_max の候補: [-log(k2/k1)/(k1 - k2)]
t_max = -log(k2/k1)/(k1 - k2)

この結果 t_{max} = - \frac{\ln(k_2/k_1)}{k_1 - k_2} = \frac{\ln(k_1/k_2)}{k_1 - k_2} は、物理化学の教科書でよく見る式ですよね。これを一瞬で導出できるのは便利でしょ。

3. グラフで可視化する

最後に、得られた式に具体的な数値を代入してプロットしてみましょう。

SymPy の式（B_expr など）は sp.lambdify を使うことで、NumPy 対応の高速な関数に変換できます。

# 数値を代入するための関数化 (lambdify)
# k1=1.0, k2=0.5, A0=1.0 とする
param_dict = {k1: 1.0, k2: 0.5, A0: 1.0}

func_A = sp.lambdify(t, sol_A.rhs.subs(param_dict), 'numpy')
func_B = sp.lambdify(t, sol_B.rhs.subs(param_dict), 'numpy')
func_C = sp.lambdify(t, sol_C.rhs.subs(param_dict), 'numpy')

# プロット用のデータ作成
t_vals = np.linspace(0, 10, 100)
A_vals = func_A(t_vals)
B_vals = func_B(t_vals)
C_vals = func_C(t_vals)

# プロット
plt.figure(figsize=(8, 5))
plt.plot(t_vals, A_vals, label='[A] (Reactant)')
plt.plot(t_vals, B_vals, label='[B] (Intermediate)')
plt.plot(t_vals, C_vals, label='[C] (Product)')

# t_max の位置を計算して表示
t_max_val = float(t_max.subs(param_dict))
plt.axvline(x=t_max_val, color='gray', linestyle='--', label=f't_max = {t_max_val:.2f}')

plt.xlabel('Time')
plt.ylabel('Concentration')
plt.title(f'Consecutive Reaction: k1={param_dict[k1]}, k2={param_dict[k2]}')
plt.legend()
plt.grid(True)
plt.show()


実行結果

グラフを見ると、A が減少し、B が一度増えてから減少し、C が増え続ける様子が分かります。また、B のピーク位置も計算通りであることが確認できますね。

まとめ

今回は SymPy を紹介しました。

  • Point 1: 面倒な微積分や方程式を Python に任せて正確に解ける。
  • Point 2: 化学反応速度式（微分方程式）の解析解を導くのに最適。
  • Point 3: 導出した式をそのままグラフ化や数値計算に利用できる。

実験データのフィッティングなどで「理論式が複雑すぎて導出できない…」と悩んだときは、ぜひ一度 SymPy に計算させてみてください。意外とあっさり解が見つかるかもしれませんよ。

ぜひ試してみてください。

参考リンク

Discussion