An Important Peer-Reviewed Paper: Part 1
We’ll now look at some of the results presented in the paper.
Introduction and Background
The authors have introduced the subject of convergence of numerical methods into the field of chaotic dynamical systems. This field is very important in many areas of current intense study and investigation. Numerical models and solution methods exhibit chaotic dynamical-system characteristics in weather and climate modeling, direct numerical and large eddy simulations of turbulent flows, as well as the classical studies of chaotic systems through nonlinear ODEs as introduced by Lorenz and others. The author’s paper seems to be the first in the literature to present results of systematic investigations of convergence into this important field of research and applications.
We will use a conventional notion that convergence means that the difference between the calculated numbers, at all locations along a trajectory, will decrease as the step size is refined. Plots analogous to those in Figure 2 in the paper, for example, will reflect only the order of the numerical solution methods. The importance of convergence has been discussed previously on this blog.
Given the overriding considerations of the importance of convergence in all other areas of numerical analysis and applications, the present situation relative to chaotic dynamical systems is very surprising.
Anyone can easily convince themselves that potentially there is a problem in demonstrating that standard numerical solution methods will converge as the step size is refined, when applied to chaotic dynamical system model equations. Simple explicit Euler or Runge-Kutta methods can be used to illustrate the problem. I have investigated the properties of several standard numerical solution methods applied to chaotic model equations. None can be shown to even begin to converge. The original 3D Lorenz system of 1963 was the primary system studied, with preliminary investigations of other model.
Hundreds, maybe even thousands, of calculations of the various Lorenz, and other, chaotic ODEs have appeared in the literature. For all the reports and papers that I have reviewed, none have investigated convergence of the numerical approximations. Texts that present the many and varied aspects of chaotic systems seem to not give consideration to the possibility that the numerical solution methods might be introducing artifacts into the calculated numbers. In fact, in some texts the numerical solution methods used to obtain the calculated results are not even specified.
While there seems to have been no previous publications of the potential problems with standard numerical solution methods, a few studies have provided hints that problems exist. Studies using very non-standard solution methods, self-validating interval and enclosure methods, present plots of results for only a limited range of time units before validation of the solution is no longer possible.
Lyapunov Exponents
As an interesting sidelight, publications of numerical calculations of the Lyapunov exponents for the original Lorenz system always show that as the step size is refined the exponents approach fixed values. Apparently this is an indication that successful calculations of the Lyapunov exponents cannot be used to judge the chaotic response of systems of nonlinear ODEs. Calculations of the Lyapunov exponents seem to have focused on convergence relative to the exponents, but not the ODEs. As accurate calculation of the Lyapunov exponents seems to be independent of convergence of the ODEs, the exponents seem to be a poor estimators for chaotic responses.
Additional discussions of the results presented in the paper are given in the following paragraphs.
Convergence and Numerical Methods
The authors have presented an excellent contribution into investigations of the status of the accuracy of standard numerical solution methods when applied to simple systems of nonlinear ODEs that exhibit chaotic responses. The additional information presented regarding numerical solutions of the systems of algebraic equations arising from discrete approximations to PDEs is also welcomed. Consistency, stability, and convergence for discrete approximations to continuous equations are of fundamental importance in all studies of mathematical models of physical phenomena and processes.
The numerical solutions of all discrete approximations to continuous equation systems are valid only when the calculated numbers have been shown to be independent of the methods used to solve the equations. For ODEs, as shown in the subject paper, independence is demonstrated by showing convergence; the calculated numbers are independent of the size of the discrete increment used in the solution methods. Until this independence and convergence has been demonstrated, the calculated discrete approximations do not satisfy the continuous equations. The displayed numbers are simply that; numbers, not solutions.
The temporal and spatial discrete increments used in numerical solution methods for PDEs are equally important. While the focus of the paper was the temporal increment in the case of PDEs, demonstration of independence and convergence is only attained when the spatial increment is also investigated. All numerical solutions are required to be shown to be independent of both the temporal and all spatial discrete increments. Convergence, an important part of Verification of numerical solution methods, is required to be demonstrated prior to investigating model error and associated ensemble-averaging methodologies.
The standard approach to demonstrate numerical convergence is to show that the calculated numbers, for all time values, approach the same exact value as the step size is refined. We could present mathematical definitions here, but the process just described is generally accepted as a working definition.
Convergence can generally be interpreted as accurate resolution of all important temporal and spatial scales contained in the mathematical models. The use of (semi- or fully-) implicit numerical solution methods is generally employed whenever some of these scales are thought to be of lesser importance so that resolution is not required nor attempted. Generally, if at later times it is decided that accurate resolution is needed, smaller discrete-increment sizes must be employed in attempts to recover the resolutions. This approach is counter productive, however, generally requiring significantly more CPU operations and associated computer resources than explicit methods. Additionally, the numerically-induced dispersion and diffusion obtained with semi- or fully-implicit numerical methods generally requires that excessively smaller increments be used.
More on Under-Resolution
As noted in the paper, there are in the literature investigations of some of the effects of under-resolution in the solutions of both linear and non-linear small systems (including a single scalar equation) of ODEs. The results of these investigations have shown that numerical solution methods can in fact result in numbers that appear to represent chaotic responses when chaos should not be present.
We will not present much at all regarding the effects of under-resolution of numerical solution methods by use of discrete step sizes that do not resolve significant scales of the equation systems. We only mention that some of the effects of under-resolution are just as of significant importance as the lack of convergence on which we focus. Lack of convergence and under-resolution, after all, are related issues. Under-resolution is the subject of this post on this blog.
Status of Numerical Solutions of the Lorenz Equations
Calculations of the Lorenz 3D equations with all standard numerical solution methods for ODEs will never indicate convergence. It is possible that this statement is holds also for all non-linear ODEs that are assumed to exhibit chaotic behavior. The same question relative to PDEs is still a very much open issue. The numerous well-founded problems associated with attaining accurate numerical solutions ot PDEs make it very difficult to separate these issues from the issues associated with chaotic response. Little of a basic theoretical nature is known about chaotic response and PDEs. It is important to note, however, that some methods of solution for PDEs depend of the same type of ODEs that are used as illustrations of chaotic response.
Regarding PDEs, it seems that chaos/chaotic response is simply invoked whenever the numerical solutions are found to be sensitive to initial values. Sensitivity to initial conditions, along with non-linearity, are not sufficient in themselves to warrant the assumption of chaos. Additionally, invoking shadowing to give credence to numerical solutions of PDEs is also not based on a solid foundation. Successful shadowing requires that other specific requirements be met, including starting near a true solution. If chaotic equations cannot be accurately integrated in a manner that ensures that true attractors have been initially determined, how can shadowing be ensured. Additionally, the breakdown in the numerical solutions, due solely to numerical errors, cannot be overcome by shadowing.
My Experiences
My own investigations with standard numerical methods and non-linear ODEs, specifically the original 3D Lorenz system, have shown us that the results presented in the paper are indeed the expected results. However, in contrast to the assumption stated in the paper, the calculations of the Lorenz equations presented in the paper are clearly not converged. Use of the smallest step size as the reference case is but a convenient rationale for not continuing the search for convergence. Any other step size could have been taken as the reference case.
Results of research by me indicate that chaotic response of models based on simple non-linear ODEs like the 3D Lorenz equations, cannot be, and have never been demonstrated to be accurate solutions of the ODEs when standard numerical solution methods are used. This conclusion holds no matter what standard numerical method is used. The numbers calculated do not satisfy the continuous form of the ODEs; never have been, never will be. The numbers are solely numerical errors, and the chaotic response is due completely to ‘numerical ODE chaos’.
In the present discussions we will expand on some of the issues addressed in the paper and present additional information regarding convergence of numerical solutions of nonlinear equation systems that exhibit chaotic responses. Our focus is also the numerical solution of systems like the Lorenz 3D ODEs used in the paper, along with a few ideas regarding numerical solutions of PDEs.
Review of the Results in the Paper
The results in Figure 2 indicate that the magnitudes of the difference between the calculated responses are of the order of the magnitude of the dependent variables themselves. It is very striking that for the case of even the next-largest step size above the reference case, the error clearly grows as the integration continues out in time. The differences even out past the saturation time show remain large. If convergence was being approached then the lines in Figure would be more-or-less parallel lines that reflected only the truncation order of numerical solution methods. My calculations are entirely consistent with the results presented in Figure 2.
Additional numerical artifacts are shown in Figure 3 of the paper. In Figure 3a, for example, the calculations indicate that the value of the asymptotic fixed points are clearly a function of the step size. (Note that the introductory text for this figure indicates that r = 17, while other text discussing the figure, and the figure caption indicate r = 19.) And for Figure 3b, the calculations for the mid-range step size remain ‘chaotic’, while step sizes both larger and smaller approach fixed points. Both these results are clear indications that there are problems in the numerical solution methods. Again, my calculations are entirely consistent with the results presented in Figure 3.
The plots of results of calculations with the QG and NOGAPS models and codes, representing solutions of systems of PDEs and PDEs plus algebraic modeling equations, respectively, also indicate lack of convergence. The discussions in the paper about these codes also indicate that only the temporal step size was investigated. A more nearly complete investigation of convergence should also include refinement of the spatial increments in addition to the temporal increment.
The calculations with the NOGAPS model/code show clearly why convergence studies are an absolute necessity. Without them it would not have been known that the asymptotic range for the solution methods for the finite-difference equations had not been reached. Equally important, calculations performed with assumed higher-order numerical methods and increased spatial resolutions might not in fact represent calculations at higher resolutions.
The behaviors and characteristics seen for the NOGAPS model/code can easily be present in all models/codes that have important physical phenomena and processes represented by algebraic equations. Careful treatment of the algebraic equations is required in order that they do not adversely impact the analysis of truncation errors. Development of higher-order numerical solution methods for the continuous PDEs will be to no avail if the algebraic equations reduce the effective order.
Summary
The calculations of the Lorenz equations presented in the paper are clearly not converged. This result is consistent with my own calculations of small systems of nonlinear ODEs that exhibit chaotic response.
The assumption that the smallest discrete step size represents convergence is not a good assumption. The results presented in the paper clearly show that convergence of the numerical solutions has not been attained.
All the calculated results given in the paper are functions of the step sizes used in the numerical solution methods. This is an unacceptable situation especially relative to models of physical systems. The step size is not a parameter of the physical systems of interest.
As noted in the paper, semi-implicit methods have been developed so as to be able to use larger step sizes. These larger step sizes will very likely lead to under-resolution and the associated problems discussed above in these notes and on this blog. Notable among the adverse effects of under-resolution are the false indications of chaotic responses. Extraction of evidence of chaotic responses by use of the displayed numbers from a calculation with a very complex mathematical model and associated computer code provides, at the very best, a measure solely of that calculation. No conclusions relative to the properties and characteristics of the continuous equations, and especially the modeled physical system, is possible and should never be accepted.
The sizes of both the temporal and spatial discrete increments are equally important. While the focus of the paper was the temporal increment in the case of PDEs, demonstration of independence and convergence is only attained when the spatial increment is also investigated. The numerical solutions are required to be shown to be independent of both the temporal and all spatial discrete increments. In the absence of convergence, the calculated numbers do not satisfy the continuous equations.
The importance of careful numerical-methods treatment of any algebraic equations used to model physical phenomena and processes has been demonstrated in the paper. These have actually reduced the order of the numerical methods relative to the theoretical order of the discrete equations.
Convergence is required to be demonstrated prior to investigating model error and associated ensemble-averaging methodologies. Use of the time-step size as a parameter for variation in the ensemble-averaging approach is not ever appropriate for any models of physical phenomena and processes.
Models Methods Software 
A comment on this paper, and a response by the authors, has been published by the JAS. The TOC for the issue is here. The comment can be read; the response is not yet (September 4, 2007) available online.