'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