- Email: [email protected]
Computation of true chaotic orbits using cubic irrationals
Physica D 268 (2014) 100–105
Contents lists available at ScienceDirect
Physica D journal homepage: www.elsevier.com/locate/physd
Computation of true chaotic orbits using cubic irrationals Asaki Saito a,b,∗ , Shunji Ito c,1 a
Future University Hakodate, 116-2 Kameda Nakano-cho, Hakodate, Hokkaido 041-8655, Japan
b
PRESTO, Japan Science and Technology Agency, 4-1-8 Honcho, Kawaguchi, Saitama 332-0012, Japan
c
Kanazawa University, Kakuma-machi, Kanazawa, Ishikawa 920-1192, Japan
highlights • • • •
We propose a method for computing true orbits of 1D piecewise linear fractional maps. The method uses cubic irrationals as numbers and involves only integer arithmetic. True orbits generated by the method display the same properties as typical orbits. We can simulate maps whose simulation has been difficult, such as the Bernoulli map.
article
info
Article history: Received 28 June 2013 Received in revised form 31 October 2013 Accepted 8 November 2013 Available online 16 November 2013 Communicated by Y. Nishiura Keywords: True orbit Piecewise linear fractional map Algebraic numbers Integer arithmetic Simulation method Typical behavior
abstract We introduce a method that enables us to generate long true orbits of discrete-time dynamical systems defined by one-dimensional piecewise linear fractional maps with integer coefficients. The method uses cubic irrationals to represent numbers and involves only integer arithmetic to compute true orbits. By applying the method to the Bernoulli map and a modified Bernoulli map, we show that it successfully generates true chaotic and intermittent orbits, respectively, in contrast with conventional simulation methods. We demonstrate through simulations concerning invariant measures and the power spectrum that the statistical properties of the true orbits generated agree well with those of typical orbits of the two maps. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Chaotic phenomena occur in diverse areas of nonlinear science, and simulations using digital computers are essential for studying these phenomena [1–3]. The accuracy of simulations is an important element in determining the validity of such studies. However, the difficulty of generating true (exact) chaotic orbits using digital computers is also well known. This difficulty involves the following factors when discrete-time dynamical systems are considered. If we use a representation of numbers with fixed precision (i.e., a representation where the total number of bits used to represent
∗ Corresponding author at: Future University Hakodate, 116-2 Kameda Nakanocho, Hakodate, Hokkaido 041-8655, Japan. Tel.: +81 138 34 6128; fax: +81 138 34 6594. E-mail address: [email protected] (A. Saito). 1 Present address: Toho University, 2-2-1 Miyama, Funabashi, Chiba 274-8510, Japan. 0167-2789/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.physd.2013.11.003
the numbers is fixed), such as the double-precision floating-point numbers that are usually used in numerical computations, then it is impossible to generate true chaotic orbits due to inevitable round-off errors. In addition, because of the fixed bit length, all orbits become eventually periodic in principle. Accordingly, we can say that number representations with fixed precision are, in the first place, incompatible with the fundamental properties of chaos: sensitive dependence on initial conditions and aperiodicity. In fact, the validity of simulations of chaotic phenomena has been a matter of serious concern for decades (see, e.g., Refs. [4–7]). There is indeed a view (see, e.g., Refs. [4,5]) that the qualitative or statistical properties of a system are preserved, even if one performs discretization using a number representation with fixed precision, especially when the total number of possible discrete values is sufficiently large. However, it is also known that in some cases the use of fixed-precision number representations can destroy even the system’s macroscopic, qualitative properties, as exemplified in spatial bifurcation phenomena in open flow systems [8,9]. For the next standard number representation other than that of fixed-precision numbers, we can cite the fractional representation
A. Saito, S. Ito / Physica D 268 (2014) 100–105
of numbers. A digital computer can perform integer arithmetic, and thus the four arithmetic operations on rational numbers, without errors. Therefore, if we consider a dynamical system described by a rational map with integer coefficients and choose rational numbers as a number representation, we can generate errorless true orbits of the map in principle (the initial value(s) are to be chosen from the rational numbers). However, this scheme is also known to present two difficulties. 1. For a general rational map with integer coefficients, the representation of a point on a true orbit grows huge very rapidly. In fact, the numerator(s) and denominator(s) of the coordinate(s) of the point grow superexponentially, making the cost of the true orbit computation extremely high [10]. For example, even with the most powerful computers available, only the first few tens of time steps of the true orbit can be generated for the logistic map with generic rational parameters and initial values. 2. In the case of piecewise linear or linear fractional maps [11,12], the cost of computing true orbits is relatively low. However, there is still a danger of generating only atypical ones [2,13]. For example, for the Bernoulli map (i.e., the 2x modulo 1 map), the orbit from the initial value given by a rational number becomes eventually periodic, which makes the Bernoulli map a wellknown map that is difficult to simulate. The fact that true orbits are not necessarily typical constitutes another difficulty with true orbit computation because it is usually intended for simulations to reproduce the typical behavior of a target system. In the above discussion, we restricted the numbers used for the true orbit generation of a rational map to the rational numbers, but we can also consider expanding them to algebraic numbers (the coefficients of a rational map can also be expanded to algebraic numbers). Even if we expand numbers in this manner, it is still true as in the case of rational numbers that, for a general rational map, the computational cost necessary for generating true orbits is extremely high. However, mathematically, such dynamics over algebraic sets has been extensively investigated in the field of arithmetic dynamics. (See, e.g., Ref. [10] for an introduction to this developing field of research.) In particular, a related study was conducted on the evolution of (roots of) polynomials under a general rational map [14]. The method provided in Ref. [14] is theoretical but of high generality. For example, it can take arbitrary algebraic numbers as the coefficients of both a rational map and a polynomial. In our study, we will focus on piecewise linear fractional maps (including piecewise linear maps), for which we can expect to generate long true orbits by using integer arithmetic. More specifically, we will deal with one-dimensional piecewise linear fractional maps with integer coefficients, whose eventually periodic points are usually rational numbers and quadratic irrationals (i.e., the roots of irreducible quadratic polynomials with integer coefficients). Consequently, if we use algebraic numbers whose degree is greater than 2, we can usually guarantee the aperiodicity of the generated true orbits. In this paper, by using real cubic irrationals having complex conjugates as numbers, we propose a true orbit generation method involving only integer arithmetic and allowing exact simulations of discrete-time dynamical systems defined by one-dimensional piecewise linear fractional maps with integer coefficients. Moreover, we simulate examples of the Bernoulli map and a modified Bernoulli map and demonstrate that long true orbits obtained by this method display the same properties as typical orbits. We also confirm at most linear growth of the number representation, which is a distinctive characteristic of our method of computation. Before we proceed to the main subject, it is worth pointing out the relation of our method to the issue of shadowing. Shadowability has been focused on as a requirement for ensuring the validity of computer-generated (pseudo-)orbits [15]. For hyperbolic
101
dynamical systems, the shadowing lemma assures the existence of true orbits that trace pseudo-orbits for arbitrarily long times [16,17]. Typical physical systems, however, are considered nonhyperbolic, and for this reason researchers have been concerned with the question of to what extent (e.g., for how long) pseudoorbits remain traced by true orbits in the case of nonhyperbolic systems (see, e.g., Refs. [18–20]). In contrast, the method that we report here can directly generate true orbits regardless of whether a system is hyperbolic or not. Consequently, unlike before, there is no need with our method to verify the shadowability of generated orbits. 2. The proposed method 2.1. Preliminaries Let us consider a cubic equation, px + qx2 + rx + s = 0, 3
(1)
where the coefficients p (̸=0), q, r, s are integers and the polynomial px3 + qx2 + rx + s is irreducible. We suppose that Eq. (1) has a single real root α . In other words, α is a real cubic irrational having complex conjugates. We denote by S the set of such real cubic irrationals having complex conjugates. A number α ∈ S satisfying Eq. (1) is denoted by a vector (p, q, r , s) ∈ Z4 . Notice that this representation of α is not unique, e.g., (1, 1, 1, −1) = (−2, −2, −2, 2). Next, we consider the following linear fractional map: f (x) =
ax + b cx + d
,
(2)
where a, b, c, d are integers satisfying ad − bc ̸= 0 (i.e., f is invertible). It is easy to see that f becomes a map from S to S.2 Moreover, a representation of the image of α ∈ S under f can be easily obtained by a matrix operation: Let α ∈ S be a number represented by (p, q, r , s). Then, α ′ = f (α) has a representation (p′ , q′ , r ′ , s′ ) with
(p′ , q′ , r ′ , s′ )T = A (p, q, r , s)T ,
(3)
where the superscript T indicates the transpose and A is a matrix given by d3 −3bd2 A= 3b2 d −b 3
−cd2 2 ad + 2bcd −2abd − b2 c ab
2
c2d −2acd − bc 2 a2 d + 2abc −a 2 b
−c 3
3ac 2 . −3a2 c a3
(4)
One can obtain Eqs. (3) and (4) by substituting α = (dα ′ − b)/ (−c α ′ + a) into Eq. (1). Notice that A is determined solely by the coefficients a, b, c, d of f . 2.2. True orbit generation In this study, we deal with discrete-time dynamical systems defined by a one-dimensional piecewise linear fractional map, which can be expressed in the following form [11]: f (x) = fi (x) =
ai x + b i ci x + di
if x ∈ (ei , ei+1 ) , i = 0, 1, . . . , N − 1. (5)
2 If α ∈ S, then α ′ = f (α) ∈ Q(α), where Q(α) denotes the cubic number field obtained by adjoining α to the rational number field Q. Hence, α ′ is an algebraic number with degree 1 or 3. From α = (dα ′ − b)/(−c α ′ + a) with a, b, c, d ∈ Z, we see that α ′ ̸∈ Q, so α ′ must be a real cubic irrational. It is obvious that α ′ has σ complex conjugates since, by using one of the conjugates of α ′ (denoted as α ′ ), σ σ that of α (denoted as α σ ) can be represented by α σ = (dα ′ − b)/(−c α ′ + a). (See, e.g., Ref. [21] for an introduction to algebraic number fields.)
102
A. Saito, S. Ito / Physica D 268 (2014) 100–105
Here, fi is the linear fractional map on the ith interval (ei , ei+1 ), where ai , bi , ci , di are integers, and ai di − bi ci ̸= 0. For simplicity, we consider the case where the number of intervals partitioning the domain, denoted as N, is finite, but our method is equally applicable to the countably infinite case. Also, it is applicable even if the intervals include semi-infinite one(s). To simplify the interval determination explained below, we let all the endpoints of the intervals be rational numbers satisfying e0 < e1 < · · · < eN . Let Ai be the matrix corresponding to fi by Eq. (4). If α ∈ S represented by (p, q, r , s) is included in the ith interval (ei , ei+1 ), we obtain Ai (p, q, r , s)T representing f (α) from Eq. (3). One can exactly determine the interval including α , and therefore the matrix to be applied to (p, q, r , s) as follows: Recalling that α ∈ S is the unique real root of px3 + qx2 + rx + s, we see that α ∈ (ei , ei+1 ) if and only if the sign of pe3i + qe2i + rei + s is different from that of pe3i+1 + qe2i+1 + rei+1 + s. Therefore, to determine the matrix that should be applied to a given (p, q, r , s), we only need to check the sign of the cubic polynomial evaluated at the endpoints of each interval sequentially, until we find the only interval where the polynomial has different signs at its endpoints. Note that such a determination can be exactly performed only by using integer arithmetic when the endpoints of the intervals are rational numbers, and that the interval including α can be found in a finite amount of time because the number of intervals is at most countable. By repeatedly applying the interval determination and the consequently selected matrix to the newly obtained representation, we can generate a true orbit of the map expressed by Eq. (5). The following should be noted. In view of Eq. (3), we see that, for piecewise linear fractional maps, the coefficients of the polynomial grow at most exponentially with time (this will be confirmed later in the section on the simulations), in contrast with the superexponential growth for general rational maps. This implies that the computational cost of generating true orbits of the piecewise linear fractional maps is much lower than the cost of general rational maps. Nevertheless, compared with conventional simulation methods for generating pseudo-orbits, our method has a much higher computational cost. In the next section, we will show that it is indeed possible with our method to generate true orbits of the piecewise linear fractional maps long enough to allow statistical analysis. Moreover, we will demonstrate that the sufficiently long true orbits generated by our method display the same statistical properties as typical orbits. 3. Simulations In the following, we report results concerning the application of our true orbit generation method to two maps: the Bernoulli map and a modified Bernoulli map. We also explain why conventional simulation methods have serious difficulties in generating orbits of these maps. 3.1. Example 1: the Bernoulli map The Bernoulli map, xn+1 = 2xn modulo 1, where xn is a state at discrete times n = 0, 1, 2, . . . , is a piecewise linear map frequently considered in chaotic dynamics, but it is also a piecewise linear fractional map because one can consider the denominator to be 1. As mentioned above, the Bernoulli map is famous as a difficult map to simulate. To explain why, the following two characteristics of the Bernoulli map should be discussed. First, for the Bernoulli map, the set of eventually fixed points (i.e., the set of initial values from which the orbit eventually maps directly onto the unstable fixed point at x = 0) coincides with the set of finite binary decimals. As a result, if we use the doubleprecision binary floating-point numbers that are usually used
a
b
Fig. 1. (a) xn versus n and (b) ln versus n for the true orbit of the Bernoulli map starting from (p0 , q0 , r0 , s0 ) = (1, 1, 1, −1).
in simulating dynamical systems, then, even without numerical errors, we can generate only those orbits that eventually reach the unstable fixed point, and we cannot generate chaotic orbits that are typical for the Bernoulli map. The second characteristic is that, for the Bernoulli map, the set of eventually periodic points (i.e., the set of initial values from which the orbit eventually maps directly onto a periodic orbit) coincides with the set of rational numbers on the unit interval. Therefore, even if one improves the number representation to the fractional representation, one can only generate atypical orbits that eventually become periodic. Now we explain the practical application of our true orbit generation method to the Bernoulli map. As mentioned already, the matrix corresponding to the linear fractional map (Eq. (2)) is given by the form of Eq. (4). Therefore, the matrices A0 and A1 corresponding to the left and right branches of the Bernoulli map, i.e., f0 (x) = 2x/1 and f1 (x) = (2x − 1)/1, respectively, are as follows:
1 0 A0 = 0 0
0 2 0 0
0 0 4 0
0 0 , 0 8
1 3 A1 = 3 1
0 2 4 2
0 0 4 4
0 0 . 0 8
We take as an initial condition (p0 , q0 , r0 , s0 ) = (1, 1, 1, −1). That is, we choose x0 = α ∈ S satisfying x30 + x20 + x0 − 1 = 0 as an initial value. The specific procedure for the true orbit generation is as follows: Let (pn , qn , rn , sn ) represent xn . The signs of the cubic polynomial pn x3 + qn x2 + rn x + sn are evaluated at x = 0 and x = 1/2. If these two signs are different (and thus xn ∈ (0, 1/2)), the matrix A0 corresponding to f0 is applied to the representation at time n in order to obtain the representation at time n+1, as follows:
(pn+1 , qn+1 , rn+1 , sn+1 )T = A0 (pn , qn , rn , sn )T . Otherwise (and thus xn ∈ (1/2, 1)), we use the matrix A1 corresponding to f1 to obtain the representation at the next time step,
(pn+1 , qn+1 , rn+1 , sn+1 )T = A1 (pn , qn , rn , sn )T . By repeating this procedure, we obtain the true orbit (pn , qn , rn , sn ) (n = 0, 1, 2, . . .). To obtain the value of xn from (pn , qn , rn , sn ), we have only to find the unique real root of pn x3 + qn x2 + rn x + sn = 0; for example, if we use the bisection method, it is possible to obtain the value of xn with arbitrary precision. Fig. 1(a) shows the true orbit of the Bernoulli map starting from the initial condition (p0 , q0 , r0 , s0 ) = (1, 1, 1, −1). This figure
A. Saito, S. Ito / Physica D 268 (2014) 100–105
103
Fig. 2. (Color online) The estimated invariant density (dots) obtained from the true orbit of the Bernoulli map starting from (p0 , q0 , r0 , s0 ) = (1, 1, 1, −1). The line represents the density function ρ(x) = 1.
clearly shows chaotic behavior, which indicates that our method succeeds in generating a true chaotic orbit of the Bernoulli map, in contrast with conventional simulation methods. Next, we illustrate a distinctive characteristic of the computation, that is, at most exponential growth of the coefficients of polynomials. The first stage of the true orbit (pn , qn , rn , sn ) (n = 0, 1, 2, . . .) in Fig. 1(a) appears as follows:
(p0 , q0 , r0 , s0 ) = (1, 1, 1, −1), (p1 , q1 , r1 , s1 ) = (1, 5, 11, −1), (p2 , q2 , r2 , s2 ) = (1, 10, 44, −8), (p3 , q3 , r3 , s3 ) = (1, 20, 176, −64), (p4 , q4 , r4 , s4 ) = (1, 40, 704, −512), (p5 , q5 , r5 , s5 ) = (1, 83, 2979, −1199).
Fig. 3. The modified Bernoulli map.
a
b
To show the growth more clearly, we define ln , the virtual length in decimal notation of the representation of a point on the true orbit, by ln = ⌊log10 (max {|pn |, 1}) + 1⌋ + ⌊log10 (max {|qn |, 1}) + 1⌋
+ ⌊log10 (max {|rn |, 1}) + 1⌋ + ⌊log10 (max {|sn |, 1}) + 1⌋ , where ⌊x⌋ is the integer part of the real number x. Then, the increase in ln should be at most linear with n in our computation. Fig. 1(b) shows a plot of ln versus n for the true orbit starting from (1, 1, 1, −1). We can see a clear linear increase in ln with n for the Bernoulli map. Now we investigate the statistical properties of the generated true orbit by estimating the invariant density of the Bernoulli map. The Bernoulli map has an invariant measure with the density function ρ(x) equal to 1 on the unit interval. Fig. 2 shows the estimated density obtained from the true orbit with length 106 , starting from the same initial value represented by (1, 1, 1, −1) as before. The estimated density takes values very close to 1 and coincides well with the density of the Bernoulli map. 3.2. Example 2: the modified Bernoulli map The second example is the modified Bernoulli map. A modified Bernoulli map was proposed as a map exhibiting intermittent behavior [22]. Here, we treat the following modified Bernoulli map [13], which is the piecewise linear fractional map version of the original one, as well as the Bernoulli version of the map associated with the mediant convergents algorithm [23]:
xn+1 =
xn −x + 1
if xn ∈ [0, 1/2)
2xn − 1
if xn ∈ [1/2, 1) .
n
xn
The graph of this map is shown in Fig. 3. As in the case of the Bernoulli map, the following characteristics of the modified Bernoulli map make the simulation very difficult. First, for the modified Bernoulli map, the set of eventually fixed points coincides with the set of rational numbers on the unit interval. This indicates that, like for the case of the Bernoulli map,
Fig. 4. (a) xn versus n and (b) ln versus n for the true orbit of the modified Bernoulli map starting from (p0 , q0 , r0 , s0 ) = (1, 1, 1, −1).
we cannot generate intermittent orbits that are typical for the modified Bernoulli map by using the fractional representation.3 A different viewpoint reveals another difficulty in simulating this map. This map displays a 1/f spectrum [13]. This property depends on the local structure around the fixed point at x = 0 (and at x = 1). However, the use of the fixed-precision number representations destroys this property, because inevitable roundoff errors break the local structure, at least in principle. The matrices needed to generate true orbits – that is, A0 and A1 corresponding to the left and right branches of the modified Bernoulli map, respectively – can be obtained from Eq. (4), like for the case of the Bernoulli map. Here, we can use the same initial condition and procedure for the true orbit generation as we used before. Fig. 4(a) shows the true orbit of the modified Bernoulli map starting from the initial condition (p0 , q0 , r0 , s0 ) = (1, 1, 1, −1). From the figure, we can clearly see intermittent behavior and that our method successfully generates a true intermittent orbit of the modified Bernoulli map, in contrast with conventional simulation methods. Here again, we examine the growth of the number representation. Like for Fig. 1(b), we consider ln , the length of the representation of a point on the true orbit, also for the modified Bernoulli
3 For the modified Bernoulli map, the set of eventually periodic points coincides with the set of rational numbers and quadratic irrationals on the unit interval.
104
a
A. Saito, S. Ito / Physica D 268 (2014) 100–105
b
Fig. 5. (Color online) Statistical properties obtained from the true orbit of the modified Bernoulli map starting from (p0 , q0 , r0 , s0 ) = (1, 1, 1, −1). (a) The estimated invariant density (dots). The curve represents the density function ρ(x) = x−1 (1 − x)−1 . (b) The estimated power spectral density (dots) plotted on double-logarithmic scales. The line shows 1/f divergence.
map and examine the dependence on n of ln in Fig. 4(b). It is clear that ln increases with n due to the aperiodicity of the true orbits, but, compared with the case for the Bernoulli map, the increase in ln is much slower. For example, l104 = 15,056 for the Bernoulli map, whereas l104 = 1872 for the modified Bernoulli map. Also, ln for the modified Bernoulli map does not show a simple linear increase such as that in Fig. 1(b). That is, ln in Fig. 4(b) has plateaus that correspond to laminar phases, during which the true orbit has an almost constant behavior near x = 0 and x = 1 (cf. Fig. 4(a)). Next, let us consider again the question of whether or not the statistical properties of the generated true orbit are good. The modified Bernoulli map has a σ -finite but infinite invariant measure µ with density dµ = x−1 (1−x)−1 dx on the unit interval [13]. Fig. 5(a) shows the estimated density obtained from the true orbit with length 106 starting from (1, 1, 1, −1) (see the Appendix for details on the density estimation). The estimated density agrees well with that of the modified Bernoulli map. Similarly, Fig. 5(b) shows the estimated power spectral density (periodogram) obtained by performing the discrete Fourier transform of the same true orbit. The estimated power spectral density is in good agreement with the 1/f spectral density of the modified Bernoulli map. 4. Conclusions By using real cubic irrationals having complex conjugates and integer arithmetic, we have introduced a true orbit generation method enabling exact simulations of discrete-time dynamical systems defined by one-dimensional piecewise linear fractional maps with integer coefficients. We have demonstrated that the true orbits generated by our method show clear chaotic and intermittent behaviors for the Bernoulli map and modified Bernoulli map, respectively, in contrast with the case for conventional simulation methods. In addition, for each case, we have confirmed at most linear growth of the number representation, which is a distinctive characteristic of our method of computation. In particular, we have shown that the statistical properties obtained by the generated long true orbits match those of typical orbits of the above two maps in simulations concerning invariant measures and the power spectrum. Expansion of the true orbit generation method to higher dimensional maps, as well as application of our method to, e.g., analysis of noise-induced order phenomena [24] and generation of pseudo-random numbers, will be reported elsewhere.
Appendix. The method for estimating the density of an infinite invariant measure It is known that, unlike the case for dynamical systems preserving a probability measure, it is very difficult to statistically treat those preserving an infinite measure (see, e.g., Refs. [25–27]); for example, time-averaged values obtained from a single orbit do not converge to actual values (i.e., values of integrals with respect to the invariant measure). Nevertheless, as seen below, we can estimate the density of the σ -finite invariant measure for a conservative ergodic map. In what follows, we explain this through the density estimation, shown in Fig. 5(a), of the modified Bernoulli map M. Let µ be the invariant measure on the unit interval X = [0, 1), whose density is ρ(x) = cx−1 (1 − x)−1 with an indefinite constant c ∈ (0, ∞), and let x ∈ X be the point at which we want to estimate the density. Consider a small interval A ⊂ X of length ∆ such that x ∈ A, and also consider some arbitrarily chosen set B ⊂ X , satisfying 0 < µ(A), µ(B) < ∞. By Hopf’s ratio ergodic theorem [28,29], we obtain, for a.e. x0 ∈ X , that
♯ {0 ≤ n < N | M n (x0 ) ∈ A} µ(A) = , N →∞ ♯ {0 ≤ n < N | M n (x0 ) ∈ B} µ(B) lim
where ♯S denotes the number of elements in the set S. Since µ(A) is approximated by ρ(x)∆, we see that if x0 is a typical initial condition, then
ρ(x) ≈
♯ {0 ≤ n < N | M n (x0 ) ∈ A} µ(B) ♯ {0 ≤ n < N | M n (x0 ) ∈ B} ∆
for sufficiently large N.4 The detailed procedure for obtaining the estimated density in Fig. 5(a) is as follows: We first divide X into 25 subintervals Ai = (i/25, (i + 1)/25) (i = 0, 1, 2, . . . , 24) of length ∆ = 1/25. Also, we choose the interval (1/3, 2/3) as B. Notice that µ(B) in Eq. (A.1) depends on the indefinite constant c, but in this study we consider 2/3 c = 1, so it has a definite value: µ(B) = 1/3 x−1 (1 − x)−1 dx =
1.386294 · · ·.5 We then generate the true orbit with length N = 106 starting from x0 represented by (1, 1, 1, −1), and count the number of times the orbit visits each Ai and B. Lastly, using Eq. (A.1), we estimate the density values for the 23 subintervals Ai (i = 1, 2, . . . , 23) satisfying µ(Ai ) < ∞. Notice that we cannot apply Eq. (A.1) to A0 and A24 , since they are sets of infinite measure and beyond the scope of Hopf’s ratio ergodic theorem. Finally, if x0 represented by (1, 1, 1, −1) is a typical initial condition, then the estimated density obtained by the above procedure should approximately coincide with ρ(x) = x−1 (1 − x)−1 . Fig. 5(a) shows a good coincidence between the two, which supports the basic premise of our method, that is, cubic irrationals are typical initial conditions, so true orbits starting from them possess good statistical properties.
4 More generally, for f , g ∈ L+ (µ), we have 1 N −1
f dµ ≈ X
n =0 N −1
Acknowledgments
f (M n (x0 )) g(
Mn
g dµ, X
(x0 ))
n=0
and we can estimate the value of X f dµ from that of X g dµ. 5 In fact, the evaluation of µ(B) is unnecessary if we merely desire to obtain one of the invariant densities: By replacing µ(B) in Eq. (A.1) with some arbitrary positive constant, say 1, we obtain an invariant density ρ(x) = cx−1 (1 − x)−1 with
We thank E. Brooks-Pollock, R.S. MacKay, J. Tamura, and S. Yasutomi for their suggestions. This research was supported by the JST PRESTO program and by a Grant-in-Aid for Young Scientists (B) 21700256 from MEXT of Japan.
(A.1)
c = 1/
2/3 1/3
x−1 (1 − x)−1 dx.
A. Saito, S. Ito / Physica D 268 (2014) 100–105
References [1] E. Ott, Chaos in Dynamical Systems, second ed., Cambridge University Press, Cambridge, 2002. [2] E. Atlee Jackson, Perspectives of Nonlinear Dynamics, Cambridge University Press, Cambridge, 1991. [3] Y. Ueda, The Road to Chaos, second ed., Aerial Press, Santa Cruz, CA, 2001. [4] F. Rannou, Numerical study of discrete plane area-preserving mappings, Astronom. Astrophys. 31 (1974) 289–301. [5] P.M. Binder, R.V. Jensen, Simulating chaotic behavior with finite-state machines, Phys. Rev. A 34 (1986) 4460–4463. [6] C. Beck, G. Roepstorff, Effects of phase space discretization on the long-time behavior of dynamical systems, Physica D 25 (1987) 173–180. [7] D.J.D. Earn, S. Tremaine, Exact numerical studies of Hamiltonian maps: iterating without roundoff error, Physica D 56 (1992) 1–22. [8] K. Kaneko, Spatial period-doubling in open flow, Phys. Lett. A 111 (1985) 321–325. [9] A. Yamaguchi, On the mechanism of spatial bifurcations in the open flow system, Int. J. Bifurcat. Chaos 7 (1997) 1529–1538. [10] J.H. Silverman, The Arithmetic of Dynamical Systems, Springer, New York, 2007. [11] F. Schweiger, Ergodic Theory of Fibred Systems and Metric Number Theory, Clarendon Press, Oxford, 1995. [12] F. Schweiger, Multidimensional Continued Fractions, Oxford University Press, Oxford, 2000. [13] A. Saito, Computational aspects of a modified Bernoulli map, Prog. Theoret. Phys. Supp. 161 (2006) 328–331. [14] F. Vivaldi, Dynamics over irreducible polynomials, Nonlinearity 5 (1992) 941–960.
105
[15] C. Grebogi, L. Poon, T. Sauer, J.A. Yorke, D. Auerbach, Shadowability of chaotic dynamical systems, in: B. Fiedler (Ed.), in: Handbook of Dynamical Systems, vol. 2, Elsevier, Amsterdam, 2002, pp. 313–344. [16] D.V. Anosov, Geodesic flows on closed Riemannian manifolds of negative curvature, Proc. Steklov Inst. Math. 90 (1967) 1–235. [17] R. Bowen, ω-limit sets for axiom A diffeomorphisms, J. Differential Equations 18 (1975) 333–339. [18] C. Grebogi, S.M. Hammel, J.A. Yorke, T. Sauer, Shadowing of physical trajectories in chaotic dynamics: containment and refinement, Phys. Rev. Lett. 65 (1990) 1527–1530. [19] S. Dawson, C. Grebogi, T. Sauer, J.A. Yorke, Obstructions to shadowing when a Lyapunov exponent fluctuates about zero, Phys. Rev. Lett. 73 (1994) 1927–1930. [20] T. Sauer, C. Grebogi, J.A. Yorke, How long do numerical chaotic solutions remain valid? Phys. Rev. Lett. 79 (1997) 59–62. [21] B.L. van der Waerden, Algebra I, Springer, Berlin, 1966. [22] Y. Aizawa, T. Kohyama, Asymptotically non-stationary chaos, Progr. Theoret. Phys. 71 (1984) 847–850. [23] S. Ito, Algorithms with mediant convergents and their metrical theory, Osaka J. Math. 26 (1989) 557–578. [24] K. Matsumoto, I. Tsuda, Noise-induced order, J. Stat. Phys. 31 (1983) 87–106. [25] J. Aaronson, An Introduction to Infinite Ergodic Theory, American Mathematical Society, 1997. [26] T. Akimoto, Y. Aizawa, New aspects of the correlation functions in nonhyperbolic chaotic systems, J. Korean Phys. Soc. 50 (2007) 254–260. [27] T. Akimoto, E. Barkai, Aging generates regular motions in weakly chaotic systems, Phys. Rev. E 87 (2013) 032915. [28] E. Hopf, Ergodentheorie, Springer, Berlin, 1937. [29] R. Zweimüller, Hopf’s ratio ergodic theorem by inducing, Colloq. Math. 101 (2004) 289–292.
Contents lists available at ScienceDirect
Physica D journal homepage: www.elsevier.com/locate/physd
Computation of true chaotic orbits using cubic irrationals Asaki Saito a,b,∗ , Shunji Ito c,1 a
Future University Hakodate, 116-2 Kameda Nakano-cho, Hakodate, Hokkaido 041-8655, Japan
b
PRESTO, Japan Science and Technology Agency, 4-1-8 Honcho, Kawaguchi, Saitama 332-0012, Japan
c
Kanazawa University, Kakuma-machi, Kanazawa, Ishikawa 920-1192, Japan
highlights • • • •
We propose a method for computing true orbits of 1D piecewise linear fractional maps. The method uses cubic irrationals as numbers and involves only integer arithmetic. True orbits generated by the method display the same properties as typical orbits. We can simulate maps whose simulation has been difficult, such as the Bernoulli map.
article
info
Article history: Received 28 June 2013 Received in revised form 31 October 2013 Accepted 8 November 2013 Available online 16 November 2013 Communicated by Y. Nishiura Keywords: True orbit Piecewise linear fractional map Algebraic numbers Integer arithmetic Simulation method Typical behavior
abstract We introduce a method that enables us to generate long true orbits of discrete-time dynamical systems defined by one-dimensional piecewise linear fractional maps with integer coefficients. The method uses cubic irrationals to represent numbers and involves only integer arithmetic to compute true orbits. By applying the method to the Bernoulli map and a modified Bernoulli map, we show that it successfully generates true chaotic and intermittent orbits, respectively, in contrast with conventional simulation methods. We demonstrate through simulations concerning invariant measures and the power spectrum that the statistical properties of the true orbits generated agree well with those of typical orbits of the two maps. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Chaotic phenomena occur in diverse areas of nonlinear science, and simulations using digital computers are essential for studying these phenomena [1–3]. The accuracy of simulations is an important element in determining the validity of such studies. However, the difficulty of generating true (exact) chaotic orbits using digital computers is also well known. This difficulty involves the following factors when discrete-time dynamical systems are considered. If we use a representation of numbers with fixed precision (i.e., a representation where the total number of bits used to represent
∗ Corresponding author at: Future University Hakodate, 116-2 Kameda Nakanocho, Hakodate, Hokkaido 041-8655, Japan. Tel.: +81 138 34 6128; fax: +81 138 34 6594. E-mail address: [email protected] (A. Saito). 1 Present address: Toho University, 2-2-1 Miyama, Funabashi, Chiba 274-8510, Japan. 0167-2789/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.physd.2013.11.003
the numbers is fixed), such as the double-precision floating-point numbers that are usually used in numerical computations, then it is impossible to generate true chaotic orbits due to inevitable round-off errors. In addition, because of the fixed bit length, all orbits become eventually periodic in principle. Accordingly, we can say that number representations with fixed precision are, in the first place, incompatible with the fundamental properties of chaos: sensitive dependence on initial conditions and aperiodicity. In fact, the validity of simulations of chaotic phenomena has been a matter of serious concern for decades (see, e.g., Refs. [4–7]). There is indeed a view (see, e.g., Refs. [4,5]) that the qualitative or statistical properties of a system are preserved, even if one performs discretization using a number representation with fixed precision, especially when the total number of possible discrete values is sufficiently large. However, it is also known that in some cases the use of fixed-precision number representations can destroy even the system’s macroscopic, qualitative properties, as exemplified in spatial bifurcation phenomena in open flow systems [8,9]. For the next standard number representation other than that of fixed-precision numbers, we can cite the fractional representation
A. Saito, S. Ito / Physica D 268 (2014) 100–105
of numbers. A digital computer can perform integer arithmetic, and thus the four arithmetic operations on rational numbers, without errors. Therefore, if we consider a dynamical system described by a rational map with integer coefficients and choose rational numbers as a number representation, we can generate errorless true orbits of the map in principle (the initial value(s) are to be chosen from the rational numbers). However, this scheme is also known to present two difficulties. 1. For a general rational map with integer coefficients, the representation of a point on a true orbit grows huge very rapidly. In fact, the numerator(s) and denominator(s) of the coordinate(s) of the point grow superexponentially, making the cost of the true orbit computation extremely high [10]. For example, even with the most powerful computers available, only the first few tens of time steps of the true orbit can be generated for the logistic map with generic rational parameters and initial values. 2. In the case of piecewise linear or linear fractional maps [11,12], the cost of computing true orbits is relatively low. However, there is still a danger of generating only atypical ones [2,13]. For example, for the Bernoulli map (i.e., the 2x modulo 1 map), the orbit from the initial value given by a rational number becomes eventually periodic, which makes the Bernoulli map a wellknown map that is difficult to simulate. The fact that true orbits are not necessarily typical constitutes another difficulty with true orbit computation because it is usually intended for simulations to reproduce the typical behavior of a target system. In the above discussion, we restricted the numbers used for the true orbit generation of a rational map to the rational numbers, but we can also consider expanding them to algebraic numbers (the coefficients of a rational map can also be expanded to algebraic numbers). Even if we expand numbers in this manner, it is still true as in the case of rational numbers that, for a general rational map, the computational cost necessary for generating true orbits is extremely high. However, mathematically, such dynamics over algebraic sets has been extensively investigated in the field of arithmetic dynamics. (See, e.g., Ref. [10] for an introduction to this developing field of research.) In particular, a related study was conducted on the evolution of (roots of) polynomials under a general rational map [14]. The method provided in Ref. [14] is theoretical but of high generality. For example, it can take arbitrary algebraic numbers as the coefficients of both a rational map and a polynomial. In our study, we will focus on piecewise linear fractional maps (including piecewise linear maps), for which we can expect to generate long true orbits by using integer arithmetic. More specifically, we will deal with one-dimensional piecewise linear fractional maps with integer coefficients, whose eventually periodic points are usually rational numbers and quadratic irrationals (i.e., the roots of irreducible quadratic polynomials with integer coefficients). Consequently, if we use algebraic numbers whose degree is greater than 2, we can usually guarantee the aperiodicity of the generated true orbits. In this paper, by using real cubic irrationals having complex conjugates as numbers, we propose a true orbit generation method involving only integer arithmetic and allowing exact simulations of discrete-time dynamical systems defined by one-dimensional piecewise linear fractional maps with integer coefficients. Moreover, we simulate examples of the Bernoulli map and a modified Bernoulli map and demonstrate that long true orbits obtained by this method display the same properties as typical orbits. We also confirm at most linear growth of the number representation, which is a distinctive characteristic of our method of computation. Before we proceed to the main subject, it is worth pointing out the relation of our method to the issue of shadowing. Shadowability has been focused on as a requirement for ensuring the validity of computer-generated (pseudo-)orbits [15]. For hyperbolic
101
dynamical systems, the shadowing lemma assures the existence of true orbits that trace pseudo-orbits for arbitrarily long times [16,17]. Typical physical systems, however, are considered nonhyperbolic, and for this reason researchers have been concerned with the question of to what extent (e.g., for how long) pseudoorbits remain traced by true orbits in the case of nonhyperbolic systems (see, e.g., Refs. [18–20]). In contrast, the method that we report here can directly generate true orbits regardless of whether a system is hyperbolic or not. Consequently, unlike before, there is no need with our method to verify the shadowability of generated orbits. 2. The proposed method 2.1. Preliminaries Let us consider a cubic equation, px + qx2 + rx + s = 0, 3
(1)
where the coefficients p (̸=0), q, r, s are integers and the polynomial px3 + qx2 + rx + s is irreducible. We suppose that Eq. (1) has a single real root α . In other words, α is a real cubic irrational having complex conjugates. We denote by S the set of such real cubic irrationals having complex conjugates. A number α ∈ S satisfying Eq. (1) is denoted by a vector (p, q, r , s) ∈ Z4 . Notice that this representation of α is not unique, e.g., (1, 1, 1, −1) = (−2, −2, −2, 2). Next, we consider the following linear fractional map: f (x) =
ax + b cx + d
,
(2)
where a, b, c, d are integers satisfying ad − bc ̸= 0 (i.e., f is invertible). It is easy to see that f becomes a map from S to S.2 Moreover, a representation of the image of α ∈ S under f can be easily obtained by a matrix operation: Let α ∈ S be a number represented by (p, q, r , s). Then, α ′ = f (α) has a representation (p′ , q′ , r ′ , s′ ) with
(p′ , q′ , r ′ , s′ )T = A (p, q, r , s)T ,
(3)
where the superscript T indicates the transpose and A is a matrix given by d3 −3bd2 A= 3b2 d −b 3
−cd2 2 ad + 2bcd −2abd − b2 c ab
2
c2d −2acd − bc 2 a2 d + 2abc −a 2 b
−c 3
3ac 2 . −3a2 c a3
(4)
One can obtain Eqs. (3) and (4) by substituting α = (dα ′ − b)/ (−c α ′ + a) into Eq. (1). Notice that A is determined solely by the coefficients a, b, c, d of f . 2.2. True orbit generation In this study, we deal with discrete-time dynamical systems defined by a one-dimensional piecewise linear fractional map, which can be expressed in the following form [11]: f (x) = fi (x) =
ai x + b i ci x + di
if x ∈ (ei , ei+1 ) , i = 0, 1, . . . , N − 1. (5)
2 If α ∈ S, then α ′ = f (α) ∈ Q(α), where Q(α) denotes the cubic number field obtained by adjoining α to the rational number field Q. Hence, α ′ is an algebraic number with degree 1 or 3. From α = (dα ′ − b)/(−c α ′ + a) with a, b, c, d ∈ Z, we see that α ′ ̸∈ Q, so α ′ must be a real cubic irrational. It is obvious that α ′ has σ complex conjugates since, by using one of the conjugates of α ′ (denoted as α ′ ), σ σ that of α (denoted as α σ ) can be represented by α σ = (dα ′ − b)/(−c α ′ + a). (See, e.g., Ref. [21] for an introduction to algebraic number fields.)
102
A. Saito, S. Ito / Physica D 268 (2014) 100–105
Here, fi is the linear fractional map on the ith interval (ei , ei+1 ), where ai , bi , ci , di are integers, and ai di − bi ci ̸= 0. For simplicity, we consider the case where the number of intervals partitioning the domain, denoted as N, is finite, but our method is equally applicable to the countably infinite case. Also, it is applicable even if the intervals include semi-infinite one(s). To simplify the interval determination explained below, we let all the endpoints of the intervals be rational numbers satisfying e0 < e1 < · · · < eN . Let Ai be the matrix corresponding to fi by Eq. (4). If α ∈ S represented by (p, q, r , s) is included in the ith interval (ei , ei+1 ), we obtain Ai (p, q, r , s)T representing f (α) from Eq. (3). One can exactly determine the interval including α , and therefore the matrix to be applied to (p, q, r , s) as follows: Recalling that α ∈ S is the unique real root of px3 + qx2 + rx + s, we see that α ∈ (ei , ei+1 ) if and only if the sign of pe3i + qe2i + rei + s is different from that of pe3i+1 + qe2i+1 + rei+1 + s. Therefore, to determine the matrix that should be applied to a given (p, q, r , s), we only need to check the sign of the cubic polynomial evaluated at the endpoints of each interval sequentially, until we find the only interval where the polynomial has different signs at its endpoints. Note that such a determination can be exactly performed only by using integer arithmetic when the endpoints of the intervals are rational numbers, and that the interval including α can be found in a finite amount of time because the number of intervals is at most countable. By repeatedly applying the interval determination and the consequently selected matrix to the newly obtained representation, we can generate a true orbit of the map expressed by Eq. (5). The following should be noted. In view of Eq. (3), we see that, for piecewise linear fractional maps, the coefficients of the polynomial grow at most exponentially with time (this will be confirmed later in the section on the simulations), in contrast with the superexponential growth for general rational maps. This implies that the computational cost of generating true orbits of the piecewise linear fractional maps is much lower than the cost of general rational maps. Nevertheless, compared with conventional simulation methods for generating pseudo-orbits, our method has a much higher computational cost. In the next section, we will show that it is indeed possible with our method to generate true orbits of the piecewise linear fractional maps long enough to allow statistical analysis. Moreover, we will demonstrate that the sufficiently long true orbits generated by our method display the same statistical properties as typical orbits. 3. Simulations In the following, we report results concerning the application of our true orbit generation method to two maps: the Bernoulli map and a modified Bernoulli map. We also explain why conventional simulation methods have serious difficulties in generating orbits of these maps. 3.1. Example 1: the Bernoulli map The Bernoulli map, xn+1 = 2xn modulo 1, where xn is a state at discrete times n = 0, 1, 2, . . . , is a piecewise linear map frequently considered in chaotic dynamics, but it is also a piecewise linear fractional map because one can consider the denominator to be 1. As mentioned above, the Bernoulli map is famous as a difficult map to simulate. To explain why, the following two characteristics of the Bernoulli map should be discussed. First, for the Bernoulli map, the set of eventually fixed points (i.e., the set of initial values from which the orbit eventually maps directly onto the unstable fixed point at x = 0) coincides with the set of finite binary decimals. As a result, if we use the doubleprecision binary floating-point numbers that are usually used
a
b
Fig. 1. (a) xn versus n and (b) ln versus n for the true orbit of the Bernoulli map starting from (p0 , q0 , r0 , s0 ) = (1, 1, 1, −1).
in simulating dynamical systems, then, even without numerical errors, we can generate only those orbits that eventually reach the unstable fixed point, and we cannot generate chaotic orbits that are typical for the Bernoulli map. The second characteristic is that, for the Bernoulli map, the set of eventually periodic points (i.e., the set of initial values from which the orbit eventually maps directly onto a periodic orbit) coincides with the set of rational numbers on the unit interval. Therefore, even if one improves the number representation to the fractional representation, one can only generate atypical orbits that eventually become periodic. Now we explain the practical application of our true orbit generation method to the Bernoulli map. As mentioned already, the matrix corresponding to the linear fractional map (Eq. (2)) is given by the form of Eq. (4). Therefore, the matrices A0 and A1 corresponding to the left and right branches of the Bernoulli map, i.e., f0 (x) = 2x/1 and f1 (x) = (2x − 1)/1, respectively, are as follows:
1 0 A0 = 0 0
0 2 0 0
0 0 4 0
0 0 , 0 8
1 3 A1 = 3 1
0 2 4 2
0 0 4 4
0 0 . 0 8
We take as an initial condition (p0 , q0 , r0 , s0 ) = (1, 1, 1, −1). That is, we choose x0 = α ∈ S satisfying x30 + x20 + x0 − 1 = 0 as an initial value. The specific procedure for the true orbit generation is as follows: Let (pn , qn , rn , sn ) represent xn . The signs of the cubic polynomial pn x3 + qn x2 + rn x + sn are evaluated at x = 0 and x = 1/2. If these two signs are different (and thus xn ∈ (0, 1/2)), the matrix A0 corresponding to f0 is applied to the representation at time n in order to obtain the representation at time n+1, as follows:
(pn+1 , qn+1 , rn+1 , sn+1 )T = A0 (pn , qn , rn , sn )T . Otherwise (and thus xn ∈ (1/2, 1)), we use the matrix A1 corresponding to f1 to obtain the representation at the next time step,
(pn+1 , qn+1 , rn+1 , sn+1 )T = A1 (pn , qn , rn , sn )T . By repeating this procedure, we obtain the true orbit (pn , qn , rn , sn ) (n = 0, 1, 2, . . .). To obtain the value of xn from (pn , qn , rn , sn ), we have only to find the unique real root of pn x3 + qn x2 + rn x + sn = 0; for example, if we use the bisection method, it is possible to obtain the value of xn with arbitrary precision. Fig. 1(a) shows the true orbit of the Bernoulli map starting from the initial condition (p0 , q0 , r0 , s0 ) = (1, 1, 1, −1). This figure
A. Saito, S. Ito / Physica D 268 (2014) 100–105
103
Fig. 2. (Color online) The estimated invariant density (dots) obtained from the true orbit of the Bernoulli map starting from (p0 , q0 , r0 , s0 ) = (1, 1, 1, −1). The line represents the density function ρ(x) = 1.
clearly shows chaotic behavior, which indicates that our method succeeds in generating a true chaotic orbit of the Bernoulli map, in contrast with conventional simulation methods. Next, we illustrate a distinctive characteristic of the computation, that is, at most exponential growth of the coefficients of polynomials. The first stage of the true orbit (pn , qn , rn , sn ) (n = 0, 1, 2, . . .) in Fig. 1(a) appears as follows:
(p0 , q0 , r0 , s0 ) = (1, 1, 1, −1), (p1 , q1 , r1 , s1 ) = (1, 5, 11, −1), (p2 , q2 , r2 , s2 ) = (1, 10, 44, −8), (p3 , q3 , r3 , s3 ) = (1, 20, 176, −64), (p4 , q4 , r4 , s4 ) = (1, 40, 704, −512), (p5 , q5 , r5 , s5 ) = (1, 83, 2979, −1199).
Fig. 3. The modified Bernoulli map.
a
b
To show the growth more clearly, we define ln , the virtual length in decimal notation of the representation of a point on the true orbit, by ln = ⌊log10 (max {|pn |, 1}) + 1⌋ + ⌊log10 (max {|qn |, 1}) + 1⌋
+ ⌊log10 (max {|rn |, 1}) + 1⌋ + ⌊log10 (max {|sn |, 1}) + 1⌋ , where ⌊x⌋ is the integer part of the real number x. Then, the increase in ln should be at most linear with n in our computation. Fig. 1(b) shows a plot of ln versus n for the true orbit starting from (1, 1, 1, −1). We can see a clear linear increase in ln with n for the Bernoulli map. Now we investigate the statistical properties of the generated true orbit by estimating the invariant density of the Bernoulli map. The Bernoulli map has an invariant measure with the density function ρ(x) equal to 1 on the unit interval. Fig. 2 shows the estimated density obtained from the true orbit with length 106 , starting from the same initial value represented by (1, 1, 1, −1) as before. The estimated density takes values very close to 1 and coincides well with the density of the Bernoulli map. 3.2. Example 2: the modified Bernoulli map The second example is the modified Bernoulli map. A modified Bernoulli map was proposed as a map exhibiting intermittent behavior [22]. Here, we treat the following modified Bernoulli map [13], which is the piecewise linear fractional map version of the original one, as well as the Bernoulli version of the map associated with the mediant convergents algorithm [23]:
xn+1 =
xn −x + 1
if xn ∈ [0, 1/2)
2xn − 1
if xn ∈ [1/2, 1) .
n
xn
The graph of this map is shown in Fig. 3. As in the case of the Bernoulli map, the following characteristics of the modified Bernoulli map make the simulation very difficult. First, for the modified Bernoulli map, the set of eventually fixed points coincides with the set of rational numbers on the unit interval. This indicates that, like for the case of the Bernoulli map,
Fig. 4. (a) xn versus n and (b) ln versus n for the true orbit of the modified Bernoulli map starting from (p0 , q0 , r0 , s0 ) = (1, 1, 1, −1).
we cannot generate intermittent orbits that are typical for the modified Bernoulli map by using the fractional representation.3 A different viewpoint reveals another difficulty in simulating this map. This map displays a 1/f spectrum [13]. This property depends on the local structure around the fixed point at x = 0 (and at x = 1). However, the use of the fixed-precision number representations destroys this property, because inevitable roundoff errors break the local structure, at least in principle. The matrices needed to generate true orbits – that is, A0 and A1 corresponding to the left and right branches of the modified Bernoulli map, respectively – can be obtained from Eq. (4), like for the case of the Bernoulli map. Here, we can use the same initial condition and procedure for the true orbit generation as we used before. Fig. 4(a) shows the true orbit of the modified Bernoulli map starting from the initial condition (p0 , q0 , r0 , s0 ) = (1, 1, 1, −1). From the figure, we can clearly see intermittent behavior and that our method successfully generates a true intermittent orbit of the modified Bernoulli map, in contrast with conventional simulation methods. Here again, we examine the growth of the number representation. Like for Fig. 1(b), we consider ln , the length of the representation of a point on the true orbit, also for the modified Bernoulli
3 For the modified Bernoulli map, the set of eventually periodic points coincides with the set of rational numbers and quadratic irrationals on the unit interval.
104
a
A. Saito, S. Ito / Physica D 268 (2014) 100–105
b
Fig. 5. (Color online) Statistical properties obtained from the true orbit of the modified Bernoulli map starting from (p0 , q0 , r0 , s0 ) = (1, 1, 1, −1). (a) The estimated invariant density (dots). The curve represents the density function ρ(x) = x−1 (1 − x)−1 . (b) The estimated power spectral density (dots) plotted on double-logarithmic scales. The line shows 1/f divergence.
map and examine the dependence on n of ln in Fig. 4(b). It is clear that ln increases with n due to the aperiodicity of the true orbits, but, compared with the case for the Bernoulli map, the increase in ln is much slower. For example, l104 = 15,056 for the Bernoulli map, whereas l104 = 1872 for the modified Bernoulli map. Also, ln for the modified Bernoulli map does not show a simple linear increase such as that in Fig. 1(b). That is, ln in Fig. 4(b) has plateaus that correspond to laminar phases, during which the true orbit has an almost constant behavior near x = 0 and x = 1 (cf. Fig. 4(a)). Next, let us consider again the question of whether or not the statistical properties of the generated true orbit are good. The modified Bernoulli map has a σ -finite but infinite invariant measure µ with density dµ = x−1 (1−x)−1 dx on the unit interval [13]. Fig. 5(a) shows the estimated density obtained from the true orbit with length 106 starting from (1, 1, 1, −1) (see the Appendix for details on the density estimation). The estimated density agrees well with that of the modified Bernoulli map. Similarly, Fig. 5(b) shows the estimated power spectral density (periodogram) obtained by performing the discrete Fourier transform of the same true orbit. The estimated power spectral density is in good agreement with the 1/f spectral density of the modified Bernoulli map. 4. Conclusions By using real cubic irrationals having complex conjugates and integer arithmetic, we have introduced a true orbit generation method enabling exact simulations of discrete-time dynamical systems defined by one-dimensional piecewise linear fractional maps with integer coefficients. We have demonstrated that the true orbits generated by our method show clear chaotic and intermittent behaviors for the Bernoulli map and modified Bernoulli map, respectively, in contrast with the case for conventional simulation methods. In addition, for each case, we have confirmed at most linear growth of the number representation, which is a distinctive characteristic of our method of computation. In particular, we have shown that the statistical properties obtained by the generated long true orbits match those of typical orbits of the above two maps in simulations concerning invariant measures and the power spectrum. Expansion of the true orbit generation method to higher dimensional maps, as well as application of our method to, e.g., analysis of noise-induced order phenomena [24] and generation of pseudo-random numbers, will be reported elsewhere.
Appendix. The method for estimating the density of an infinite invariant measure It is known that, unlike the case for dynamical systems preserving a probability measure, it is very difficult to statistically treat those preserving an infinite measure (see, e.g., Refs. [25–27]); for example, time-averaged values obtained from a single orbit do not converge to actual values (i.e., values of integrals with respect to the invariant measure). Nevertheless, as seen below, we can estimate the density of the σ -finite invariant measure for a conservative ergodic map. In what follows, we explain this through the density estimation, shown in Fig. 5(a), of the modified Bernoulli map M. Let µ be the invariant measure on the unit interval X = [0, 1), whose density is ρ(x) = cx−1 (1 − x)−1 with an indefinite constant c ∈ (0, ∞), and let x ∈ X be the point at which we want to estimate the density. Consider a small interval A ⊂ X of length ∆ such that x ∈ A, and also consider some arbitrarily chosen set B ⊂ X , satisfying 0 < µ(A), µ(B) < ∞. By Hopf’s ratio ergodic theorem [28,29], we obtain, for a.e. x0 ∈ X , that
♯ {0 ≤ n < N | M n (x0 ) ∈ A} µ(A) = , N →∞ ♯ {0 ≤ n < N | M n (x0 ) ∈ B} µ(B) lim
where ♯S denotes the number of elements in the set S. Since µ(A) is approximated by ρ(x)∆, we see that if x0 is a typical initial condition, then
ρ(x) ≈
♯ {0 ≤ n < N | M n (x0 ) ∈ A} µ(B) ♯ {0 ≤ n < N | M n (x0 ) ∈ B} ∆
for sufficiently large N.4 The detailed procedure for obtaining the estimated density in Fig. 5(a) is as follows: We first divide X into 25 subintervals Ai = (i/25, (i + 1)/25) (i = 0, 1, 2, . . . , 24) of length ∆ = 1/25. Also, we choose the interval (1/3, 2/3) as B. Notice that µ(B) in Eq. (A.1) depends on the indefinite constant c, but in this study we consider 2/3 c = 1, so it has a definite value: µ(B) = 1/3 x−1 (1 − x)−1 dx =
1.386294 · · ·.5 We then generate the true orbit with length N = 106 starting from x0 represented by (1, 1, 1, −1), and count the number of times the orbit visits each Ai and B. Lastly, using Eq. (A.1), we estimate the density values for the 23 subintervals Ai (i = 1, 2, . . . , 23) satisfying µ(Ai ) < ∞. Notice that we cannot apply Eq. (A.1) to A0 and A24 , since they are sets of infinite measure and beyond the scope of Hopf’s ratio ergodic theorem. Finally, if x0 represented by (1, 1, 1, −1) is a typical initial condition, then the estimated density obtained by the above procedure should approximately coincide with ρ(x) = x−1 (1 − x)−1 . Fig. 5(a) shows a good coincidence between the two, which supports the basic premise of our method, that is, cubic irrationals are typical initial conditions, so true orbits starting from them possess good statistical properties.
4 More generally, for f , g ∈ L+ (µ), we have 1 N −1
f dµ ≈ X
n =0 N −1
Acknowledgments
f (M n (x0 )) g(
Mn
g dµ, X
(x0 ))
n=0
and we can estimate the value of X f dµ from that of X g dµ. 5 In fact, the evaluation of µ(B) is unnecessary if we merely desire to obtain one of the invariant densities: By replacing µ(B) in Eq. (A.1) with some arbitrary positive constant, say 1, we obtain an invariant density ρ(x) = cx−1 (1 − x)−1 with
We thank E. Brooks-Pollock, R.S. MacKay, J. Tamura, and S. Yasutomi for their suggestions. This research was supported by the JST PRESTO program and by a Grant-in-Aid for Young Scientists (B) 21700256 from MEXT of Japan.
(A.1)
c = 1/
2/3 1/3
x−1 (1 − x)−1 dx.
A. Saito, S. Ito / Physica D 268 (2014) 100–105
References [1] E. Ott, Chaos in Dynamical Systems, second ed., Cambridge University Press, Cambridge, 2002. [2] E. Atlee Jackson, Perspectives of Nonlinear Dynamics, Cambridge University Press, Cambridge, 1991. [3] Y. Ueda, The Road to Chaos, second ed., Aerial Press, Santa Cruz, CA, 2001. [4] F. Rannou, Numerical study of discrete plane area-preserving mappings, Astronom. Astrophys. 31 (1974) 289–301. [5] P.M. Binder, R.V. Jensen, Simulating chaotic behavior with finite-state machines, Phys. Rev. A 34 (1986) 4460–4463. [6] C. Beck, G. Roepstorff, Effects of phase space discretization on the long-time behavior of dynamical systems, Physica D 25 (1987) 173–180. [7] D.J.D. Earn, S. Tremaine, Exact numerical studies of Hamiltonian maps: iterating without roundoff error, Physica D 56 (1992) 1–22. [8] K. Kaneko, Spatial period-doubling in open flow, Phys. Lett. A 111 (1985) 321–325. [9] A. Yamaguchi, On the mechanism of spatial bifurcations in the open flow system, Int. J. Bifurcat. Chaos 7 (1997) 1529–1538. [10] J.H. Silverman, The Arithmetic of Dynamical Systems, Springer, New York, 2007. [11] F. Schweiger, Ergodic Theory of Fibred Systems and Metric Number Theory, Clarendon Press, Oxford, 1995. [12] F. Schweiger, Multidimensional Continued Fractions, Oxford University Press, Oxford, 2000. [13] A. Saito, Computational aspects of a modified Bernoulli map, Prog. Theoret. Phys. Supp. 161 (2006) 328–331. [14] F. Vivaldi, Dynamics over irreducible polynomials, Nonlinearity 5 (1992) 941–960.
105
[15] C. Grebogi, L. Poon, T. Sauer, J.A. Yorke, D. Auerbach, Shadowability of chaotic dynamical systems, in: B. Fiedler (Ed.), in: Handbook of Dynamical Systems, vol. 2, Elsevier, Amsterdam, 2002, pp. 313–344. [16] D.V. Anosov, Geodesic flows on closed Riemannian manifolds of negative curvature, Proc. Steklov Inst. Math. 90 (1967) 1–235. [17] R. Bowen, ω-limit sets for axiom A diffeomorphisms, J. Differential Equations 18 (1975) 333–339. [18] C. Grebogi, S.M. Hammel, J.A. Yorke, T. Sauer, Shadowing of physical trajectories in chaotic dynamics: containment and refinement, Phys. Rev. Lett. 65 (1990) 1527–1530. [19] S. Dawson, C. Grebogi, T. Sauer, J.A. Yorke, Obstructions to shadowing when a Lyapunov exponent fluctuates about zero, Phys. Rev. Lett. 73 (1994) 1927–1930. [20] T. Sauer, C. Grebogi, J.A. Yorke, How long do numerical chaotic solutions remain valid? Phys. Rev. Lett. 79 (1997) 59–62. [21] B.L. van der Waerden, Algebra I, Springer, Berlin, 1966. [22] Y. Aizawa, T. Kohyama, Asymptotically non-stationary chaos, Progr. Theoret. Phys. 71 (1984) 847–850. [23] S. Ito, Algorithms with mediant convergents and their metrical theory, Osaka J. Math. 26 (1989) 557–578. [24] K. Matsumoto, I. Tsuda, Noise-induced order, J. Stat. Phys. 31 (1983) 87–106. [25] J. Aaronson, An Introduction to Infinite Ergodic Theory, American Mathematical Society, 1997. [26] T. Akimoto, Y. Aizawa, New aspects of the correlation functions in nonhyperbolic chaotic systems, J. Korean Phys. Soc. 50 (2007) 254–260. [27] T. Akimoto, E. Barkai, Aging generates regular motions in weakly chaotic systems, Phys. Rev. E 87 (2013) 032915. [28] E. Hopf, Ergodentheorie, Springer, Berlin, 1937. [29] R. Zweimüller, Hopf’s ratio ergodic theorem by inducing, Colloq. Math. 101 (2004) 289–292.