On the lower bound error for discrete maps using associative property

ABSTRACT This paper introduces a class of pseudo-orbits which guarantees the same lower bound error (LBE) for two different natural interval extensions of discrete maps. In previous work, the LBE was investigated along with a simple technique to evaluate numerical accuracy of free-run simulations of polynomial NARMAX or similar discrete maps. Here we prove that it is possible to calculate the LBE for two pseudo-orbits, extending so the results of previous work in which the LBE is valid for only one of the two pseudo-orbits. The main application of this technique is to provide a simple estimation of the LBE. We illustrate our approach with the Logistic Map and Hénon Map. Using double precision, our results show that we ought simulate the Logistic Map and Hénon Map with less than 100 iterations, which is, for instance, far less than the number usually considered as transient to build bifurcation diagrams.


Introduction
Discrete maps, also called recursive functions, are used to account for many types of systems (Feigenbaum, 1978), in particular nonlinear dynamics systems (Ferreira, Nepomuceno, & Cerqueira, 2006;Hammel, Yorke, & Grebogi, 1987;Lozi, 2013;Pereira, Kurcbart, & Nepomuceno, 2005). To investigate these functions in nonlinear dynamical systems, numerical computation plays a key role (Galias, 2013;Lozi, 2013), and a myriad of numerical experiments (Peck, 2004) have been performed since the work of Lorenz (1963) to understand the behaviour of such nonlinear dynamical systems (Hammel et al., 1987). They are used as direct tool to model many types of dynamical systems (Feigenbaum, 1978) and many types of discretization schemes for continuous nonlinear dynamical systems make use of discrete maps (Letellier, Mendes, & Mickens, 2007).
Numerous researchers are confident in numerical solutions of nonlinear dynamical systems. All sort of dynamics have been classified, such as chaotic (Ott, 2002) or intermittent (Hirsch, Huberman, & Scalapino, 1982) dynamics, whose conclusions are based on numerical simulations. However, there exist doubts about these dynamics. The work of Ford (1986), with his intriguing title: 'Chaos: Solving the unsolvable, predicting the unpredictable', is one of the first attempts to articulate some discussions on the discovered evidences by computer simulations. Later Lozi (2013) asks if 'In the simple case of a dynamic discrete CONTACT E.G. Nepomuceno nepomuceno@ufsj.edu.br system (of Hénon map) there are doubts as to the nature of the computational results: long unstable pseudo-orbits or strange attractors?'. Similarly, Galias (2013) expresses the importance of developing methods to prove the existence of chaotic attractors; the existence for Lorenz system have been carried successfully (Tucker, 1999a(Tucker, , 1999b, while existence for Chua's circuit is still open. Nevertheless, it is important to state, that the proof in Tucker (1999aTucker ( , 1999b is based on 'combination of normal form theory and rigorous numerics'. As any discretization scheme may be seen as a recursive function, difficulties have been reported to keep reliability on the results of numerical simulation (Amigó, Kocarev, & Szczepanski, 2007;Berthé, 2012;Elabbasy, Elsadany, & Zhang, 2014;Letellier & Mendes, 2005;Lima, Claro, Ribeiro, Xavier, & López-Castillo, 2013;Mendes & Letellier, 2004;Nepomuceno & Mendes, 2017;Rodríguez & Barrio, 2012;Teixeira, Reynolds, & Judd, 2007). Hence, the error propagation in computer simulations requires a careful analysis (Lozi, 2013). Based on a preliminary work about convergence of recursive functions on a computer (Nepomuceno, 2014), Nepomuceno and Martins (2016) proposed a method to calculate the lower bound error (LBE) for free-run simulations of recursive functions, with particular attention to polynomial NARMAX (Billings, 2013). In that work, a simulation is performed by means of two natural interval extensions, derived from mathematical properties, such as commutativity, distributivity or associativity. These two natural interval extensions produce two different pseudo-orbits. The LBE has already been shown useful in works such as Mendes and Nepomuceno (2016), where the authors have developed a fast and robust method to calculate the positive largest Lyapunov exponent, and Nepomuceno and Mendes (2017), where the authors analyse the consequences of error propagation in discretization schemes to solve nonlinear differential equations.
In Nepomuceno and Martins (2016), the results on the LBE have been proved to be valid for only one of the two pseudo-orbits. In this work, we present a step further in the consolidation of the LBE. We found a class of natural interval extensions, which derives pseudoorbits with equal error bound when performing the rules of float point based on the IEEE 754-2008 standard (Goldberg, 1991;Institute of Electrical and Electronics Engineers (IEEE), 2008;Overton, 2001). The investigation of different behaviours from different expressions for the Logistic Map has also been studied in Yabuki and Tsuchiya (2013). We may think that an arithmetical operation in computer is not exact but exact within a factor of (1 + ε), where ε is the machine epsilon (Overton, 2001). The associative property of multiplication allows us to produce different pseudo-orbits, while keeping the same error bound. In this way, we prove a theorem that states the lower bound error for this kind of natural interval extension, called arithmetic interval extension, and for two pseudo-orbits. We also focus our attention on discrete maps, using the Logistic Map and Hénon Map as examples of application. Another contribution regards the notation of the LBE to clarify a connection with absolute error and relative error in floating-point arithmetic.
The remainder is organized as follows. First, preliminary concepts are presented, then we propose a new theorem on the lower bound of two pseudo-orbits. We apply the methodology for the Logistic Map and Hénon Map, and the results are validated by means of a comparison with an arbitrary precision.

Preliminary concepts
In this section, we present the concepts of natural interval extensions, orbits and pseudo-orbits and a brief overview of floating point arithmetic as specified by IEEE floating point standard 754-2008 Institute of Electrical and Electronics Engineers (IEEE) (2008). These concepts are used to establish the main contribution of this paper, as the lower bound error is derived from two different natural interval extensions and their pseudo-orbits under the arithmetic carried out by IEEE 754-2008.

Natural interval extension
An interval is set of real numbers such that any number that lies between two numbers in the set is also included in the set. As for notation, an interval X is denoted [X,X], i.e. X = {x : X ≤ x ≤X}. In an degenerated interval, we have X =X and such an interval amounts to a real number x = X =X. In this context, Moore, Kearfott, and Cloud (2009) give the definition of natural interval extensions of a function.  An interval extension of f is thus an interval valued function which has real values when the arguments are real (degenerate intervals) and coincides with f in this case. The natural interval extension is achieved by combining the function rule f (x) with the equivalents of the basic arithmetic and elementary functions.
There are natural interval extensions that are equivalent in terms of intervals. This leads us to the following definition.  Here, only G(X) and H(X) are equivalent interval extensions.
It is imperative to stress that from the perspective of classic mathematics on real numbers, diverse syntactic partners of two equivalent natural interval extensions yield a similar outcome. However, the outcome might be distinctive in interval arithmetic, specifically when the arithmetic operations are performed on computers using float point operations.

IEEE floating point representation
Finite IEEE floating point numbers can be expressed as in the following form Overton (2001): where p is the precision of the floating point system; b 0 to b p−1 is the binary expansion of significand S; E is the exponent.
Example 2.4: Let us consider the single precision, that is, 32 bits, where 1 bit is for the signal, 8 is for the exponent (represented in a bias format) and 23 is reserved for the significand. The decimal 0.1 is approximated by a binary string 00111101110011001100110011001101 where the first 0 represents the positive signal of Equation (4), the following binary string 01111011 represents the biased exponential, E = 123 − 127 = −4, and finally the significand is given by 1001100110011001100 1101. Recall that b 0 = 1 is not represented and it is denominated hidden bit. Using Equation (4), we have: Let x ∈ R be a real in a computer using the IEEE standard. We define x − to be the floating point number nearest to x that is less than or equal to x, and define x + to be the floating point number nearest to x that is greater than or equal to x. The IEEE standard defines the correctly rounded value of x, denoted here as round(x). If x is a floating point number, then round(x) = x. Otherwise, round(x) depends on which rounding modes is in effect, as follows: whichever is nearer to x.
The absolute rounding error associated with x is defined as and the relative rounding error is where Therefore, for a binary floating point system with precision p, we have: for some δ and κ satisfying where ε = 2 −p and γ = ε × 2 E if round mode in effect is rounding to nearest and considering x within the normalized range (Overton, 2001).
A key feature of the IEEE standard is that it requires correctly rounded arithmetic operations. Let x and y denote float point numbers, let +, −, × and ÷ the standard arithmetic operations, and let ⊕, , ⊗ and the equivalent operations as they are actually implemented on computers. The rule is as follows: if x and y are floating point numbers, then (Overton, 2001) x We can write that the result of an arithmetic operation relies in the interval given by the multiplication of (1 − ε) and (1 + ε). Taking the addition as example, we have which is equivalent to say We may continue likewise for subtraction, multiplication and division. Henceforth we may revise Equations (14) to (17) using a description of intervals as follows.
It is known that certain properties of classic arithmetic, such as associativity of addition, do not hold in the context of floating point arithmetic (Institute of Electrical and Electronics Engineers (IEEE), 2008; Overton, 2001). Consequently, two mathematical equivalent sequences of arithmetic operations can yield two different results, as in the case of orbits and pseudo-orbits of maps.

Orbits and pseudo-orbits
Let n ∈ N, a metric space M ⊂ R, the relation where f : M → M, is a recursive function or a map of a state space into itself and x n denotes the state at the discrete time n. The sequence {x n } obtained by iterating Equation (24) starting from an initial condition x 0 is called the orbit of x 0 (Gilmore, Lefranc, & Tufillaro, 2012).

Definition 2.5 (Orbit):
Given a map x n+1 = f (x n ), an orbit of the map is a sequence of values of the map, represented by The calculation of an orbit is usually carried out by a finite-precision computer, leading to a pseudo-orbit. A pseudo-orbit of a map is an approximation of a mathematical orbit in a specific hardware or software. There is no unique pseudo-orbit, as there are different hardware, software and numerical precision standards, such as IEEE 754-2008, which may yield different outputs for each extension interval. Here it is important to emphasize that some pseudo-orbits may produce better results than others. For this reason, there has been different investigations to obtain the most precise pseudo-orbits; for example by manipulating natural interval extensions, see e.g. Horner's method (Lambers & Sumner, 2016;Muller et al., 2010;Rodríguez & Barrio, 2012;Stahl, 1997), or by using specific interval extensions, such as mean value forms, which usually produce narrower widths (Caprani & Madsen, 1980;Dymowa, 2011;Rall, 1983). In this paper we do not focus on the reduction of the propagation error or the width of the error bounds.

Definition 2.6 (Pseudo-orbit): Given a map
with the absolute error given by where The relative error is given by where Thus, for a floating point system, we have: for some δ i,n and κ i,n satisfying It is clear that for n = 0, Equations (29) and (32) are equivalent to Equations (10) and (13), that is, δ i,0 = δ and ε i,0 = ε, wherein the index i could be dropped as there is no pseudo-orbit in action, but the process of rounding. It is important to stress out that the IEEE standard does not mention any general approach to calculate ε i,n or γ i,n for n ≥ 1. This has been the focus of many works, particularly those related to interval arithmetic (Moore et al., 2009).
Here we may define an interval associated with each value of a pseudo-orbit, I i,n such that Let {x a,n } and {x b,n } be two pseudo-orbits, here we denote the LBE as ,n , where is the LBE calculated at a point n using a set of pseudo-orbits indicated by . In this paper, is defined with just two interval natural extensions, denoted by = [{x a,n }, {x b,n }]. A preliminary result in Nepomuceno and Martins (2016) indicates that γ a,n ≥ ,n or γ b,n ≥ ,n . In this paper we provide a class of pseudo-orbit which allows to state that γ a,n = γ b,n ≥ ,n . It is important to stress that γ a,n = γ b,n does not imply κ a,n = κ b,n . In other words, the class of natural interval extensions we are going to introduce here presents the equivalent bounds for the absolute error, but the absolute error for a specific instant n may be different (see Figure 1(a)).

Lower bound error
The key point of this paper is the analysis of a specific class of natural interval extensions based on the associative property of multiplication. This leads us to the following definition of equivalent arithmetic intervals, cf. Definition 2.2.

Definition 3.1 (Arithmetic interval extension):
Two equivalent natural interval extensions G and H of a function f are arithmetic interval equivalent if they produce the same resulting interval following Equation (20) to (23).
Manipulations on terms of a function by means of associative property of multiplication will keep the same bounds of error. We indicate this key aspect in the following Lemma. Proof: Let us consider a generic mathematical operation a × b × c, its float point equivalent a ⊗ b ⊗ c, and the machine epsilon ε. We have  (1 − x)) are arithmetic interval equivalent, since we just apply the associative property abc = a(bc). Indeed the resulting interval for G(x) is as follows whereas for H(x) we have Thus As we can see, the intervals in (35) and (36) are the same. Although, the computational outcome for G(x) and H(x) is likely unique, they present the same error bounds due to round off.
Another important point to emphasize in this case is that f (·) is a map, and the error produced by round off is the same for G(x) and H(x) for any iteration of a map. This is the key aspect that we are going to apply to develop a stronger result on lower bound error.
Arithmetic equivalent extensions can be used to refine the following theorem proved in Nepomuceno and Martins (2016). Theorem 3.4 is limited in that at least one of the two pseudo-orbits has an error bound equal or greater than the lower bound error ,n . We refine here this result for a class of natural interval extensions, namely the class of equivalent arithmetic interval extensions (Definition 3.1), for which error bound are the same. So a stronger result for the lower bound error is given by the following theorem. Proof: The two interval natural extensions are arithmetic interval extensions. By Lemma 3.2, it holds that ε a,0 = ε b,0 = ε. As γ = ε × 2 E , we also have γ a,0 = γ b,0 . By induction, this result can be further extended to γ a,n = γ b,n . Now we need to show the definition of ,n implies γ a,n = γ b,n ≥ ,n . It is known that I a,n ∩ I b,n = ∅ if and only ifx b,n − γ b,n ≤x a,n + γ a,n andx a,n − γ a,n ≤x b,n + γ b,n .
Theorem 3.4 establishes that at least one of the two pseudo-orbits must have a bound error greater or equal to the lower bound error. Theorem 3.5 extends this result and states that two pseudo-orbits, which are derived from two arithmetic interval extensions, have their bound error greater than the lower bound error. This difference has a practical meaning: if this lower bound error is greater than the required precision, we must interrupt our simulation, as both pseudo-orbits are no longer reliable. Table 1. Summary of the main notations for error of floating point representation and computation of a pseudo-orbit. For n = 0, we have κ i,0 = κ, γ i,0 = γ , δ i,0 = δ and ε i,0 = ε, where the index i can be just dropped as we have only the round mode to the nearest in effect.

Definition
Floating Point Pseudo-orbit

Summary of notation
Notation is recapitulated in Table 1. We can notice that for n = 0, the floating point is the initial condition of the pseudo-orbit, that is, κ i,0 = κ, γ i,0 = γ , δ i,0 = δ and ε i,0 = ε, where the index i can be dropped without any damage.
The main contribution of this paper is to state that γ a,n = γ b,n ≥ ,n .

Numerical experiments
Numerical experiments were conducted by means of an arbitrary precision technique. We used the software Maxima to compute a pseudo-orbit with a precision much higher than common computations. We used 1000 digits in Maxima and compare our results with 16 digits precision in Scilab 5.5.2. We are interested in noticing whether the error calculated by the absolute difference between Scilab and Maxima results are equal or greater than the lower bound error. We illustrate our results with the Logistic Map and Hénon Map.

Logistic map
Let the logistic map (May, 1976) defined by with r = 3.9 and x 0 = 0.1. Let two arithmetic interval extensions be: and We compared the first ten iterations using double precision, i.e around 16 digits of precision, with another simulation with 1000 digits of precision. In both cases, our lower bound is correct, that is, the module of the difference between the Maxima result and the two natural interval extensions G(X n ) and H(X n ) are equal or greater than the lower bound error. Table 2 shows these ten first values. Appendix presents the code executed by Maxima. The code line with true as result indicates that the error of each pseudo-orbit is greater than the lower bound. The value obtained by Maxima after the 10th iteration is given by 537948645688298809066288385414421252457  904090752732170600455542715349546559676  182749500840788738608085567373632215672  744170192816742919433742996537129905709  746766127038898588764092437192223368127  212213527694597407931652212238049148048  511344308788342883935422845687587707841  315189287670577453218653237764283713863  665017413010842537037054504348421366464  616846822403037740290136707985403129494  055555373233114708689612217293241135643  356324713734411729802482278651521288550  738883540304381914688553216548441436963  390397750773329223657066700828287398478  411010681662922523469885906359952736570  361460429260964913757169005885846182835  324439815759501634227118527512730337696  588933677857091039906043305882918361403  093098361698215523311209139941041229444  005656944236824812007534572812618929386  310410242541746697089279560798912764244  285365756004530014850331177395907198257  858824610321778535193528841461304953684  816117265437059953849384086883856191280  584004759890502596051919722861142653632  405227843041774850104642355298188103284  865207401   10 1023 Since ,10 = 2.275957200481570908 · 10 −15 , we can see that the lower bound error is valid by proceeding with the following basic operations: κ G,10 = |x 10 − G(X 10 )| = 1.110223024625157 · 10 −14 > ,10 and κ G,10 = |x 10 − H(X 10 )| = 6.550315845288424 · 10 −15 > ,10 .
Figure 1(a) shows the simulation of the two pseudoorbits (39) and (40) from iteration 70-101. Figure 1(b) shows the evolution of lower bound error in logarithm scale. When n = 83, the lower bound error is superior than 0.1. It means that all the digits have been lost significance. From this results, it is clear that we cannot trust the simulation of Logistic Map with r = 3.9 and x 0 = 0.1 using Scilab and double precision beyond n = 83, which is a contrary view of works, such as Ott (2002), that suggests a transient of 500 points to build the bifurcation diagram.

Hénon map
Let the Hénon map (Hénon, 1976) be defined by where a = 1.4, b = 0.3 with x 0 = y 0 = 0.1. Let two arithmetic interval extensions be: and Observe that in this case only associative property is applied, and Equation (44) is kept the same for both extensions. Following Lemma 3.2 it is clear that Z G = Z H , in a similar way as Equation (37). Figure 2(a) shows the simulation of the two pseudoorbits (45) and (46) from iteration 70-101. Figure 2(b) shows the evolution of lower bound error in logarithm scale. When n = 92, the lower bound error is superior than 0.1. It means that all the digits have lost significance.
To check the validity of the lower bound error, we proceed in a similar way that the logistic map. The value Table 2. Simulation of the first ten values of the algorithm presented in Section A.1, where x n and y n are arithmetic interval extension extension G(x n ) and H(x n ), respectively, and presented in Equation (39) and (40) for the logistic map. The third column is the lower bound error of these two pseudo-orbits. obtained by Maxima after 10 iterates is given below.

Symbolic validation
We checked our results using ten iterations of symbolic computations by means of Maxima. In both cases, when n = 10 the error of the used pseudo-orbits appeared greater than the lower bound error. The code is presented in Section A.3 and A.4. A close look on Table 2 and 3 indicates some values equal to zero. One of the reasons is that the difference caused by two different pseudo-orbit may not be observed in the first few iterations, as the rounding mode may produce equal float point numbers for both. Table 3. Simulation of the first ten values of the algorithm presented in Section A.2, where x n and y n are arithmetic interval extensions G(x n ) and H(x n ) of Equations (45) and (46) for the Hénon map. The third column is the lower bound error.

Conservativeness of the result
Although it seems that the LBE presents a deep conservative position, that is, the LBE produces an overestimation of the error, there are some evidences that show the opposite. First, the numerical examples show the LBE quite close to a very rigorous computation of the precise error (see eqs. (41)- (47) and (42)- (48)). Besides, the fact that it has been possible to calculate the largest positive Lyapunov Exponent from the LBE in other work (see Mendes & Nepomuceno, 2016) is an evidence that computation is very close to the error propagation. As explained in Mendes and Nepomuceno (2016), starting two chaotic maps from the same initial condition, two different pseudo-orbits usually produce the same divergence of those observed due the sensitivity of initial conditions which are often observed for chaotic systems.

Conclusion
This paper investigates a technique to evaluate the LBE of discrete maps. We found a class of natural interval extensions that produce an equivalent bound error in recursive functions or discrete maps. This step allows us to compute the LBE for two pseudo-orbits, while previous results were applied for only one of the two pseudo-orbits. The focus of this work is on discrete maps. The relevance of the method relies on its simplicity to calculate the LBE and its multiple applications, for example in discretization schemes of differential equations, optimization techniques, polynomial NARMAX, Lyapunov Exponent computation to cite a few computational techniques. The propagation of the error is performed by means of a simple code. The results agree in general with the literature that shows an exponential growing of the error for chaotic systems (Adler, Kneusel, & Younger, 2001). For instance, Figure 2(b) is very similar with results presented in Figure 1 in Galias (2013). Both results indicate a maximum of around 100 iterations for Hénon Map using double precision before losing the significance of simulation, that is, when the error reaches the magnitude of the variable. The comparison of the proposed method and the technique described in Galias (2013) seems an interesting direction for further research.
Much of efforts and results in the the study of nonlinear dynamics relies on the use of numerical simulations, to validate (Billings, 2013), to verify accuracy (Hammel et al., 1987) and even to make mathematical proofs of existence of strange attractors (Tucker, 1999a(Tucker, , 1999b. Therefore our method may be useful as an additional tool to falsify non-reliable simulations. We verified our results for the first ten iterates using an arbitrary precision by means of the software Maxima. As mentioned by Adler et al. (2001), a mathematical package cannot be used as a black box. Accordingly we would like to stress the importance of having access to the code and full details of numerical simulations, otherwise, it may be impossible to reproduce results from only given equations.
We hope that in a near future we can investigate other classes of pseudo-orbits with similar results. Another interesting research is the identification of classes of pseudo-orbits for which the LBE grows or tends to zero. In the same line, we intend to investigate LBE methods to reduce the propagation of error in recursive functions, as well as to narrow the 1 width of interval arithmetic methods.

Note
1. Scilab is a free software which may be downloaded from www.scilab.org.