Models Methods Software

Chaos Part 0: ODE Chaotic Response is Numerical Noise

In this post we look at a every important issue relative to the foundations of chaotic response as inferred from numerical calculations with systems of ordinary differential equations (ODEs). The relationship of the contents of these notes to calculations with numerical weather prediction (NWP) and general climate models (GCMs) might be discussed in future posts.

Summary Statement
Converged numerical solutions of systems of ordinary differential equations that exhibit chaotic response have never been published. Converged numerical solutions have yet to be obtained for such equation systems. The preceding statements need to be qualified by limiting them to long ranges of values of the independent variable. And of course the statements apply only to the situations for which a given equation system exhibits chaotic response. Additionally, the statements are based on what I will denote as standard numerical solution methods. Standard numerical solution methods include just about all methods that are widely used. Solution methods other than standard, denoted as self-validating, will be specified later in these notes.

Given that chaotic response of complex dynamical systems at its foundation is based almost entirely on numerical solutions of the ODE systems, this situation is almost beyond comprehension. The recently-published paper reviewed here and here on this blog clearly shows the lack of convergence as the step size is refined.

My investigations have for the most part focused on the original 1963 system given by Lorenz 1963. However I have preliminary indications that the statements above will obtain for all systems of ODEs that exhibit chaotic response. Anyone can easily conduct a convergence study by use of off-the-shelf ODE solver software, or by constructing their own solver. The explicit Euler method, for example, is very simple and readily implemented into a computer program. You will find that the calculated numbers will change as the step size is changed and will not converge as the step size is refined. Additionally, the numbers will change if different solution methods are used and as different options within off-the-shelf software are used. This is a very disturbing situation. I have FORTRAN source codes for a few equation systems using a few solution methods that are available to those interested.

Background
Let’s first review some of the definitions of chaotic response of complex dynamical systems. Very frequently non-linearity, complexity, and sensitive dependence of specified initial conditions are invoked as characteristics of ODE equation systems that exhibit chaotic response. However, none of these are sufficient conditions for systems to show chaotic response. Each may be a necessary condition, but even taken together they are not sufficient conditions. Unpredictability and an appearance of randomness are usually taken as one summation of these characteristics whenever numerical methods are used to calculate chaotic systems. Non-linear, complex dynamical systems is also frequently used as a summary description.

The properties of non-linearity, complexity of the ‘modeled’ phenomena and processes, and sensitive dependence on specified initial conditions are generally invoked as a characterization of chaotic response of ODEs. In many cases complexity and chaos are often used to mean the same thing.

However, solutions of mathematical models of complex physical phenomena and processes are not always sensitive to initial conditions. Additionally, chaotic response can also be obtained from simple mathematical constructs; especially, algebraic-map equations. Equally important, chaotic response is sometimes obtained when linear ODEs, even a single one, are not correctly solved by numerical methods. This situation has been discussed previously on this blog along with citations to examples in the literature.

Wikipedia has a few short discussions of chaos; Chaos (physics) and Chaos are two examples. The latter states; “Mathematically, chaos means an aperiodic deterministic behavior which is very sensitive to its initial conditions, i.e., infinitesimal perturbations of boundary conditions for a chaotic dynamic system originate finite variations of the orbit in the phase space; see chaos theory.”

Again note the somewhat subjective nature of these definitions and especially note the complete lack of mathematical statements of necessary and sufficient (necessity and sufficiency) conditions.

The resources of the Mathematica Web site provide this information. Note that this description practically states that a mathematical definition in terms of necessary and sufficient conditions is not yet available. Also note the following taken directly from the note: “Finally, Rasband (1990, p. 1) says, “The very use of the word ‘chaos’ implies some observation of a system, perhaps through measurement, and that these observations or measurements vary unpredictably. We often say observations are chaotic when there is no discernible regularity or order.”” Additional mathematical conditions are given in that document.

For the case of systems of ODEs that exhibit chaotic response, the observations noted in the previous quotation are in fact numbers calculated by numerical methods applied to the system of ODEs. That is, relative to systems of ODEs, the definitions of chaos are somewhat circular. It is only after the numbers have been calculated that observations and analysis of the results can be used to infer chaotic response. I hope we don’t get stuck on trying to pin down a definition of chaos. It is clear to me that chaotic response depends on observation and analysis of calculated response and that a priori necessary and sufficient conditions have yet to be determined.

The necessity for using the calculated numbers to infer chaotic response means that accurate solutions to the ODEs are vitally important. If the numbers that represent the data to be observed are not correct, the conclusions reached from the observations cannot be correct. Unfortunately, as will be indicated by these notes, the numbers are not correct.

“Pretty pictures, wrong numbers.” Anonymous.

Lorenz-like ODE Systems
We’ll first review some basic properties and characteristics of the Lorenz-like equation systems. These systems of ODEs seem to form the basis of much of the research and publication about systems that exhibit chaotic responses. While chaotic response of complex dynamical systems is basically a mathematical enterprise, there should be some connection with physical reality if the concept of ‘chaotic phenomenology’ is to be attached to physical phenomena and processes of interest.

Firstly I think it is important to note that the original 1963 Lorenz equations are not accurate descriptions of any fluid motions. Lorenz-like ODEs are said to represent complex physical phenomena and processes associated with heated fluid flows, although in an extremely simplified manner. The equations, especially the original system of 1963, are severely truncated approximations to an already overly-simplified ‘model’ of a fluid heated from below. Generally, however, the systems, as a model of fluid flow, are not valid beyond the onset of fluid motions. Additionally, almost all examples of calculations with the equations are conducted using a value of the Prandtl number parameter that is outside the range applicable for the fluids of interest. The Prandtl number for air, for example, is about, but less than, 1.0, while for liquid water at atmospheric pressure and temperature conditions is about 7.0. The Prandtl number value of 10.0 is usually taken for most calculations.

Importantly, the values of the Rayleigh number usually used in calculations is far removed from values that would represent the onset of fluid motions. The parameter in the equations is the ratio of the Rayleigh number at the operating conditions to the critical value of the Rayleigh number. As such, numerical values near unity represent the range of applicability of the equations. Tritton 1988 has given calculated results for smaller values of the Rayleigh number. The large values of the Rayleigh number represent conditions for which significant fluid motions would be present in the flow field and the equations are not valid under this condition

There are no physical phenomena or processes relative to fluid motions contained in the equation system that can be invoked as an analogy to any real-world applications at all. More nearly analyses of the fluid motions that the Lorenz-like equations are said to represent have been investigated by Burroughs and colleagues; Burroughs 2003 , Burroughs et al. 2002, Burroughs et al., and Burroughs et al. 2005.

For certain ranges of the parameters, the equations have the property that calculations with the systems show a sensitive dependence on the initial values specified for the solutions. Small changes in the specified initial values lead to large changes as the integration proceeds. The original 1963 Lorenz equations do not exhibit chaotic response for all values of the parameters in the equations. This situation is additional evidence that physical complexity and non-linearity are not sufficient conditions for calculated chaotic response. Either alone, or as in this case, even when both are present.

No analytical solutions are known for any of the Lorenz-like ODE systems; non-linearity prevents obtaining analytical solutions. As a result, all calculations of the dependent variables are preformed by use of numerical solution methods. Generally, this is true for all systems of ODEs that exhibit chaotic response.

We next look into some of the peer-reviewed published literature of results from numerical calculations of Lorenz-like ODE systems.

Literature Review
The focus of the literature review will be on the numerical methods used and the results of convergence studies as the step size is refined. The review has shown that the first published peer-reviewed paper reporting a convergence study for the original Lorenz 1963 equations appeared only in early 2007 [Teixeira, et al. 2007]. The calculations in the paper indicate that convergence is not possible when using standard numerical solution methods. Yao 2003 and 2005 has presented results of convergence studies, some of these papers have not yet appeared in peer-reviewed publications. My own investigations were started late in 2006 and early in 2007 I ran across the papers by Yao. The results of my work are in complete agreement with those of Teixeira et al. and Yao. Convergence of the solutions cannot be obtained. Yao has determined a cause for the lack of convergence.

Lorenz 1963
The original system of three autonomous equations was published by Lorenz in 1963. The numerical method used to integrate the equations, a modified Euler method, was presented in the paper. The step size used for the integrations was 0.01 Lorenz time units (LTU). The equations were integrated from 0 to 3000 LTU, so that 300,000 steps were needed to complete the integrations.

The primary objective of this paper was to introduce the original concepts into the literature. However, as the focus of the paper was the numbers obtained by numerical integration, it is very surprising that convergence of the numerical method was not mentioned. This paper apparently set the tone for almost all subsequent publications of numerical investigations of the chaotic response of complex dynamical systems. A very unfortunate precedent that set the standard for over 40 years.

Absence of a demonstration that integration of the equations has converged, the numbers are known not to be solutions of the continuous ODEs. History has in fact shown that the numbers are not solutions of the continuous equations.

Lorenz [Tellus, Vol. 42A, pp. 378-389, 1990]
A system of three non-autonomous equations was presented in which forcing functions were present in the equations. The equations were first given in Lorenz [Tellus, Vol. 36A, pp. 98-110, 1984]. The numerical method used to integrate the equations is reported to be a fourth-order Taylor-series method; whether implicit or explicit is not given. The step size reported to be 0.025 LTUs, which is reported to correspond to 3 hours. The integrations were carried out for 12 months of time, corresponding to a full year. Convergence testing is not discussed in the paper.

Pielke and Zeng 1994
These authors used the Lorenz 1984 equations for a long term integration study. A fourth-order Runge-Kutta method was used with a step size of 0.025 LTUs, the same as above. (I assume that an explicit Runge-Kutta method was used.) Numerical integrations were carried out for 1101 years, a little over 401,000 days (at 365 days per year). Convergence testing is not discussed in the paper. The integrations required over 385 million step increments (representing an enormous number of operations using a fourth-order Runge-Kutta method). I suspect round-off and global error should have been monitored for these calculations.

Tritton 1988
Tritton 1988 presents discussions of the original Lorenz 1963 system in Chapter 24 of this excellent text. There is no mention of numerical solution methods, step size, or of convergence studies. The author does note that periodic behavior has been reported by Sparrow [1982] for values of the Rayleigh number higher than the chaotic-transition value for some ranges of the other two parameters in the equations (the Prandtl number and the geometric ratio). This finding is another indication that complexity and non-linearity are not sufficient conditions for determination of chaotic response. Calculation with the Lorenz system is given in Chapter 24. The numerical method is not noted and convergence is not mentioned.

Thompson and Stewart 2002
The book by Thompson and Stewart 2002 contains discussions of several systems of ODEs said to exhibit chaotic response. While I have not reviewed in detail all the numerical results presented in the book, it appears that in many cases not even the method used to integrate the equation is mentioned. The step size used receives equally little attention and convergence studies are not mentioned.

Frisch
Dissipative dynamical systems have more than one attractor. Frisch states, “Each attractor has an associated basin. The statistical properties of the solution will then depend on to which basin the initial condition belongs. Thus, not only may the detailed behavior of the orbits be unpredictable (because of the sensitivity to initial conditions), but even their statistical properties may be unpredictable , insofar as it may be impossible to determine to which basin the initial condition belongs. Translated into meteorological vocabulary, this is equivalent to stating that not only the weather but also the climate may be unpredictable.”

Roy and Musielak 2006
These authors have derived Lorenz-like equations systems having more that the three equations of the original system. The systems are determined by requiring them to conserve energy in the dissipationless limit. The authors state that some previous generalized Lorenz systems do not satisfy this requirement. The first paper is based on selecting vertical modes that conserve energy in the dissipationless limit and lead to systems that have bounded solutions. Calculated results are given for the 9-equation system obtained by this procedure. The solution method, step size, and convergence are not mentioned. Lyapunov exponents and power spectra are also calculated and stated to indicate that the response from the equation system is chaotic.

In the second part of the paper horizontal modes are selected under the requirements listed above and a 5-equation system is obtained. The comments relative to the numerical methods given above apply here. Convergence is not mentioned.

Lyapunov Exponents
An interesting situation has developed relative to calculation of the Lyapunov exponents for ODE systems exhibiting chaotic response. As indicated above convergence studies have not been a part of almost all investigations of chaotic response of ODE systems. On the other hand, however, numerical calculations of the Lyapunov exponents associated with such systems seems to almost always include a convergence study. It is kind of ironic that while those calculating the Lyapunov exponents have apparently focused completely on the convergence of the numerical values of the exponents those calculating the solutions to the equation systems have seldom considered convergence. A couple of papers that and resources related to calculation of the exponents are listed here.

Viswanath 1998 has written a PhD dissertation on the subject. The results for the original 1963 Lorenz system are discussed in Chapter 5 and the numerical results for the exponents are given in Table 5.1, 5.2, and 5.3. A ‘standard’ 4th-order Runge-Kutta method was used for most of the numerical integrations. I do not know if ‘standard’ refers to the explicit or implicit version. The step size used ranged from 1/16 to 1/256, with most calculations done with 1/128. An interesting discussion of the possibility of numerical artifacts being present in some of the numerical integrations is also given in Chapter 5. As noted in the previous paragraph above, the author focused on convergence of the exponents while apparently convergence of the numerical integration for the dependent variables was not studied.

Dieci and colleagues have developed theory and computer software for calculation of the Lyapunov exponents. Calculation of the exponents for the Lorenz 3D system is given in . A variety of integration methods were used in the work. Software for calculating the exponents is available online from here. The LESNLL version of the non-linear code available here is set up for calculating the exponents for the original Lorenz system. Interestingly, the values of the dependent variables at the end of a 100 LTU startup transient are displayed as a part of the printed output from the code. If you modify the source code to use smaller step sizes you will see that the values of the dependent variables change at the end of the startup transient. Reprints of publications by this group are available online.

I’m sure that many other reports and papers and methods and software are available online relative to calculating the Lyapunov exponents. I think that the examples given here show that investigations into calculating the exponents have focused on the effects of the numerical methods in more detail than have those investigating numerical solutions for the dependent variables.

The results of investigations into calculating the exponents indicate that the numerical values of the exponents can not be taken as reliable indications of chaotic response. The lack of convergence for the dependent variables means that the calculated exponents have measured numerical noise and not solutions of the continuous equations.

Yao 2003 2005
Yao seems to be among the first to have conducted the classical convergence study for systems of ODEs that exhibit chaotic response. Two papers have been written and posted here in 2003 and here in 2005. The results show that convergence cannot be obtained.

Teixeira et al. 2007
This paper seems to have been the first peer-reviewed and published paper addressing the convergence problem for the original 1963 Lorenz equations. The results given in the paper clearly indicate that convergence is not attained when using a 4th-order Runge-Kutta method. Calculations with NWP codes illustrate that these codes also do not indicate convergence. They recommend that the temporal step size be a parameter in the ensemble-forecasting methodology. This paper has been discussed here and here on this blog. The authors did not investigate the effects of the spatial increment on calculations with the NWP codes. Some of the expected effects of under-resolution when calculating numerical solutions of partial differential equations have been discussed on this blog. It is interesting that Lorenz [Physica D, Vol. 35, pp. 299-317, 1989] has also discussed spurious chaotic response when too-large step sizes are used in integrations of systems of ODEs.

Self-Validating Integrations of ODEs
Recent developments in numerical solution methods for ODEs have focused on self-validating methodologies. An introduction to the methodology by Stauning is available. Self-validating methods have been implemented into the COSY Infinity Site and VNODE software codes. A publications list associated with the COSY software is also available.

A description of the VNODE software is here. Professor Nedialkov’s Web site, from which reprints of reports and papers are available. His thesis is available from here.

The Web site for the VALidation of state ENClosure using Interval Arithmetic for Initial Value Problems (ValEncIA-IVP) project is is here.

I have not yet done any calculations with any of these methods and codes.

It is very striking that the results from self-validating integration methods show the same general trend as than obtained with standard numerical methods. That is, the standard methods seem to indicate for the case of the original 1963 Lorenz system, the calculated numbers begin to show divergence at about 10-to-20 Lorenz time units. The self-validating methods show approximately the same results.

Nonstandard Finite Difference Methods

Summary and Conclusion
Many publications investigating numerical calculations with Lorenz-like systems of ODEs, and other complex dynamical systems, do not include a convergence study as a part of the research. For mathematical investigations of an inferred phenomena founded upon numerical calculations of systems of non-linear ODEs, this would seem to be absolutely unacceptable. The very ‘chaotic phenomenology’ that is the subject of the research is available only from numerical calculations of ODE systems. One would think that a mathematical phenomenon that is based entirely on numerical results would require that those results be demonstrated to be absolutely correct. The results of this literature review shows that this is not the case.

The papers by Yao 2005 and 2003 and Teixeira et al. 2007 are exceptions to the above statement. Yao has presented convergence studies and his results indicate that convergence is not possible. I have independently performed the same kinds of calculations with a variety of standard numerical methods for ODEs. Convergence cannot be demonstrated. Both Yao and Teixeira et al. have presented results of convergence studies and those results show that convergence was not attained. Self-validating numerical methods for solutions of systems of ODEs have not yet successfully integrated the original Lorenz system beyond about 20 Lorenz time units.

Many reports and papers in which calculations of chaotic response have been reported have not addressed the issue of convergence. The few reports and papers that have address this very important issue have shown that convergence cannot be attained.

The calculated numbers that are the data to be observed to determine chaotic response are not solutions of the continuous equations. Chaotic response calculated by use of Lorenz-like systems of ODEs is very likely to be nothing more than numerical noise. It is very likely that chaotic response from all systems of ODEs exhibiting chaotic response is in fact nothing more than numerical noise.

It is known that ‘numerical chaos’ can result from inaccurate numerical integrations of ODEs that should not exhibit chaotic response. This numerical noise does not represent solutions to the continuous equations. The overwhelming evidence from numerical integrations of systems of non-linear ODEs for which chaotic response has been reported also indicates that the calculated numbers are not solutions for the continuous equations.

Because the numbers are numerical noise and not solutions to the continuous equations, the numbers cannot be used to infer anything about the mathematical model used to generate the numbers. And most certainly and importantly, the numbers cannot be used to infer anything about any physical systems.

REFERENCES
The references have been posted in this post.

April 28, 2007 - Posted by Dan Hughes | Uncategorized | Chaos, Numerical methods Verification, ODEs