问题描述
我正在尝试建立一种将ODE系统适合我的模拟化学反应数据的方法。我正在尝试找到一种构建工具来从实验数据中查找某些参数(速率常数,例如k1_f,k1_r,k2_f等)的方法。
据我了解,这些速率常数/参数是针对所追踪的每种物质/化学物质的化学反应速率/ ODE定义的。因此,如果我们将原始数据适合ODE,则可以提取这些参数。而且,我必须能够在没有每种物质/化学物质数据的情况下进行此操作(因为有时无法通过实验确定这些数据)。
我在下面尽了最大的努力,从Kitchin组网站上构建了一个简单的示例,但无法弄清楚如何使用curve_fit将ODE系统拟合到数据中。
我正在尝试解决此问题,部分是为了找到可行的方法,另一方面是为了更好地了解合适的ODE的工作原理。
下面是我正在处理的代码:
#Defining initial Data
tspan = [0.,0.6,1.2,1.8,2.4,3.,3.6,4.2,4.8]
Prod_data = [0.,0.60420736,0.73895125,0.7591921,0.76204478,0.76244296,0.76249842,0.76250616,0.76250722]
#Defining Function to get fit
def fitfunc1(t,k):
k1f = k[0]
k1r = k[1]
k2f = k[2]
k2r = k[3]
def oderxn(C,t):
#Defining differential form for each species
Ce = C[0]
Cs = C[1]
Ces = C[2]
Cp = C[3]
#Defining rate laws for forward and reverse reaction
frxn1 = k1f * Ce * Cs
rrxn1 = k1r * Ces
frxn2 = k2f * Ces
rrxn2 = k2r * Ce * Cp
dEdt = -frxn1+rrxn1+frxn2-rrxn2
dSdt = -frxn1+rrxn1
dESdt = frxn1-rrxn1-frxn2+rrxn2
dPdt = frxn2-rrxn2
return [dEdt,dSdt,dESdt,dPdt]
C0 = [0.4,1,0]
Csol = odeint(oderxn,C0,t)
return Csol[:,0]
C_int = (0.4,Prod_data)
print(C_int)
k_fit_a = curve_fit(fitfunc1,tspan,C_int)
print (k_fit_a)
tfit = np.linspace(0,5);
fit = fitfunc1(tfit,k_fit_a)
import matplotlib.pyplot as plt
plt.plot(tspan,A_data,'ro',label='data')
plt.plot(tfit,fit,'b-',label='fit')
plt.legend(loc='best')
以下是使用Symfit的示例(我必须开始工作),以帮助解决问题:
from symfit import variables,parameters,ODEModel,D,Fit
from symfit.core.support import key2str
import numpy as np
import matplotlib.pyplot as plt
#Defining initial Data
tspan = [0.,0.76250722]
# Define our ODE model
t,E,S,ES,P = variables('t,P')
k1f,k1r,k2f,k2r = parameters('k1f,k2r')
S_0 = 1 #Inital [S]0 from above
E_0 = 0.4 #Initial [E]0 from above
frxn1 = k1f * E * S
rrxn1 = k1r * ES
frxn2 = k2f * ES
rrxn2 = k2r * E * P
model_dict = {
D(E,t): -frxn1+rrxn1+frxn2-rrxn2,D(S,t): -frxn1+rrxn1,D(ES,t): frxn1-rrxn1-frxn2+rrxn2,D(P,t): frxn2-rrxn2,}
model = ODEModel(
model_dict,initial={t: 0.0,E: E_0,S: S_0,ES: 0.0,P: 0.0}
)
print(model)
#Fitting the Data
#k1_f.min,k1_f.max = 1,20
#k1_r.min,k1_r.max = 1,20
#k2_f.min,k2_f.max = 0,20
#k2_r.min,k2_r.max = 0,20
fit = Fit(model,t=tspan,E= None,S=None,ES=None,P=Prod_data)
fit_result = fit.execute()
print(fit_result)
#Graphing Results
taxis = np.linspace(0,5,1000)
model_fit = model(t=taxis,**fit_result.params)._asdict()
plt.scatter(tspan,Prod_data,label='[P] Data',color='blue')
colors = ('red','blue','green','yellow')
for var in model:
plt.plot(taxis,model_fit[var],label='[{}]'.format(var.name),linewidth=2.0)
plt.legend()
plt.show()
解决方法
暂无找到可以解决该程序问题的有效方法,小编努力寻找整理中!
如果你已经找到好的解决方法,欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@)