Scientific standards for numerical errors

Question

In my field of research the specification of experimental errors is commonly accepted and publications which fail to provide them are highly criticized. At the same time I often find that results of numerical computations are provided without any account of numerical errors, even though (or maybe because) often questionable numerical methods are at work. I am talking about errors which result from discretization and finite precision of numerical computations etc. Sure, these error estimates are not always easy to obtain, such as in the case of hydro-dynamical equations but often it seems to result from laziness while I believe that the specification of numerical error estimates should be standard just as much as they are for experimental results. Hence my question: Are there resources which discuss in some detail the treatment of numerical errors or propose scientific standards for the specification of numerical errors which result from typical approximations such as discretization?

Answer 1

Your question is asking about model Verification. You can find numerous resources on methods and standards by searching for Verification and Validation (Roache 1997, 2002, 2004, Oberkampf & Trucano 2002, Salari & Knupp 2000, Babuska & Oden 2004), as well as the broader topic of Uncertainty Quantification. Rather than elaborate on methods, I would like to highlight a community that took a firm stand on the issue.

In 1986, Roache, Ghia, and White established the Journal of Fluids Engineering Editorial Policy Statement on the Control of Numerical Accuracy which opens with

A professional problem exists in the computational fluid dynamics community and also in the broader area of computational physics. Namely, there is a need for higher standards on the control of numerical accuracy.

[...] The problem is certainly not unique to the JFE and came into even sharper focus at the 1980-81 AFOSRHTTM-Stanford Conference on Complex Turbulent Flows. It was a conclusion of that conference's Evaluation Committee that, in most of the submissions to that conference, it was impossible to evaluate and compare the accuracy of different turbulence models, since one could not distinguish physical modeling errors from numerical errors related to the algorithm and grid. This is especially the case for first-order accurate methods and hybrid methods.

They conclude with very direct guidelines:

The Journal of Fluids Engineering will not accept for publication any paper reporting the numerical solution of a fluids engineering problem that fails to address the task of systematic truncation error testing and accuracy estimation.

[...] we must make it clear that a single calculation in a fixed grid will not be acceptable, since it is impossible to infer an accuracy estimate from such a calculation. Also, the editors will not consider a reasonable agreement with experimental data to be sufficient proof of accuracy, especially if any adjustable parameters are involved, as in turbulence modeling.

The current version contains a comprehensive set of criteria and represents a standard that, in my opinion, other fields should aspire to match. It is shameful that even today, awareness about the importance of model verification is absent in so many fields.

Answer 2

No such standards exist, as reliable error estimates often cost much more than the approximate calculations.

Basically there are four kinds of error estimates:

(i) Theoretical analyses proving that a numerical method is numerically stable. This doesn't really give an error bar as the analysis only guarantees that the error made is not worse than a quantified error in the input arguments. It suffices for most scientific calculations as the inputs are also approximate only, so the error made with a numerically stable method is not worse than having used a slightly different (but unknown) input. Most highly regarded numerical methods are accompanied by a numerical stabitly analysis, though one finds hardly any implementation that reports upon request the resulting so-called backward error.

(ii) Asymptotic error estimates. These assume that products of all errors (input errors, rounding errors or discretization errors being the most common sources) can be neglected (questionable if functions are very nonlinear), and use sensitivity analysis to propagate input errors. Together with a numerical stability analysis, this can also capture the effect of rounding errors or discretization errors. The resulting error bars are as realiable as the validity of the assumptions on which they are based. Using automatic differentiation tools, the cost of the error estimate is typically a factor of 1 or 2 in addition to the cost for the approximation. Thus this kind of error estimate is fairly frequent in practice.

[Edit] For example, the Oettli-Prager theorem gives easily computable backward error estimates for the solution of linear systems. Sensitivity analysis tells that these errors must be multiplied by the norm of the matrix inverse, which can be estimated using Hager's estimator (built in into modern condition number estimators).

(iii) Stochastic error analysis: (CESTAC, http://www.sciencedirect.com/science/article/pii/0378475488900705) This is done by overloading all operations with a corresponding stochastic variant that evaluates three sets of arguments and afterwards adds an artificial random rounding error. the final three results are used to compute a mean and a standard deviation the square root of (sum of squares of deviations from the mean divided by 2 = 3-1). This gives a fairly useful accuracy estimate of the rounding error part. This does not account of discretization error, however, which is typically the dominating error in ODE and PDE computations. The cost depends on the programming language because of the overhead in executing overloaded operations. Assuming (which is almost never the case) the overloading carries no time penalty, the cost for result plus error estimate is a factor of 3 compared to computing the approximation only. With a special hardware chip for the stochastic operation, this method could be made as fast as evaluating the approximation itself.

(iv) Interval analysis: This gives rigorous bounds for all error sources if done properly, but except in simple cases it requires lots of experience (or software embodying it) to do it in a way that the bounds do not severely overestimate the true errors. Good interval software is available among others for linear algebra (e.g., IntLab http://www.ti3.tu-harburg.de/rump/intlab/; cost a factor of about 6 if the dimension is large) and global optimization (e.g., COCONUT http://www.mat.univie.ac.at/~coconut/coconut-environment/; may be far more expensive or even be cheaper than approximate global optimization, depending on the problem features). But many other classes of problems easy to treat accurrately approximately (e.g., enclosing the trajectories of the big planets of the solar system over 10 years) are completely out of reach for the current generation of interval methods.

Answer 3

Sort of. There are theoretical error bounds that have been derived by numerical analysts that are usually overestimates, and may not be as useful in practice, because they may involve information that is difficult to obtain for problems in practice. A good example would be the bounds on numerical errors in the solution of ordinary equations, which you can find in Hairer and Wanner's books. Nick Higham's book, Accuracy and Stability of Numerical Algorithms (I may be slightly off about the title) also provides some error bounds on common numerical operations and linear algebra algorithms. The numerical analysis literature is rife with such bounds.

Interval analysis methods have also been used to calculate error bounds; these methods are rigorous, and tend to provide stronger error bounds than theoretical error bounds, but these methods may still grossly overestimate the error in a numerical calculation. These methods have been best exploited (to my knowledge) in global optimization, but are also finding use in uncertainty quantification. Arnold Neumaier has written at least one book on interval analysis methods, and is better qualified to comment on this topic in detail. In addition to potential overestimation issues, interval analysis methods suffer from requiring additional computational infrastructure that requires a retrofitting of existing large numerical simulation packages (such as PETSc, Trilinos, CLAWPACK/PyClaw, etc.) to include interval arithmetic and automatic differentiation (for Taylor-based methods). From what I've seen, there aren't many permissively licensed interval arithmetic and automatic differentiation packages out there, although there are some. Even then, sometimes, these libraries have limited functionality; it's been hard to find a permissively-licensed (LGPL, or BSD-like) interval arithmetic library with BLAS-like functionality.

A posteriori error estimates can be obtained more easily, but aren't rigorous. I'm most familiar with these estimates from work on ordinary differential equations, but they also exist for many methods used to calculate solutions to partial differential equations.

More broadly, methods from uncertainty quantification, such as the use of polynomial chaos expansions, Monte Carlo methods, or other sampling methods can be used to quantify the uncertainty in calculations due to variations in input parameters. These methods should be able to provide some sort of heuristic "error bar" due to variations in parameters, but won't give rigorous bounds.

I believe that you're absolutely right on when it comes to specification of numerical errors: computational science should be just as rigorous about presenting its results as experiment- based physical sciences. There is a great deal of work being done in this area (under the umbrella terms "uncertainty quantification" and "numerical analysis"), and it is my hope that error bars will be included when discussing most computational results at some point in the future.

Answer 4

In addition to the other answers, there are a few additional points to consider.

1. Numerical discretization errors, or at least the order of the schemes, can be determined analytically. The discussion of these errors may be omitted from papers if they use a commonly known scheme.
2. Grid refinement studies where the same problem, usually something simple, is run on progressively finer grids. These are compared with an exact solution, or the solution on a ridiculously fine grid, to find the L-norm, typically L2. The slope of this error estimate gives the order of accuracy.
3. In problems where different numerical schemes are available but grid refinement or exact solutions are not, another method called Richardson Extrapolation will give bounds on the error terms. A good review describing these methods can be found in this paper
4. Lastly, each journal sets its own standards for acceptance. Some are strict, others are not. For example, the AIAA set out its standards here. Other journals have similar information for authors.