Skip to contents

Getting started with odelia

Source: vignettes/odelia.Rmd
odelia.Rmd

This vignette walks through the minimal setup for solving an ODE system with odelia, using the classic Lorenz system as the worked example. The Lorenz system is built into the package, so everything below runs against library(odelia) directly — no C++ compilation needed.

The core abstractions

A simulation in odelia is assembled from a few pieces:

  • a system — the ODE model itself (parameters, state, and a rule for the rates dy/dt);
  • a control object (OdeControl) — the solver’s tolerances and step sizes;
  • a solver (or runner) — which drives the adaptive RK4-5 stepper forward over the times you request and collects the solution history.

Define the system

The Lorenz system has three state variables and three parameters (sigma, R, b). We create it and set an initial state. The Lorenz system is autonomous (its rates do not depend on time directly), so no time is required.

lz <- LorenzSystem$new(sigma = 10, R = 28, b = 8 / 3)

# Inspect stored parameters
lz$pars()
#> [1] 10.000000 28.000000  2.666667

# Set the initial state (x, y, z)
lz$set_state(c(1, 1, 1))

# We can query the current state and rates directly
lz$state()
#> [1] 1 1 1
lz$rates()
#> [1]  0.000000 26.000000 -1.666667

Solve the system

Create a control object (the defaults are sensible) and a solver, then advance with adaptive time stepping over the times we want output at.

ctrl <- OdeControl$new()

runner <- Lorenz_Solver$new(lz$ptr, ctrl$ptr)
runner$advance_adaptive(seq(0, 100, by = 0.01))

out <- runner$history()
head(out)
#> # A tibble: 6 × 7
#>    time     x     y     z  dxdt  dydt    dzdt
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>
#> 1  0     1     1    1      0     26   -1.67  
#> 2  0.01  1.01  1.26 0.985  2.47  26.1 -1.35  
#> 3  0.02  1.05  1.52 0.973  4.75  26.8 -0.997 
#> 4  0.03  1.11  1.80 0.965  6.91  28.1 -0.583 
#> 5  0.04  1.19  2.09 0.962  9.02  30.0 -0.0858
#> 6  0.05  1.29  2.40 0.964 11.1   32.4  0.520

history() returns a tibble with one row per requested output time. Because the stepper is adaptive, the solver may take more internal steps than there are output rows — the internal evaluation times are available via runner$times().

Plot the attractor 

plot(out$x, out$z,
  type = "l",
  xlab = "x", ylab = "z",
  main = "Lorenz attractor (via odelia)"
)

Controlling the solver

OdeControl exposes the solver’s tuning knobs — absolute/relative tolerances, state and derivative scaling, and minimum/maximum/initial step sizes. You can inspect the current settings or tighten them:

ctrl$get_controls()
#> $tol_abs
#> [1] 1e-08
#> 
#> $tol_rel
#> [1] 1e-08
#> 
#> $a_y
#> [1] 1
#> 
#> $a_dydt
#> [1] 0
#> 
#> $step_size_min
#> [1] 1e-08
#> 
#> $step_size_max
#> [1] 10
#> 
#> $step_size_initial
#> [1] 1e-06

# Tighten tolerances for a more accurate (but slower) solve
ctrl$set_tol_abs(1e-8)
ctrl$set_tol_rel(1e-8)

Choosing an integration method

odelia ships two steppers, chosen with the method argument when you create a solver:

  • "rkck" (default) — an explicit adaptive Cash-Karp Runge-Kutta 4(5) method. Fast and accurate for non-stiff problems like Lorenz.
  • "rodas" — an implicit, adaptive RODAS4(3) Rosenbrock method for stiff problems. Each step does one exact-Jacobian evaluation (via automatic differentiation) plus a handful of linear solves, so the step size is limited by accuracy rather than stability. On a stiff problem the explicit method is forced into tiny steps (or fails), while RODAS keeps large steps.
runner_stiff <- Lorenz_Solver$new(lz$ptr, ctrl$ptr, method = "rodas")

Both steppers share the same adaptive step-size controller and OdeControl tolerances, so you can switch between them without changing anything else.

To use method = "rodas" with your own system, the system’s C++ class must provide a template<class U> System<U> rebind() const method (a one-liner that copies the system with its scalar type swapped) so its right-hand side can be differentiated for the Jacobian; see the Lorenz example. The implicit stepper currently runs on the passive (non-AD) solver; differentiating a fit through RODAS is planned (odelia issue #35).

Where next