【化学でPython】SciPy：複雑な積分や最適化問題を解く

はじめに

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

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

SciPy とは？

SciPy は、NumPy を基盤とした科学技術計算のための総合ライブラリです。

積分・最適化・統計・線形代数など、高度な数値計算機能を提供します。

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

主な特徴

• 豊富な数値計算機能: 積分、微分方程式、最適化、補間、統計など幅広い科学計算をカバー
• 高速な計算: 内部的にFortranやC言語で書かれた最適化されたコードを使用
• NumPyとの親和性: NumPy配列をそのまま扱えるため、データ処理がスムーズ
• 化学への応用: 反応速度論やスペクトル解析など、多様な用途に対応

インストール

pip で簡単にインストールできます。

terminal

pip install scipy

Google Colab での利用
Google Colab には SciPy がプリインストールされているため、通常はインストール不要です。

基本的な使い方

まずは、基本的な使い方を見てみましょう。

SciPy は機能ごとにサブモジュールに分かれていて、必要なものだけをインポートして使います。

import numpy as np
from scipy import optimize, integrate, stats

# 例1: 簡単な関数の最小値を求める
def f(x):
    return x**2 + 5*x + 6  # 2次関数

# 最小値を探索(初期値 x=0 から開始)
result = optimize.minimize(f, x0=0)
print(f"最小値: {result.fun:.3f}")
print(f"最小値を与えるx: {result.x[0]:.3f}")

# 例2: 定積分を計算する
# ∫[0,1] x^2 dx を計算
integral_result, error = integrate.quad(lambda x: x**2, 0, 1)
print(f"\n積分結果: {integral_result:.3f}")
print(f"誤差推定: {error:.2e}")

実行結果

最小値: -0.250
最小値を与えるx: -2.500

積分結果: 0.333
誤差推定: 3.70e-15

このように、SciPy を使えば複雑な数値計算を数行のコードで実行できます。

実践例①: 反応速度論のシミュレーション

実践的な例として、「一次反応の速度論シミュレーション」 をやってみます。

化学反応 A → B において、反応物 A の濃度変化を常微分方程式(ODE)で解き、時間変化をグラフ化します。

1. 反応速度式の定義

一次反応の速度式は以下のように表されます:

$$\frac{d[A]}{dt} = -k[A]$$

ここで、$k$ は反応速度定数です。

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# 反応速度式を定義(微分方程式)
def reaction_rate(C, t, k):
    """
    一次反応 A → B の速度式
    C: 濃度 [A]
    t: 時間
    k: 反応速度定数
    """
    dCdt = -k * C  # d[A]/dt = -k[A]
    return dCdt

# パラメータ設定
k = 0.5  # 反応速度定数 [1/s]
C0 = 1.0  # 初期濃度 [mol/L]
t = np.linspace(0, 10, 100)  # 0〜10秒を100点で分割

# 常微分方程式を数値的に解く
C = odeint(reaction_rate, C0, t, args=(k,))

# 結果をプロット
plt.figure(figsize=(8, 5))
plt.plot(t, C, 'b-', linewidth=2, label='[A] (Numerical Solution)')

# 理論解(解析解)と比較
C_analytical = C0 * np.exp(-k * t)
plt.plot(t, C_analytical, 'r--', linewidth=2, label='[A] (Analytical Solution)')

plt.xlabel('Time (s)', fontsize=12)
plt.ylabel('Concentration [A] (mol/L)', fontsize=12)
plt.title('Concentration Change of First-Order Reaction A → B', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

print(f"半減期: {np.log(2)/k:.2f} 秒")

実行結果

半減期: 1.39 秒

数値解(青線)と解析解(赤破線)がぴったり重なることが確認できますね。

これにより、scipy.integrate.odeint が正確に微分方程式を解いていることがわかります。

化学的な読み解き

• 反応速度定数 $k = 0.5 \, \text{s}^{-1}$ の場合、半減期は約1.39秒です
• 実験データから $k$ を求めたい場合は、後述のカーブフィッティングを使います

実践例②: スペクトルデータのピークフィッティング

次に、「UV-Visスペクトルのガウシアンピークフィッティング」 を行います。

実験で得られたスペクトルデータから、ピークの位置（波長）、高さ、幅を定量的に求めます。

1. 模擬スペクトルデータの作成

まず、ノイズを含む模擬的なスペクトルデータを作成します。

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# ガウシアン関数の定義
def gaussian(x, amplitude, center, width):
    """
    ガウシアン関数
    amplitude: ピークの高さ
    center: ピークの中心位置(波長)
    width: ピークの幅(標準偏差)
    """
    return amplitude * np.exp(-(x - center)**2 / (2 * width**2))

# 模擬スペクトルデータの生成
np.random.seed(42)
wavelength = np.linspace(400, 700, 300)  # 波長 400-700 nm

# 真のパラメータ(これを推定するのが目標)
true_amplitude = 1.5
true_center = 550  # nm
true_width = 30    # nm

# ノイズを含むスペクトルデータ
absorbance = gaussian(wavelength, true_amplitude, true_center, true_width)
absorbance += np.random.normal(0, 0.05, wavelength.shape)  # ノイズ追加

print("=== 模擬スペクトルデータを生成しました ===")
print(f"真のピーク位置: {true_center} nm")
print(f"真のピーク高さ: {true_amplitude:.2f}")
print(f"真のピーク幅: {true_width} nm")

実行結果

=== 模擬スペクトルデータを生成しました ===
真のピーク位置: 550 nm
真のピーク高さ: 1.50
真のピーク幅: 30 nm

2. カーブフィッティングの実行

scipy.optimize.curve_fit を使って、実験データからパラメータを推定します。

# 初期推定値(適当な値でOK)
initial_guess = [1.0, 500, 50]

# カーブフィッティング実行
popt, pcov = curve_fit(gaussian, wavelength, absorbance, p0=initial_guess)

# フィッティング結果
fitted_amplitude, fitted_center, fitted_width = popt

# 標準誤差の計算
perr = np.sqrt(np.diag(pcov))

print("\n=== フィッティング結果 ===")
print(f"推定ピーク位置: {fitted_center:.1f} ± {perr[1]:.1f} nm")
print(f"推定ピーク高さ: {fitted_amplitude:.3f} ± {perr[0]:.3f}")
print(f"推定ピーク幅: {fitted_width:.1f} ± {perr[2]:.1f} nm")

# フィッティング曲線の生成
fitted_curve = gaussian(wavelength, *popt)

# プロット
plt.figure(figsize=(10, 6))
plt.plot(wavelength, absorbance, 'o', markersize=3, alpha=0.5, label='Experimental Data (with noise)')
plt.plot(wavelength, fitted_curve, 'r-', linewidth=2, label='Fitted Curve')
plt.axvline(fitted_center, color='green', linestyle='--', alpha=0.7, label=f'Peak Position: {fitted_center:.1f} nm')
plt.xlabel('Wavelength (nm)', fontsize=12)
plt.ylabel('Absorbance', fontsize=12)
plt.title('Gaussian Fitting of UV-Vis Spectrum', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

実行結果

=== フィッティング結果 ===
推定ピーク位置: 549.8 ± 0.7 nm
推定ピーク高さ: 1.497 ± 0.012
推定ピーク幅: 30.1 ± 0.7 nm

ノイズを含むデータからでも、真の値にかなり近いパラメータを推定できました。

化学的な読み解き

• 定量分析: ピークの高さ(吸光度)から物質の濃度を算出できます
• 物質同定: ピーク位置(波長)は物質固有の値なので、未知物質の同定に使えます

実践例③: 吸着等温線のフィッティング

最後に、「Langmuir吸着等温線のフィッティング」 を行います。

吸着実験のデータから、最大吸着量と吸着平衡定数を求めます。

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# Langmuir吸着等温線の式
def langmuir(C, qmax, K):
    """
    Langmuir吸着等温線
    C: 平衡濃度 [mg/L]
    qmax: 最大吸着量 [mg/g]
    K: Langmuir定数 [L/mg]
    """
    return qmax * K * C / (1 + K * C)

# 実験データ(模擬)
C_exp = np.array([0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0])  # 平衡濃度 [mg/L]
q_exp = np.array([2.3, 4.1, 6.8, 12.5, 18.2, 24.1, 31.5, 35.8])  # 吸着量 [mg/g]

# カーブフィッティング
# 適切な初期推定値を与える (qmaxは最大吸着量付近、Kは正の値)
initial_guess_langmuir = [40.0, 0.1] # qmax=40 mg/g, K=0.1 L/mg として初期推定
popt, pcov = curve_fit(langmuir, C_exp, q_exp, p0=initial_guess_langmuir)
qmax_fit, K_fit = popt

# フィッティング曲線
C_fit = np.linspace(0, 100, 200)
q_fit = langmuir(C_fit, qmax_fit, K_fit)

# 結果表示
print("=== Langmuir吸着等温線フィッティング結果 ===")
print(f"最大吸着量 qmax: {qmax_fit:.2f} mg/g")
print(f"Langmuir定数 K: {K_fit:.4f} L/mg")

# 決定係数(R²)の計算
residuals = q_exp - langmuir(C_exp, qmax_fit, K_fit)
ss_res = np.sum(residuals**2)
ss_tot = np.sum((q_exp - np.mean(q_exp))**2)
r_squared = 1 - (ss_res / ss_tot)
print(f"Coefficient of Determination R²: {r_squared:.4f}")

# プロット
plt.figure(figsize=(10, 6))
plt.plot(C_exp, q_exp, 'o', markersize=8, label='Experimental Data')
plt.plot(C_fit, q_fit, 'r-', linewidth=2, label=f'Langmuir Fit (qmax={qmax_fit:.1f}, K={K_fit:.3f})')
plt.xlabel('Equilibrium Concentration C (mg/L)', fontsize=12)
plt.ylabel('Adsorption Amount q (mg/g)', fontsize=12)
plt.title('Langmuir Adsorption Isotherm', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

実行結果

=== Langmuir吸着等温線フィッティング結果 ===
最大吸着量 qmax: 38.86 mg/g
Langmuir定数 K: 0.0908 L/mg
決定係数 R²: 0.9967

決定係数 R² = 0.9967 と非常に高い値が得られ、Langmuirモデルがこのデータによく適合していることがわかります。

化学的な意義

• 吸着材の性能評価: 最大吸着量は、吸着材の性能を示す重要な指標です
• 吸着親和性: Langmuir定数が大きいほど、吸着材と吸着質の親和性が高いことを示します

まとめ

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

• Point 1: 積分・最適化・統計など、化学研究に必要な高度な数値計算を簡単に実行できます
• Point 2: 反応速度論のシミュレーション、スペクトル解析、吸着等温線フィッティングなど、実践的な化学計算に幅広く応用できます
• Point 3: NumPyとの組み合わせで、データ処理から高度な解析まで一貫して行えます

SciPy は「数値計算の万能ツール」。実験データの解析、理論モデルの検証、反応機構の解明など、あらゆる場面で活躍するはず。

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

参考リンク

• SciPy 公式ドキュメント
• GitHub Repository
• SciPy Tutorial
• SciPy Lecture Notes