Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

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?

share|improve this question

4 Answers 4

up vote 20 down vote accepted

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.

share|improve this answer
1  
The final point about tunable parameters reminds me of a quote by John von Neumann, "With four parameters I can fit an elephant, and with five I can make him wiggle his trunk." –  Jed Brown Oct 9 '12 at 22:48
    
This only addresses the effect of discretization error, not that of rounding errors, which are almost universally ignored in fluid simulations as they typically contribute less to the total error. But in applications not involving differential equations or integrals, there is no discretization error. In these problems, rounding errors and iteration truncation errors are the dominant source of errors. –  Arnold Neumaier Oct 10 '12 at 6:49
    
Do you think a rigorous a posteriori estimate could replace the grid independence requirement? As stated it is a little vague, how much do you refine the grid to ensure independence? A good a posteriori estimator on the other hand should leave no room for ambiguity. –  Reid.Atcheson Nov 24 '12 at 22:57
2  
@Reid.Atcheson "Yo dawg, I heard you like error estimates so I got you an error estimate for your error estimate so you can estimate the error in your ...", but unfortunately, practical error estimates tend not be come with rigorous worst-case bounds. What bounds there are available are frequently pessimistic and rely on assumptions that are not true for most engineering problems. Even with the supreme error estimator, you need to verify that it was implemented correctly. By all means, use an error estimator if you have one, but the error estimator is not stand-alone verification. –  Jed Brown Nov 25 '12 at 0:47

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.

share|improve this answer
    
Thanks. I would like to see (ii) at work. Can you give an example where authors combine sensitivity analysis and numerical stability analysis to determine global error estimates. –  highsciguy Nov 20 '12 at 14:06
    
@highsciguy: see my edit to (ii) –  Arnold Neumaier Nov 20 '12 at 16:00
    
Thanks. Jacobian matrix inverse? –  highsciguy Nov 20 '12 at 16:46
    
@highsciguy: don't understand your last question. The example was solving a linear system of equations, no Jacobian figures. –  Arnold Neumaier Nov 21 '12 at 9:07
1  
I overlooked linear. Then it's clear. Are there attempts to generalize the Oettli-Prager theorem to non-linear systems? –  highsciguy Nov 21 '12 at 14:11

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.

share|improve this answer
    
do you have a good reference of an overview article on these polynomial chaos expansions? I've seen the term popping up regularly and would like to learn a bit more about them. Thanks. –  GertVdE Oct 9 '12 at 17:43
2  
Dongbin Xiu generally writes accessible papers on polynomial chaos expansions. Here's one of the more general overview papers he's written: dam.brown.edu/scicomp/media/report_files/BrownSC-2003-07.pdf –  Geoff Oxberry Oct 9 '12 at 17:45

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.
share|improve this answer
    
I just want to comment on point 2. I feel that in numerically solving a real world problem, it is much more likely to see the norm comparing the difference between two successive refinements, rather than with "a ridiculously fine grid". If you can solve the very fine grid, why are you bother with much coarser ones? –  Godric Seer Oct 9 '12 at 23:30
    
The grid refinement studies are impractical on real problems typically. That's where point 3 comes in, it's much easier to determine error bounds by changing scheme order than progressively refine the grid. For instance, we validate our codes using an inviscid vortex which is still non-linear but we can run very fine "exact" solutions to verify accuracy. But on a full swirl combustor, we can't practically do that so we use different schemes. –  tpg2114 Oct 9 '12 at 23:34
    
Also note that refining the grid and checking the difference in the answer shows grid independence, which is a different thing than checking error bounds. Grid refinement is still required even in practical cases, although grid refinement for something like Large Eddy Simulations opens a whole other can of worms. –  tpg2114 Oct 9 '12 at 23:35
    
Thank you, I had incorrectly related grid independence with error bounds. I guess I had always assumed that grid independence implied some qualitative level of guarantee on discretization error. –  Godric Seer Oct 9 '12 at 23:37
3  
Grid independence implies a finer grid will not improve the answer but it doesn't tell you how accurate your answer is or how rapidly refining the grid improved the accuracy. –  tpg2114 Oct 9 '12 at 23:40

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.