R help - Jun 2009 - External signal in ODE written in C (using deSolve and approx1?)

Dear list

The deSolve package allows you to specify the model code in C or Fortran.
Thanks to the excellent vignette this works fine. However I have not yet
managed to use forcing functions in C code.

In pure R code this works very well with approxfun() specified outside the
model:

###############################################
#Model
lvml <- function(t, x, parms) {
       with(as.list(c(parms, x)), {
     import <- sigimp(t)
         dS <- import - b*S*P + g*K     #substrate
         dP <- c*S*P  - d*K*P           #producer
         dK <- e*P*K  - f*K             #consumer
         res <- c(dS, dP, dK)
         list(res)
       })
     }

## Parameters
     parms  <- c(b = 0.0, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.0)
## vector of timesteps
     times  <- seq(0, 100, length = 101)
## external signal with rectangle impulse
     signal <- as.data.frame(list(times = times,
                                 import = rep(0,length(times))))
signal$import[signal$times >= 10 & signal$times <= 21] <- 0.2
sigimp <- approxfun(signal$times, signal$import, rule = 2)
 ## Start values for steady state
y <- xstart <- c(S = 1, P = 1, K = 1)
## Solving
out2 <- as.data.frame(lsoda(xstart, times, lvml, parms))
###############################################

I would like to do the same thing in C, and my guess is that the approx1
function has to be used in some way. So far did not figure out how. If
anyone has managed to do so, or know how to approach this problem, please
let me know.

Best wishes
Glenn

	[[alternative HTML version deleted]]

Glenn Woodart wrote:
> Dear list
> 
> The deSolve package allows you to specify the model code in C or Fortran.
> Thanks to the excellent vignette this works fine. However I have not yet
> managed to use forcing functions in C code.
> 
> In pure R code this works very well with approxfun() specified outside the
> model:
[... example deleted, see original post ...]
> I would like to do the same thing in C, and my guess is that the approx1
> function has to be used in some way. So far did not figure out how. If
> anyone has managed to do so, or know how to approach this problem, please
> let me know.
> 
> Best wishes
> Glenn
Hi Glen,

this problem is on our radar for a while, however, we have not found the 
time yet to implement this in a general way. There are, in principal, 
several different possibilities:

1) use a solver with a fixed internal step size, e.g. rk4 and supply the 
data in the appropriate discretization (note also the intermediate 
steps). The upcoming deSolve 1.3 on R-Forge is shortly before to be 
released and has now all solvers in C resp. Fortran.

2) make a back-call from your C model to the R function approx. This is 
the most flexible method, but makes your code considerably slower. See 
"Writing R Extensions" how to do this.

3) write your own linear interpolation function in C. This is quite 
simple as even R uses the bisection method (in approx1, as you correctly 
identified). So take approx1, make a copy and simplify it to your needs.

Solution (1) needs the least effort in C programming and, in my opinion, 
well behaving ODE systems that require considerable interpolation effort 
for forcing data are one of the very few cases where rk4 may have 
retained its ecological niche. But even in such cases there is still the 
problem of unknown numerical accuracy.

The best method is certainly (3) and we would be very glad to integrate 
such a function (or a documented example) into the deSolve or simecol 
package.

Thomas Petzoldt


-- 
Dr. Thomas Petzoldt
Technische Universitaet Dresden
Institut fuer Hydrobiologie        thomas.petzoldt at tu-dresden.de
01062 Dresden                      http://tu-dresden.de/hydrobiologie/
GERMANY

http://tu-dresden.de/Members/thomas.petzoldt

Dear Glenn, dear list,

this is just a short notice, that a new version 1.5 of package deSolve 
was released yesterday. It now supports the feature requested below. 
Details are documented in the package vignette "Writing Code in Compiled 
Languages" that comes with the package and is also available online:

http://cran.r-project.org/web/packages/deSolve/vignettes/compiledCode.pdf

See section 6. Using forcing functions.

Thomas Petzoldt (co-author of deSolve)


Glenn Woodart wrote:
> Dear list
> 
> The deSolve package allows you to specify the model code in C or Fortran.
> Thanks to the excellent vignette this works fine. However I have not yet
> managed to use forcing functions in C code.
> 
> In pure R code this works very well with approxfun() specified outside the
> model:
[...]