R: зависящие от времени параметры в дифференциальных уравнениях для расчета химии (рабочие) - PullRequest
Новая
       75

R: зависящие от времени параметры в дифференциальных уравнениях для расчета химии (рабочие)

4 голосов
/ 20 сентября 2019

edit: в формуле произошла ошибка, и теперь все работает как положено.Спасибо за предложения, которые заставили меня пересмотреть код

Я пытаюсь изучить довольно известное уравнение химии Шестак – Берггрен, которое представляет собой мощный инструмент для описания кинетических данных методами подбора моделейвзяты из DOI: 10.1016 / j.tca.2011.03.030.Нет свободного / открытого пакета для моделирования, поэтому я пытаюсь добавить его в свой кинетический пакет takos https://cran.r -project.org / web / packages / takos / index.html .Поскольку это «вариация» логистической функции, я подумал использовать пакет deSolve.Единственная проблема в том, что у меня есть параметр, который зависит от времени.Я думал, что это может быть решено, используя приблизительно забаву, но я застрял, есть код [теперь работает] 

 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!

Я не получаю сигмоид, как ожидалось, но просто плоская линия.Любая помощь приветствуется!

...