rodeint_integrate: Integrate a System of ODEs
In richfitz/rodeint: R interface to odeint

Description Usage Arguments Details Author(s) See Also Examples

Integrate a system of ordinary differential equations (ODEs). This is the low-level interface that corresponds exactly to the ideint interface. An easier to use, more R-ish, interface is available using make_integrate.

integrate_const(stepper, ode_system, y, t0, t1, dt, save_state = FALSE)

integrate_n_steps(stepper, ode_system, y, t0, dt, n, save_state = FALSE)

integrate_adaptive(stepper, ode_system, y, t0, t1, dt, save_state = FALSE)

integrate_times(stepper, ode_system, y, times, dt)

`stepper`	A `stepper` object, created by `make_stepper`
`ode_system`	A system of ODEs, created by `ode_system`
`y`	Initial conditions
`t0`	Time to start the integration
`t1`	Time to finish the integration
`dt`	Step size
`n`	Number of steps to take (`integrate_n_steps` only)
`times`	Vector of times (`integrate_times` only)
`save_state`	Logical: return information about intermediate points as an attribute. Not applicable for `integrate_times`, which always returns information about intermediate points.

There are four different integration functions here. See the odeint documentation for more detail about the precice differences. The summary, though, is:

integrate_const "Equidistant observer calls"
integrate_n_steps "Integrate a given number of steps"
integrate_adaptive "Observer calls at each step"
integrate_times "Observer calls at given time points"

This page provides much more detail. The functions here are direct wrappers of the odeint functions, so the interpretation is identical.

If you just want the end-points, integrate_adaptive is probably the function to use. If you want to know intermediate values at particular times too, use integrate_times. The functions integrate_const and integrate_n_steps are provided more for completeness than for utility, but they may be useful in specific applications.

Rich FitzJohn

make_stepper for building steppers and for information about possible algorithms, ode_system for building a sytem of ODEs, and make_integrate for a higher-level interface to these functions.

The odeint documentation also has useful reference information on differences between the four integrate functions, and how dense output steppers work.

## Picking up with the harmonic oscillator from the ode_system
## example:
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)

## We also need a stepper (see the ?stepper help page)
s <- make_stepper("dense", "runge_kutta_dopri5")

## Suppose we want output at these times:
times <- seq(0, 10, length=101)

## Starting from this initial state:
y0 <- c(0, 1)

## Initial step size guess:
dt <- 0.01

## This intergrates the system (all integrate_ functions take stepper,
## system, y as the first theree arguments).
y <- integrate_times(s, sys, y0, times, dt)

## Here are the variables changing over time:
matplot(times, y, type="l")

## The result also has useful attributes.  The "t" attribute is the
## times (this is the same as 'times' here)
attr(y, "t")

## The "steps" attributes contains the number of steps the system took
attr(y, "steps")

## Here, the number is lower than the number of times!  This is
## because we used a "dense output stepper" which can interpolate back
## over times that have passed.  See the "See also" section for a link
## to the odeint help about this.

## The "y" attribute contains the final system state - here it will be
## the same as the last row of 'y' itself
attr(y, "y")
y[nrow(y),]

## If you don't care about intermediate times, the
## "integrate_adaptive" function is probably the right one to use.
y_a <- integrate_adaptive(s, sys, y0, min(times), max(times), dt)
y_a

## Here we jump all the way to the final time as quickly as possible.
y_a
attr(y, "y")

## It's possible to remember the steps taken by usin the "save_state"
## argument.
y_a <- integrate_adaptive(s, sys, y0, min(times), max(times), dt,
                          save_state=TRUE)

## Now y_a contains three attributes, the same as returned by
## integrate_times.  "t" is the vector of times that the stepper
## stopped at:
attr(y_a, "t")

## "y" contains the y values at these points, wherever they are:
attr(y_a, "y")

## Here is the dense output (thick lines) with the actual observed
## spots overlaid.  You can see the step size change at the beginning
## of the problem until the stepper works out how fast it can go while
## achieving the required accuracy.
matplot(times, y, type="l", col=c("grey", "pink"), lty=1, lwd=10,
        xlim=c(0, 2))
matlines(attr(y_a, "t"), attr(y_a, "y"), type="o", pch=19)

## Manually specifying the stepper and step size and passing around 6
## or so parameters all the time gets old, so there is a function
## "make_integrate" to help collect everything together.  See
## ?make_integrate for more information.