using a time series of parameters to solve ODE in R

Question

I am trying to solve a simple ODE in R using deSolve: `dQ/dt = f(Q)*(P - E).`The whole thing is a time series of Q. The trick is that P and E are fixed time series of data themselves, so the diff eq is effectively in Q alone. It's relatively straightforward to solve this iteratively with a fixed time step, but I am trying to find a way to use an adaptive time step in R. Because P and E are discrete, successive timesteps may have the same value of P and E, which is fine. I have been playing around with deSolve, but have not been able to work this out. Ideally I would like to use standard 4th order Runge-Kutta. Any ideas? Do it in MATLAB? Edited for reproduceable example. I would like to be able to do this calc with a variable time step Runge-Kutta 4 method. I could program a fixed time step rk4 pretty easily, not so much adaptive. ``` working <- structure(list(datetime = structure(c(1185915600, 1185919200, 1185922800, 1185926400, 1185930000, 1185933600, 1185937200, 1185940800, 11859444…

Accepted Answer

Maybe not the most sophisticated approach, but the quick and dirty approach is to keep the time-dependent forcing variables as external (global) variables.

I use `pmax(1,ceiling(t))` to convert from the time index to the row index of the data frame (the `pmax` is necessary if you want to start from `t=0` because `ceiling(0)` is 0, and `x[0]` in R is generally an empty vector, which then breaks stuff). There are probably other ways to do the indexing.

```
library(deSolve)
gradfun <- function(t,y,parms) {
   pcp <- working$p[pmax(1,ceiling(t))]
   pet <- working$e[pmax(1,ceiling(t))]
   list(y^2*((pcp-pet)/y),NULL)
}

gradfun(0,working$qsim[1],1)   ## test
ode1 <- ode(y=c(qsim=working$qsim[1]),func=gradfun,
                time=seq(0,nrow(working),length.out=101),
                parms=NULL,method="rk4")
plot(ode1)
```

using a time series of parameters to solve ODE in R

Solution

Problem