Convergence Criterium of Numerical Chaotic Solutions Based on Statistical Measures

Solutions of most nonlinear differential equations have to be obtained numerically. The time series obtained by numerical integration will be a solution of the differential equation only if it is independent of the integration step h. A numerical solution is considered to have converged, when the difference between the time series for steps h and h/2 becomes smaller as h decreases. Unfortunately, this convergence criterium is unsuitable in the case of a chaotic solution, due to the extreme sensitivity to initial conditions that is characteristic of this kind of solution. We present here a criterium of convergence that involves the comparison of the attractors associated to the time series for integration time steps h and h/2. We show that the probability that the chaotic attractors associated to these time series are the same increases monotonically as the integration step h is decreased. The comparison of attractors is made using 1) the method of correlation integral, and 2) the method of statistical distance of probability distributions.


Introduction
When using as nonlinear models ordinary differential equations (ODE's) or delay differential equations (DDE's), it is necessary, in general, to resort to numerical simulations for analyzing their dynamical behaviour.Numerical solutions of differential equations are obtained by discretizing them, and it should be recalled that solutions of the continuous equation are always solutions of the discretized equation but not vice-versa.The time series obtained by numerical integration will be a solution of the continuous differential equation if it is independent of the integration step h, indicating that convergence of the numerical integration has been reached.Usually the criterium of convergence of the numerical integration is based on the direct comparison of the numerical time series calculated using steps of integration h and h/2: the difference between these time series should decrease as h is decreased, and convergence is achieved when this difference is smaller than a given precision ε.It should be remarked that the comparison has to be made after stationarity is reached, that is, the transient must be removed before comparison is made.
Nonlinear systems can exhibit complex behaviour and, in the case of chaotic solutions, the direct comparison of the calculated time series, as the integration step is halved, is not appropriate for establishing numerical convergence, due to their characteristic sensitivity to initial conditions.In fact, it is very common the assertion that chaotic solutions do not converge (see, for instance, Teixeira et al. (2007) [1], Yao and Hughes (2007) and (2008) [2,3].The sensitivity to initial conditions causes the blow up of the difference between the time series calculated using integration steps h and h/2, respectively, regardless of it being initially extremely small (Teixeira et al. [1] state that "for chaotic systems, numerical convergence cannot be guaranteed forever").Sauer (2002) [4] has shown that numerical integration cannot produce approximate correct individual (or even average) chaotic trajectories due to shadowing breakdown, even though it may produce an approximately correct chaotic attractor.So, in the case of chaotic solutions, as convergence criterium of the numerical integration, we propose the comparison of the chaotic attractors associated with the aperiodic time series calculated using integration steps h and h/2.The chaotic attractor is the object that is invariant, thus our proposal of a convergence criterium of the numerical integration based on the comparison of the attractors as the integration step is halved.Of course this criterium is applicable also to the case of time series corresponding to periodic or quasiperiodic solutions.
DDEs, due to the presence of delays, have infinite dimension, but the corresponding discretized equations are finite dimensional so that, in this case, numerical convergence is a problem even more delicate, so we shall illustrate our method using the multilooped negative feedback equations that have been proposed by Glass et al. [5] for modelling physiological control systems.It should be stressed that our method can be applied for ODE's also.In Section 2 we briefly present the multilooped negative feedback equations for three delays, and introduce the notation for the time series solution obtained numerically.We exhibit the blowing up of the difference between the time series calculated using steps h and h/2, in the case of aperiodic solution.In Section 3 we present the criterium of convergence based on the comparison of the attractors corresponding to the time series calculated using h and h/2.We then apply it to the nonlinear model presented in Section 2, and show that the probability of the attractors being the same increases monotonically as the step of integration h decreases.The comparison of the attractors is made using 1) the method of correlation integral, and 2) the method of statistical distance of probability distributions (divergence measure).Final comments are presented in Section 4.

The Multilooped Negative Feedback Equation with Three Delays
The DDE model, proposed by Glass et al. [5] for modelling physiological control system, consists of N loops described by variables x i , the variable of primary interest being the average The x i are controlled by feedback where τ i is the delay associated to the loop corresponding to the variable x i , and The function F i is a monotonically decreasing function, the parameters n and θ i governing the steepness and the threshold, respectively.The equation for the variable x is easily obtained by summing the equations for variables x i (2) and dividing by N: The initial condition will be a continuous function, We shall consider the case of N = 3 that has been studied in detail by Glass and Malta [6] and Bastos de Figueiredo et al. [7] with the following parameter values: The initial condition used was It was shown that by increasing n (increasing the steepness of F i (x) (3)) a period-doubling cascade was obtained.
For n = 45 the solution had period 2, and for n = 75 the solution was shown to be aperiodic.The numerical integration of Equation ( 4) is performed using the three-step integration method of Gear.This method has been tested by Malta and Teles (2000) [8], and was found to be more efficient than the 4th order Runge-Kutta method for the above equation.The step h is such that the delays (τ i = k i h, with k i integer.Given an initial condition, the numerical integration of equation ( 4), with step h, after discarding the transient (t ≤ t 0 ), produces the stationary series x h (t 0 + mh), m =1, 2,  , M. Comparison of the stationary series obtained for steps h and h/2 is made by calculating the following quantity For a given precision ε, the usual convergence criterium is that convergence is achieved if If the parameter values used correspond to a periodic or quasi periodic solution, the quantity (6) can be kept within the required precision for any time interval (M as large as possible).This is illustrated in Figures 1 and 2 that display the time series difference ∆ h for n = 45 (period 2) and n = 75 (aperiodic), respectively.
Figure 1 shows that the condition ( 7) is a good convergence criterium in the case of periodic solution, while Figure 2 shows that, in the case of aperiodic solution, after some time the condition (7) is no longer satisfied (this has also been pointed out in [1][2][3]).The smaller h is, the longer it takes for this to take place.This amplification of ∆ h as time evolves could be due to lack of stationarity of the time series (transient not removed completely).To check the stationarity of the time series we  have used the correlation integral.Grassberger and Procaccia [9] have observed that the probability that two points of the same attractor fall inside a box of size r may be approximated by the correlation integral that gives the probability that two points of an attractor are separated by a distance ≤ r.The correlation integral is related to the correlation dimension (that is an invariant) of an attractor, thus it will not depend on the time interval considered, unless the series is not stationary.The correlation integral is defined as where m is the dimension of the embedding space,  is the Heaviside function, and x i (m) is the m-dimensional vector constructed from the time series, with time lag p: We have used m = 2 and p = max(τ i ) = 2.00 to calculate the correlation integral (8) for points of the series contained in the time interval [kδt 1 , kδt 1 + δt 2 ], (k = 1, 2, 3,  ).This correlation integral will be denoted by C(2, r, k), and stationarity will be achieved at t = kδt 1 if increasing k causes no significant change.We have taken δt 1 = 25 and δt 2 = 200, and in the Figure 3 we show the curves

Comparison of Attractors
Since the chaotic attractor is known to be invariant, we propose that the criterium of convergence be based on the comparison of the attractors associated with the time series calculated using h and h/2 as integration step.We shall present two procedures for comparing the attractors, one based on the Kolmogorov-Smirnov test [10,11], and another one based on the calculation of the attractors divergence measure proposed by Diambra (2001) [12].It should be stressed that given any pair of time series, these procedures can be used for comparing them.

Kolmogorov-Smirnov Test
The Kolmogorov-Smirnov test provides the probability that two sets of measurements correspond to the same , , , N

  
 be two sets of measurements of the same observable.The corresponding cumulative distribution functions (not exceeding a given value η) are given by The Kolmogorov-Smirnov statistics is defined as Smirnov [11] has demonstrated that, given λ > 0, the probability gives the significance of acceptance of the null hypothesis that the two sets of data derive from the same distribution.
Albano, Rapp and Passamante [13] proposed the use of the correlation integral for calculating the Kolmogorov-Smirnov statistics to be used for comparing the attractors.The correlation for the 2-dimensional embedding of the attractor corresponding to the time series where the integer k 3 = τ 3 /h, and (2)   , i h x is the ith vector of the two-dimensional embedding of the attractor: , .
The Kolmogorov-Smirnov statistics related to the correlation integral is given by We have calculated C h (2, r) for r in the interval [0.01, 1.0] (N 1 = N 2 = 1000), using the time series resulting from the numerical integration of the three delay equation with h = 0.01/2 k , k = 0, 1, 2, •••, 9 (with parameter values (5), t 0 = 100).The results are displayed in Figure 4(a) and Figure 4(c) for n = 75, and n = 45, respectively.In Figures 4(b) and 4(d) we display the corresponding Q KS , calculated using the statistics (15).
Q KS , gives the probability that the attractors related to the series calculated with time steps h and h/2 are the same.Figure 4(b) shows that for the chaotic case (n = 75) this probability gets close to 100% for h ≤ 0.01/2 5  (probability is greater than 99.5%) while for the period-2 case (n = 45) (Figure 4(d)) this occurs for h = 0.01 (probability 100%).In the Figure 5 we display ∆h (6) for the chaotic case with h = 0.01/2 8 , and we can see the amplification of ∆ h for t > 300.

Convergence Criterium Based on the Statistical Distance
This method uses the divergence measure proposed by  Diambra [12] for comparing attractors.The distribution of points of the attractor in phase space is analogous to a probability distribution, therefore the statistical distance of probability distributions can be used to evaluate the attractors similarity.The statistical distance of the distribution probabilities p and g (it measures the amount of information associated with p relative to g) is given by (see [14]) where f is a functional operator for the system entropy.Diambra [12] uses the generalized entropy function [15] given by Substituting (17) in Equation ( 16) (q = 2) gives To calculate (18) numerically, it is necessary to discretize it.The two-dimensional phase space attractor region is divided in B squares of side c (see illustration in the Figure 6), and then the points of the attractors inside each square cell is counted (p ij and g ij , respectively).
Normalizing the distribution of points, In the Figure 7 we exhibit the statistical distance D 2 (p: g) where p (g) corresponds to the attractor associated with the time series calculated with step h (h/2).The number of boxes used was B = 2000 (c = 0.0005).Comparison of Figures 4(b), (d) and 7(a), (b) shows that the convergence of D 2 (p:g) to zero is slower than the convergence of Q KS to 100%.According to the convergence criterium based on the attractors statistical distance, if we require D 2 (p:g) ~ 0, in the period-2 case (n = 45) convergence is obtained for h = 0.01/2 2 , compared to h = 0.01 of the Kolmogorov-Smirnov test for which Q KS is 100%.In the chaotic case (n = 75), D 2 (p:g) ~ 0 for h = 0.01/2 8 , and Q KS is 99.5% for h = 0.01/2 5 .

Final Comments
It is well known that if N = 1 Equation ( 4) does not exhibit complex dynamics.In the case of N ≥ 3, chaotic solutions were found in the case of the feeedback function F i being piecewise constant [5] or in the case of the smooth sigmoidal function (3) [6].In the two feedback loops case [7] (N = 2) the comparison of the attractors corresponding to the time series for integration step h and h/2 was fundamental to establish the existence of chaotic solutions.One of the chaotic solutions was found for n = 45 with the parameters value.
To confirm that this time series corresponded to a chaotic solution, we constructed the bifurcation diagram, displayed in the Figure 8, by varying n in the interval [39,46] (computing 300 points in the Poincaré section for each value of n).
The bifurcation diagram in Figure 8 shows the characteristic route of period-doubling bifurcation cascade to chaos, indicating the existence of chaos in the two-loops negative feedback system studied by Bastos de Figueiredo et al. [7].Our numerical investigation has shown that the criterium of numerical convergence based on the comparison of attractors corresponding to the time series calculated using time steps h and h/2 (using either (15) or ( 19)) works well for all types of solutions: periodic, quasi-periodic or aperiodic.Figures 4 (Kolmogorov-Smirnov test) and 7 (statistical distance) show that the probability of the attractors being the same increases (monotonically) as the step h is decreased.
In the period-2 case (n = 45) the Kolmogorov-Smirnov test (12) (using statistics (15)) indicates convergence for h = 0.01 (see Figures 4(c) and 4(d)), and the test (19) of attractors statistical distance gives the same result if we require that D 2 (p:g) < 2.2 (see Figure 7(b)).If we require that D 2 (p:g) < 1.0, it is necessary to use time step h = 0.005.If we require D 2 (p:g) < 0.1, it is necessary to use h = 0.0025.In the case of periodic or quasi-periodic solutions it is possible to use (7) as convergence criterium for the numerical integration.However, in the chao-tic case (n = 75), due to the dynamical properties of this type of solution, the difference (6) blows up as shown in Figures 2 and 5. Nevertheless, the comparison of the phase space chaotic attractors corresponding to the time series calculated with time steps h and h/2 indicates that the probability of the corresponding chaotic attractors being the same increases monotonically as the time step h decreases (see Figures 4(a  and 7(a)).Therefore, the fact that the time series calculated by Teixeira et al. [1] satisfies condition (7) during a certain time interval does not necessarily mean that convergence has been achieved.It could well be that their numerical integration has not converged as already pointed out by Yao and Hughes [2].
Our results indicate that D 2 (p:g) (Equation ( 19)) is more sensitive than Q KS (Equation ( 12)): the value of h that gives probability ≥ 99.5% of the attractors being the same corresponds to D 2 (p:g) < 1.0 (see Figures 4(b Our method is very useful for detecting deviations from a given reference time series that corresponds to a behaviour that is desirable in a process being monitored as in Stienstra et al. [16].

Figure 3 .
Figure 3. Curves log C(2, r, k) × r, k = 1, 2, •••, 6 for the time series with n = 75 calculated with step h = 0.0001 (same time series used in the Figure 2(c)).logC(2, r, k) × r for k = 1, 2, •••, 6. Stationarity is achieved for k = 4 (t = 100) and this rules out the possibility of the amplification of ∆ h observed in Figure 2 being due to nonstationarity effects.This indicates that the amplification of ∆ h observed in Figure 2 for the case n = 75, is a consequence of the system dynamics.

Figure 6 .
Figure 6.(a) Illustrating the square division of the plane region occupied by the attractor; (b) Detail of the insert in (a) to illustrate the counting of points in a square cell.

Figure 8 .
Figure 8. Bifurcation diagram for two feedback loops with parameters (20).The Poincaré section contains 300 points for each value of n; We considered 100 values of n in the interval [39,46].