# Models Methods Software

## Lack of Convergence, Under-Resolution, and Numerical Errors

The Basic Hypotheses

The following is well known but because it is the focus of this discussion I list it for handy reference.

Convergence is Paramount; Nothing else Counts

The fundamental objective of numerical solution of systems of algebraic equations, ordinary differential equations (ODEs) and partial differential equations (PDEs) is to ensure that the approximations made in order to solve the equations do not in fact influence the solutions. In the case of systems of algebraic equations, it must be shown that the stopping criteria applied to iterative solution methods does not influence the accepted solutions. The solutions are independent, to an acceptable level, of the stopping criteria, and the calculated numbers satisfy the equations.

In the case of numerical solutions of ODEs and PDEs, the discrete approximations and solution methods must fulfill the following requirements; (1) the discrete approximations must be consistent with the continuous equations, (2) the solution procedure for the discrete equations must be shown to be stable, and (3) the calculated results must be shown to be independent of the size of the discrete increments used in the solution. The last requirement is generally denoted to be convergence. Convergence means that the solutions of the discrete approximations approach the solution of the continuous equations. Generally this is demonstrated as the sizes of the discrete representation of the independent variables are made smaller; usually at some fixed ratio that measures important time and space scales in the problem.

Generally it is accepted that consistency plus stability ensures convergence for both ODEs and PDEs. That all three are obtained for any given mathematical model, numerical solution methods, and computer code must be demonstrated by analysis and calculations. Lack of demonstration of convergence of the numerical solution methods as coded should automatically indicate that the numbers being generated are very likely numerical noise and are most definitely not solutions of the continuous equations. Several professional engineering societies will no longer accept papers for publications for which Verification and convergence of the numerical solution methods as coded is not demonstrated. One example of such editorial policies is online and discussed here and here.

A Very Important Fact

Here is a very important fact. While we almost always treat consistency, stability, and convergence as equals, convergence is what really matters. In the absence of convergence of the solution of the discrete approximations to the solution of the continuous equations, all we have are numbers. These numbers in fact do not represent solutions of the continuous equations. See detailed discussions of convergence in all textbooks on analysis of numerical solution methods, and any specialized application area, ODEs as well as PDEs.

If we do not calculate solutions to the continuous equations, of what use are the continuous equations?

Numerical Solutions of ODEs
Numerical solution methods for ODEs are a well-established area of mathematics, science, and engineering, yet fruitful research continues to this day. Recent research has focused on numerical solutions for systems of nonlinear ODEs that exhibit chaotic responses. It is becoming well-known that the dynamics of discrete approximations of ordinary differential equations can differ significantly from the dynamics of the underlying continuous equations. Numerical chaos, spurious steady-state solutions, spurious fixed points, and spurious periodicity. Some recent references are listed below in these notes. You can demonstrate this for yourself by looking at a single ODE and the corresponding discrete numerical approximations. The numerical approximation can transform the continuous equation into an algebraic mapping in which the step size is an additional parameter and for which the algebraic map has dynamical-system responses that do not correspond to those of the continuous equation.

Numerical Solutions of PDEs
The numerical solution method requirements given above apply also to numerical solution methods applied to PDEs. However, the presence of more than a single independent variable greatly complicates the matter, and by an enormous amount. There remains a certain amount of art in the science of numerical solution methods for PDEs.

This latter statement is especially applicability to the case of real-world models of inherently complex physical phenomena and processes that involve multi time scales, multi material constitutes, and complex geometries, and when very extensive temporal and spatial scales occur in the applications of interest. These physical aspects combined with the presence of several independent variables make the numerical solution methods exceedingly complex in their own right. Given the size of the algebraic equation systems that arise from the discrete approximations, and the extent of the temporal and spatial regions of interest, computing power is frequently the limiting process relative to demonstration of convergence to the solution of the continuous equations. The discrete approximations probably result in several millions of unknowns that need to be solved at several millions of spatial and temporal locations; a massive computing problem.

The numerical problems mentioned above for ODEs apply for PDEs as well. The short bibliography given below has citations to papers in which the numerous problems that come with spatial and temporal under-resolution within the framework of PDEs.

Under Resolution in Space and Time is the Norm
Available computing power frequently limits the extend to which rigorous investigations into convergence can be conducted. Additionally it is frequently the case, for practical considerations alone, that the time spent on a given calculation is necessarily limited. Some form of semi-implicit time-differencing methods are generally used so that time-step sizes larger than those imposed by the time constants dictated by the physical phenomena and processes can be employed. That is, the objective of semi-implicit approximations is to be able use larger time-step sizes and the solution methods remain stable. Accurate calculations of the physical phenomena and processes requires that all important time constants be resolved. Under these conditions less than ideal spatial and temporal resolutions are used in modeling calculations of inherently complex phenomena and processes.

This situation does not exempt developers and users of such models and codes from providing detailed information about the limitations of the presented results. The presented results are numbers, not solutions of the underlying continuous equations.

A (very) Short list of References about Under-Resolution
There are a significant number of investigations into the effects of under-resolution in a variety of equation and solution-methods settings. Cases of temporal and spatial under-resolution have been investigated. Both ODEs and PDEs have been studied with both stand-alone relatively simple models and methods as well as within the framework of complex models, methods and computer codes. Summaries of some of the effects of under-resolution are given in the following publications.

L. D. Cloutman, “A Note on the Stability and Accuracy of Finite Difference Approximations to Differential Equations”, Cloutman 1996

L. D. Cloutman, “Chaos and Instabilities in Finite Difference Approximations to Nonlinear Differential Equations”, Cloutman 1998

H. C Yee, P. K. Sweby, and D. F. Griffiths, “Dynamical Approach Study of Spurious Steady-State Numerical Solutions of Nonlinear Differential Equations. Part I”, H. C Yee, P. K. Sweby, and D. F. Griffiths 1990

H. C. Yee and P. K. Sweby, “Dynamical Approach Study of Spurious Steady-State
Numerical Solutions of Nonlinear Differential Equations II”, H. C. Yee and P. K. Sweby 1992

H. C. Yee, J. R. Torczynski, S. A. Morton, M. R. Visbal, and P. K. Sweby, “On Spurious Behavior of CFD Simulations”, H. C. Yee, J. R. Torczynski, S. A. Morton, M. R. Visbal, and P. K. Sweby 1997

P. K. Sweby, H. C. Yee and D. F. Griffiths “On Spurious Steady-State Solutions of Explicit Runge-Kutta Schemes”, P. K. Sweby, H. C. Yee and D. F. Griffiths 1990

H. C. Yee and P. K. Sweby, “Nonlinear Dynamics and Numerical Uncertainties in CFD”, H. C. Yee and P. K. Sweby 1996

H. C. Yee and P. K. Sweby, “Some Aspects of Numerical Uncertainties in Time-Marching to Steady-State numerical Solutions”, H. C. Yee and P. K. Sweby 1996

Timothy D. Sauer, “Shadowing Breakdown and Large Errors in Dynamical Simulations of Physical Systems”, Timothy D. Sauer

D. Orrell, “Estimating Error Growth and Shadow Behavior in Nonlinear Dynamical Systems”, D. Orrell 2004

Jon Lee, “Chaos and Direct Numerical Simulation in Turbulence”, Lee 1995

Some of the problems identified in these investigations include:

(2) Chaotic responses from non-chaotic equations.

(3) Spurious solutions to the discrete equations that do not correspond to solutions of the continuous equations.

(4) Spurious fixed points that do not correspond to those of the continuous equations and spurious periodicity.

Very importantly, all of these are strictly numerical problems. The numbers generated do not satisfy the continuous equations.

AOLGCM Codes
The AOLGCM codes are known to not be capable of attaining calculated results that are independent of the sizes of the discrete increments used in the numerical integration of the model equations. So far as I can determine, there is no attempt at all to demonstrate that the numerical solutions converge to the solutions of the continuous equations. I think this situation is unique in all of science and engineering; maybe unprecedented in that hundreds (thousands?) of papers and reports have been published based on the calculated results from AOLGCM models/codes. The spatial increments used in the integrations are enormous; many spatially significant weather/climate processes are not resolved by the calculations. Many smaller-scale phenomena and processes are represented by algebraic parameterizations. I assume that semi-implicit numerical approximations are used in order to be able to use larger time step sizes. The papers and reports listed above in these notes have demonstrated the adverse consequences of under-resolution of both spatial and temporal step sizes in numerical solution methods.

Identification, or elimination, of the presence of any of the problems associated with under-resolution should be associated with Verification of the properties and characteristics of the numerical solution methods as coded. Excellent starting references for information on these procedures include:

(1) Patrick Roache, “Verification and Validation in Computational Science and Engineering,” published by Hermosa Press.

(2) William L. Oberkampf, Timothy G. Trucano, and Charles Hirsch, “Verification, Validation, and Predictive Capability in Computational Engineering and Physics,” Sandia National Laboratories Report SAND 2003-3769, 2003. See also.

These publications contain extensive reference citations to the published literature on these subjects. A Google search will find many more resources. In particular, the group at Sandia has remained very active in these areas for several years now.

The AOLGCM models/codes have not presented documentation that these Verification processes have been conducted. The numbers presented as calculated by AOLGCM codes are very likely numerical noise that do not satisfy the underlying continuous equations.

February 13, 2007 -

1. I have attempted to upload a pdf file that has the reference citations given in the post along with the abstracts. I hope it is here UnderResolveRefs Comment by Dan Hughes | February 13, 2007 | Reply

2. Excellent post Dan. What you say has been known for decades yet AOLGCM models are accepted by many respected scientists. The willful ignorance is stunning.

Dan,
Is it possible that if you code an AOLGCM on different software development platforms or run it on different hardware processor platforms you will get divergent runs of the same model? Same data, same algorithms, different software development and/or operating system, different hardware processors. I would bet that the same exact AOLGCM algorithms and data would produce divergent model runs on different platforms. Comment by Reid | February 13, 2007 | Reply

3. Dan, this all sounds very familiar from the sort of modelling I am used to. Can you point me at a simple worked example of one of these problems, so I can try to recast it in the way I am used to dealing with things ? Comment by fFreddy | February 14, 2007 | Reply

4. Several professional engineering societies

I’d fudge this by inserting that the societies will no longer accept papers for publications for which Verification and convergence

I think it’s more accurate to say: societies have published policies stating they will no longer accept will no longer accept papers for publications.

The sad fact is, scholarly journals use a peer review process that involves sending papers to mostly univerisity professors, who, let me assure you, do not check to find out whether tests for grid independence etc. were done. They all agree that it should be done. Most perform these checks when doing their own computations (though I’m sure many don’t). However, the peer review system is appallingly bad at *actually checking* *this particular* sort of thing and you’ll almost never see any proof these checks were done appearing in a scholarly article.

I think in the bit on Numerical Solutions of PDEs, you might be able to say something rather stronger in regards to GCM’s. I suspect that GCMs are known to give different results using different grid sizes. Since all computations involve coarse grids, this could lead to systematic errors or biases if we decide to try to determine what might be “real” by averaging the results of a variety of codes all of which are known to use coarse grids. (Basically, if they all share the same type of flaw, the average will retain this flaw. It won’t average out.)

On the semi-implicit time step: Do GCM codes use semi-implicit time steps? (My impression from reading the text of the GISS papers is they run into the courant condition, and overcome it by inserting false numerical diffusion at the poles. Maybe they use explicit time stepping? You could look at the actual code for that– or check the wording in the papers. I don’t happen to know.)

You might have to dig to find the references that tell use we know the solutions are not grid independent, but I’m pretty sure I’ve read that somewhere.

This situation does not exempt developers and users of such models and codes from providing detailed information about the limitations of the presented results.

Absolutely true.

The presented results are numbers, not solutions of the underlying continuous equations.

I guess I’d be a bit less severe here. I think they may not be solutions. (But then, I haven’t looked at these things all that long. ) Comment by Margo | February 16, 2007 | Reply

5. fFreddy, the reports by Cloutman contain relatively simple examples that illustrate some of the problems. Here is a short online note that shows numerical chaos based on numerical integration problems. The 3-equation Lorenz system is not very complex either.

Let me know if these are what you’re looking for Comment by Dan Hughes | February 17, 2007 | Reply

6. I updated the formatting of this post to aid readability. For some reason my eyes glaze over if there are not enough paragraphs or differentiation between paragraph titles and paragraphs.

I hope you don’t mind. Comment by John A | February 17, 2007 | Reply