4. How accurate is the result?

fobis verify says whether a number changed. fobis assess says how far it can be trusted: it runs an ensemble of variants of the case and reads the spread. Three families here:

[assess:area]
families = refinement precision rounding
refine   = n: arg 1 values 8 16 32 64 128 inverse
budget   = simpson: order >= 3.9; trapezoid: order >= 1.9

• refinement: the case at n = 8 ... 128 intervals (arg 1: the program's first argument), inverse because the interval width is h = 1/n. The observed order of convergence comes from successive triples.
• precision: the case rebuilt with its reals promoted (gfortran's -freal-8-real-16): how many digits of the base run survive in a more precise run.
• rounding: the case rebuilt with compiler-level perturbations (optimisation level, FMA contraction, -ffast-math, vectorisation): how many digits every variant shares.

The budget makes the assessment a check: the orders the two rules must reach.

$ fobis assess area
assess: build base
assess: build precision-promoted
assess: build rounding-O0
assess: build rounding-O3
assess: build rounding-fma
assess: build rounding-no-fma
assess: build rounding-fast-math
assess: build rounding-no-vectorize
assess: run 13 member(s)
assess area (mode release): 8 build(s), 13 run(s)
refinement (5 member(s); method: fobis kb show numerics:assessment#refinement-family)
  note: n = 8 16 32 64 128, h = 1/n
  n: the refined parameter itself (its values are the ladder's): not assessed
  trapezoid: observed order 2.000 (triples: 2.003, 2.001, 2.000), monotonic, asymptotic
    extrapolated 1.999999984, error estimate 0.0001, GCI 6.27e-05
  simpson: observed order 4.001 (triples: 4.021, 4.005, 4.001), monotonic, asymptotic
    extrapolated 2, error estimate 4.03e-09, GCI 2.52e-09
precision (1 member(s); method: fobis kb show numerics:assessment#precision-family)
  note: default promotion for gnu: -freal-4-real-8 -freal-8-real-16
  n: every digit agrees with precision-promoted: 16 vs 16
  trapezoid: 16.0 digit(s) agree with precision-promoted: 1.993570344 vs 1.993570344
  simpson: 15.7 digit(s) agree with precision-promoted: 2.000016591 vs 2.000016591
rounding (6 member(s); method: fobis kb show numerics:assessment#rounding-family)
  note: compiler-level perturbations (Verificarlo's Monte Carlo arithmetic is not used)
  n: every member gives the same value; range 16 .. 16
  trapezoid: 15.9 significant digit(s) (>= 15.6 at 95%); range 1.993570344 .. 1.993570344
  simpson: 16.0 significant digit(s) (>= 15.8 at 95%); range 2.000016591 .. 2.000016591
budgets
  simpson: order >= 3.9: 4.001 ok
  trapezoid: order >= 1.9: 2 ok

Reading it:

• refinement: the trapezoidal rule converges with order 2, Simpson's with order 4, both in the asymptotic range (the last triples agree). The Richardson extrapolation and the error estimate say how far the 128-interval result is from the limit; the GCI is the grid convergence index, a conservative error band.
• precision and rounding: about 16 digits agree. In this program the error that matters is the discretisation error, not round-off.
• budgets: both orders met. fobis assess exits with 1 when one is not.

Each family names the section of the numerics knowledge pack that describes its method and its assumptions:

$ fobis kb show numerics:assessment#refinement-family | head -n 12
## Assessing the numerical error of a result > Refinement family  [numerics:assessment#refinement-family]

Run the case at three or more levels of a resolution parameter `h` (grid spacing, time step, tolerance) and
measure how a scalar result converges. With levels 3 (coarse), 2, 1 (fine), refinement ratios `r21 = h2/h1`,
`r32 = h3/h2` and differences `e21 = f2 - f1`, `e32 = f3 - f2`:

- ratio `R = e21/e32`: `0 < R < 1` monotonic convergence, `-1 < R < 0` oscillatory, `|R| >= 1` divergent
  (Stern et al. 2001);
- observed order `p` by fixed-point iteration on `p = |ln|e32/e21| + q(p)| / ln r21`,
  `q(p) = ln((r21^p - s)/(r32^p - s))`, `s = sign(e32/e21)` (Celik et al. 2008; with a constant ratio `q = 0`);
- Richardson extrapolation `f_ext = (r21^p f1 - f2) / (r21^p - 1)` (Richardson 1911) and the error estimate
  `|f_ext - f1|`;

What you learned

[assess:NAME], the refinement, precision and rounding families, budgets. Other families vary the compiler, the thread count and the inputs. Reference: error and uncertainty assessment.

Next: 5. Is a refactoring harmless?.