16  Midpoint Rule vs. Gauss-Legendre Quadrature

The release of v1.0.1 of the {forestIPM} R package introduces a new integration approach: the Gauss-Legendre Quadrature. In this chapter, we compare this new approach agains the original midpoint rule.

Method Description Grid Meshpoints
Midpoint rule Uniform bins of width \(h = 1\) mm \([L,\, L_{\max}]\) \((L_{max} - 127)/h\)
Gauss-Legendre (GL) Optimal quadrature nodes and weights \(n_{\text{gl}}\) nodes on \([L,\, L_{\max}]\) \(n_{gl}\) nodes

The midpoint rule was used in the original cluster simulations (covariates_perturbation) and serves as the historical reference throughout this report. GL quadrature is the new implementation: it places nodes at positions that exactly integrate polynomials of degree \(\leq 2n_{\text{gl}} - 1\), achieving faster convergence for smooth integrands.

In practice, we run this simulation to addresses two questions:

  1. Accuracy — does GL recover the same \(\lambda\) as the midpoint rule, and at what n_gl?
  2. Speed — is GL faster than the midpoint rule regardless of species size range?

Simulation design

Results were pre-computed by compare_versions_gl.R:

  • 15 stratified environmental scenarios drawn from 216 027 combinations of species × plot × climate.
  • 100 Bayesian posterior parameter draws per scenario.
  • n_gl tested: 20, 50, 100, 200, 500.
  • Reference: midpoint rule (\(h = 1\) mm) from final_output.RDS.
  • Metric: dominant eigenvalue \(\lambda\) of the IPM kernel K.
  • Tolerance: \(|\lambda_{\text{GL}} - \lambda_{\text{midpoint}}| < 10^{-4}\).

GL accuracy vs. midpoint rule

Absolute deviation from midpoint

Absolute deviation of GL(n) lambda from the midpoint reference, by n_gl. Boxes show median and IQR across all 15 scenarios × 100 draws. Orange dashed line = GL tolerance (1e-4).

Relative deviation from midpoint

Relative deviation of GL(n) lambda from the midpoint reference. Red dashed = 0.1%, orange dotted = 1% tolerance.

Pass/fail at absolute tolerance 1e-4

GL lambda convergence vs. midpoint rule (absolute deviation)
n_gl N draws Pass (< 1e-4) % Pass Max |dev| Median |dev|
20 1500 0 0.0% 3.28e+00 1.28e+00
50 1500 13 0.9% 6.86e-01 6.60e-03
100 1500 529 35.3% 8.35e-03 2.16e-04
200 1500 840 56.0% 7.19e-03 7.69e-05
500 1500 1073 71.5% 6.58e-03 3.27e-05

Per-species deviation heatmap

Median absolute deviation of lambda from the midpoint reference, per species and n_gl.

Per-species convergence profiles

Each panel shows how \(\lambda\) deviation from the midpoint rule evolves as n_gl increases. The ribbon spans the IQR across the 15 scenarios × 100 parameter draws.

Per-species absolute deviation of GL lambda from midpoint. Line = median; ribbon = IQR. Orange dashed = tolerance 1e-4.

GL internal convergence

This section examines convergence within the GL family — how smaller n_gl values compare to GL(500), independent of the midpoint reference.

Absolute deviation of GL(n) from GL(500). Convergence within the GL family is much faster than convergence toward the midpoint rule.
GL(n) vs. GL(500): absolute deviation in lambda
n_gl Median |dev| Max |dev|
20 1.28e+00 3.28e+00
50 4.37e-03 6.86e-01
100 2.41e-04 5.83e-03
200 7.82e-05 2.65e-03

Structural difference: GL vs. midpoint

Even at GL(500) — which is fully converged within the GL family — \(\lambda\) differs systematically from the midpoint rule. The figure below shows this plateau: as n_gl increases, deviation from the midpoint rule does not approach zero.

Mean absolute deviation of GL(n) from midpoint (blue) vs. GL internal deviation from GL(500) (grey). The GL-vs-midpoint curve plateaus, revealing a structural difference between the two integration schemes.

The plateau arises because GL and midpoint are computing fundamentally different discrete approximations of the same continuous integral. GL uses non-uniform nodes and weights optimised for smooth integrands; the midpoint rule uses uniform spacing. At finite resolution, neither converges to the other — they converge to the same continuous-limit value at different rates.

In practice, switching from midpoint to GL will always produce a step-change in \(\lambda\), regardless of n_gl. This is expected and not a numerical error; it reflects the change in integration scheme.

Computational benchmark

Timing by species and method

Median computation time per lambda evaluation (mesh setup + kernel assembly + eigenvalue). Each box shows 30 independent repetitions. Species chosen to bracket the midpoint mesh-size range: PINBAN (fewest bins) and QUERUB (most bins).

Speedup relative to midpoint rule

GL timing relative to midpoint rule (median over 30 repetitions)
Species Midpt bins Midpt median(s) n_gl GL median(s) GL min(s) GL max(s) Speedup
PINBAN 453 0.1750 50 0.001 0.001 0.002 175.0
PINBAN 453 0.1750 100 0.006 0.005 0.009 29.2
PINBAN 453 0.1750 150 0.014 0.012 0.016 12.5
QUERUB 2269 24.3975 50 0.001 0.001 0.002 24397.5
QUERUB 2269 24.3975 100 0.007 0.007 0.010 3485.4
QUERUB 2269 24.3975 150 0.018 0.015 0.025 1355.4

Kernel matrix assembly and eigenvalue computation scale as \(\mathcal{O}(n^2)\) and \(\mathcal{O}(n^3)\) respectively, so the reduction in node count from midpoint to GL compounds quickly for large-domain species.

Recommendation

Based on the accuracy analysis (15 scenarios × 100 parameter draws, midpoint h = 1 as reference):

n_gl Pass ( dev < 1e-4) Pass (rel < 1%)
20 0.0% 0.5% 3.28e+00 Exploration only
50 0.9% 59.5% 6.86e-01 Fast screening
100 35.3% 100.0% 8.35e-03 good balance
200 56.0% 100.0% 7.19e-03 default Production
500 71.5% 100.0% 6.58e-03 Reference / validation

Note on comparing outputs across methods

Lambda values obtained with GL quadrature cannot be directly compared with historical midpoint-rule results without accounting for the structural step-change between schemes (see Structural difference). For continuity with final_output.RDS, use the midpoint-rule lambda as the reference when reporting changes relative to historical runs.

Methods summary

Item Value
GL algorithm Golub-Welsch (Jacobi tridiagonal eigendecomposition)
C++ solver Eigen::SelfAdjointEigenSolver via RcppEigen
Size domain \([127\,\text{mm},\, L_{\max}]\) (species-specific \(L_{\max}\))
Vital rate models von Bertalanffy growth × logistic survival + log-linear recruitment
Metric Dominant eigenvalue \(\lambda\) of the IPM kernel K
Reference Midpoint rule (\(h = 1\) mm) from final_output.RDS
n_gl tested 20, 50, 100, 200, 500
Scenarios 15 stratified rows from 216 027 combinations
Draws per scenario 100 Bayesian posterior samples
Absolute tolerance \(\|\lambda_{\text{GL}} - \lambda_{\text{midpoint}}\| < 10^{-4}\)
Benchmark species PINBAN (smallest domain), QUERUB (largest domain)
Benchmark repetitions 30 per method per species