'The number of derivatives returned by func() (2) must Equal the length of the initial conditions vector (4)' in R

Question

I’m using ode, package deSolve to solve an R ED system, but I get the error:

Error in checkFunc(Func2, times, y, rho) : The number of derivatives returned by func() (2) must equal the length of the initial conditions vector (4)

My return vector has size 4 as the initial condition, so I don’t know what the problem is.

monod = function(t, c, parametros){
  
  # Parametros
  um = parametros$um 
  Ks = parametros$ks 
  Yx = parametros$Yx 
  Yp = parametros$Yp 
  up = parametros$up 
  Kp = parametros$Kp 
  
  # Concentração das espécies
  X = c["X"]
  P = c["P"]
  S = c["S"]
  V = c["V"]
  
  # derivadas d/dt 
  dX = um*S*X/(S+Ks)-(0.7/V)*X                 #dS/dt
  dP = up*S*X/(S+Kp)-(0.7/V)*P                 #dE/dt
  dS = -(1/Yx)*dX-(1/Yp)*dP-(0.7/V)*(S-30)     #dC/dt
  dV = 0.7                                     #dV/dt
  
  return( list( c(dX, dP, dS, dV) ) )
}

#Teste de saída
c_inicial = c(X = 1.7, P = 0, S = 40, V = 1)

# tempo
t     = seq(0, 96, by = 1)

# Parametros
parametros = list(um = 0.1, Ks = 0.18, Yx = 0.25, Yp = 0.68, up = 0.05, Kp = 0.0002)

# solver (método pardrão lsoda)
out = ode(y = c_inicial, times = t, func = monod, parms = parametros)
head(out)

If anyone can help me, I’d appreciate it!

Answer 1

I think the problem is how you define the initial values and parameters within the model function. Looking at the package documentation deSolve and following the guidelines on how to build the model with the function with(), will work:

library(deSolve)

monod = function(t, c, parametros){
    with(as.list( c(c, parametros)), {

        # derivadas d/dt 
        dX = um*S*X/(S+Ks)-(0.7/V)*X                 #dS/dt
        dP = up*S*X/(S+Kp)-(0.7/V)*P                 #dE/dt
        dS = -(1/Yx)*dX-(1/Yp)*dP-(0.7/V)*(S-30)     #dC/dt
        dV = 0.7                                     #dV/dt
        
        return( list( c(dX, dP, dS, dV) ) )
    })
}

#Teste de saída
c_inicial = c(X = 1.7, P = 0, S = 40, V = 1)

# tempo
t     = seq(0, 96, by = 1)

# Parametros
parametros = list(um = 0.1, Ks = 0.18, Yx = 0.25, Yp = 0.68, up = 0.05, Kp = 0.0002)

# solver (método pardrão lsoda)
out = ode(y = c_inicial, times = t, func = monod, parms = parametros)
head(out)