The characteristic method
implementation
This page builds on The size-structured PDE. Read it first if the notation here is unfamiliar.
The idea first
Imagine trying to track how a population of differently sized plants changes through time. One option is to place the population on a fixed size grid and continually update the density at every grid point. This is difficult because growth constantly moves plants from one size class to another, coupling size and time in a partial differential equation (PDE).
The method of characteristics takes a different route. Instead of watching the grid, you follow individual plants—or cohorts of plants that started life together—as they grow. Each one traces a path through time as its size changes, and along that path you can also track how many individuals remain and how their local density evolves. Because you move with the plant rather than watching it pass through fixed points, the partial differential equation reduces to an ordinary differential equation for each trajectory. These equations are much easier to solve one step at a time with a standard numerical stepper.
plant uses this characteristic-method approach to solve the size-density distribution (O. Angulo & López-Marcos, 2004; O. Angulo, López-Marcos, & López-Marcos, 2014, 2016). Earlier versions used the related Escalator Boxcar Train technique (De Roos, 1988; De Roos, Diekmann, & Metz, 1992; De Roos, Taljapurkar, & Caswell, 1997).
When simulating an individual plant or the development of a patch, we need to solve for individual plants’ size, survival, and seed output. When solving for the size-density distribution in a large patch, we also need to estimate the average abundance of individuals. Each problem is formulated as an initial-value problem (IVP)—an ODE with specified starting values—which can be solved using an ODE stepper.
All these IVPs must be advanced through time. plant uses an embedded Runge-Kutta Cash-Karp 4-5 algorithm with an adaptive time step (Cash & Karp, 1990), based on code ported from the GNU Scientific Library (Galassi, 2009). The controls ode_tol_rel and ode_tol_abs set the relative and absolute error tolerances. See ODE stepping for a full explanation.
Written as a set of ODEs, the equations to be solved are as follows.
ODEs along each characteristic
Size
We obtain an individual’s size from the size equation (defined on the size-structured PDE page), which is solved via the IVP
\[ \frac{dy}{dt} = g(x, y, t) , \]
\[ y(0) = H_0(x). \]
Survival
The probability that an individual survives from patch age \(a_0\) to patch age \(a\) is given by the individual-survival equation (see the size-structured PDE page), which is solved via the IVP
\[ \frac{dy}{dt} = d(x, H_i(t) , E_t), \]
\[ y(0) = - \ln\left(S_{\rm G} (x, H_0, E_{a0})\right) . \]
The survival probability is then
\[ S_{\rm I} (x, a_0, a) = \exp\left( - y(a) \right). \]
Seed production
We obtain lifetime seed production from the lifetime seed-production equation, which is solved via the IVP
\[ \frac{dy}{dt} = S_{\rm D} \, f(x, H_i(t), E_t) \, S_{\rm I} (x, a_0,t) \, S_{\rm P} (a_0, t), \]
\[ y(0) = 0, \]
where \(S_{\rm I}\) is calculated as described above and \(S_{\rm P}\) is calculated as in the patch-survival equation.
Size-density of individuals
By integrating along the characteristics of the governing PDE (defined on the size-structured PDE page), the size-density of individuals born with height \(H_0\) and traits \(x\) at patch age \(a_{0}\) is given by (De Roos et al., 1997; O. Angulo & López-Marcos, 2004)
\[ N(H | x, a) = N(H_0 | x, a_0) \exp \left( - \int _{a_0}^{a} \left[\frac{\partial g(x, H(x, a_0, a^\prime), E_{a^\prime})}{\partial H} + d(x, H(x, a_0, a^\prime), E_{a^\prime})\right] {\rm d} a^\prime \right). \tag{1}\]
Equation 1 says that size-density \(N\) at patch age \(a\) starts from the density at age \(a_0\) and is then adjusted for growth and mortality. Mortality lowers density in the usual way. Growth changes density by spreading or compressing the characteristic trajectories. If growth slows with size (\(\partial g / \partial H < 0\)), nearby trajectories converge, creating a compression effect that raises density relative to what mortality alone would produce. Total density increases only if that compression is stronger than the mortality loss. If growth accelerates with size (\(\partial g / \partial H > 0\)), the trajectories spread apart, so growth and mortality both reduce density.
Let \(\left[H_0, H_{ + } \right)\) denote the range of heights attainable by any individual. To solve metacommunity dynamics, our algorithm subdivides this interval into a series of nodes with heights \(H_0 < H_1 < \ldots < H_k\) at the initial points of the characteristic curves. These nodes are then transported along the characteristics of the PDE. Their placement is controlled indirectly by the schedule of patch ages at which new nodes are introduced into the metacommunity. We then track the demography of each node.
The integral in Equation 1 is solved via the IVP
\[ \frac{dy}{dt} = \frac{\partial g(x, H_i(t), E_t)}{\partial H} + d(x, H_i(t), E_t), \tag{2}\]
\[ y(0) = - \ln\left(N(H_0 | x, a_0) \, S_{\rm G} (x, H_0, E_{{\rm a}0}) \right), \tag{3}\]
from which we obtain the size-density
\[ N(H_0 | x, a_0) = \exp( - y(a)). \tag{4}\]
Controls on approximation error
The equations above cannot usually be solved exactly, so several parts of the calculation introduce numerical approximation:
- estimating light at a given height requires integration over the patch’s size-density distribution;
- calculating assimilation requires integration of leaf-level photosynthesis over a plant’s crown;
- simulating patch dynamics requires choosing when to introduce nodes and then stepping their size, survival, and fecundity forward through time; and
- finding initial height from seed mass, or equilibrium birth rates for a metacommunity, requires numerical root finding.
Each method is accurate only to a chosen tolerance. The following sections explain the main sources of error and the plant controls that govern them. A worked example of changing numerical controls is also available in Appendix S3 of Falster et al. (2016).
Initial plant heights
When a seed germinates, it produces a seedling of a given height. We assume that seedling height varies with seed mass. Because there is no analytical solution relating seedling height to seed mass – at least when using the default FF16 physiological model – we must solve for this height numerically. The calculation is performed by the function height_seed within the physiological model, using the Boost library’s one-dimensional bisect routine (Schäling, 2014; Eddelbuettel, Emerson, & Kane, 2015). The accuracy of the solution is controlled by the parameter offspring_production_tol.
Approximation of size-density distribution via the characteristic method
Two sources contribute to errors in the approximation of the size-density distribution: (i) coarse stepping of nodes through time and (ii) poor spacing of nodes across the attainable size range.
As described above, the stepping of the ODE solver is controlled by two control parameters for relative and absolute accuracy, ode_tol_rel and ode_tol_abs (see ODE stepping).
A second source of error is the derivative in Equation 1, which plant estimates using standard finite differences (Abramowitz & Stegun, 2012). Setting node_gradient_richardson to TRUE adds Richardson extrapolation (Stoer & Bulirsch, 2002); node_gradient_richardson_depth controls how far that refinement continues. The finite-difference scale is controlled by node_gradient_eps.
The primary factor controlling node spacing is the schedule of node introduction times. Because the system of equations to be integrated is deterministic, the schedule of node introduction times determines the spacing of nodes throughout the entire development of a patch. Poor node spacing introduces error because various emergent properties – such as total leaf area, biomass, or seed output – are estimated by integrating over the size-density distribution. The accuracy of these integrations declines with inappropriate spacing of nodes. Our algorithm therefore aims to build an appropriately refined schedule, which allows the required integrations to achieve the desired accuracy at every time point. For computational efficiency, however, we want to use as few nodes as possible. The general idea of adaptively refining node spacing times was first applied by Falster, Brännström, Dieckmann, & Westoby (2011), and was described further by Falster, Brännström, Westoby, & Dieckmann (2015).
For a worked example illustrating schedule refinement (run_scm(refine_schedule = TRUE), backed by SCM::refine_schedule()), see node spacing.
Calculation of light environment and influence on assimilation
To solve the system of ODEs, we must calculate the amount of shading on each node caused by all the other plants in the patch.
Calculating the canopy openness \(E_a(z)\) at a given height \(z\) in a patch of age \(a\) requires that we integrate over the size-density distribution (the light equation on the size-structured PDE page). This integration is performed using the trapezium rule, within the function area_leaf_above in species.h. Node spacing is the main factor controlling the accuracy of this integration. The node introduction times determining the spacing of nodes are adaptively refined as described above. This implies that also the trapezium integration within the area_leaf_above function is adaptively refined via SCM::refine_schedule() (exposed in R as run_scm(refine_schedule = TRUE)). See node spacing for detail.
The cost of calculating \(E_a(z)\) increases linearly with the number of nodes in the metacommunity. Since the same calculation must be repeated for every node, the overall computational cost of a step increases as \(O(k^2 )\), where \(k\) is the total number of nodes across all species. This disproportionate increase in computational cost with the number of nodes is highly undesirable.
We reduce the computational cost from \(O(k^2)\) to \(O(k)\) by approximating \(E_a(z)\) with a spline. Canopy openness increases monotonically with height, so a piecewise-continuous spline can represent it using a limited number of points. Once fitted, the spline provides cheap additional evaluations of the competitive environment. A new spline is constructed at every time step.
The accuracy of the spline interpolation depends on the number of points used in its construction and on their placement along the size axis. We select the number and locations of points via an adaptive algorithm. Starting with an initial set of light_availability_spline_nbase points (17 by default), we assess how much each point contributes to the accuracy of the spline fit at the location of each node, first via exact calculation, and second by linearly interpolating from adjacent nodes. The absolute difference in these values is compared to the spline tolerance light_availability_spline_tol (default 1e-4). Note that this is a property of the light spline configured on the environment (e.g. FF16_Environment), not a field of the Control object. If the error is greater than this tolerance, the interval is bisected and the process repeated, up to light_availability_spline_max_depth (16 by default) levels (see adaptive_interpolator.h for details).
Integration over light environment
Plants have leaf area distributed over a range of heights. Estimating a plant’s assimilation at each time step therefore requires integrating leaf-level rates over the plant. We perform this integration using Gauss-Kronrod quadrature; a rule controls its level of detail.
Solving for offspring production
For a single species, solving for \(Y_x\) is a straightforward one-dimensional root-finding problem, which we solve with a simple bisection algorithm.
Solving for offspring production in metacommunities with multiple species is significantly harder because there is no generally applicable method for multi-dimensional root finding. plant therefore implements several approaches, available through the package {regnans}.