integrate_sys: Integrate an ODE system using ODEINT

In odeintr: C++ ODE Solvers Compiled on-Demand

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

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.

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.

Author(s)

Timothy H. Keitt

See Also

compile_sys

Examples

## 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")

## End(Not run)

Example output

   user  system elapsed 
  0.705   0.016   0.735 
   user  system elapsed 
  2.804   0.004   2.841 

odeintr documentation built on May 2, 2019, 2:08 a.m.