R for solving differential equations: deSolve - Number of derivatives greater than initial conditions

Question

I am trying to solve a differential equation that describes the change in dissolved oxygen (DO) concentration in a pond between a specific set of days. The equation I am solving is: dc/dt=-kC+p where k is a first order DO decay constant (1/day) and p is a zero order DO production rate. I have an initial concentration and an equilibrium that is ultimately reached and have successfully solved this equation for daily average DO concentration using 'ode' in the deSolve package. I am now attempting to solve this differential for hourly concentration where decay occurs constantly but DO production only occurs between the hours of 6AM and 6PM. I have implemented the function as:

dc.dt<-function(time,y,parms) {
  with(as.list(c(y,parms)), {
    dC<-ifelse(times%%1>=0.25 & times%%1<=0.75, (-k)*C+(p), (-k)*C)
    list(dC)
  })
}
y<-c(C=6.2) 
parms<-c(p=(1.752),k=0.1)
times<-as.numeric((seq(from=as.POSIXct("2016-06-01 09:00",format="%Y-%m-%d 
%H:%M"),to=as.POSIXct("2016-06-30 23:00",format="%Y-%m-%d %H:%M"),
        by="hours")))
times<-(times-as.numeric(as.POSIXct("2016-06-01",format="%Y-%m-%d")))/86400

out<-ode(y,times,dc.dt,parms)

Times is in decimal day format beginning at 2016-06-01 9:00 and ending 2016-06-30 23:00

> head(times)
[1] 0.3750000 0.4166667 0.4583333 0.5000000 0.5416667 0.5833333

I am attempting to use modulus 1 of the times vector in my ifelsestatement where times %% 1 = 0.25 is 6AM and times %% 1 = 0.75 is 6PM.

Solving using ode seems to want to return the correct number of results but I am getting this error message:

> out<-ode(y,times,dc.dt3,parms)
Error in checkFunc(Func2, times, y, rho) : 
  The number of derivatives returned by func() (711) must equal the length of    
  the initial conditions vector (1)

The initial conditions vector that this is referring to if my vector for y which is simply an initial concentration. I'm not new to R but I am new to differential equations and this is my first time using R to solve these types of equations. Using a daily timestep, C starts at the initial value and reaches an equilibrium value. Using an hourly timestep I expect to see a diurnal pattern where DO starts at the initial value, exceeds the equilibrium value, and then goes below the equilibrium value again for each day. Any help that someone can offer on where I may have something wrong would be greatly appreciated!

J_F J_F · Accepted Answer · 2016-04-27T12:39:23

You have to change only a litte bit:

dc.dt<-function(time,y,parms) {
  with(as.list(c(y,parms)), {

    ifelse(times%%1>=0.25 & times%%1<=0.75, dC <- (-k)*C+(p), dC <-   (-k)*C)
    return(list(dC))

  })
}

Now it works. Because of the "ifelse" command the program thinks, that dC is of length 711. You have to allocate dC in the ifelse command.

Best wishes,

Johannes

R for solving differential equations: deSolve - Number of derivatives greater than initial conditions

1 Answers