Chaotic Numerical Instabilities Arising in the Transition from Differential to Difference Nonlinear Equations

For computational purposes, a numerical algorithm maps a differential equation into an often complex difference equation whose structure and stability depends on the scheme used. When considering nonlinear models, standard and nonstandard integration routines can act invasively and numerical chaotic instabilities may arise. However, because nonstandard schemes offer a direct and generally simpler finite-difference representations, in this work nonstandard constructions were tested over three different systems: a photoconductor model, the Lorenz equations and the Van der Pol equations. Results showed that although some nonstandard constructions created a chaotic dynamics of their own, there was found a construction in every case that greatly reduced or successfully removed numerical chaotic instabilities. These improvements represent a valuable development to incorporate into more sophisticated algorithms.


I. INTRODUCTION
Realistic modeling of a physical or social dynamical process involves the formulation of a nonlinear system of differential equations, which in general has to be solved by numerical techniques.This implies to express the original differential equa- tions into finite-difference representations.For a "well behaved" dynamics, the numerical method may play a secondary role.However, there could be a set of parameters of interest where the dynamics may enter an irregular or even chaotic regime, due to the nonlinear nature of the model.In this case the choice of the numerical algorithm is crucial as the scheme itself could importantly interact with the model, producing a dynamics of its own.Such interactions may arise mainly in the transition from a continuous-time dynamic model to a discrete one, as the numerical algorithm maps a differential equation into a difference equation.Thus for different algorithms there will be different finite representations, and the goal is that the resulting discrete-time system reproduces accurately the dynamics of the continuous-time system.Otherwise the finite representation is said to be numerically unstable.
In particular, in a nonlinear discrete-time repre- sentation, the iteration Xk+l F(xk) may become artificially chaotic.These results could be regarded as chaotic numerical instabilities arising from the imposed iterative structure and/or as a consequence of the expansion of the space parameter with the step size as a new "parameter" in the finite- difference construction.Every scheme has its own impact on the model it solves, but the effect could be much more pronounced in the chaotic regime, either because the algorithm interferes with chaos or because chaos is a numerical consequence.
There is an ample variety of numerical methods available, and development of new schemes is an important research.Now, many of the commonly used schemes and even adaptive sophisticated schemes use or include, at some level, very basic standard algorithms, most remarkably the forward Euler.And it is at this level where elementary chaotic numerical instabilities could emerge.The purpose of this work is to handle these instabilities through nonstandard simple finite representations [1], with the goal that if significant improvement is achieved, then nonstandard constructions could replace standard building blocks of far more advanced algorithms.
To this end, nonstandard finite-difference schemes, to solve a system of ordinary differential equations (ODEs) of well-known examples for which chaotic and nonchaotic dynamics are well established, are constructed.These are a photoconductor model [2,3], the Lorenz equations [4] and the Van der Pol equations [5].As there is no unique representation, the "right" finite scheme could be determined (or discarded) by requesting the same stability properties of the continuos-time system.

II. NONSTANDARD RULES
A formal introduction and development to these nonstandard rules for ordinary as well for partial differential equations are given in Ref. [1].For ODEs, the most relevant rules are: (i) The order of the difference equations should be the order of the differential equations.
(ii) (iv) Special conditions that the solutions must satisfy should also be satisfied by the discrete solutions of the finite-difference models.A very important case is for physical systems with positive solutions, thus the computed solution should be nonnegative.

III. NUMERICAL EXPERIMENTS
In this work three systems of ordinary nonlinear differential equations have been studied, namely a photoconductor model [2,3], the Lorenz equations [4] and the Van der Pol equations [5].Several nonstandard constructions were used, and in particular, the forward Euler scheme was a first choice because of its simplicity to represent the derivative terms.For comparison purposes, a conventional forward Euler algorithm was used in all cases.

III.1. The Photoconductor Model
The dynamics of photoconductor equations has been discussed earlier [2,3].By means of three nonlinear coupled ODEs, this system models the kinetic of charge carriers in a semiconductor when illuminated with the proper wavelength, and is given by

FIGURE
Solution for the electron density from the photoconductor equations according to a nonstandard forward Euler scheme.
It has been shown that numerical integration through a fixed step fourth order Runge-Kutta produced a complex transient solution where intermittent periodic stages are alternated with chaotic ones, until a breakdown is reached and solu- tion converges to the asymptotic state.In Refs.[2,3] is shown that this chaotic time series is numerically induced, as integration with an adaptive Gear algorithm replaces all chaotic behavior by a "well behaved" solution.
As an alternative technique, the photoconductor model was integrated through several nonstandard schemes, using the same (realistic) parameters that make this system considerably stiff (that is why Gear is a suitable integrator).It was found [6] that a nonstandard forward Euler scheme with a proper denominator function successfully resolves the model, giving exactly the same solution found by the Gear integrator, see Fig. 1.The difference equations for this system according to a non- standard forward Euler construction are given as (2) and q(h)-T(1-exp ()), r-v-', (3)   where the denominator function (h) takes into account the (smallest) time scales of the system.The use of denominator functions illustrates in the present case the significance of using these non- standard step functions, as they allow the scheme to handle a stiff system as efficiently as a conventional Gear stiffly-stable algorithm.
Besides the nonstandard forward Euler scheme (2), other schemes were tested.First of all, accord- ing to rule (iv) of Section 2, a given scheme must satisfy the condition that the charge densities are positive.In order to help determine the best scheme, the linear stability properties of each finite con- structions were investigated and checked against those obtained for the original system (1).The photoconductor solution of interest has eigenvalues with real negative parts, and thus one stable equilibrium point (n,,m,,p) [2].The nonstandard construction (2) indeed produced positive solutions with negative real part eigenvalues.On the other hand, a standard forward Euler and a nonstandard central difference schemes, even with different choices of denominator functions, invariably yielded one eigenvalue with a positive real part.This was consistent with the numerical overflow observed.It is worth to mention that a central dif- ference scheme produces a second order difference equation, and thus, according to rule (i) of Section 2, a numerical instability is to be expected [6].
In any case, here is clearly shown how a stability analysis can be a helpful tool for defining a suitable finite construction.

III.2. The Lorenz Equations
The Lorenz equations have been extensively studied, have a well-established dynamics [4] and thus represents, a suitable system to test nonstan- dard schemes.The equations are

(4c)
The parameters are cr l0 and b-8/3, with r as a control parameter Two sets of initial conditions (IC) (x0, Y0, z0) were considered: IC1 (-7.004, 10.778, 18.083) and IC2 (0, 1,0).Differently from the physical meaning of the photoconductor equations, there is no positiveness condition to enforce in the solution to indicate which nonlocal substitutions are to be made.Thus, because the Lorenz equations is an artificial system, is somewhat arbitrary the nonlocal modelling to be made.Therefore, we have constructed three nonstandard forward Euler models, depending on the nonlocal substitutions considered.In addition, a standard forward Euler representation was made and integration with an adaptive Adams-Bashfort integrator was used as well.
For all the nonstandard constructions the sub- stitution Xk --+ Xk + and zk zk + in Eqs.(4a) and (4c), respectively was performed.These nonstan- dard equations are given as Xk+ Zk+ Xk qhcryk +her (5) zk + hXkYk l+bh Equations (5) remains fixed and only variations for the construction on the second equation (4b) will be considered.
The first nonstandard construction (NS0) is made by taking y -+ y + in Eq. (4b), similarly to Eqs. (5) and is given as: The second nonstandard construction (NS1) is made by taking Xk Xk + and Yk ---> Yk + in Eq. (4b) for the first and second nonderivative terms respectively, so the new equation is given as y + hrxk hXkZk (7) Y+ + h The third nonstandard construction (NS2) is made by taking Yk --+ Yk + in Eq. (4b) in the second nonderivative term and x -+ xk + in the nonlinear term, so the new equation is given as The first nonstandard set (that is, Eqs. ( 5) with ( 6)) was completely unstable, as it produces numer- ical overflow for any value of the parameter r, with any step size and for both sets of initial conditions.Thus it is discarded.
For the second nonstandard set (Eqs. ( 5) with ( 7)), the results depended strongly on the size of the integration step h, in particular for the chaotic regime.For r=22.4 and IC1, where a chaotic transient is expected [4,7] it was necessary to reduce substantially the integration step to a computation- ally expensive value of h 10 -4 in order to get a breakdown of chaos.As the origin of this transient is due to the formation of homoclinic orbits and the onset of a boundary crisis [7], that means that even for small step sizes the second nonstandard con- struction (NS1) is openly interfering with the Lorenz dynamics, see Fig. 2(a).For comparison purposes, Fig. 2(b) shows the expected results with the help of an Adam-Bashfort adaptive integrator.Similar results were obtained for r= 28 and IC1,  5) and ( 7), for h= 0.01.(b) Solution according to an adaptive Adam-Bashfort integrator.(c) Solution with a nonstandard scheme NS2, see Eqs. ( 5) and ( 8), for h=0.01.
where the Lorenz chaos is fully developed, and the time series shows intermittent quasiperiodic stages.With the NS1 scheme such structure was not displayed unless the step was again sensibly reduced, see Fig. 3(a) and (b).Therefore, this scheme shows a great variation depending sensi- tively on the integration step, and is discarded due to these chaotic numerical instabilities.Interestingly, though, the Adams-Bashfort integrator was useless for IC2 because of a numerical overflow, which was not present for the NS1 scheme (which ).(a) Solution with the nonstan- dard scheme NS1 for h 0.01.(b) Solution with the nonstan- dard scheme NS1 for h=0.001.(c) Solution with the nonstandard scheme NS2, for h 0.01.yielded a transient for r 22.4 and chaos for r 28 if h < 10 4 with IC2).
The third nonstandard set formed with Eqs. ( 5) and (8) gave very good results.For the chaotic regime, the NS2 scheme converged to the expected outcome for the very conventional step size value of h=0.01, see Figs. 2(c) and 3(c).Here the main difference with the previous schemes is in the delocalization of the nonlinear term of Eq. (4b).The effect is that this scheme does not interfere with the boundary crisis at r 22.4 or the intermittent quasi- periodic structure (interior crisis?) at r=28.The NS2 construction is robust to changes of step size, and supports the two initial condition sets con- sidered.
Finally, a construction with a completely stan- dard forward Euler construction was made for Eqs.(4).In this case none of the time series showed a chaotic behavior despite the value of r, or the initial conditions or how small (and impractical) the step size was.Although this standard scheme yielded unrealistic results, this does not mean to exclude others based on forward Euler standard schemes, such as the Runge-Kutta method, which indeed accurately follows the Lorenz dynamics.The point here is that direct and relatively simpler nonstandard constructions can also successfully resolve a model.
In conclusion, these results indicate first, that a nonstandard scheme can interfere with the (chaotic) dynamics of the model.Secondly, for a nonlinear model, in which is not evident an outcome with physical meaning to enforce, and with the possibi- lity of chaotic dynamics, the nonlinear terms should be modeled nonlocally, as nonlinear terms are responsible, after all, for chaotic dynamic.

III.3. The Van der Pol Equations
The Van der Pol equations model a nonlinear electronic circuit with self-sustained oscillations [5].Solution of these equations includes a well- known limit cycle in the phase space representation [7] and a stable fixed point at the origin.Therefore, a proper nonstaldard scheme must comply with these well-established stability properties.These equations are given as where it is easily seen that (0, O) is a fixed point with eigenvalues given as C hv@ 2 4 /1,2 ---]-2 (10) such that e re<0 stableFP, Re(A) 5 => ], (11) e > 0 unstable FP.
According to (11), for e < 0 the fixed point is stable and thus orbits converge to the origin.For e >0, the fixed point is unstable and a stable limit cycle is born [7].For lel < 2, there is a pair of complex conjugates eigenvalues, for which at e 0, the real part of (10) is zero.In addition, dRe/ dell=0>0; all these conditions imply a Hopf bifurcation at the bifurcation parameter c0=0.
Thus a limit cycle is present, and whose stability has been investigated for periodic perturbations.
It has been shown [8] that there is a supercritical Hopf bifurcation, and thus a stable limit cycle is expected.These are the stability requirements that a proper integration scheme is expected to satisfy.
First of all, a nonstandard finite exact scheme for Eqs.(9) and thus with no numerical instabilities for all step sizes have been found, see Ref. [1].This exact scheme is obtained by working "in reverse" with the linear part of the nonlinear system and using rather complex numerator and denominator functions of h and e.The result is semi-explicit equations given as where (h,e) and (I)(h,c) are the numerator and denominator functions.
The set (12) agrees with the stability properties outlined above for system (9), see Fig. 4(a) for e--0.1 and all step sizes.For e 0.1 a limit cycle is observed up to h 1.22, see Fig. 4(b); for h > 1.22,  o c' =-0.1 = 0.1 FIGURE 4 Solutions for the Van der Pol equations according to an exact nonstandard scheme.(a) For c -0.1, a stable fixed point at the origin is obtained for all step sizes.(b) For c =0.1, a stable limit cycle is observed up to h= 1.22.
the limit cycle simply disappears without any chaotic transition (described as follows) following thereafter numerical overflow.Now, for practical reasons, simpler representa- tions than (12) were constructed, hoping that if not exact these constructions could be useful schemes if numerical instabilities are eliminated.Thus, a forward Euler standard and nonstandard schemes were investigated.

(14)
where it is seen that in a discrete representation the stability criterion for the fixed point is equivalent to the real part be inside the unit circle, that is, Re < 1.For e < 0, in particular --0.1, such situation is satisfied only if h <_ 0.08; for h-0.1 it is observed a limit cycle but highly unstable, as overflow is observed for h > 0.1.Therefore, this scheme has a threshold numerical instability, since stability for (0, 0) should hold for all step sizes and not cease at h0.08.
For e--0.1, where the fixed point is not stable and a limit cycle should be present, a limit cycle up to h _< 0.1 was observed; for 0.1 < h < 0.5 the limit cycle begins to become distorted.For h > 0.5 a chaotic transition starts and for h-0.7 there is a fully developed chaotic attractor, see Fig. 5(a).For h > 0.7 numerical overflow is present.Thus, the standard scheme has a chaotic numerical instability with a threshold instability at h > 0. The eigenvalues are exactly as those in Eq. ( 14).
The results for c=0.1 are equal to those of the standard scheme.For e =-0.1 both standard and nonstandard schemes are similar, except that the threshold step for chaos is substantially larger in the latter.That is, for e -0.1, the limit cycle is seen up to h 1.2 (as in the exact nonstandard case, but with distorsion); thereafter, a chaotic transition starts with a fully developed attractor at a quite large step size h 7, see Fig. 5(b).For h > 7 there is a progressively less complex attractor, and, differ- ently from the standard scheme, no overflow is a h=0.6 h=0.7 h=5 h=7 FIGURE 5 Solutions for the Van der Pol equations for c =0.1.(a) According to a standard forward Euler scheme for h=0.6 and h =0.7.(b) According to a nonstandard forward Euler scheme, see Eq. ( 14), for h 5 and h 7.
observed.Of course, at these large steps the numerical significance is lost, and h is more of an ordinary parameter to the difference equations than an integration step.These results indicate that neither the standard (13) nor the nonstandard (15) constructions are appropriate schemes because of the threshold and chaotic instabilities observed.However, the non- standard scheme is comparatively much better for a "reasonable" range of step sizes.These findings encourage the search for other different (and simpler than system (12)) nonstandard construc- tions, including, for instance, nonlocal modeling for both Eq. ( 9) and/or the use of denominator functions.

IV. CONCLUSIONS
In solving a system of differential equations the transition to difference equations via integration algorithm may result in numerical instabilities.For a nonlinear system of equations these instabilities could induce chaos as seen in the photoconductor or in the Van der Pol equations; or interfere with a chaotic dynamics as in the Lorenz equations.It is shown that simple appropriate nonstandard schemes greatly improve or eliminate these instabil- ities.Thus nonstandard constructions are a poten- tially valuable substitute of the fundamental standard schemes present at the core of the more advanced integrators.

7 .
The forward Euler nonstandard construction of Eqs.(9) is based on modeling nonlocally only the nonlinear term of the second equation.The non- standard set is given as xk+l Xk / hyk, Yk + h(eyk xk) The discrete derivative takes the form, for example in nonstandard Forward Euler, as