Detection of the Onset of Numerical Chaotic Instabilities by Lyapunov Exponents

It is commonly found in the fixed-step numerical integration of nonlinear differential equations that the size of the integration step is opposite related to the numerical stability of the scheme and to the speed of computation. We present a procedure that establishes a criterion to select the largest possible step size before the onset of chaotic numerical instabilities, based upon the observation that computational chaos does not occur in a smooth, continuous way, but rather abruptly, as detected by examining the largest Lyapunov exponent as a function of the step size. For completeness, examination of the bifurcation diagrams with the step reveals the complexity imposed by the algorithmic discretization, showing the robustness of a scheme to numerical instabilities, illustrated here for explicit and implicit Euler schemes. An example of numerical suppression of chaos is also provided.


I. INTRODUCTION
The numerical solutions of a nonlinear set of differential equations are discrete approximations of closed form solutions, which are usually difficult, impractical or impossible to obtain.The reliability of such approximations will depend largely on the integration algorithm and, accord- ingly, the first effort is devoted to the search of an accurate scheme.The stability of the algorithm employed may become the next important issue, followed, finally, by considerations regarding the 121 computational efficiency in terms of computa- tional speed.The unbalance of these three aspects of a numerical scheme is mostly due in that the interest in a system with a complex or even chaotic dynamics is more of an object of study than a matter of practical importance.Notoriously, con- siderations on speed have been traditionally relegated, as long as the scheme is accurate and stable.These requirements are inconsistent with the dynamical nature of the systems being solved, and many problems may demand a fast as possible computation of trajectories, because they are occurring in real time and/or because a simulation is being modeled by a large number of equations, a situation often found in engineering.A direct way to explore the balance between accuracy, stability and speed is by means of numerical integration with conventional fixed-step schemes, in terms, specifically, of the size of the integration step.Typically, very "small" step-sizes usually guaran- tee the fidelity of the solution (even though too small steps may have the opposite effect [1]), but to the expense of longer calculations.
Widely used fixed-step schemes are the second and fourth order Runge-Kutta, and the explicit and implicit Euler schemes.In the implementation of these schemes a continuous system of differ- ential equations is typically mapped into a discrete representation of difference equations, which hopefully will share most of the properties of its continuous counterpart.Therefore, the diver- gences produced by this algorithmic discretization will be regarded as numerical instabilities, and the algorithm is said to be numerical unstable.This work presents a procedure to estimate the largest possible integration step size of elementary standard fixed-step algorithms, before the onset of certain numerically instabilities arising in non- linear systems of differential equations.The instabilities considered here are the numerically induction or suppression of chaos, because the nonlinear nature of the system being integrated allows the possibility of masking a chaotic or periodic regime.Therefore, the method reported is based upon a standard evaluation of the largest Lyapunov exponent [2,3], on account of the observation that this number changes distinctively with the step size.
To this end we worked some numerical exam- ples of physical importance known to have a limit cycle, and/or a chaotic attractor, depending on the parameter settings.Algorithmic discretization is performed employing a standard forward Euler scheme.This elementary scheme is important because despite its simplicity, most sophisticated methods are based on simple Euler routines and used to solve very complex problems, such as parabolic and hyperbolic partial differential equa- tions [4].Moreover, in many other numerous cases combination with/from Euler schemes (like several forms of predictor-corrector schemes) are com- monly used to solve systems of ordinary differ- ential equations.For comparison purposes, a nonstandard backward Euler scheme is employed.In Section II are presented the examples with the Euler discretization.Section III presents results illustrating the procedure and a discussion on the relation of certain numerical effects to computational errors.Finally, in Section IV are given some conclusion remarks.

II. NUMERICAL EXAMPLES
We employ three well-known examples of non- linear oscillators, namely, the Lewis oscillator, the Van der Pol oscillator and the two-well Dulling oscillator.The two first systems undergo a Hopf bifurcation, where a fixed point is losing stability giving birth to a limit cycle; this type of bifurcation is analytical found by examination of the eigen- values when varying a system parameter.Likewise, the forced Duffing oscillator is able to display a periodic cyclic behavior in addition to a well defined chaotic dynamics.The equations for these systems are as follows: (i) The Lewis oscillator is given by, dt y (la) dt with e being a parameter.(ii) The Van der Pol oscillator is given by, d--y (2a) --x + ey-ex2y (2b) where e is the parameter of the system.
In system form, the Duffing system is represented by, dt dz dt where x is position, y is velocity and z and F are the phase and amplitude respectively of the periodic driving force.Parameters values were taken as a 1, b 1, c-0.5, co= and initial conditions x0 Y0 z0 0.01.
We employ the Duffing oscillator to illustrate the algorithmic discretization.According to a standard forward Euler scheme (SFE), the first order differential equations (4) are mapped into difference equations given by, x+ x / hy (5a) y/+l y/(1 hc)+ h(axl bx3 + Fcos (z/)) (5b) where the derivative terms were modeled following the basic calculus definitions xk + xk/h y, with the step size h--tg+l-t being the integration time interval; x, yk and z are the discretized variables.The discretization yields an iterative structure (5) for the variables x+ 1, y + and z+ l, which is implemented straightforwardly.Analogously, it was constructed a backward Euler discrete representation for the Duffing oscillator (4).Backward schemes are, in general, based on the implicit representations Xk+a= Xg + hf (Xk + 1, tk + 1), in contrast to the explicit constructions as in Eqs. ( 5).In addition, a nonstandard Euler construction is simply made by replacing nonlinear terms by nonlocal repre- sentations, such as yZ--+y+ly, or y3__+y2 (Yk + +Yk)/2, etc.The derivative denominator, given as h, can be replaced by a more elaborated function ff(h) of the integration step; in this work it was simply kept as (h)=h.For details on nonstandard techniques, see Ref. [5].
An example of a nonstandard backward con- struction is given by the Van der Pol equations (2); here, it was only considered the replacement Y --+ Yg-instead of y -Yk in the nonlinear term of Eq. (2b), yielding the iterative equation, Yk+I Yk (1 + he heXk+) hxk+ (6)   It will become manifest in the following discussion of results that these apparently "innocuous" changes provided great flexibility to the numerical integration, by reducing considerably many of the numerical instabilities founded in the standard counterparts, including the usually more robust implicit versions.

III. RESULTS AND DISCUSSION
Numerical instabilities are regarded as those dissimilarities between the dynamical behavior of the continuous system and the numerical solution produced by a given scheme.In fixed-step schemes there are a variety of numerical instabilities: the most common, by far, is the threshold instability, where beyond a critical step size, numerical solutions began to differ [5]; creation instabilities, because of the appearance of spurious additional points and typically found in higher order integrator [6]; chaotic instabilities [7-11, 14], were chaotic output is generated by the numerical discretization; numerical overflow, highly undesir- able because of the dramatic outbreak of the computation when a critical value of the step size is reached [9]; and computational alterations of Hopf bifurcations [10].
In the present work we examine the onset of unwanted chaotic/periodic behavior, which will be monitored by the evaluation of the Lyapunov characteristic exponents (LCE).These Lyapunov exponents provide a powerful dynamical diagnos- tic of the chaotic status of a system, as they are related to the exponentially divergence or conver- gence of nearby orbits in phase space [2,3].
At this point is convenient to mention that in one or higher dimensional maps is enough to have one positive Lyapunov exponent for chaos, while for continuous dissipative system chaos is present for one or more positive Lyapunov exponents, provided we have no less than a 3-dimensional system [1-3].Therefore, the presence of a chaotic attractor, that is, one LCE > 0 in the 2-dimen- sional Lewis (1) or Van der  Pol (2) oscillators, will be considered a numerical breaking resulting from a discretization effect.

III.1. Detection of Numerical Chaos
The 2-D Lewis oscillator provides a typical example of computational chaos.The real part of the eigenvalues A [1,3] of the continuous sys- tem (1) is given by Re[A] =e/2.The system has a fixed point at the origin which is stable for e < 0 and a stable limit cycle for e > 0. That means that the computed largest LCE for this system should be negative for e < 0 and equal to zero for >0.
Figure shows the largest Lyapunov exponent for the Lewis oscillator (1) as a function of the step size when integrated according to a standard forward Euler algorithm ( =0.1 and initial con- ditions x0=Y0=0.01).This figure shows unequi- vocally that up to a step size h*0.45 the computed LCE remains approximately zero, as expected for the limit cycle.Beyond that value begins an increasingly intricate dynamics, due to an assortment of numerically induced instabilities, that includes numerous n-periodicity windows, chaotic regions and finally a breakdown at h 1.178 due to numerical overflow.This complex dynamics is easily visualized with the help of the corresponding bifurcation diagram, which com- plements the information provided by the LCE vs. h plot, see Figure 1.
Figure shows the start of instabilities, clearly indicated by evident changes in the LCEs for h > h*.Furthermore, is important to notice that this step size h* is more than a mark signal to limit computations; these limiting steps values were found to be usually larger than the "typical" setting h 0.0001 0.01.This means that if the step size could be safely increased by one or more orders of magnitude, there will be a substantial reduction on the CPU time, which may be a crucial factor of practical importance.
The Duffing oscillator offers an interesting example of deterministic chaos appearing in systems with a nonlinear potential.In the present example, the amplitude F of the driving force will determine a behavior ranging from periodic, at F <_ 0.1, to fully chaotic, as in F--0.4.In Figure 2 are plotted the largest Lyapunov characteristic exponent as a function of the step size for various values of F. For F-0.1 the computed LCEs are approximately zero, as expected for the limit cycle.At a step size h* 0.16 the LCEs becomes positive, explicitly indicating the onset of computational chaos.Again, a complex dynamics sets in and finally, at h-0.391, numerical overflow stops further calculation.
A chaotic dynamics is expected for F--0.4 and, in agreement, the computed LCEs are positive.
However, for h > 0.03, the shape of the computed attractor in phase space becomes different, and a more broaden and/or "scattered" chaotic attractor takes place as the step size approaches the value that overflows calculations.In fact, the scattered appearance of the attractor appears to be a frequent feature of numerically induced chaos, as observed in other systems in our work [9][10][11].This increased degree of chaoticity (reflected by an increase of the LCEs) in which a chaotic attractor is changing its shape, may be caused by the occurrence of an internal crisis [12] forced by round-off and/or truncation errors [1, 4, 13, 14].
For the same reason, a (internal) crisis may

FIGURE
Largest Lyapunov characteristic exponent and bifurcation plots as a function of the step size of the standard forward Euler scheme used to integrate the Lewis oscillator ( 1) for e --0.1.The onset of numerical instabilities is marked distinctively about h* 0.45.0 0.1 0.2 0.3 h FIGURE 2 Largest Lyapunov characteristic exponent as a function of the step size of the standard forward Euler scheme used to integrate the two-well Duffing oscillator (3) for various values of the amplitude of the driving force.For F=0.1 (where a limit cycle is expected) is clearly observed the numerical induction of chaos at h* 0.16.In contrast, for F=0.7 (where a chaotic regime is expected) it is shown a numerically suppression of chaos at h* 0.05, followed at h** 0.14 by a numerical transition periodicity/ chaos similar to that for F--0.1.
produce other type of attractors, like the very pronounced quasiperiodic regime shown in Figure 2 for the valley in F= 0.4.
All these dynamical behavior could be conven- tionally understood as derived from the new system created by the algorithmic discretization, with the "parameter" h as our control parameter.But more interesting is the connection of this dynamics with the ideas of shadowing [15,16].Formally, the critical step h* that signals the onset of the new dynamics, could be related to the maximum step value defined in the shadowing theories of Refs.[15][16][17].According to these ideas, a true or suitable solution fi of a dynamical system du/dt-F(u) follows closely or "shadows" a computed solution, or more specifically, a pseudo- orbit Pn, which is the computed solution together with the computational errors (round-off and local errors).In single-step and multi-step discretiza- tion, the numerical method will reproduce, within a tolerance, the true orbit for h _< [17], such that It(nh)-pn] <_ ChI':, where C and are positive constants and K is the order of the numerical method.In the context of the present work, clearly h*.For h > h*, the computed solutions are no longer reliable and a complete new dynamics takes place, as shown distinctively in Figures and 2.
Therefore, this h* represents a limit between (approximately) continuous and discrete dynamics.

III.2. The Numerical Suppression of Chaos
This procedure of examining the LCEs as a function of the step size also helps in the detection of the numerical suppression of a chaotic dynamics.And like the numerical induction of chaos on the previous examples, the actual suppression occurs in a well-located range of step-sizes.This situation is illustrated in Figure 2 for the thick curve F 0.7.Here is possible to see that for a step size lower than h* 0.05, the computed LCEs are positive, as expected, but thereafter they sharply become negative for an extended range of steps sizes, where a numerically induced periodicity takes place [13,18].And again, towards the values of numerical overflow (h 0.168 for F= 0.7), a numerical crisis occurs, and the characteristic scattered numerical chaotic attractor takes place.In Figure 3 is shown the bifurcation diagrams for the Duffing oscillator for F E [0.1, 1] (eliminating 5000 first points), for h 0.001 and  3), integrated with a standard forward Euler scheme with step size h 0.001 and h 0.147.Arrows indicate the suppression of chaos by the advance of a numerically induced periodic front.h 0.147.Notice the initial region for F < 0.3 at the left of the bifurcation plots, where the expected periodic behavior is sustained.From then on, the rough borders in Figure 3 indicate chaotic regions and the subsequent smooth borders corresponds to periodic regions.The effect of the computa- tional suppression of chaos in the Euler scheme is strikingly clear in these plots: in the upper plot, at F 0.7 the system is well within a chaotic region, which becomes artificially periodic with h increasing, because a periodic front beginning at F 0.82 literally swallows the chaotic region up to F 0.490 (so then F 0.4 still remains chaotic) as indicated by the arrows pointing to the left.
From Figure 3, is observed that at F the Duffing system will be unaffected by the advance of this numerical periodicity and will remain in the periodic region for the range of step-sizes con- sidered, including for the value h 0.147, very close to the step that overflows F= 1.In consequence, the LCEs for F 1, are approxi- mately zero, see the gray curve in Figure 2.

Nonstandard Scheme
In Figure 4 is shown both the LCEs vs. h and the bifurcations plots for the Van de Pol system (2); integration is done according to the nonstandard implicit scheme (6) described in Section II, for e 0.1 and initial conditions x0 y0 0.01.The same integration carried out with a forward Euler standard algorithm produced a much more com- plex bifurcation diagram (similar to Fig. 1) than that shown in Figure 4 for the nonstandard backward Euler scheme.Compared to the bifurca- tion plots obtained for an explicit SFE, Figure 4 shows the improvements of the nonstandard scheme, because the limit cycle is sustained for larger values of the step, and the chaotic and bifurcations regions are much more reduced, preserving a more uniform dynamics.In the integration of the Van der Pol system with the nonstandard backward Euler scheme (6), the deviation from the limit cycle was observed at about h* 1.20 (omitting a small region at h 0.9) compared to h* 0. Van der Pol oscillator (2) for 0.1.Expected limit cycle and a mostly instability-free dynamics is sustained for an extended range of steps, about h* 1.20.explicit version.And the nonstandard implicit scheme overflows the calculations at a step h 1.71, which is comparatively much larger than h 0.73, of the SFE integrator.Although is well known that implicit schemes are usually more robust [1], superior results were obtained with the simple nonlocal replacement already mentioned.In our work, the effect of nonstandard replace- ments was almost always to enhance stability and to allow larger steps before the onset of instabil- ities, which usually were greatly reduced or even eliminated, see Refs.[5, 7, 9, 19].Accordingly, this example confirms such properties by means of the LCEs vs. (h) and the bifurcation plots technique, see Figure 4.

III.4. Effect of the Computational Errors
Although there is an increasing evidence of the incidence of numerical chaotic instabilities due to computational errors [10, 13, 14, 18, 19], the exact mechanisms are still poorly understood.In part, because it is not well defined the build-up of such errors, and how the interaction model/scheme affects the numerical output [9].The relation of computational errors (mainly truncation and round-off) to the step size has been regarded mainly as deterministic, like 0(hK), where K is an integer, usually the order of the scheme [1,4]; however, some evidence points out to a chaotic buildup of errors with the step size [10], probably a more consistent approach when dealing with nonlinear dynamics.
The algorithmic discretization of a continuous nonlinear system introduces the step size as an extra parameter, affecting the dynamic of the "equivalent" difference system, which has now the dynamical properties of a nonlinear map, see Eqs. ( 5) and (6).For certain values of the h parameter, the spurious apparent chaos is presumably trig- gered by round-off errors in a mechanism similar to sensitivity to initial conditions.From this per- spective, it may be possible to explain the manifest crisis events, that is, collisions with saddle-type objects being artificially generated and producing the many bifurcations shown in the bifurcation plots of Figures and 4.
In the other hand, the forced periodicity ob- served in Figure 2 for the Duffing system with F 0.7, as well as the numerous periodic windows shown in the bifurcation plots of Figures and 4, could be explained in terms of the finite arithmetic precision of the computation, which forces chaotic trajectories to became periodic [18].This could be because the truncation and round-off errors ex- cludes the possibility of an aperiodic (infinite digits) dynamical evolution and after some time of computing, the numerical trajectories may begin to repeat themselves.
In short, the combined effects of computational errors, enhanced by larger values of our parameter h in the discrete system, could affect the computed trajectories by washing out the correlation (shadowing) with a "true" trajectory after some time.

IV. CONCLUDING REMARKS
The manifest changes experienced by the LCEs with the step size provide explicit information about the numerical effects over periodic and chaotic dynamics.Construction of the LCEs vs. h plots allows establishing a range of workable values for step-sizes before instabilities, and shows a characteristic way to undergo numerical chaos.
To complement the LCEs plots, the examina- tion of the bifurcation diagrams facilitate the visualization of numerically induced instabilities as a function of the step size and thus help to assess the robustness and/or quality of the integration scheme.
In summary, the procedure outlined in this work provides a simple and direct criterion for the selection of much-larger-than-usual step-sizes of commonly used fixed step algorithms, under the premises of minimum instabilities for the shortest computation time possible.
The driven two-well Duffing oscillator is given by, d2x dt 2 ax + bx + C--z-F cos cot dt

1 FIGURE 3
FIGURE3 Bifurcation plots as a function of the amplitude of the driving force for the Dutting oscillator (3), integrated with a