5

编辑：公式中有错误，现在一切正常。感谢让我查看代码的建议

我正在尝试研究一个非常著名的化学方程式 Šesták–Berggren，它代表了一个强大的工具，用于通过取自 DOI 的模型拟合方法描述动力学数据：10.1016/j.tca.2011.03.030。没有可用于模拟它的免费/开放包，所以我试图将它添加到我自己的动力学包 takos https://cran.r-project.org/web/packages/takos/index.html中。由于它是逻辑函数的“变体”，我正在考虑使用 deSolve 包。唯一的问题是我有一个取决于时间的参数。我认为它可以使用 approx fun 来解决，但我被困在那里的代码 [现在工作]

 rm(list=ls())
library(deSolve)
A=10^6.3 #parameter needed by the function which is "fixed"
Ea=80000 #an example of activation energy
R=8.314  #gas constant
npoints=100 #just 1000 points to start
qqq=5  #ratio
T0=0   #starting temperature
T.end=1200 #end temperature
Ts=273.15+T0 #transformation in K
time.e=(T.end-T0)/(qqq/60) #estimated time for the analysis
time.s=seq(0.1,time.e, length.out=npoints) #vector with all the times
Temp=Ts+(time.s*(qqq/60)) #temperatures calculated at each time
tm=time.s
m=1 #parameter in the Sestak-Berggen Equation da/dt=A*exp(-Ea/RT)*a^m*(1-a)^n
n=2#secon parameter

eqRG = function(tm, state, parms)  #mod of a working example of https://www.researchgate.net/post/How_can_I_solve_a_system_of_ODEs_with_time_dependent_parameters_in_R_or_in_Python
{
with(as.list(c(tm, state, parms)),
{
   a1 = parms[["a1"]]

   dy1 =   A*  exp((-Ea)/(R*a1(tm)))  *(y1)*(1-y1)^2 #this should do the trick
   #print("exp") #to check
   #print (exp((-Ea)/(R*a1(tm))))
   #print("exp2")
   #print((y1)*(1-y1)^2)
   #print("exp3")
   #print(dy1)
   return(list(c(dy1)))
   })
}
#tm = seq(0, 10, len = 100)
state = c(y1 = 0.01) #starting value
a1 = approxfun(tm,Temp) #function that changes in time
P = list(a1 = a1)
sol = ode(y = state, times = tm, parms = P, func = eqRG) #it works but gives a flat result
plot(sol) #correct!

我没有像预期的那样获得 sigmoid，而只是一条平线。任何帮助表示赞赏！

4

0 回答 0