在Python中将一组数据拟合到二阶动力学速率方案

Question

我正在处理一组数据

时间 = [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速率常数。因此，我希望生成一个速率常数，显示 R 和 P 转化为 AB 的速率，仅包含 R 和 P 的初始浓度以及整个时间范围内 AB 的浓度。

我首先尝试使用以下代码将此动力学方案视为曲线拟合模型：

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()

得到这个结果：

我认为这没有产生我所希望的结果，因为产生的“拟合”与数据不太匹配。此外，我认为这不是最好的方法（甚至不是正确的方法），因为我正在努力了解时间进入上面的模型函数（kinet）的位置，并且我觉得它应该包含在内我们正在将某些东西映射为时间的函数。所以无论如何，然后我尝试在网上寻找可以帮助我的资源。然后我尝试用这行代码将速率方程重写为微分方程：

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

但我正在努力弄清楚如何从那里开始。是否有人有将数据拟合到仅可用产品浓度值的二阶动力学数据方案的经验？或者有任何有用的提示可以帮助我弄清楚如何做到这一点，或者即使这是一个合理的要求？谢谢！

Answer 1

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize

解决动态

首先让我们看看如何用

solve_ivp

求解动态系统。 我们为每种物质写下比率：

def kinetic(t, y, k):
    return np.array([
        - k * y[0] * y[1],  # [P]
        - k * y[0] * y[1],  # [R]
        + k * y[0] * y[1]   # [AB]
    ])

我们解决它的一些常数：

k = 1.345e-1
tlin = np.linspace(0, 10, 200)

sol = integrate.solve_ivp(kinetic, [tlin.min(), tlin.max()], y0=[1., 1., 0.], t_eval=tlin, args=(k,))

#  message: The solver successfully reached the end of the integration interval.
#  success: True
#   status: 0
#        t: [ 0.000e+00  5.025e-02 ...  9.950e+00  1.000e+01]
#        y: [[ 1.000e+00  9.933e-01 ...  4.280e-01  4.267e-01]
#            [ 1.000e+00  9.933e-01 ...  4.280e-01  4.267e-01]
#            [ 0.000e+00  6.713e-03 ...  5.720e-01  5.733e-01]]
#      sol: None
# t_events: None
# y_events: None
#     nfev: 38
#     njev: 0
#      nlu: 0

解决方案如下：

解决适配问题

现在，如果我们不想费心数学来推导我们的底层模型，我们可以在 IVP 解决方案之上构建：

def model(t, k):
    sol = integrate.solve_ivp(kinetic, [t.min(), t.max()], y0=[1., 1., 0.], t_eval=t, args=(k,))
    return sol.y[-1,:]

我们记得我们的假设是单一反应物浓度和零初始产物常数。

我们将

(0, 0)

点添加到您的数据集中：

t = np.array([0, 1, 1.25, 2.5, 3.75, 5, 6.25, 7.5, 8.75, 10])
C = np.array([0, 0.041355887, 0.228856274, 0.283712222, 0.401528071, 0.450842768, 0.514348728, 0.550876642, 0.61845291, 0.663312161])

然后简单地拟合函数：

popt, pcov = optimize.curve_fit(model, t, C)
# (array([0.17119482]), array([[0.00011828]]))

渲染效果如下：