validation

Numerical Simulation

Methodology for evaluating the P(t) master formula numerically: discretising the integral over energy and angle, supplying inputs, and checking convergence.

Last updated: 2026-07-20

Because the Schubart Master Formula is a double integral over energy and angle, evaluating it for any concrete case requires numerical methods. This chapter describes the general methodology for doing so — how to discretise the integral, how to supply the inputs, and how to confirm that a computed result has converged. Throughout, the emphasis is on method rather than on specific numbers; any figures mentioned are illustrative placeholders framed as such, not measured results.

Discretising the integral

The starting point is the expanded physical representation,

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

The most direct numerical approach replaces the continuous integrals with sums over a grid. The energy axis is partitioned into $N_E$ intervals and the angular axis into $N_\theta$ intervals, and the integrand is evaluated at representative points in each cell. Writing the integrand as $f(E,\theta,t)$ , a simple quadrature reads

P(t) \approx \sum_{i=1}^{N_E}\sum_{j=1}^{N_\theta} w_{ij}\,f(E_i,\theta_j,t)\,\Delta E_i\,\Delta\theta_j ,

where $\Delta E_i$ and $\Delta\theta_j$ are the widths of the cells and $w_{ij}$ are quadrature weights. For the midpoint rule the weights are unity and the sample points sit at cell centres. More accurate schemes — trapezoidal, Simpson, or Gauss–Legendre quadrature — replace the midpoint evaluation with weighted sums that reduce discretisation error for a given number of points. The choice of scheme is a standard trade-off between accuracy and evaluation cost, and nothing about the framework favours one particular quadrature over another.

Discretisation strategy and convergence order

Because the domain is a product of the energy and angle axes, the two-dimensional rule is most simply built as a tensor product of one-dimensional rules, one for $E$ and one for $\theta$ . This factorisation matters for both cost and accuracy. On the cost side, a tensor grid uses $N_E \times N_\theta$ integrand evaluations, so refining both axes together is quadratically expensive; when one axis carries all the structure — a sharp resonance in $\sigma(E)$ along the energy axis, say — it is often far cheaper to refine that axis alone. On the accuracy side, each one-dimensional rule carries its own convergence order, the exponent $p$ in the asymptotic error scaling

\left| P_{\text{exact}} - P_{N} \right| \;\sim\; C\,N^{-p} ,

where $N$ is the number of points along the axis and $C$ depends on the smoothness of the integrand. For the trapezoidal rule $p = 2$ , so halving the step size cuts the error roughly fourfold; for Simpson’s rule $p = 4$ on a smooth integrand; and for Gauss–Legendre quadrature the convergence is faster still, effectively exponential for an analytic integrand with no singularities in the domain. These orders are exactly the classical results tabulated in standard references such as Press et al., Numerical Recipes, and they carry over to the double integral through the tensor construction. The practical consequence is that the smoothness of the integrand, not the nominal order of the scheme, sets the achievable accuracy: a fourth-order rule applied across a near-discontinuity in $\sigma(E)$ degrades to first-order behaviour locally, and the global error is then dominated by that one badly-resolved region.

Supplying the inputs

A simulation is only as meaningful as its inputs. Four kinds of input must be provided. First, the flux model $\Phi(E,\theta,t)$ : a description of the ambient spectrum and its angular distribution, either measured or drawn from an assumed distribution. Second, the cross-section $\sigma(E)$ : the energy-dependent coupling, typically supplied as a tabulated or parameterised function; its analytic behaviour is discussed in sigma analysis. Third, the momentum-transfer, transport, and efficiency parameters $\Delta p$ , $v_{\text{ph}}$ , and $\eta_{\text{conv}}$ , which may be constants or may themselves depend on $E$ and $\theta$ . The definitions and units of all of these are fixed in variable definitions. Because several of these inputs carry uncertainty, a responsible simulation propagates that uncertainty rather than reporting a single deterministic number — a point taken up in Monte Carlo validation.

Checking convergence

The central quality-control question for a deterministic quadrature is whether the result has converged with respect to grid resolution. The standard procedure is to compute $P(t)$ on a sequence of successively finer grids — doubling $N_E$ and $N_\theta$ at each step — and to observe how the estimate changes. When successive refinements alter the result by less than a chosen tolerance, the grid is deemed sufficient. A more quantitative version of this test uses the assumed convergence order directly: if the error scales as $N^{-p}$ , then computing the integral on three grids with spacings $h$ , $h/2$ , and $h/4$ and comparing the successive differences reveals both the observed order and, by Richardson extrapolation, an improved estimate. Writing $P_h$ , $P_{h/2}$ , $P_{h/4}$ for the three results, the ratio

\frac{P_{h} - P_{h/2}}{P_{h/2} - P_{h/4}} \;\approx\; 2^{p}

should recover the expected order; a value far from $2^p$ is a warning that the integrand is not as smooth as assumed and that the nominal error estimate cannot be trusted. For a smooth integrand a well-chosen quadrature converges rapidly; sharp features in $\sigma(E)$ or in the flux, such as a narrow resonance, demand finer sampling in the region of the feature and are a common cause of slow convergence. Adaptive schemes, which subdivide cells wherever a local error estimate exceeds a threshold and thereby place more points where the integrand varies quickly, are the natural remedy.

As an illustrative example only, one might find that increasing the energy grid from a coarse setting to a moderately fine one shifts a placeholder output estimate by a few percent, after which further refinement changes it by well under a percent — the signature of a converged calculation. Such figures are purely notional and would need to be established anew for any real case; they are given here only to show the shape of a convergence study, not to assert a result.

Error estimation and budgeting

A defensible numerical result comes with an explicit error budget that separates the distinct contributions. The discretisation error is the $C\,N^{-p}$ term above, estimated in practice by the successive-refinement difference $|P_{h/2} - P_{h/4}|$ , which for a converging sequence bounds the remaining error to within a small factor. The truncation error arises from cutting the energy domain at a finite upper limit; it is bounded by the integral of the neglected tail, and the domain should be extended until that bound falls below the discretisation error. The input error, propagated from uncertainty in $\Phi$ , $\sigma(E)$ , $\Delta p$ , $v_{\text{ph}}$ , and $\eta_{\text{conv}}$ , is usually the largest of the three and is handled by the stochastic methods of Monte Carlo validation rather than by grid refinement. A common and instructive mistake is to refine the grid to machine precision while ignoring an input uncertainty of tens of percent — polishing a number whose true error is dominated by a factor the quadrature cannot touch. The discipline is to report all three contributions and to let the largest one govern how much effort each further refinement deserves.

Sources of numerical error

Beyond discretisation error, a careful practitioner watches for a few recurring pitfalls. Truncating the energy integral too early discards flux in the tails; the domain should extend until the integrand is genuinely negligible, as noted in mathematical architecture. Using an angular grid that ignores the correct measure — for instance omitting a $\sin\theta$ factor where a solid-angle integration requires it — introduces a systematic bias that no amount of refinement will remove, because it is an error in the integrand, not in its sampling. And interpolating a coarsely tabulated $\sigma(E)$ across a sharp feature can under- or over-estimate the overlap with the flux. Each of these is avoidable with standard numerical hygiene, and each should be documented alongside any reported result.

From deterministic to statistical evaluation

Deterministic quadrature answers the question “given these exact inputs, what is $P(t)$ ?” It does not, by itself, answer “given the uncertainty in these inputs, how confident should we be in $P(t)$ ?” That second question calls for stochastic methods, in which the inputs are sampled from their distributions and the integral is estimated statistically. Those methods, and the confidence intervals they produce, are the subject of Monte Carlo validation and sigma analysis. A converged deterministic quadrature also serves a second role there: run on a single fixed set of inputs, it is the reference against which a Monte Carlo estimator is cross-checked, so the two methods are complementary rather than alternative.