[R-sig-ME] [FORGED] Subjecting data to differential equations deSolve

[Rolf Turner]

Maybe I'm obtuse (some would say that there's no "maybe" about it! :-) ) 
but I cannot for the life of me see what this question has to do with 
mixed models.

Why not ask it on plain vanilla r-help???

cheers,

Rolf Turner

On 8/11/19 2:12 AM, Torsten Hauffe wrote:
> Dear list,
> 
> My apologies if this question doesn't fit the scope of the mailing list.
> 
> I try to solve some differential equations with the deSolve::ode and need
> to feed some values at certain time-points into the equations.
> 
> Anyone has a guess why the results are off by factor 10?
> 
> Thanks for any input and cheers,
> Torsten
> 
> Here is a simplified reproducible example:
> 
> # Time of observations
> ObsTime <- c(1.9, 2.0, 2.1, 3.9, 4.0, 4.1, 9.9, 10.0, 10.1, 23.9, 24.0,
> 24.1)
> 
> # Value of the observations
> Obs <- c(0, 6, 0, 0, 3, 0, 0, 3, 0, 0, 3, 0)
> 
> # following the second example of ode(), it needs a function from where
> values could be obtained as function of time
> obs3Func <- approxfun(ObsTime, Obs, method = "constant", rule = 2)
> 
> # Subjecting the time to the approx-function works
> obs3Func(ObsTime)
> 
> # A simple function that I expect to return the cumulative sum of Obs (i.e.
> 15)
> Test_rhs <- function(t, x, parms, obsFunc)
> {
>    with(as.list(c(parms, x)), {
>      Sa <- obsFunc(t)
>      list(Sa * s)
>    })
> }
> 
> Test <- ode(y = 0,
>              times = seq(0, 65, by = 0.1),
>              func = Test_rhs,
>              parms = list(s = 1),
>              obsFunc = obs3Func)
> 
> # Using s = 1, the results are off (my expectation) by factor 10
> plot(Test)
> 
> # With s = 10 the result is as expected but I don't know why
> Test2 <- ode(y = 0,
>              times = seq(0, 65, by = 0.1),
>              func = Test_rhs,
>              parms = list(s = 1),
>              obsFunc = obs3Func)
> plot(Test2)
> # Changepoints & values
> Time <- c(2, 4, 10, 24)
> abline(v = Time, lty = 2, col = "grey")
> points(ObsTime, cumsum(Obs), pch = 19)