[R-sig-dyn-mod] Skinny square waves (John Harrold)

[Thomas Petzoldt]

Hello,

thanks for this interesting discussion. A model that skips external 
forcings was indeed my first problem report regarding "odesolve" to 
Woody Setzer long time ago and a beginning of a great adventure. We know 
now different ways how to handle this, all having their pros (+) and 
cons (-). Here an incomplete list of what we do in our lab:

1) set times to match external forcings
- slows down simulation,
- produce lots of output data

2) set hmax to limit the internal time step
- slows down simulation
+ much better than (1)
+ output data still under control (better than 1)

3) use another solver, e.g. Runge-Kutta
- less efficient than any of the odepack solvers
+ steps are known before (incl. intermediate steps),
    so that data can be supplied in correct resolution

4) make the external signal more smooth, e.g. by using (steep) sigmoidal 
instead of rectangular signals.

In all these cases, forcing data needs to be interpolated, either in the 
model function (1-4) or beforehand (3). It is also possible, to pass a 
"hook function" to the derivative function for interpolation. This can 
enhance efficiency for simulations with big amount of external data, 
especially if the interpolator has "knowledge" about the data structure.

It is also good to know that most solvers are able to interpolate the 
states. In case of doubt, we set defaults towards robustness instead of 
maximum efficiency. There are still some options for fine-tuning, and in 
case of doubt, set hmax.

5) always check precision (atol, rtol) especially if your states occur 
at very small or large scales.

6) use events!
+ elegant to implement (thanks to Karline)
- slows down simulation (but not as much as (1) or (2)), because the 
solver does "internal stop-and-go" between integration steps.
+ can be used in combination with 1-3

Events are good if external manipulations are very short, to be 
considered as "instantaneous change of state" while processes are quite 
slow during the rest of the time, e.g. if you want to follow a patient 
over months that gets an 1min injection every morning.

7) use "external stop-and-go" mode by chaining calls to ode()
+ external forcing data can be kept constant like parameters during each 
step
+ optimum performance of the solver (no intermediate steps)
- little bit more programming (e.g. loops)
- some computational overhead, due to repeated invocation of ode()

... see also the vignette of package "rodeo" about this: 
https://cran.r-project.org/web/packages/rodeo/vignettes/rodeoVignette.pdf#14

To sum up: R is a full programming language and many roads lead to Rome.

Thomas


Am 09.03.2017 um 14:45 schrieb Setzer, Woodrow:
> Can you arrange to code this using deSolve's "events"? You would have to code the model so that it is state variables that change abruptly.
>
> Woody
>
> -----Original Message-----
> From: R-sig-dynamic-models [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of Karline Soetaert
> Sent: Thursday, March 09, 2017 2:22 AM
> To: Special Interest Group for Dynamic Simulation Models in R <r-sig-dynamic-models at r-project.org>
> Subject: Re: [R-sig-dyn-mod] Skinny square waves (John Harrold)
>
> Hi John,
>
> Usually in such cases, the problem is that the default solver is too efficient and thus takes too large timesteps. You could either try a less efficient solver that by default takes smaller time steps or restrict the maximal time step (argument hmax).
>
> Karline
>
> -----Original Message-----
> From: R-sig-dynamic-models [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of John Harrold
> Sent: donderdag 9 maart 2017 2:56
> To: Special Interest Group for Dynamic Simulation Models in R <r-sig-dynamic-models at r-project.org>
> Subject: Re: [R-sig-dyn-mod] Skinny square waves (John Harrold)
>
> Howdy Steve,
>
> I've run into this problem in the past in Matlab. I tested it in R at the time and I didn't seem to have the same problem. But it seems to be somewhat situational :).
>
> I had a similar solution in Matlab to what you did. The R-equivalent of adding events at each of the switching times. So as a test I added events at each of the switching times. Basically I told ode to add 0 to the first state at each of those times, and that seemed to do it. These events should force the ode solver to stop and restart at each time point. Now I just need to automate this in my simulation routines.
>
> I think it would be good if we could just provide the ode function a list of times and have it do this automatically. This way it would make it easier for me :). And it would also allow this to be optimized better than my poor coding skills will probably allow. Is this the appropriate list to make these suggestions?
>
> John
>
> On Wed, Mar 8, 2017 at 8:33 PM, Stephen Paul Ellner <spe2 at cornell.edu>
> wrote:
>
>> John,
>>
>> the efficient deSolve solvers can't deal with such rapid on-off
>> changes (even if you make the changes continuous). And even with start
>> and stop times in the output times, the solvers may well leap over the
>> 30 minute periods when the forcing is 'on' and give you interpolated
>> output values for output times during the 'on' periods.
>>
>> The way I've always dealt with such situations is to solve the 'on'
>> and 'off' periods separately.
>>
>> for(j in 1:whatever) {
>>         solve from stop[j] to start[j+1] with forcing=0
>>         solve from start[j+1] to stop[j+1] with forcing = value }
>> storing the results along the way in one big data frame or matrix.
>>
>> Steve
>>
>> Stephen P. Ellner
>> Department of Ecology & Evolutionary Biology, E339 Corson Hall Cornell
>> University, Ithaca NY 14853-2701

>
> --
> -------------------------------------
> John M. Harrold
> john.m.harrold _at_gmail
> -------------------------------------

-- 
Dr. Thomas Petzoldt
Technische Universitaet Dresden
Faculty of Environmental Sciences
Institute of Hydrobiology
01062 Dresden, Germany

E-Mail: thomas.petzoldt at tu-dresden.de
http://tu-dresden.de/Members/thomas.petzoldt