[R] External signal in ODE written in C (using deSolve and approx1?)

[Thomas Petzoldt]

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