Note on improvement precision of recursive function simulation in floating point standard

DINCON 2017 CONFER ˆ ENCIA BRASILEIRA DE DIN ˆ AMICA, CONTR OLE E APLICAC ¸ ˜ OES 30 de outubro a 01 de no v em bro de 2017 – S˜ ao Jos ´ e do Rio Preto/SP Note on impro v emen t precision of recursiv e function sim ulation in floating p oin t standard Melanie Ro drigues e Silv a 1 Group of Con trol and Mo deling, UFSJ, S˜ ao Jo˜ ao del-Rei, MG, Brazil Eriv elton Geraldo Nep om uceno 2 Departmen t of Electrical Engineering, UFSJ, S˜ ao Jo˜ ao del-Rei, MG, Brazil Samir Angelo Milani Martins 3 Departmen t of Electrical Engineering, UFSJ, S˜ ao Jo˜ ao del-Rei, MG, Brazil Abstract . An improv e men t on precision of recursive function simulation in IEEE floating p oin t standard is presen ted. It is shown that the av erage of rounding to w ards negative infinite and rounding tow ards positive infinite yields a b etter result than the usual standard rounding to the nearest in the simulation of recursive functions. In general, the metho d impro v es one digit of precision and it has also b een useful to av oid divergence from a correct stationary regime in the logistic map. Numerical studies are presen ted to illustrate the metho d. Keyw ords . IEEE floating p oint, rounding mo de, numerical sim ulation, recursive functions and Logistic map. 1 In tro duction The mathematical implementation of functions and systems in computers enables the scien tific in v estigation in many areas. According to [19], interferences and noises of real systems are not incorp orated in to the simulations, though it is known that it is hardware as a set of real physical systems then the noises are alwa ys presen t. How ever, limitations of soft w are and hardw are are obstacles to the reliability of results [14, 15, 18, 22, 24, 26]. Nev ertheless, alternatives hav e b een developed to circumv ent hardw are shortcomings and impro v e results, suc h as in [3, 6, 7, 17, 25]. These alternativ es that restricts the effects of hardware limitations are usually called as rigorous computing, whic h is based on the dev elopmen t of refined metho ds for the implementation of algorithms [11, 14]. Chaotic systems suc h as discrete maps and problems in v olving recursive functions ha ve a differentiated degree of complexit y , since the curren t iteration dep ends on the previous 1 me ro drigues silv a@hotmail.com 2 nep om uceno@ufsj.edu.br 3 martins@ufsj.edu.br 2 ones [8]. The study of c haotic systems still receive atten tion and many of scien tific conclu- sions in this area relies up on computer sim ulation. Chaotic b ehaviour has b een asso ciated to many real applications, from electronic circuits [2] to ultrasonic cutting system, whic h dynamic b eha vior is analyzed using a t w o-degree-of freedom Duffing oscillator mo del [21]. Moreo v er, these iterations are usually ev aluated as n umerical ill-conditioned [5]. In gene- ral, the simulation pro cess of these types of functions inv olves errors and its reliability has b een questioned in many w orks [14, 15]. In this sense, it is w ell kno wn that in a floating p oin t environmen t, a computer simulation of a c haotic system ma y b e limited to a short n um ber of iterates [26]. T o deal with this problem, man y to ols hav e b een in v estigated. W e may summarize these to ols in tw o categories. The first is fo cused on hardware p erfor- mance using parallel and cluster computation [16]. The second ma jor categories is dev oted to improv e the algorithm, using in terv al analysis, significance arithmetic or noisy-mo de computation [11]. In this pap er, w e fo cus our atten tion on the noisy-mode computation, but instead of adding pseudo-random digits, as suggested in [11], w e exploit the random nature of rounding error [4, 20]. Considering that the rounding errors are uniformly dis- tributed, we presen t the background necessary for the understanding of the guidelines and the algorithm of a metho d that improv es the precision of the computational sim ulation of recursive functions. In order to pro v e the efficiency of the proposed strategy , we bring examples of the logistic map case. 2 Bac kground This section contains some definitions ab out recursive functions, orbits, and interv al pseudo-orbits. Let n ∈ N , b e a metric space such that I ⊆ R , and f : I → R . Thus, recursiv e function can b e defined as [18]: x n = f ( x n − 1 ) , (1) or b y comp osite functions: x n = f 1 ( x n − 1 ) = f 2 ( x n − 2 ) = . . . = f n ( x 0 ) , (2) where f : I → I is a recursiv e function or a map of a state space into itself and x n denotes the state at the discrete time n . As suggested by [13], the sequence x n obtained b y iterating Eq. (1) starting from an initial condition x 0 is called the orbit of x 0 . Definition 1. A n orbit is a se quenc e of values of a map, r epr esente d by x n = [ x 0 ,x 1 ,...,x n ] . The Definition 2. is suggested by [24] and highligh ts the presence of the difference b et w een real and computational pseudo-orbits. Definition 2. L et i ∈ N r epr esents a pseudo-orbit, which is define d by an initial c ondition and a c ombination of softwar e and har dwar e. A pseudo-orbit is an appr oximation of an orbit and c an b e r epr esente d as { ˆ x i,n } = [ ˆ x i, 0 , ˆ x i, 1 , . . . , ˆ x i,n ] (3) 3 such that, | x n − ˆ x i,n |≤ ξ i,n (4) wher e ξ i,n ∈ R is the limit of err or and ξ i,n ≥ 0 . There is not an unique pseudo-orbit, as there are different hardw are, soft w are, numeri- cal precision standard and discretization schemes, which can pro duce differen t results and consequen tly differen t errors ξ i,n . 3 Metho d for error reduction It is reasonable to assume that round-off errors are uniformly distributed [4, 9, 20, 23]. According to IEEE 754-2008 [10], the arithmetic op erations are also rounded. In this sense, w e also assume here that a finite set of arithmetic function presen ts an error randomly distributed. Rounding to the nearest is the most usual for computational arithmetic. How ever, there are soft w are that enables the user to set the rounding mo de tow ards p ositive infinite or negativ e infinite. With these considerations in mind, w e establish the main con tribution of this letter in the follo wing Lemma. Lemma 1. L et ˆ x − i,n and ˆ x + i,n b e the c alculate d value by r ound towar ds p ositive infinite and r ound towar ds ne gative infinite, r esp e ctively. The arithmetic aver age given by ˆ x j,n = ˆ x + i,n + ˆ x − i,n 2 (5) such as ˆ x j,n = f ( ˆ x j,n − 1 ) + δ j,n , pr esents an r ound-off err or smal ler than than the r ound-off err or due the r ound to ne ar est as n → ∞ , ther efor e δ j,n < ξ i,n . Pr o of. Assuming Eq. (5), considering that ˆ x − i,n = f ( ˆ x j,n − 1 )+ δ − i,n and ˆ x + i,n = f ( ˆ x j,n − 1 )+ δ + i,n , then: ˆ x j,n = f ( ˆ x j,n − 1 ) + δ − i,n + f ( ˆ x j,n − 1 ) + δ + i,n 2 ⇒ ˆ x j,n = 2 f ( ˆ x j,n − 1 ) + δ − i,n + δ + i,n 2 (6) and ˆ x j,n = f ( ˆ x j,n − 1 ) + δ − i,n + δ + i,n 2 (7) But, as w e hav e b een considered δ − i,n and δ + i,n uniformly distributed and it is w ell kno wn that av eraging a random v ariable leads to a reduction of the noise p o w er in n , which is in this case is 2. And that completes the pro of.  The Algorithm 1 presents a pseudo-co de based on the Lemma 1. The round to the nearest is the standard mo de. The operators r ound − and r ound + stands for round tow ards negativ e infinite and round tow ards p ositive infinite. Where there is no indication of the round mo de, one should consider as the round to the nearest. In this pap er, w e implemen ted this algorithm in Matlab R2016a . The comparison is made with the high precision v alues pro vided b y the VP A to olb ox of Matlab R2016a . 4 Algorithm 1: Simulation of recursive function based on Lemma 1. 1 ˆ x i, 0 ; %Initial c ondition 2 N ; %Numb er of iter ates 3 ˆ x − i, 0 ← r ound − ( ˆ x i, 0 ); 4 ˆ x + i, 0 ← r ound + ( ˆ x i, 0 ); 5 ˆ x j, 0 ← ( ˆ x − i, 0 + ˆ x + i, 0 ) / 2; %Aver age of the two r ound mo des 6 ˆ x − i, 0 ← ˆ x j, 0 ; 7 ˆ x + i, 0 ← ˆ x j, 0 ; 8 F or n ← 1 to N %Main lo op 9 ˆ x − i,n ← r ound − ( f ( ˆ x − j,n − 1 )); 10 ˆ x + i,n ← r ound + ( f ( ˆ x + j,n − 1 )); 11 ˆ x j,n ← ( ˆ x − i,n + ˆ x + i,n ) / 2; 12 ˆ x − i,n ← ˆ x j,n ; 13 ˆ x + i,n ← ˆ x j,n ; 14 EndF or 4 Numerical Examples This sections presen ts three n umerical examples using the metho d prop osed in this pap er. All the examples are based on the logistic map [1] x n +1 = r x n (1 − x n ) , (8) where the initial condition and the v alue of bifurcation parameter r are indicated in each example. Example 4.1. L o gistic map with x 0 = 1 / 3 . 9 and r = 3 . 9 . Only using the round to nearest, the pseudo-orbit presen ts a chaotic behaviour. How e- v er, using the Algorithm 1, it is observed that the logistic map function conv erges to a fixed point as shown in T able 1. As one can see, the ˆ x j,n reac hes a fixed p oint, whic h is the exact answ er for this example, as v erified in the follo wing equations: x 1 = 3 . 9(1 / 3 . 9)(1 − (1 / 3 . 9)) = (1 − 10 / 39) = 29 / 39 (9) x 2 = 3 . 9(29 / 39)(1 − 29 / 39) = 29 / 39 (10) and 29 / 39 ≈ 0 . 743589743589744 is a fixed p oint, whic h is the v alue that ˆ x j,n con v erges to. On the other hand, the pseudo-orbit ˆ x i,n div erges and presen ts a c haotic b ehaviour. Example 4.2. L o gistic map with r = 3 . 9 and x 0 = 0 . 01 . The Figure 1 sho ws the logarithm (base 10) of the error for the first 20 iterates of the logistic map, using the Algorithm 1 and a traditional metho d based on the round to the nearest. The Algorithm 1 starts with an error greater than the traditional metho d, but after few iterates the error becomes smaller. The error presented in Figure 1 has b een calculated b y means of VP A to olb ox of Matlab with 1000-digit precision. 5 T ab ela 1: Iterates of the logistic map with x 0 = 1 / 3 . 9 and r = 3 . 9. ˆ x i,n is follo wing the rounding to the nearest and ˆ x j,n is following the Algorithm 1. The v alues are shown in decimal and hexadecimal notation. Notice that the third iterate is equal to the second for the Algorithm 1. This can b e seen only in the hex format. n ˆ x i,n (hex) ˆ x i,n (dec) ˆ x j,n (hex) ˆ x j,n (dec) 1 3fd0690690690691 0.256410256410256 3fd0690690690691 0.256410256410256 2 3fe7cb7cb7cb7cb8 0.743589743589744 3fe7cb7cb7cb7cb8 0.743589743589744 3 3fe7cb7cb7cb7cb7 0.743589743589744 3fe7cb7cb7cb7cb8 0.743589743589744 4 3fe7cb7cb7cb7cb9 0.743589743589744 3fe7cb7cb7cb7cb8 0.743589743589744 5 3fe7cb7cb7cb7cb5 0.743589743589743 3fe7cb7cb7cb7cb8 0.743589743589744 6 3fe7cb7cb7cb7cb d 0.743589743589744 3fe7cb7cb7cb7cb8 0.743589743589744 7 3fe7cb7cb7cb7cae 0.743589743589743 3fe7cb7cb7cb7cb8 0.743589743589744 8 3fe7cb7cb7cb7ccb 0.743589743589746 3fe7cb7cb7cb7cb8 0.743589743589744 9 3fe7cb7cb7cb7c93 0.743589743589740 3fe7cb7cb7cb7cb8 0.743589743589744 10 3fe7cb7cb7cb7cfd 0.743589743589751 3fe7cb7cb7cb7cb8 0.743589743589744 0 5 10 15 20 25 n -18 -17 -16 -15 -14 -13 -12 -11 Error Figura 1: Propagation error for the logistic map with r = 3 . 9 and x 0 = 0 . 01. T raditional metho d (blue). Algorithm 1 (red) Example 4.3. L o gistic map with thr e e sets of initial c onditions and bifur c ation p ar ameter, given by ( x 0 , r ) = [(0 . 1 , 4 . 2); (0 . 2 , 4); (0 . 41 , 3 . 85)] . T able 2 shows results of ξ i,n and δ j,n for the first ten iterations of the logistic map case, considering these three differen t set of initial conditions and bifurcation parameter. It is p ossible to note a general decrease in the error that keeps along the iteration pro cess. W e arbitrary chosen these parameters, but we ha v e p erformed our tests with many other sets with similar results. 6 T ab ela 2: Results of ξ i,n and δ j,n of the logistic map case, considering different initial conditions, x 0 and r v alues. x 0 = 0 . 1 and r = 4 . 2 x 0 = 0 . 2 and r = 4 x 0 = 0 . 41 and r = 3 . 85 n ξ i,n δ j,n ξ i,n δ j,n ξ i,n δ j,n 1 9.33254 e-18 9.33254 e-18 4.57088 e-16 1.20226 e-16 3.63078 e-16 2.69153 e-17 2 3.98107 e-17 3.31131 e-17 7.58577 e-16 9.12010 e-17 7.76247 e-16 6.76082 e-18 3 1.38038 e-16 8.31763 e-17 1.99526 e-15 2.08929 e-16 1.25893 e-15 4.67735 e-17 4 3.38844 e-16 1.17489 e-17 1.41253 e-15 6.76082 e-17 2.81838 e-15 1.28824 e-16 5 1.86208 e-15 7.58577 e-16 5.24807 e-15 2.69153 e-16 3.23593 e-15 2.51188 e-16 6 1.90546 e-15 6.91830 e-16 1.62181 e-14 8.51138 e-16 9.33254 e-15 7.07945 e-16 7 7.76247 e-15 2.81838 e-15 1.28824 e-14 6.76082 e-16 6.02559 e-15 4.67735 e-16 8 2.88403 e-14 1.04712 e-14 4.78630 e-14 2.51188 e-15 1.99526 e-14 1.54881 e-15 9 6.30957 e-14 2.29086 e-14 1.34896 e-13 7.07945 e-15 4.16869 e-14 3.23593 e-15 10 1.41253 e-13 5.12861 e-14 4.16869 e-15 1.94984 e-16 5.88843 e-14 4.46683 e-15 5 Conclusions This pap er presents a metho d based on the a v erage of the round mo des which promo- tes an b etter p erformance in the simulation of recursiv e function. The metho d exploits the sto chastic nature of rounding. It consists on the av erage of tw o round mo de whic h promotes a reduction of the p ow er noise by a factor of 2. This is a go o d example of, although a simulation can b e seen as a theoretical exp erimen t, in fact, it presents some real world dimension, as sto c hasticit y , and therefore, usual to ols presented in Engineering can b e applied. W e applied the method in three n umerical examples. In the first case, w e sho w that the metho d is able to k eep the correct fixed p oin t, while the rounding to the nearest div erge to c haotic b eha viour. In the second and third w e show the general b ehaviour of the metho d in reducing the error of the sim ulation. It is imp ortant to sa y that there is no significant increase of computational time and this metho d can b e easily extended to man y other applications, suc h as n umerical solution of differen tial equations, as seen applied in [21]. In the future, we intend to inv estigate the use of in terv al extensions and other rounding mo des as further discussed in [24, 26]. Ac kno wledgements This w ork has b een supp orted b y the CNPq and CAPES. A sp ecial thanks to Con trol and Mo delling Group (GCOM). Refer ˆ encias [1] R. May Simple mathematical mo dels with very complicated dynamics, Nature, vol. 261, 459–467, (1976). [2] T. Matsumoto, A chaotic attractor from Chua’s circuit, Circuits and Systems, IEEE T ransactions on, vol. 31, 1055–1058, (1984). 7 [3] K. Sugihara and M. Iri, Tw o design principles of geometric algorithms in finite- precision arithmetic, Applied Mathematics Letters, v ol. 2, 203–206, (1989). [4] K. F uruya, H. Date and K. F ukudome, Effects of finite register length in the fast hartley transform, Electronics and Comm unications in Japan (Part I I I: F undamental Electronic Science, v ol. 75, 76–86, (1992). [5] R. M. Corless, What go o d are n umerical simulations of chaotic dynamical systems?, Computers & Mathematics with Applications, v ol. 28, 107–121, (1994). [6] P . V. Mcdonough and J. P . No onan and G. H. Hall, A new c haos detector, Computers & Electrical Engineering, v ol. 21, 417–431, (1995). [7] G. S. Maharana and P . K. Meher, Algorithm for efficien t interpolation of real-v alued signals using discrete Hartley transform, Computers & Electrical Engineering, vol. 23, 129–134, (1997). [8] E. Ott, Chaos in dynamical systems, Cam bridge: Cambridge Univ ersit y Press, (2002). [9] R. Landau and M. J. Paez and C. C. Bordeianu, Computational physics : problem solving with computers, Wiley-V CH John Wiley distributor, (2007). [10] IEEE. Institute of Electr. and Electr. Eng., 754-2008 - IEEE standard for floating- p oin t arithmetic, IEEE, (2008). [11] B. Parhami, Computer arithmetic algorithms and hardware arc hitectures, Oxford Univ. Press, New Y ork, (2012). [12] R. Gilmore and M. Lefranc, The top ology of chaos: Alice in stretch and squeezeland, John Wiley & Sons, (2012). [13] R. Gilmore and M. Lefranc, The top ology of chaos: Alice in stretch and squeezeland, John Wiley & Sons, (2012). [14] Z. Galias, The Dangers of Rounding Errors for Simulations and Analysis of Nonli- near Circuits and Systems - and Ho w to Av oid Them, IEEE Circuits and Systems Magazine, v ol. 13, 35–52, (2013). [15] R. Lozi, Can w e trust in numerical computations of chaotic solutions of dynamical systems?, W orld Scientific Publishing Co. Pte. Ltd., vol. 1, 63–98, (2013). [16] S.J. Liao and P .F. W ang, On the mathematically reliable long-term simulation of c haotic solutions of Lorenz equation in the interv al [0,10000], Science China: Ph ysics , Mec hanics and Astronom y , v ol. 57, 330-335, (2014). [17] C. W ang and T. T ang, Recursive least squares estimation algorithm applied to a class of linear-in-parameters output error moving av erage systems, Applied Mathematics Letters, v ol. 29, 35–41, (2014). 8 [18] E. G. Nep omuceno, Conv ergence of recursive functions on computers, The Journal of Engineering, 1–3, (2014). [19] G. M. V. Ladeira and G. V. Lima and C. A. Ara ˜ A o jo, Analysis of nonlinear dy- namics using electronic circuits, 23rd ABCM In ternational Congress of Mec hanical Engineering - COBEM2015, (2015). [20] Y. A. Gadzhiev, On the v ariance of a centered random v alue roundoff error, Signal Pro cessing, Elsevier, v ol. 106, 30–40, (2015). [21] K. Y u and C. W ang and H. Y au and J. Chen, Numerical computation and nonlinear dynamic analysis of ultrasonic cutting system, Computers & Electrical Engineering, v ol. 51, 270–283, (2016). [22] E. G. Nep omuceno, and S. A. M. Martins, A lo w er bound error for free-run sim ulation of the p olynomial NARMAX, Systems Science & Control Engineering, vol. 4, 50–58, (2016). [23] R. Li and S. Nadara jah, Mean and v ariance of round off error, Signal Pro cessing, v ol. 127, 185–190, (2016). [24] E. M. A. Mendes and E. G. Nepomuceno, A V ery Simple Metho d to Calculate the (P o- sitiv e) Largest Lyapuno v Exp onent Using In terv al Extensions, In ternational Journal of Bifurcation and Chaos, v ol. 26, 1650226, (2016). [25] W. Shi and J. W u and S. Lam and T. Srik anthan, Algorithms for bi-ob jectiv e multiple- c hoice hardw are/soft w are partitioning, Computers & Electrical Engineering, v ol. 50, 127–142, (2016). [26] E. G. Nep omuceno and E. M. A. Mendes, On the analysis of pseudo-orbits of con- tin uous chaotic nonlinear systems simulated using discretization schemes in a digital computer, Chaos, Solitons & F ractals, vol. 95, 21–32, (2017).