'Differential Equation in R and Python with one parameter that changes in time (odeint and deSolve) [solved]
I'm trying to "translate" a script I have in R that can solve a DE of a function that has a parameter changing in time. I've compared with the results that I obtained from a commercial software and the results are comparable. While translating the code I took note of the different tolerance used by LSODA in deSolve and in Python (or at least I tried to using the Differences between Rs deSolve and Pythons odeint)
there's the two codes PYTHON
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
npoints=30
qqq=10
T0=0
TempEnd=500
timeStart=0
R=8.314
Ts = 273.15+T0
timeEnd = (TempEnd - T0)/(qqq/60) #estimated end time for the analysis depending on the rate
t = np.linspace(start=timeStart,stop=timeEnd, num=npoints)
Temperature = Ts + t * qqq/60
A=1995262
m=1
n=2
y0=5.00082970E-5
Ea=80000
#here we define the function that changes according to time
def SB(y1, t, qqq, Ts, Ea, m, n, A):
R=8.314
Temperature = Ts + t*qqq/60
coeff = np.exp(-Ea/(R*Temperature))
dydt = A * coeff * pow(y1,m) * pow((1-y1),n)
return dydt
yt = odeint(SB, y0, t,args=(qqq,Ts,Ea,m,n,A))
R
npoints=30
qqq=10
T0=0
TempEnd=500
time.start=0
R <- 8.314
Ts <- 273.15+ T0 #transformation in K
timeEnd=(TempEnd - T0)/(qqq/60) #estimated end time for the analysis depending on the rate
timeSeconds=seq(time.start,timeEnd, length.out=npoints) #vector with all the times for solving the equations time in seconds KARLINE: THIS HAS NO IMPACT ON THE SOLUTION...
a1_tm = Ts + timeSeconds*qqq/60
A=1995262
m=1
n=2
a0= 5.00082970E-5
Ea=80000
tm = timeSeconds
Temp = Ts+(timeSeconds*(qqq/60)) #temperatures calculated at each time
#function to solve
SBf = function(tm, state, parms)
{
with(as.list(c(state, parms)),
{
a1_tm = Ts+(tm*qqq/60)
dy1 = A* exp(-Ea/(R*a1_tm)) * y1^m*(1-y1)^n
return(list(dy1, # first return value = derivative (used to integrate)
Temp = a1_tm, # second value = calculated temperature
fi = dy1)) # the derivative again - now used for output
})
}
rootfun <- function(tm, y1, parms) {return(y1 - 1)}
state = c(y1 = a0) #
P <- list(Ts = Ts, qqq = qqq, A = A)
sol <- ode(y = state, times = tm, parms = P, func = SBf, rootfun=rootfun,atol = 1.49012e-8, rtol = 1.49012e-8)
sol[,2]
[1] 5.000830e-05 5.000834e-05 5.000861e-05 5.001002e-05 5.001632e-05 5.003160e-05 5.008512e-05 5.025359e-05 5.071842e-05 5.195406e-05 5.508693e-05
[12] 6.297770e-05 8.403674e-05 1.516641e-04 4.795905e-04 4.091840e-03 1.424349e-01 7.889817e-01 9.317398e-01 9.688010e-01 9.834006e-01 9.904015e-01
[23] 9.941441e-01 9.962867e-01 9.975735e-01 9.983746e-01 9.988878e-01 9.992246e-01 9.994502e-01 9.996042e-01
Sources
This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.
Source: Stack Overflow
| Solution | Source |
|---|
