Implementation
implementation
The biological equations describe what the plant model should do. The implementation determines how a computer solves those equations accurately and efficiently. This section explains the main numerical ideas without assuming that you already know the structure of the C++ code.
How the pieces fit together
A patch simulation can be understood as a sequence of five connected tasks:
- The characteristic method turns the population-level partial differential equation into trajectories for representative cohorts.
- An ODE solver advances the size, survival, and reproduction of those cohorts through time.
- The node-spacing algorithm decides when another cohort is needed to represent the population accurately.
- The canopy model combines the cohorts into a vertical light profile, which then affects their growth.
- Numerical tools such as quadrature, finite differences, and splines make the repeated integrations and derivative calculations practical.
The first four topics each have a dedicated page. The final part of this page introduces the smaller numerical tools used throughout the solver.
Characteristic method
The size-structured population model is written as a partial differential equation (PDE). Rather than solving that PDE on a fixed size grid, plant follows representative cohorts as they grow. Each cohort traces a characteristic through time, turning the original problem into a collection of ordinary differential equations (ODEs).
Read The characteristic method for the equations and their connection to growth, survival, reproduction, and population density.
ODE stepping
Once the model has been expressed as ODEs, a numerical stepper moves their state forward through time. The default adaptive solver takes small steps when the dynamics change rapidly and larger steps when they are smooth. Its error tolerances control the balance between accuracy and computation time.
Read ODE stepping and control for the adaptive RK45 solver, its C++ interface, and the fixed-step alternative.
Adaptive node spacing
Each node represents a cohort of plants born at roughly the same time. Too few nodes give a poor approximation to the population; too many make the simulation unnecessarily expensive. An ordinary run uses the supplied introduction schedule. When greater numerical accuracy is needed, users can opt in to a refinement step that adds nodes only where removing a node would make an important population integral too inaccurate.
Read The node-spacing algorithm for the refinement test and worked diagnostics.
Canopy and shading models
Every cohort both experiences and contributes to a shared light environment. Different canopy models trade biological detail for speed, so the appropriate choice depends on whether the goal is stand dynamics, demographic equilibrium, or eco-evolutionary analysis.
Read Canopy models for the available shading models, their computational cost, and the reason PPA requires a smoothed light profile under adaptive ODE stepping.
Numerical tools inside the solver
The methods above rely on several smaller numerical operations. Most users do not need to call these directly, but knowing their roles makes the larger algorithm easier to follow.
Gauss-Kronrod quadrature
Gauss–Kronrod quadrature estimates the integral of a function from values at a chosen set of points. plant uses it, for example, to integrate leaf-level assimilation over crown depth.
The QK class performs the quadrature. The control variable function_integration_rule selects the number and arrangement of evaluation points from rules defined in qk_rules.cpp, including QK15, QK21, QK31, QK41, and QK51. The default rule uses 21 points. This code was ported from the GNU Scientific Library (GSL) by Rich FitzJohn.
The QAG class adds adaptive refinement around QK. It compares a Gauss rule with its Kronrod extension to estimate approximation error, then adds points where needed. The implementation is adapted from the QAG algorithm in QUADPACK. For general background, see Gaussian quadrature.
Numerical derivatives
The characteristic method needs the derivative of growth rate with respect to plant size when it updates population density. plant estimates this derivative with finite differences and can optionally improve the estimate using Richardson extrapolation. The relevant controls are described on The characteristic method.
Trapezium integration
The trapezium rule approximates the area under a curve by joining neighbouring points with straight lines. plant uses this rule to integrate quantities such as leaf area over the discrete cohort distribution. The internal helper is trapezium(xx, yy); its accuracy therefore depends strongly on node spacing.
Adaptive splines for the light environment
Calculating canopy openness directly at every cohort would make each solver step scale poorly as more nodes are added. Instead, plant calculates openness at an adaptively chosen set of heights, fits a spline, and then evaluates that much cheaper approximation wherever it is needed. The ResourceSpline class manages the adaptive sampling and stores the result in an odelia::interpolator::Interpolator.
Readers working on the C++ implementation can follow the calculation in two directions.
To construct the environment, the patch calls environment.compute_environment(f, height_max()), with f = compute_competition. That calculation passes from patch to species to node to individual to strategy. At the node level, density weights the individual competitive effect. The strategy ultimately evaluates the leaf area above height \(z\) from an expression of the form k_I * area_leaf(height) * Q(z, height). ResourceSpline samples this result across height and builds the light-availability spline.
To use the environment, compute_rates passes it from patch to species to node to individual to strategy. In FF16_Strategy, net_mass_production_dt calls assimilation. A quadrature::QK integrator then evaluates assimilation_leaf(environment.get_environment_at_height(z)) * q(z, height) over crown depth \(z\). Earlier versions used the adaptive QAG layer here, but it could take about twice as long without noticeably changing results, so the current calculation uses the fixed QK rule selected by function_integration_rule.