reference

Methods: Discretisation, Sampling and Convergence

A methods-section treatment with numbered equations: the master balance, the representative-slice reduction, quadrature discretisation of the double integral, the Monte-Carlo estimator with 1/√N error, and the relative-variance propagation law, with a discussion of convergence.

Last updated: 2026-07-26

This chapter sets out, in the style of a methods section, the computational procedures by which a number is extracted from the Schubart Master Formula. It proceeds from the continuous master balance to its practical evaluation in five numbered steps: the exact integral, its reduction to a representative slice, its discretisation as a quadrature sum, its estimation by Monte Carlo, and the propagation of input uncertainty. Each key relation is numbered so that the surrounding text and later chapters can refer to it unambiguously as “Eq. (n)”. The full analytic derivation that motivates these steps is given in the full derivation; here the emphasis is on how the integral is actually computed and on the errors each method carries.

Equation (1) — the expanded balance

The starting point is the expanded physical representation of the Schubart Master Formula — the accounting identity for the instantaneous converted power written out in its spectral–transport form,

P(t) = \iint \Phi(E,\theta,t)\,\sigma(E)\,\Delta p\,v_{\text{ph}}\,\eta_{\text{conv}}\;dE\,d\theta \tag{1}

where the integration runs over the incident energy $E$ and the angular coordinate $\theta$ , and the five factors are the flux density, the scattering cross-section, the effective momentum transfer, the phonon velocity, and the conversion efficiency. Eq. (1) is a double integral over a two-dimensional domain; every subsequent method in this chapter is a way of evaluating it. The dimensional bookkeeping that shows the right-hand side reduces to watts is carried out in the dimensional analysis, and is assumed here.

Equation (2) — representative slice

When the integrand is slowly varying over a well-chosen window $\Delta E \times \Delta\theta$ , or when a first estimate is wanted, Eq. (1) reduces to a single representative slice by replacing the integral with the integrand evaluated at representative values times the window size,

P \approx \Phi\,\sigma\,\Delta p\,v_{\text{ph}}\,\eta_{\text{conv}}\;\Delta E\,\Delta\theta . \tag{2}

Eq. (2) is the mean-value form of Eq. (1): it is exact if the factors are constant across the window and a controlled approximation otherwise. Its validity is exactly the representative-averaging condition of Assumption C in the assumptions and limitations — the window must be wide enough to be representative yet narrow enough that the integrand does not swing across it. Eq. (2) is invaluable for order-of-magnitude estimates and for checking that a fuller calculation has not gone astray, but it must not be trusted where the integrand is sharply peaked, for instance across a cross-section resonance.

Equation (3) — quadrature discretisation

To evaluate Eq. (1) without assuming a single representative slice, the domain is partitioned into a grid of cells indexed by $j$ (in energy) and $k$ (in angle), and the integral is approximated by a Riemann or quadrature sum over those cells,

P \approx \sum_{j,k} \Phi_{jk}\,\sigma_j\,\Delta p\,v_{\text{ph}}\,\eta_{\text{conv}}\;\Delta E_j\,\Delta\theta_k , \tag{3}

where $\Phi_{jk}$ and $\sigma_j$ are the factor values sampled at the cell centres and $\Delta E_j$ , $\Delta\theta_k$ are the cell widths. Eq. (3) is simply Eq. (2) summed over many slices; as the cells shrink it converges to Eq. (1). A uniform midpoint rule as written has a discretisation error per cell of order $(\Delta E)^3$ and $(\Delta\theta)^3$ in the second derivative of the integrand, so the global error of the midpoint or trapezoidal rule scales as $O(h^2)$ in the cell size $h$ ; Simpson’s rule improves this to $O(h^4)$ on smooth integrands. Non-uniform grids that place cells more densely where the integrand varies fastest — around a spectral peak, say — recover accuracy far more cheaply than uniform refinement. This grid-based evaluation is the backbone of the deterministic numerical simulation.

Equation (4) — the Monte Carlo estimator

Grid quadrature becomes expensive when the integrand is high-dimensional, sharply peaked, or awkwardly shaped, because the number of cells grows as a power of the dimension. The Monte Carlo alternative rewrites the integral as an expectation and estimates it by random sampling. Drawing $N$ points $x_n = (E_n,\theta_n)$ from the domain and writing $f(x)$ for the full integrand of Eq. (1) times the domain volume, the estimator is the sample mean

\hat P = \frac{1}{N}\sum_{n=1}^{N} f(x_n), \qquad \mathrm{SE}\big[\hat P\big] \propto \frac{1}{\sqrt{N}} . \tag{4}

Eq. (4) is unbiased — its expectation equals the true integral — and its standard error falls as $1/\sqrt{N}$ , independent of the dimension of the domain. That dimension-independence is the reason Monte Carlo overtakes tensor-grid quadrature for complicated integrands, and it is also why the method extends naturally to uncertainty propagation: by drawing the input factors themselves from their uncertainty distributions in an outer loop, each Monte Carlo pass yields one realisation of $P$ , and the ensemble reproduces the output distribution. The full nested procedure is set out in the Monte Carlo validation.

Equation (5) — uncertainty propagation

The final method concerns not the integral’s value but its uncertainty. Because Eq. (1) is a product of factors, the relative variances of the factors add. For $P = \prod_i x_i$ with statistically independent factors, a first-order propagation gives

\left(\frac{\sigma_P}{P}\right)^{2} = \sum_i \left(\frac{\sigma_{x_i}}{x_i}\right)^{2} . \tag{5}

Eq. (5) follows from the log-differential $\mathrm{d}\ln P = \sum_i \mathrm{d}\ln x_i$ and is derived in full, together with the correlated-input covariance terms, in the uncertainty analysis. Its role among the methods is to tell us, before any expensive sampling, which input dominates the output error: because the terms add in quadrature, the largest relative uncertainty controls the result, and refinement effort should be directed there. In the small-uncertainty limit the Monte Carlo spread from Eq. (4) reproduces Eq. (5), so the two serve as mutual consistency checks.

Convergence

The three evaluation methods — Eqs. (2), (3) and (4) — carry qualitatively different convergence behaviour, and choosing between them is a matter of matching the method to the integrand.

The representative slice, Eq. (2), does not converge in any refinement sense; it is a fixed approximation whose error is set entirely by how representative the chosen window is. It is a first estimate, not a convergent scheme.

Quadrature, Eq. (3), converges as a power of the cell size: $O(h^2)$ for the midpoint and trapezoidal rules, $O(h^4)$ for Simpson’s rule, on smooth integrands. Halving the cell size therefore cuts the trapezoidal error fourfold. This rapid, high-order convergence makes quadrature the method of choice for smooth, low-dimensional integrands — precisely the regime of the two-dimensional Eq. (1) when the factors are well-behaved. Its weakness is that the cell count grows as a power of the dimension, and that sharp peaks force very fine grids.

Monte Carlo, Eq. (4), converges as $1/\sqrt{N}$ : reducing the error tenfold demands a hundredfold increase in samples. This is slow compared with quadrature in low dimensions, but the exponent $-\tfrac12$ is independent of the dimension, so Monte Carlo wins for high-dimensional or sharply structured integrands where quadrature’s grid explodes. Its slow rate is routinely improved by variance reduction: importance sampling, which concentrates draws where the integrand is largest and reweights by $f/g$ ; stratified sampling, which guarantees coverage of important sub-regions; and control variates, which subtract a correlated quantity of known integral. Each of these lowers the effective variance $\mathrm{Var}[f]$ in the numerator of the standard error without changing the $1/\sqrt{N}$ scaling, so that a modest $N$ achieves what would otherwise require orders of magnitude more uniform samples.

In practice a robust evaluation of Eq. (1) uses all three: Eq. (2) for a sanity-check order of magnitude, Eq. (3) for the deterministic value when the integrand is smooth, and Eq. (4) with variance reduction for the uncertainty propagation and for any regime where the grid becomes impractical. The propagation law Eq. (5) then supplies the analytic error budget that the sampling is expected to reproduce.