I include the 'time' paramter in my differential equation solver and it messes it up somehow

Question

library(deSolve)
require(deSolve)

delta_conc <- function(time, current_state, params) {

  with(as.list(c(current_state, params)),{

    dX <- Y
    dY <- X - X^3 - 0.25*Y + A * sin(t)

    return(list(c(dX, dY)))
  })
}

params <- c(
  A <- 0.2645
)

initial_state <- c(
  X <- 0.9,
  Y <- 0.4
)

times <- 1:10

model <- ode(initial_state, times, delta_conc, params)

summary(model)

matplot(model, type="l",lty=1, main="Enzyme model", xlab="Time")

I get this error message when I try to run it:

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

When I exclude the 'sin(t)' part it works, so the problem is with that part, but I'm very much a beginner so I have no idea how to approach this problem

Answer 1

You should consistently use einer t or time for the actual time step. In your case t is not defined as variable, so tis interpreted as transpose-function.

The following should work:

require(deSolve)

delta_conc <- function(time, current_state, params) {
  with(as.list(c(current_state, params)),{
    dX <- Y
    dY <- X - X^3 - 0.25*Y + A * sin(time)
    return(list(c(dX, dY)))
  })
}

params <- c(
  A = 0.2645
)

initial_state <- c(
  X = 0.9,
  Y = 0.4
)

times <- 1:10
model <- ode(initial_state, times, delta_conc, params)
summary(model)

matplot.0D(model, type="l",lty=1, main="Enzyme model", xlab="Time")

In addition, the code had also some other issues:

• use either require or library and not both
• use = within c(). It is parameter matching and not assignment

Two additional suggestions:

• you can use the deSolve-built in plot function matplot.0D
• I would recommend to use times <- seq(0, 10, length.out = 100) instead of 1:10. This way the plot will get smooth. Starting time with 1 (or another value) may be ok, but is often more convenient to start time with zero.