Dense Output in rxode2: Continuous Interpolation for Fast Solving

What is Dense Output and what does `dense = TRUE` do?

When solving differential equations, adaptive-step solvers choose their step sizes dynamically to satisfy error tolerances. Normally, to evaluate the system states at your exact observation grid, the solver must adjust its step sizes and stop precisely at each observation time. This constant stopping and restarting restricts step sizes and degrades performance–especially on fine observation grids.

Dense Output resolves this by allowing the solver to step forward at its own optimal pace. Instead of stopping at observation times, the solver records polynomial coefficients for each accepted step. We then use these cheap, continuous polynomials to interpolate the system states at any observation times falling within that step.

In rxode2, dense output is supported across five solvers:

Method	Type	Stiff-Safe	Dense Polynomial / Method
`"dop853"`	Explicit Runge-Kutta 8(5,3)	No	7th-order `contd8` polynomial (Hairer native)
`"dop5"`	Explicit Runge-Kutta 5(4)	No	4th-order Shampine continuous extension
`"bs"`	Bulirsch-Stoer	No	Extrapolation polynomial dense output (odeint)
`"ros4"`	Implicit Rosenbrock 4	Yes	Implicit Rosenbrock dense output (odeint)
`"cvode"`	SUNDIALS CVODE (BDF/Adams)	Yes	CVODE native `CVodeGetDky` interpolation

The standard path (`dense = FALSE`) vs. the dense path (`dense = TRUE`)

dense = FALSE (default): The solver is stopped at each observation time, advancing from one event to the next. The adaptive step-size control re-enters from scratch at every observation point.
dense = TRUE (dop853): The native Hairer C solver runs the entire segment from one dose to the next in a single dop853() call. A solout callback fires after each accepted step and calls contd8() to evaluate the 7th-order dense polynomial at every observation time falling within that step. This is not Boost.Odeint-based; it uses Dormand & Prince’s original Fortran algorithm ported to C.
dense = TRUE (odeint-based solvers: dop5, bs, ros4): An odeint dense stepper is initialized at the start of each inter-dose interval. The event loop visits each observation individually; for any observation already within the stepper’s last completed step, calc_state() recovers the solution by polynomial interpolation without new derivative evaluations. The stepper state (current time, step size, history polynomial) is preserved across consecutive observations so the solver never restarts within an inter-dose interval.
dense = TRUE (cvode): Uses a different but complementary mechanism. CVODE advances in CV_ONE_STEP mode toward the simulation end, taking natural adaptive steps that can span many consecutive observation times. At each observation, CVodeGetDky interpolates the solution using CVODE’s internal Nordsieck polynomial – the same polynomial CVODE already maintains for its own error control. This avoids the CVodeReInit call that the standard path requires at every observation, preserving the BDF order-selection state, the step-size history, and any cached Jacobian factorization across the entire inter-dose interval. CVodeReInit is called only when a dose event, reset event, or steady-state calculation introduces a true state discontinuity.

The Critical Role of `hmax` – Handled Automatically

The default hmax (maximum step size) in rxode2 is computed from the average spacing between events in your event table. If your grid has 0.1 h spacing, hmax ~ 0.1 h. Under this default, the solver is prevented from taking steps larger than 0.1 h, which defeats the performance benefit of dense output.

To realize the full speedup of dense = TRUE, the solver must be allowed to take steps larger than the observation spacing. This requires setting hmax = NULL (unlimited step size).

Because hmax = NULL is almost always desired when dense = TRUE is requested, rxode2 configures this automatically, but you can override this by passing an explicit numeric value if you want to restrict the maximum step size:

# Explicit hmax overrides the automatic NULL
rxSolve(mod, ev, method = "dop853", dense = TRUE, hmax = 1.0)

Benchmarking Dense Output Solvers

We will compare the standard (dense = FALSE) and dense (dense = TRUE) paths across the different supported solvers using a typical non-stiff oral absorption model.

Setup for the Benchmark

library(rxode2)
library(microbenchmark)

We define a standard two-compartment oral absorption model:

mod <- function() {
  ini({
    tka <- log(0.5)
    tcl <- log(4)
    tv  <- log(70)
    tq  <- log(2)
    tvp <- log(40)
  })
  model({
    ka  <- exp(tka)
    cl  <- exp(tcl)
    v   <- exp(tv)
    q   <- exp(tq)
    vp  <- exp(tvp)
    d/dt(depot)   <- -ka * depot
    d/dt(central) <-  ka * depot - (cl + q) / v * central + q / vp * periph
    d/dt(periph)  <-  q / v * central - q / vp * periph
    cp <- central / v
  })
}

# Warm up compilations
invisible(rxSolve(mod, et(amt = 100) |> et(0:5), method = "dop853"))
invisible(rxSolve(mod, et(amt = 100) |> et(0:5), method = "dop5"))
invisible(rxSolve(mod, et(amt = 100) |> et(0:5), method = "bs"))
invisible(rxSolve(mod, et(amt = 100) |> et(0:5), method = "ros4"))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
invisible(rxSolve(mod, et(amt = 100) |> et(0:5), method = "cvode"))

Correctness Check

Let’s verify that the dense interpolation yields numerically equivalent results to standard solving:

doses <- et(amt = 100, time = 0, ii = 24, addl = 6)
ev    <- doses |> et(seq(0, 168, by = 0.1))

# Compare standard vs. dense for dop853
f_std   <- rxSolve(mod, ev, method = "dop853", dense = FALSE, hmax = NULL)
f_dense <- rxSolve(mod, ev, method = "dop853", dense = TRUE)

max(abs(f_std$cp - f_dense$cp))
#> [1] 3.774758e-15

The maximum difference is at the level of solver tolerance, confirming correctness.

Comparison Across Sampling Grids and Solvers

Below we run a benchmark comparing the median execution times (in milliseconds) for the standard and dense paths across three different sampling grid densities: sparse (1 h), medium (0.1 h), and dense (0.01 h).

run_benchmark <- function(method, ev, dense) {
  bm <- microbenchmark(
    rxSolve(mod, ev, method = method, dense = dense, cores = 1),
    times = 20L
  )
  median(bm$time) / 1e6 # nanoseconds -> milliseconds
}

for (m in c("dop853", "dop5", "bs", "ros4", "cvode")) {
  cat(sprintf("\n=== Solver: %s ===\n", m))
  for (by in c(1.0, 0.1, 0.01)) {
    ev_grid <- doses |> et(seq(0, 168, by = by))
    nobs    <- nrow(ev_grid[ev_grid$evid == 0, ])

    std_t <- run_benchmark(m, ev_grid, FALSE)
    den_t <- run_benchmark(m, ev_grid, TRUE)

    cat(sprintf("  Obs spacing = %4.2f h (n = %5d): Std = %5.1f ms | Dense = %5.1f ms | Speedup = %.2fx\n",
                by, nobs, std_t, den_t, std_t / den_t))
  }
}

Key Solver Observations:

"dop853" & "dop5": Highly optimized explicit Runge-Kutta dense output methods. They exhibit consistent speedups of 10% to 25% on dense grids (by = 0.01 h) because Runge-Kutta polynomial interpolation is extremely fast compared to full derivative evaluations.
"bs" (Bulirsch-Stoer): Neville-Aitken interpolation dense output. It is highly robust and performs beautifully on smooth, high-precision non-stiff systems.
"ros4" (Rosenbrock 4): The implicit stiff-safe dense solver. Each step requires a Jacobian evaluation and an LU-factorized linear-system solve (one per stage), but Rosenbrock is linearly implicit – there is no Newton iteration. On stiff systems with dense observation grids, "ros4" with dense = TRUE yields speedups by taking large natural steps and interpolating, rather than being forced to stop and restart step-size selection at every observation.
"cvode" (SUNDIALS CVODE BDF): The stiff-optimized dense solver using CVODE’s native CVodeGetDky interpolation. Unlike the odeint-based dense solvers, CVODE advances in CV_ONE_STEP mode, taking natural adaptive steps across the entire inter-dose interval. The speedup over the standard path comes from avoiding CVodeReInit at every observation, which preserves the BDF order-selection state, step-size history, and cached Jacobian factorization. Speedups are most pronounced on very fine grids with stiff systems, where the BDF multi-step history carries most of the convergence benefit.

Guidelines: When to Use `dense = TRUE`

Use dense = TRUE when all of the following hold:

Fine-grained observation grid (rich sampling): The break-even point is typically around 5-10 observations per dose interval. If you only have 1-3 observations, the overhead of setting up the dense output callbacks exceeds the savings.
No analytical linCmt() models: Dense output is not applicable for analytical linear-compartment solving. If the model uses linCmt(), rxode2 will print an informational message and silently fall back to the standard path.

Choosing the Right Dense Solver

For general non-stiff systems: Use "dop853" or "dop5". "dop5" is faster for moderate tolerance requirements, while "dop853" is optimal for high accuracy.
For very smooth non-stiff systems requiring high precision: Use "bs".
For stiff systems with moderate sampling grids: Use "ros4". It is linearly implicit, evaluates the Jacobian once per step, and the dense output avoids restarting step-size selection at each observation.
For stiff systems with very fine sampling grids (e.g., TMDD, QSP models with rich PK sampling, or mechanistic models where observations are spaced much closer than the solver’s natural step size): Use "cvode". The BDF method preserves its multi-step history polynomial and cached Jacobian factorization across observation points, so the dense path avoids repeated CVodeReInit calls that would otherwise discard accumulated convergence state.

Session Information

sessionInfo()
#> R version 4.6.0 (2026-04-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] microbenchmark_1.5.0 rxode2_5.1.3        
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.4.10           generics_0.1.4        lattice_0.22-9       
#>  [4] digest_0.6.39         magrittr_2.0.5        evaluate_1.0.5       
#>  [7] grid_4.6.0            RColorBrewer_1.1-3    fastmap_1.2.0        
#> [10] lotri_1.0.5           qs2_0.2.2             jsonlite_2.0.0       
#> [13] rxode2ll_2.0.14       backports_1.5.1       scales_1.4.0         
#> [16] textshaping_1.0.5     jquerylib_0.1.4       cli_3.6.6            
#> [19] symengine_0.2.11      crayon_1.5.3          rlang_1.2.0          
#> [22] units_1.0-1           cachem_1.1.0          yaml_2.3.12          
#> [25] otel_0.2.0            tools_4.6.0           memoise_2.0.1        
#> [28] checkmate_2.3.4       dplyr_1.2.1           ggplot2_4.0.3        
#> [31] vctrs_0.7.3           R6_2.6.1              lifecycle_1.0.5      
#> [34] fs_2.1.0              stringfish_0.19.0     htmlwidgets_1.6.4    
#> [37] ragg_1.5.2            PreciseSums_0.7       pkgconfig_2.0.3      
#> [40] desc_1.4.3            rex_1.2.2             pkgdown_2.2.0        
#> [43] RcppParallel_5.1.11-2 bslib_0.11.0          pillar_1.11.1        
#> [46] gtable_0.3.6          data.table_1.18.4     glue_1.8.1           
#> [49] Rcpp_1.1.1-1.1        systemfonts_1.3.2     xfun_0.58            
#> [52] tibble_3.3.1          tidyselect_1.2.1      sys_3.4.3            
#> [55] knitr_1.51            farver_2.1.2          dparser_1.3.1-13     
#> [58] htmltools_0.5.9       nlme_3.1-169          rmarkdown_2.31       
#> [61] compiler_4.6.0        S7_0.2.2

2026-06-09

What is Dense Output and what does dense = TRUE do?

The standard path (dense = FALSE) vs. the dense path (dense = TRUE)

The Critical Role of hmax – Handled Automatically