Solve a user-written ordinary differential equation

Solves dy/dt = f(t, y) forward from start_time with a ported Numerics Runge-Kutta method (P2 "math extras"): fixed-step second- or fourth-order Runge-Kutta over an equally-spaced grid (method = "rk2"/"rk4" with end_time/time_steps), the fourth-order method's single-step form (method = "rk4" with dt alone), or one of the two adaptive-step-size methods, Runge-Kutta-Fehlberg or Runge-Kutta-Cash-Karp (method = "rkf"/"cash_karp" with dt/dt_min).

Usage

ode_solve(
  f,
  initial_value,
  start_time,
  end_time = NULL,
  time_steps = NULL,
  dt = NULL,
  dt_min = NULL,
  method = c("rk4", "rk2", "rkf", "cash_karp"),
  tolerance = NULL
)

Arguments

f: a function taking two numbers (t, y) and returning dy/dt, one number.
initial_value: the value of y at start_time.
start_time: the time to start integrating from.
end_time, time_steps: the end time and the number of equally-spaced points between start_time and end_time (inclusive of both ends). Required together for method = "rk2"; for method = "rk4", supplying them selects this array form over the single-step form (see dt below) – their PRESENCE, not a separate flag.
dt: the step size. For method = "rk4" without end_time/time_steps, the single step to advance by (the ODE is solved once, at start_time + dt). For method = "rkf"/ "cash_karp", the maximum internal step size, required together with dt_min.
dt_min: the minimum internal step size for method = "rkf"/"cash_karp", required together with dt.
method: one of "rk4" (the default), "rk2", "rkf", or "cash_karp".
tolerance: the absolute error tolerance for method = "rkf"/"cash_karp" alone. NULL, the default, leaves the ported routine's own default (1e-3) in force.

Value

for method = "rk2", or "rk4" with end_time/time_steps: a numeric vector of length time_steps, the solution at each grid point (y[1] is initial_value). For every other case: a single number, the solution at start_time + dt ("rk4"'s single-step form) or at start_time + dt under the adaptive step control ("rkf"/"cash_karp").

Examples

# dy/dt = y, y(0) = 1: the exact solution is y(t) = exp(t).
y <- ode_solve(function(t, y) y, initial_value = 1, start_time = 0, end_time = 1,
               time_steps = 100)
y[100]
#> [1] 2.718282