integrate_sys: Integrate an ODE system using ODEINT

Description

Numerically integrates an ODE system defined in R

Usage

integrate_sys(sys, init, duration, step_size = 1, start = 0,
  adaptive_obs = FALSE, observer = function(x, t) x, atol = 1e-06,
  rtol = 1e-06)

Arguments

sys

a function with signature function(x, t)

init

the initial conditions

duration

time-span of the integration

step_size

the initial step size (adjusted internally)

start

the starting time

adaptive_obs

if true, call observer after each adaptive step

observer

a function with signature function(x, t) returning values to store in output

atol

absolute error tolerance

rtol

relative error tolerance

Value

A data frame, NULL if no samples were recorded and a very complicated list-of-lists if the observer returned objects of different length.

Details

The system will be integrated from start to start + duration. The method is an error controlled 5th-order Dormand-Prince. The time step will be adjusted to within error tolerances (1e-6 absolute and relative). The observer can return arbitrary data in any form that can be coerced to a list. This could be a single scalar value (no need to wrap the return with list!) or a list containing heterogeneous types. These will be inserted into the columns of the returned data frame. If the observer function returns a zero-length object (NULL or list()), then nothing will be recorded. You can use the t argument to selectively sample the output.

Examples

Run this code

## Not run: ------------------------------------
# # Lotka-Volterra predator-prey equations
# LV.sys = function(x, t)
# {
#    c(x[1] - 0.1 * x[1] * x[2],
#      0.05 * x[1] * x[2] - 0.5 * x[2])
# }
# null_rec = function(x, t) NULL
# system.time(integrate_sys(LV.sys, rep(1, 2), 1e3, observer = null_rec))
# named_rec = function(x, t) c(Prey = x[1], Predator = x[2])
# x = integrate_sys(LV.sys, rep(1, 2), 100, 0.01, observer = named_rec)
# plot(x[, 2:3], type = "l", lwd = 3, col = "steelblue")
# Sys.sleep(0.5)
# 
# # Lorenz model from odeint examples
# Lorenz.sys = function(x, t)
# {
#  c(10 * (x[2] - x[1]),
#    28 * x[1] - x[2] - x[1] * x[3],
#    -8/3 * x[3] + x[1] * x[2])
# }
# system.time(integrate_sys(Lorenz.sys, rep(1, 3), 1e2, obs = null_rec))
# x = integrate_sys(Lorenz.sys, rep(1, 3), 100, 0.01)
# plot(x[, c(2, 4)], type = 'l', col = "steelblue")
## ---------------------------------------------