implementation
plant casts the entire patch as a system of ordinary differential equations (ODEs) and steps it through time with classic ODE methods. (For details on the equations themselves see the size-structured PDE page.) This page first describes the Runge–Kutta solver and the Control() parameters that govern it, then works through the trade-off between the default adaptive stepper and a fixed-step forward-Euler mode using a work–precision diagram.
The plant model uses an explicit Runge-Kutta ODE stepper to integrate the ordinary differential equations that describe the system. Specifically, we use the Runge-Kutta Cash-Karp 4-5 (RK45) algorithm, with adaptive time step.
In brief, this means that the solver takes a step through time, and then evaluates the error in the step. Error is estimated by comparing two possible steps: one using 4 sub-intervals, and one using 5 sub-intervals. By comparing the step made with more sub-intervals, the algorithm can estimate the non-linearity in the system and error in the step. If the error is too high, the solver will attempt to take a smaller step. If the error is well below a target value, the solver will attempt to take a larger step. Importantly, the 5 step reuses information from the 4 step and thereby only requires a single extra call on rates of change, so the RK45 algorithm is much more efficient than alternatives like the RK48. The math of RK45 stepping has long been established, and the method is available in most solvers for ODE systems. You can read more about Runge-Kutta methods and adaptive time steps on Wikipedia.
The RK45 is an explicit method, as opposed to an implicit method. As wikipedia explains “Explicit methods calculate the state of a system at a later time from the state of the system at the current time, while implicit methods find a solution by solving an equation involving both the current state of the system and the later one”. The advantage of explicit methods are that they are simple to implement and fast, particularly when allowing adaptive time steps. The disadvantage is that they are less stable than implicit methods when they encounter stiff systems. As described in wikipedia
In practice, we’ve found the RK45 to be extremely stable and fast for solving the size-density distribution. However, there are circumstances where the system becomes stiff and the solver struggles. E.g.
Any extensions should therefore be wary of stiff behaviours.
The characteristic method used to solve density requires regular introduction of new nodes, which entails a resizing of the ODE system as the system progresses. To make this play nicely with the ODE stepping, Rich FitzJohn designed a custom ODE solver class based on the RK45 method. The code is built into the model, with some details ported from the GNU Scientific Library (Galassi, 2009).
The main interface is described in the class plant::ode::Solver described in the file ode_solver.h. (See folder inst/include/plant/ode_solver where this and related class are implemented.) The class is established using a template System class, meaning you can effectively cast any variety of classes as an ODE Solver. (In a similar way, you can create a vector of any variety of data types.)
From ode_solver.h:
template <class System>
class Solver {
public:
Solver(const System& system, OdeControl control_);
...
};The main solver for our system is the class SCM, (in the file sch.h). This creates an ode solver based on patch type. Simultaneously, the class SCM owns an actual instance of the patch, which contains species and environment.
patch_type patch; # Instance of patch
ode::Solver<patch_type> solver; # ODE solverAll we need to do to advance the system is make calls like the following. These advance the system, using the patch to estimate rates of change. Moreover, we can take steps of a specified fixed size, or an adaptive step where the solver determines the spacing.
solver.advance_fixed(patch, e.times);
solver.advance_adaptive(patch, e.time_end());In the background, the solver then passes information (rates of change, states) between the actual patch (which contains all the biological information informing rates of change) and ode solver (which uses the RK45 math to calculate future states from rates of change).
Additionally, we have the following function, which resizes the ODE system after a node has been introduced.
solver.set_state_from_system(patch)For the above to work, the patch class (see patch.h) must have some essential functions used by the ODE system.
// * ODE interface
// Caluclate size of ode system (number of equations).
size_t ode_size() const;
// How many auxiallary variables are we tracking. These are being collected but
// are not part of core ode system
size_t aux_size() const;
double ode_time() const;
// Retrieve ode state from patch and save into the ode solver
ode::iterator ode_state(ode::iterator it) const;
// Retrieve ode rates from patch and save into the ode solver
ode::iterator ode_rates(ode::iterator it) const;
// Retrieve auxillary variables and save into the ode solver
ode::iterator ode_aux(ode::iterator it) const;
// Set state of patch, based on estimate of future state estimated by the solver
ode::const_iterator set_ode_state(ode::const_iterator it, double time);Note the use of iterators. These enable state from a complicated data structure like a patch, to be converted into a vector, as needed by ode solver. And vice versa.
As templated C++, we can use the ode::Solver class to step any “system” that has the requisite functions. In addition to patches, we also use the ode solver to run individual plants, given an environment. This is achieved by introducing a wrapper class ode::Runner (in ode_runner.h), which according to its author simplifies “the difficulty of ownership semantics around the solver and object.”
Like the Solver class, the Runner class is templated, meaning we can setup a variety of systems.
template <typename T>
class Runner {
public:
Runner(T obj_, OdeControl control) : obj(obj_), solver(obj, control) {}
};The class plant::tools::IndividualRunner uses this approach to create a system for stepping individual plants to a given size. This is available via the R functions grow_individual_..., which in turn uses the OdeRunner.
The test suite also shows how to use the OdeRunner to step the classic Lorenz system.
The RK45 ode solver is controlled by a set of parameters that determine how the solver steps through time.
As a simple example, we can step a system using the default values for ODE stepping:
library(plant)
# configure a patch containing a single species
p0 <- scm_base_parameters("FF16")
p0$max_patch_lifetime <- 5
p <- expand_parameters(trait_matrix(0.0825, "lma"), p0)
# Controls on numerical methods
ctrl <- Control()
# run
out <- run_scm(p, ctrl = ctrl)
# number of steps
length(out$ode_times)[1] 97The precise times used by the solver are stored in out$ode_times, so the length of this object indicates the number of times used.
The number and timing of steps taken by the solver will vary with the controls on the ode solver. Looking at the ctrl object, we can extract the items beginning with ode_, which control the ODE solver:
ctrl %>% subset(grepl("(?<!n)ode_", names(.), perl = TRUE))$ode_step_size_initial
[1] 1e-06
$ode_step_size_min
[1] 1e-06
$ode_step_size_max
[1] 5
$ode_tol_rel
[1] 1e-04
$ode_tol_abs
[1] 1e-04
$ode_a_y
[1] 1
$ode_a_dydt
[1] 0These have the following effect
ode_step_size_initial: The initial step size used by the solver in the first step.ode_step_size_min, ode_step_size_max: The minimum step size permitted. Any attempt by the solver to use a smaller value will be curtailed at this value.ode_tol_abs, ode_tol_rel: When taking a step, the absolute and relative error for any single equation that is allowed. If the error is too high, the solver will attempt to take a smaller step. If the error is well below this value, the solver attempts to expand the step size in the next step.ode_a_y, ode_a_dydt: Weightings on y and dy/dt, in calculation of allowable relative error (see OdeControl::errlevel).Varying these will potentially change the number of steps taken by the solver.
However, note that unlike a regular ode solver, the number of steps taken is also influenced by the timing of node introduction. The ode stepper will only expand the time step up until the next node introduction event, for any species. So the node scheduling sets an upper limit on the size of time steps. In the default setup, the node schedule primarily determines the timesteps in the first months of the simulation. After that, the ode stepper takes a bigger control.
Let’s alter the maximum step size and compare to what we had above.
ctrl2 <- ctrl
ctrl2$ode_step_size_max <- 0.1
# run
out2 <- run_scm(p, ctrl = ctrl2)
# compare
length(out$ode_times)
length(out2$ode_times)
# maximum
f <- function(ode_times, title = "
") {
diff(ode_times) %>% hist(main = title, xlim = c(0,2), breaks =seq(0, 2, length.out=20), freq = FALSE)}
f(out$ode_times, "original")
f(out2$ode_times, "max 0.1")Now let’s alter the minimum step size and compare to what we had above.
ctrl3 <- ctrl
ctrl3$ode_step_size_min <- 0.01
ctrl3$ode_step_size_initial <- 0.01
# run
out3 <- run_scm(p, ctrl = ctrl3)
# compare
length(out$ode_times)
length(out3$ode_times)
f(out$ode_times, "original")
f(out3$ode_times, "min 0.01")We set exact step sizes by specifying max and min step sizes. But this could result in other times, depending on the cohort schedule.
ctrl <- Control()
ctrl$ode_step_size_min <- 0.1
ctrl$ode_step_size_max <- 0.1
ctrl$ode_step_size_initial <- 0.1We can also specify some exact steps, by passing in a vector of ode_times.
pc <- p
pc$node_schedule_times[[1]] <- c(0, 0.1, 0.4, 0.8, 1.2)
pc$ode_times <- seq(0, p$max_patch_lifetime, by=0.1)
result2 <- run_scm(pc, ctrl = ctrl, use_ode_times = TRUE)
result2$ode_times [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 4.6 4.7
[50] 4.8 4.9 5.0The options above still drive the full adaptive RK45 machinery at each step. For operational-realism comparisons with daily-step DGVMs, plant can instead use plain forward (explicit) Euler on a uniform grid — one derivative evaluation per step, no error estimate — by setting fixed_time_step (units: years) on the Control:
ctrl <- Control(fixed_time_step = 1 / 365) # daily step; 0 = adaptive (default)
run_scm(p, ctrl = ctrl)This is not the way to get an accurate or cheap solution: forward Euler is first-order, and matching the adaptive solver’s accuracy needs far more evaluations. The next part works this trade-off through with a work–precision diagram. Note that fixed_time_step is incompatible with the mutant-fitness replay path (save_RK45_cache / use_ode_times) and errors clearly if combined.
Would it help to have option to pass in ode_times, use if supplied?
By default, plant integrates the patch with the adaptive, error-controlled Cash–Karp Runge–Kutta (RK45) stepper described above. The solver chooses each step size to keep the local error within a tolerance: tiny steps early on, when newly introduced cohorts are diverging rapidly, and much coarser steps later, once the stand settles.
Many widely used dynamic global vegetation models (DGVMs) instead march the state forward on a fixed grid — most commonly a daily time step — using plain forward (explicit) Euler. plant can be run this way too, via the fixed_time_step control (option 3 above). This part shows how to enable it and uses a work–precision diagram to make the central trade-off concrete: at matched accuracy, adaptive stepping does far less work.
library(plant)
library(dplyr)
library(ggplot2)The integration mode is selected by a single control field, fixed_time_step (units: years):
fixed_time_step = 0 (the default) → adaptive RK45.fixed_time_step > 0 → forward Euler on a uniform grid of that spacing.A daily step is therefore fixed_time_step = 1 / 365.
p0 <- scm_base_parameters("FF16")
p0$max_patch_lifetime <- 60
p <- expand_parameters(trait_matrix(0.0825, "lma"), p0, birth_rate_list = 1.0)
env <- Environment("FF16")
ctrl_adaptive <- Control() # fixed_time_step = 0
ctrl_daily <- Control(fixed_time_step = 1 / 365)
out_adaptive <- run_scm(p, env, ctrl_adaptive)
out_daily <- run_scm(p, env, ctrl_daily)Forward Euler does one derivative evaluation per step (y <- y + h * dy/dt), whereas each RK45 step evaluates the derivatives six times. The two runs report slightly different fitness, and use a very different number of steps:
tibble::tibble(
mode = c("adaptive RK45", "daily Euler"),
R0 = c(out_adaptive$net_reproduction_ratios,
out_daily$net_reproduction_ratios),
n_steps = c(length(out_adaptive$ode_times),
length(out_daily$ode_times))
)# A tibble: 2 × 3
mode R0 n_steps
<chr> <dbl> <int>
1 adaptive RK45 4.48 166
2 daily Euler 4.49 21959Caution — don’t read this as “Euler is faster”. A single daily step is cheaper than a single adaptive step (1 derivative evaluation versus 6), so a naive timing race makes Euler look quick. But that comparison holds the step count fixed, not the accuracy. Daily Euler here is visibly less accurate than the adaptive run. The honest question is: how much work does each method need to reach the same accuracy? That is what the work–precision diagram below answers — and there adaptive wins comfortably.
We sweep the adaptive solver over a range of error tolerances, and forward Euler over a range of step sizes, and for each run record:
6 × n_steps for RK45, 1 × n_steps for Euler (this is the dominant cost of the inner loop);## Tight-tolerance adaptive reference.
ctrl_ref <- control_accurate()
ctrl_ref$ode_tol_rel <- 1e-8
ctrl_ref$ode_tol_abs <- 1e-8
R0_ref <- run_scm(p, env, ctrl_ref)$net_reproduction_ratios
measure <- function(ctrl, mode) {
stages <- if (mode == "adaptive") 6 else 1 # derivative evals per step
scm <- run_scm(p, env, ctrl)
n <- length(scm$ode_times)
tibble::tibble(
mode = mode,
work = n * stages,
error = abs(scm$net_reproduction_ratios - R0_ref) / R0_ref
)
}
adaptive_sweep <- purrr::map_dfr(c(1e-2, 1e-3, 1e-4, 1e-5, 1e-6), function(tol) {
ctrl <- Control()
ctrl$ode_tol_rel <- tol
ctrl$ode_tol_abs <- tol
measure(ctrl, "adaptive")
})
euler_sweep <- purrr::map_dfr(c(2, 1, 0.5, 0.25, 0.125, 0.0625), function(dt) {
measure(Control(fixed_time_step = dt), "euler")
})
work_precision <- bind_rows(adaptive_sweep, euler_sweep)ggplot(work_precision, aes(work, error, colour = mode)) +
geom_line() +
geom_point(size = 2) +
scale_x_log10() +
scale_y_log10() +
labs(x = "Derivative evaluations (work)",
y = expression("Relative error in " * R[0]),
colour = NULL) +
theme_bw() +
theme(legend.position = "top")
The two curves tell a consistent story:
To get daily-Euler accuracy down to the adaptive frontier you would need a sub-daily step everywhere, including the long, smooth tail of the simulation where the adaptive solver is happily taking steps of years. That is wasted work, and it is exactly the inefficiency adaptive stepping removes.
The practical conclusion: fixed-step Euler is valuable for operational-realism comparisons — running plant on the same footing as a daily-step DGVM — but it is not the way to get an accurate or cheap solution. For that, keep the default adaptive solver.
Fixed-step integration removes the solver’s safety net: forward Euler will march across anything, including a discontinuous right-hand side that the adaptive solver rightly refuses. When that happens the integrator still returns numbers — but they are an artefact of where the step grid happens to fall, not an approximation to a solution that gets better as you refine.
plant’s experimental PPA hard-step shading variant (ppa_layer_smoothing = 0, the literal field discretisation) is exactly such a case: its light profile is discontinuous in height, so the adaptive solver fails with “Cannot achieve the desired accuracy”. Under a fixed step it runs. The test of whether that result means anything is convergence: does it settle as the step shrinks? We compare it with the smoothed variant (ppa_layer_smoothing = 0.3), which is a well-posed, continuous problem.
seed_output <- function(steps_per_year, smoothing) {
ctrl <- Control(fixed_time_step = 1 / steps_per_year)
ctrl$shading_model <- "ppa"
ctrl$ppa_layer_smoothing <- smoothing
run_scm(p, env, ctrl)$offspring_production
}
convergence <- purrr::map_dfr(c(91, 182, 365, 730, 1460), function(spy) {
tibble::tibble(
steps_per_year = spy,
`hard step (s = 0)` = seed_output(spy, 0),
`smoothed (s = 0.3)` = seed_output(spy, 0.3))
}) |>
tidyr::pivot_longer(-steps_per_year, names_to = "variant", values_to = "seed")ggplot(convergence, aes(steps_per_year, seed, colour = variant)) +
geom_line() + geom_point(size = 2) +
facet_wrap(~variant, scales = "free_y") +
scale_x_log10() +
labs(x = "Steps per year (1 / fixed_time_step)", y = "Resident seed output",
colour = NULL) +
theme_bw() + theme(legend.position = "none")
The smoothed variant marches monotonically toward a limit — refining the step buys accuracy, as it should. The hard step does not converge to anything: each step size gives a different answer, so no step size gives the answer. This is why such variants must not be read as fitness or as a solution — the fixed-step solver makes them executable, not correct.
fixed_time_step. The mutant pathway replays a resident run by reusing the six cached RK sub-step environments per step; forward Euler has a single stage and no such cache. Combining fixed_time_step > 0 with save_RK45_cache (or with a pinned use_ode_times replay) raises a clear error rather than returning a wrong fitness.