edit: there was an error in the formula, and now everything works as expected. Thanks for the suggestions that made me review the code
I'm trying to study a quite famous chemistry equation Šesták–Berggren that represent a powerful tool for the description of kinetic data by the model-fitting methods as taken from DOI: 10.1016/j.tca.2011.03.030. There are no free/open package available for simulating it and so I'm trying to add it to my own kinetic package takos https://cran.r-project.org/web/packages/takos/index.html. Since it is a "variation" of a logistic function I was thinking to use the deSolve package. The only problem is that I've a parameter that depends in time. I thought it could be solved using approx fun but I'm stuck there's the code [now working]
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!
I do not obtain a sigmoid as expected but just a flat line. Any help appreciated !
