如何用Python将实验数据拟合至二级动力学速率模型？

实验数据 #

• 时间序列：Time = [1, 1.25, 2.5, 3.75, 5, 6.25, 7.5, 8.75, 10]
• 产物浓度序列：AB = [0.041355887, 0.228856274, 0.283712222, 0.401528071, 0.450842768, 0.514348728, 0.550876642, 0.61845291, 0.663312161]

拟合目标 #

将数据拟合至二级动力学速率模型，已知反应物R和P的初始浓度均为1，提取速率常数k。

错误尝试 #

最初的拟合代码未考虑时间变量，模型完全错误，拟合结果与数据匹配度极差：

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

with open('binding.csv', 'r') as d:
    reader = csv.reader(d)
    no_comments = (line for line in d if not line.lstrip().startswith('#'))
    next(no_comments, None)
    data = list(reader)
    data_array = np.array(data, dtype='float64')

Time = data_array[:,0]
AB = data_array[:,1]

def kinet(k, P, R):
    return k * (P * R)

popt, pcov = curve_fit(kinet, Time, AB)
print(f"Fitted k = {(popt[1])}")

plt.plot(Time, kinet(Time, *popt), "k--", label='Curve Fit')
plt.plot(Time, AB, 'o', label='Raw Data')
plt.xlabel("Time")
plt.ylabel("AB")
plt.legend()

之后尝试改写为微分方程，但不知如何继续：

def rate(y, t):
    c = {"AB" : AB}
    dc = dict()
    dc["AB"] = k*c["R"]*c["P"]
    dy = dc["AB"]
    return dy

模型推导 #

对于二级反应 ( R + P \rightarrow AB )，初始浓度 ( [R]_0 = [P]_0 = 1 )，设产物浓度 ( [AB] = x )，则时刻t的反应物浓度为 ( [R] = [P] = 1 - x )。

反应速率的积分形式为：
[
\frac{x}{1 - x} = k t
]
变形后得到关于x（即AB浓度）的表达式：
[
x(t) = \frac{k t}{1 + k t}
]

拟合代码实现 #

利用scipy.optimize.curve_fit直接拟合上述解析解，无需解微分方程：

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

# 读取数据
with open('binding.csv', 'r') as d:
    no_comments = (line for line in d if not line.lstrip().startswith('#'))
    next(no_comments, None)  # 跳过表头
    data = list(csv.reader(no_comments))
    data_array = np.array(data, dtype='float64')

Time = data_array[:, 0]
AB = data_array[:, 1]

# 定义二级动力学模型的解析表达式
def second_order_model(t, k):
    return (k * t) / (1 + k * t)

# 拟合曲线
popt, pcov = curve_fit(second_order_model, Time, AB, p0=[0.1])  # p0为k的初始猜测值
k_fitted = popt[0]
print(f"拟合得到的速率常数k = {k_fitted:.6f}")

# 绘制结果
plt.figure(figsize=(8, 6))
plt.plot(Time, AB, 'o', label='原始数据')
plt.plot(Time, second_order_model(Time, k_fitted), 'k--', label=f'拟合曲线 (k={k_fitted:.6f})')
plt.xlabel('时间')
plt.ylabel('产物AB浓度')
plt.legend()
plt.grid(True)
plt.show()

关键说明 #

• 直接使用积分后的解析解避免了数值求解微分方程的复杂度，拟合效率更高
• p0参数可以提供k的初始猜测值，帮助拟合算法更快收敛
• 可通过计算决定系数 ( R^2 ) 验证拟合结果可靠性：

from sklearn.metrics import r2_score
y_pred = second_order_model(Time, k_fitted)
r2 = r2_score(AB, y_pred)
print(f"决定系数R² = {r2:.6f}")

内容的提问来源于stack exchange，提问作者HAILEY RUDE