make_integrate: Make Integration Function
In richfitz/rodeint: R interface to odeint

Description Usage Arguments Details Author(s) See Also Examples

Helper function for binding ode_systems, steppers and integration functions together. This can be used to create a function f(y0, t) from your system of f\'(y0, t). These functions reurn new functions that have all the arguments that won't change set. For those interested, this form of programming is nicely described here, (though you do not need to understand that to use these functions).

make_integrate(ode_system, ..., stepper = NULL,
  integrate = integrate_adaptive)

make_integrate_pars(ode_system, ...)

`ode_system`	A ode_system function (`ode_system`).
`...`	Additional named arguments to bind. Setting `t0` when a system is time independent means `t1` will be a function of elapsed time, which can be useful. All integrate functions take a `dt` argument, so that's useful to bind too. You can also pass in `set_as_defaults=TRUE` and the arguments, including `stepper` and `ode_system` will simply be set as defaults allowing some tuning later.
`stepper`	A stepper object. By default the controlled `runge_kutta_dopri` stepper is used with default tolerances. If you want to change the tolerance, you must provide a different stepper object.
`integrate`	One of the integration functions. The default is `integrate_adaptive`.

The motivation here is that often we just want to specify a set of ODE solving parameters once and then solve a system many times - at different initial conditions, or over different parameter sets. This function simplifies this approach by remembering arguments passed in.

Note the opposite ordering of the ode_system and stepper arguments here compared with the rest of the package (following odeint. This is because a reasonable default stepper will be chosen if none is provided, but a system of ODEs is always required!

The function make_integrate_pars is a higher-higher order function. It returns a function that takes parameters of the system as an argument. When that function is run it returns a function with arguments bound as for make_integrate. This is a lot simpler than it sounds - see the final example.

Rich FitzJohn

integrate_adaptive, which this function wraps around, and ode_system for building a system of ODEs to integrate.

## A system of differential equations describing a harmonic oscillator:
derivs <- function(y, t, pars) {
  c(y[2],
    -y[1] - pars * y[2])
}

## Parameters of the system:
pars <- 0.5

## Build the system itself
sys <- ode_system(derivs, pars)

## To integrate this system, we need to specify a stepper, initial
## conditions, begining and end times and an initial step size:
s <- make_stepper("dense", "runge_kutta_dopri5")
y0 <- c(0, 1)
t0 <- 0
t1 <- 5
dt <- 0.01
integrate_adaptive(s, sys, y0, t0, t1, dt)

## If we wanted to run this from a couple of different starting points
## it gets quite tedious:
y1 <- c(1, 1)
y2 <- c(2, 1)
y3 <- c(3, 1)
integrate_adaptive(s, sys, y0, t0, t1, dt)
integrate_adaptive(s, sys, y1, t0, t1, dt)
integrate_adaptive(s, sys, y2, t0, t1, dt)
integrate_adaptive(s, sys, y3, t0, t1, dt)

## There is a lot of repetition here, and it's not actually that clear
## what is changing.  The make_integrate functions binds arguments
## that won't change (all arguments must be named!)
f <- make_integrate(sys, t0=t0, t1=t1, dt=dt, stepper=s)

## Our function now has only a single required argument (y) and an
## optional argument "save_state":
args(f)

## This is because all the arguments in integrate_adaptive that
## matched the list given were set to the values provided.

## This new function can then be run like so:
f(y0)
f(y1)
f(y2)
f(y3)

## Or, to get intermediate values out:
f(y0, save_state=TRUE)

## If we also wanted 't1' to be left free, just don't provide it to
## make_integrate:
f <- make_integrate(sys, t0=t0, dt=dt, stepper=s)
args(f)

## We now have a function f(y, t1).  Note that this approach has
## convered a function dy/dt (derivs) into a function simply of time
## and initial conditions.

## This approach works with integrate_times, too, by specifying that
## as the "integrate" argument:
f_at_times <- make_integrate(sys, dt=dt, stepper=s,
                             integrate=integrate_times)

## This returns a function taking "y" and "times" as arguments (no
## save_state, because integrate_times does not accept that argument.
args(f_at_times)

## It can be run like so:
tt <- seq(0, 5, length=21)
f_at_times(y0, tt)

## Finally, sometimes it is useful to easily apply the function over a
## series of parameters.  To do this, use make_integrate_pars.  The
## calling sequence is the same as above.
g <- make_integrate_pars(sys, t0=t0, dt=dt, stepper=s)

## This returns a function of 'pars' only:
args(g)

## When run with parameters, this function returns a function of y and t1:
args(g(pars))

## So
g(pars)(y0, t1)

## sets parametes to pars, and then integrates with initial conditions
## y0 up until time t1, using t0, dt and the stepper provided above.

## This makes it easy to explore how the system changes over parameter
## space
g <- make_integrate_pars(sys, y=y0, t0=t0, t1=1, dt=dt, stepper=s)

g(pars * 0.9)()
g(pars)()
g(pars * 1.1)()