CODEWITHSUNDEEP
pythonrscipyode
user20286118
user20286118

Reputation:

How to return ordinary variables from scipy odeint?

I have an R code, solving an ordinary differential equation using ode function from DeSolve package.

f <- function(t, y, parms) {
  with(as.list(c(y, parms)),{
    rate <- M * A/(A + K) * A
    
    dA <- -rate + D * (IN - A)
    dB <- rate - D * B
    
    Total <- A + B
    
    Conserved<- dA + dB - D * (IN - A) + D * B
    
    return (list(c(dA, dB),
                 Total = Total, rate = rate, Conserved = Conserved))
  })
}

I want to replicate this code in python, however, I cannot figure out how can I return also values of ordinary variables such as Total, rate, and Conserved. Right now, I'm recalculating everything after the integration, as in the code below, but I'm not sure how to recalculate the variable Conserved.

def f(C0, t):
    A, B = C0
    rate = M*A/(A+K)*B
    dA = - rate + D*(IN-A)
    dB = rate - D*B
    return [dA, dB]

times = np.arange(0, 0.05*365, 0.005)

C = odeint(f, C0, times)

DF = pd.DataFrame(C, columns=['A', 'B'])
DF['time'] = times
DF['Total'] = DF['A'] + DF['B']
DF['rate'] = M*DF['A']/(DF['A']+K)*DF['B']

Could you help me with converting this R code?

Upvotes: 1

Views: 316

Answers (1)

tpetzoldt
tpetzoldt

Reputation: 5813

Python's odeint does not support the output of internal variables directly, but it is possible to emulate the behavior.

Python

First create a function f2 with all desired variables. For sake of simplicity, let's also output time t and both the derivatives dA, dB and the states A, B:

import numpy as np
import scipy as sp
from scipy.integrate import odeint

def f2(C0, t):
    A, B = C0
    rate = M*A/(A+K) * A
    dA = - rate + D * (IN - A)
    dB =   rate - D * B
    Total = A + B
    Conserved = dA + dB - D * (IN - A) + D * B
    return [t, dA, dB, A, B, Total, rate, Conserved]

Now create a function f that returns only the derivatives:

def f(C0, t):
  x = f2(C0, t)
  return(x[1], x[2])

We need also some time steps, initial variables and parameters. As the OP did not provide values, let's set it to arbitrary values, e.g. 1.0:

times = np.arange(0, 0.05*365, 0.005)
C0 = 1, 1
D  = 1
K  = 1
IN = 1
M  = 1

Now we solve the model f with odeint:

C  = odeint(f, C0, times)

Finally, we run the full function with all internal variables in a loop with the results of odeint instead of C0.

ret = [None] * len(C)

for i in np.arange(1, len(C)):
  ret[i] = f2(C[i], times[i])

# last row of the result
ret[len(C)-1]

The procedure is essentially how R's deSolve::ode handles this internally. First it integrates the ode system and then it runs a second pass with the integrated states to get the internal variables.

R

Here for comparison a full R code:

library(deSolve)

f <- function(t, y, parms) {
  with(as.list(c(y, parms)),{
    rate <- M * A/(A + K) * A
    dA <- -rate + D * (IN - A)
    dB <- rate - D * B
    Total <- A + B
    Conserved <- dA + dB - D * (IN - A) + D * B
    return (list(c(dA, dB),
                 Total = Total, rate = rate, Conserved = Conserved))
  })
}

times <- seq(0, 0.05*365, 0.005)
parms <- c(D=1, K=1, IN=1, M=1)
y0 <- c(A=1, B=1)

out <- ode(y0, times = times, f, parms)
tail(out, 1)

Upvotes: 1

Related Questions