Ultradiscrete hard-spring equation and its phase plane analysis

Ultradiscretization enables us to construct a piecewise linear equation which approximates a given subtraction-free difference equation. Recently proposed “ultradiscretization with parity variables” (pUD) can treat an equation with subtraction. However, its solution may have an infinite number of branches under some specific conditions. In this paper, solutions of pUD for the hard-spring equation is investigated. The pUD equation is reinterpreted as the mapping which maps a set on the phase plane to another one, and the behaviour of the solutions are analyzed through approximative transition diagrams.As a result, the infinite branches are translated into a few finite branches and allowed easier understanding.


Introduction
Ultradiscretization [1] is a limiting procedure which constructs a cellular automaton [2] from a given difference equation. Its research arose in the field of integrable systems. It has been reported in many papers that integrable structure in the original difference equation inherits to the resulting CA. Even for some non-integrable systems, it is illustrated that ultradiscretization works well in the sense that a mathematical structure survives through this procedure [3], [4]. Hence, ultradiscretization is expected to be a new and useful tool to analyze nonlinear systems.
However, ultradiscretization has a serious restriction called 'negative difficulty'. That is, the difference equation must be subtraction-free. In addition, its solution should be positive definite. Some methods have been attempted for overcoming this difficulty [5,6]. Ultradiscretization with parity variables (pUD) [7] is one of the such methods. We explain its detail later. This procedure can treat a difference equation with subtraction or non-definite signed solutions. Instead of this advantage, a solution may lose uniqueness under a specific condition. Actually, the solution X n for the specific n is given by the inequality such as X n ≤ A (referred to as an indeterminate solution). This ambiguity is translated as branches to an infinite number of values, or as an 'infinite multi-valued mapping'. The procedure of pUD was applied to some q-Painlevé equations and its special solutions were discussed [8][9][10]. For verifying whether pUD is also useful to study other nonlinear systems, the pUD analogue of the hard-spring equation was investigated in [11]. In this previous research, the solutions for wide range of initial values, and the solutions were classified into some types. For example, the exactly four-periodic solutions, the four-periodic solutions after some time steps, the indeterminate solutions summarizable by bounded diagrams, and the indeterminate solutions with unbounded diagrams are found. These are important results but we need some kind of broad perspective for understanding global structure of these solutions.
In this paper, we again investigate the same pUD equation studied in [11]. We actively utilize the phase plane analysis [12], and reinterpret the pUD equation as the mapping which maps a set on the phase plane to another one. This "approximative" viewpoint enables us to draw the transition diagrams which have such sets as its nodes, in which we see the infinite multi-valued mapping as a finite multi-valued mapping. By these diagram, the structure of indeterminate solution is visually understood. This paper is organized as follows. We review the construction of the pUD equation in Sect. 2. In Sect. 3 following the method developed in [13], we summarize its time evolution to some tables and illustrate it by some figures of the phase plane. Then, we draw the transition diagrams and analyze the structure of the solution orbits. Finally, we give concluding remarks in Sect. 4.

Ultradiscretization with parity variables
We explain the procedure of pUD. We consider an explicit difference equation which does not include subtraction. We write its dependent variable as x n . The subtractionfree setting implies that x n takes positive values only. After transformation upon introduction of a parameter ε > 0, we apply ε log to both sides of the equation and take the limit ε → +0. Then, multiplication, division and addition for the original variable x n are replaced by addition, subtraction and the max operation among new dependent variables X n , respectively.
Here, the formula is used. This is the procedure of ultradiscretization. However, this procedure cannot handle an explicit difference equation with subtraction, whose solution may take negative values. For the sake of solving this issue, ultradiscretization with parity variables (pUD) was proposed.
In the procedure of pUD, we introduce the sign variable ξ n and the "amplitude" variable X n for the dependent variable x n by Moreover, using the function s, we perform the following transformation: After transposing the negative terms to the other side of the equation, we apply ε log to both sides and take the limit ε → +0. Here, we use the formula This is the procedure of pUD. Subtraction is handled by transposed terms. As well as the original ultradiscretization, multiplication, division and addition are replaced by addition, subtraction and the max operation, respectively. The sign is included in the resulting expression as the terms S(±ξ n ).
In general, the obtained equation by pUD takes an implicit form, max(. . . ) = max(. . . ). Therefore, its solution may lose uniqueness. This indefiniteness makes the solution space richer, but may do analysis complicated.

Ultradiscretization with parity variables of the hard-spring equation
The hard-spring equation is given by where x = x(t) and k, l > 0 are constants. This equation is an integrable system since it has the conserved quantity We give an integrable discretization of (8) given in [14]. We write t = nδ, x(t) = x n . Our discrete equation is written as Since (10) has a conserved quantity this equation is an integrable equation. For the purpose of pUD, we introduce the sign variable ξ n , the amplitude variable X n for x n by (3) and new parameters by We substitute (5) into (10), we obtain For simplicity of notation, we write X n+1 , X n , X n−1 as X + , X , X − , respectively. Moreover, we put α i + 2 = α i . If we transfer the negative terms and take the limit ε → +0, we have By considering the four types of combinations of signs ξ n+1 , ξ n , ξ n−1 , (14) is written into some equations without the function S. We use, for simplicity, the notations The rewritten equations are as follows. For ξ n+1 = ξ n = ξ n−1 , for −ξ n+1 = ξ n = ξ n−1 , for ξ n+1 = ξ n = −ξ n−1 , and for ξ n+1 = −ξ n = ξ n−1 , If we consider magnitude relationship between b and d (or b and d ), (20)-(23) are further rewritten as explicit forms summarized in Table 1 with 5 × 5 cells ("NS" means "no solution"). We explain Table 1. Assume that (ξ n , X n ) and (ξ n−1 , X n−1 ) are given. We study how (ξ n+1 , X n+1 ) is determined. We first notice the relationship of the signs ξ n and ξ n−1 . Here, we consider the case where ξ n = ξ n−1 . This case is located at the second and third columns in Table 1. However, we do not know ξ n+1 . Therefore, we try both cases of (I) ξ n+1 = ξ n = ξ n−1 and (II) −ξ n+1 = ξ n = ξ n−1 , which corresponds to (20) and (21), respectively. We further consider three cases for amplitude, that is, Let us firstly study the case (I) and (i). The equation for the amplitude X n+1 is (20) and it has no solution from max X + + a, b = max d . Hence, the (2, 2) cell in Table 1 is filled by "NS". Next, we study the case (II) and (i). The equation is (21) and it is rewritten as X + = b − a from max b = max X + + a, d . Hence, the (2, 3) cell in Table 1 is filled by X + = b − a. Secondly, we consider another relation between amplitudes, that is, (ii) b = d. In this case, both of (20) and (21) are satisfied by the inequality X + max[b, d] − a, which means X + becomes indeterminate. At this time, ξ n+1 is also indeterminate since two cases ξ n+1 = ±1 give the same result. These correspond to the (3, 2) and (3, 3) cells in Table 1. Finally, we consider the case of (iii) b < d. For (I) and (iii), our Eq. (20) is rewritten as X + = d − a from max X + + a, b = max d . This case corresponds to the (4, 2) cell in Table 1. For (II) and (iii), our Eq. (21) has no solution from max b = max X + + a, d . This case corresponds to the (4, 3) cell in Table 1. By studying the cases ξ n = −ξ n−1 and magnitude relation between b and d in a similar manner, we obtain Table 1. Note that the four expressions (20)-(23) can be represented as We go to the next step. Table 1 is written in terms of (16)-(19), which include max operations (except for d of (17)). However, for later analysis, we prefer the expressions without max operations. Hence, we further introduce the cases whether α 2 is positive or 0 or negative. Noting we find that further cases X + α 2 X − + a or X X − + a are necessary for our purpose. We use the following notations for shorter expression: By introducing these cases, we refine Table 1 to Table 2 with 11×6 cells. For example, we notice the (2, 3) cell in Table 1 filled by In the case of α 2 > 0, it becomes which is located at the (2, 4) cell in Table 2. In the case of α 2 = 0, it becomes which is located at the (5, 4) cell in Table 2. In the case of α 2 < 0, we have and this result is located at the (8, 4) cell of Table 2. By continuing similar operations, Table 2 can be created. Finally, we rewrite a of (15) according to the magnitude relation between the terms with X and other terms. Namely, we obtain Therefore, it is found that if we see the boundary of the above classification in Table 2 on the phase plane (X − , X ), it is not a straight line but a broken line at the point If we introduce the above classification and rewrite our conditions C g , C e , C and so on (see (25)), we obtain the following expressions without max operations. That is, if α 2 > 0, and if α 2 0,

Time evolution of the amplitude on the phase plane
In the previous subsection, we have studied the time evolution of (ξ n+1 , X + ) with conditions for (ξ n−1 , X − ) and (ξ n , X ). The time evolution of a discrete dynamical system is often studied by plotting the solution (x n , x n+1 ), n ≥ 0 on one phase plane.
In this subsection, we apply the same method to our pUD system. For this purpose, we define some sets on the phase plane of (X − , X ) (or (X , X + )) as follows. If These sets are shown in the left figure in Fig. 1. Note that the axes in these figure are drawn with its origin at P1 instead of (0, 0). In a similar way, we define the sets for These sets are shown in the left figure in Fig. 2 (The origin is again P1). We refer to a set P1 ∪ L2 ∪ S2 ∪ L3 ∪ S3 ∪ L4 as "cone" in this study.
Similarly, we define sets on the phase plane (X , X + ). The results are shown in the right figure in Figs. 1 and 2, respectively. Note that the axes are drawn with the origin P2. We also note that the sets in the right figure are drawn by displace those in the left figure symmetrically with respect to L1. We refer to the image of the cone mentioned above as the same word "cone." Now, we study where a point on the phase plane of the left figure in Fig. 2 (or Fig. 1) is mapped on the right figure. As the first example, we consider the case of ξ n = ξ n−1 with α 2 > 0 and summarize the result to Table 3. This case corresponds Fig. 1 Sets on phase plane with α 2 ≤ 0 Fig. 2 Sets on phase plane with α 2 > 0 to the cells from (2, 3) to (4, 4) of Table 2. Only the (2, 4) cell has a solution and the three conditions C g , C e and C (see (25)) should be considered. Firstly, under the condition C , which corresponds to the fourth column of Table 3, the time evolution is X + = X − and therefore an original point (X − , X ) is mapped to (X , X + ) = (X , X − ) being symmetry with respect to the line L1. Secondly, under the condition C e , which corresponds to the third column of Table 3, we have the same result as C . Finally, under the condition C g , which corresponds to the second column of Table 3, the time evolution is X + = α 2 + X − a. Hence, we find that the image of this mapping is this break line X + = α 2 + X − a. The upper row in Table 3 illustrates each area corresponding to 3, respectively, and the lower row shows where the original area is mapped. This is how to create Table 3.
Similarly, Table 4, 5, 6, 7 and 8 are obtained. As an important example, we also explain the case of ξ n = ξ n−1 and C e with α 2 > 0 corresponding to the third column of Table 4. In this case, the time evolution is the indeterminate type, X + ≤ X − . We can choose an infinite number of points (X + , X ) as the image of the original point (X − , X ). This situation is illustrated in the lower row, the third column by the "filled" Table 3 Mapping of amplitude for α 2 > 0, ξ n = ξ n−1 Table 4 Mapping of amplitude for α 2 > 0, ξ n = ξ n−1 Table 5 Mapping of amplitude for α 2 = 0, ξ n = ξ n−1 Next, in order to trace the time evolution in one phase plane, we draw our areas for (X − , X ) plane and (X , X + ) plane into one coordinate plane. Then, we should consider "positional relationship" between two cones, one has the edge L2 ∪ L4 ∪ P1 and the other has L5 ∪ L7 ∪ P2. We have three types of relationship: "overlapping," "sticking," and "separating." Fig. 3, 4, 5, 6 and 7 show all positions of two cones. Note that the "overlapping" type appears only for α 2 > 0. The other types appear for any α 2 .
We note that these three types of "separating," "sticking," and "overlapping" are inevitably determined by the magnitude relationship between α 1 and α 2 , which is summarized in Table 9. For example, "overlapping" happens under the condition of α 2 > 0 and α 2 > α 1 and this is given on the second column of Table 9. The other cases, "sticking" and "separating" appear on the third and fourth column of Table 9, respectively. Table 6 Mapping of amplitude for α 2 = 0, ξ n = ξ n−1 Table 7 Mapping of amplitude for α 2 < 0, ξ n = ξ n−1 Table 8 Mapping of amplitude for α 2 < 0, ξ n = ξ n−1 Fig. 4 "sticking" cones with Fig. 6 "sticking" cones with Here, for simplicity, we write +1 and −1 as + and −, respectively. We use the (ξ n−1 e X − ) v.s. (ξ n e X ) coordinates and four copies are placed on each quadrant as shown in Fig. 8. Our cones are also drawn in each quadrant. Note that we here draw only the case of "sticking" for simplicity. Table 9 Relationship between parameters and cone-position "overlapping" "sticking" "separating" For distinguishing the pair of the sign, we represent our sets defined by (27) or (28) with black-and-white shape. That is, in Fig. 8, a pair of signs (+, +) is represented by a white circle, (+, −) by a black circle, (−, −) by a black-painted diamond, and (−, +) by a white-painted diamond. Now, we generally put Z = (ξ n−1 , ξ n , X n−1 , X n ) and investigate the time evolution for a set of assigned initial values. For example, we start from a set of initial values Z 0 = (+, +, X − , X ), (X − , X ) ∈ S1. We point out in advance that the result is summarized by one diagram in Fig. 9. The evolution of the amplitude and the sign can be traced by referring to Tables 3, 4, 5, 6, 7, 8 and 2, respectively. From this initial values, we find the following time evolution for any α 2 : The sign-inverted solution is also obtained easily. These solutions are trivially four-periodic. Here, we reinterpret this situation as follows. This system maps a set S1 with the sign (+, +) to another set S6 with (+, −) (namely, S1 with (+, +) → S6 with (+, −), not S1 with (+, +) → S6 with (+, −)), and then S1 with (−, −) and finally S6 with (−, +). We call this reinterpretation "approximative", which is useful for complicated indeterminate solutions as we see later. We illustrate the above solutions as the approximative transition diagram shown in Fig. 9. We see a simple cycle among four nodes. It is noted that this solution are uniquely determined in any step, in other words, we do not encounter an indeterminate solution. This situation is expressed in Fig. 9 by which two nodes are connected with a single arrow only. This transition is common to any value of α 2 and positioning cases. Looking at Fig. 8, we can see that this transition rotates S1 or S6 or L1 with a pair of sign by 90 degrees clockwise at one time. Figure 10 shows a transition that is common to all positional cases of α 2 < 0. For example, let us start from a point Z Z 0 = (+, +, X − , X ), (X − , X ) ∈ L2. In Table 2, we can choose the (9, 3) or (9,4) cell. This means that the indeterminate solution has appeared here. If we choose the (9, 4) cell and the equal case X + = X − , the sign is determined by the (11, 4) cell as ξ n+1 = − and we have the next point Z Z 1 = (+, −, X , X − ), (X , X − ) ∈ L7. Moreover, since Z Z 1 corresponds to the C case of the (10, 5) cell, we find the after the next point is uniquely determined as Z Z 2 = (−, −, X − , X ), (X − , X ) ∈ L2. However, since the next point Z Z 3 becomes indeterminate again, we need a similar case consideration. From the above (and further calculation), we, for example, obtain the development such as This solution orbit is understood approximatively as L2 with (+, If, for the same initial point Z Z 0 , we choose another cell (9, 3) and the equal case X + = X − , we have a different point Z Z 1 = (+, +, X , X − ), (X , X − ) ∈ L7. After the next point becomes Z Z 2 = (+, −, X − , X ), (X − , X ) ∈ L2 from the case C e of the (8,4) cell. Further the next point uniquely becomes Z Z 3 = (−, −, X , X − ), (X , X − ) ∈ L7. From these and further calculation, we find the solution such as This solution orbit is understood approximatively as L2 with (+, Furthermore, for the same initial point Z Z 0 , we can choose the next point Z Z 1 = (+, +, X , < X − ), (X , < X − ) ∈ S5. Then, after the next point uniquely becomes Z Z 2 = (+, −, < X − , X ), ( < X − , X ) ∈ S2 from the cell (8,4) with C , where < X − means any value Y satisfying Y < X − . Further the next point is uniquely determined as Z Z 3 = (−, −, X , X − ), (X , X − ) ∈ L7 from the cell (10, 5) with C g , which has already appeared above. This solution behaves as which is unique except for Z Z 1 . This solution orbit is understood approximatively as L2 with (+, Investigating all (but finite number of) cases, we obtain the "approximative" transition diagram shown in Fig. 10. In usual meaning, we have encountered infinitely multi-valued mapping when we determine the next point from the initial point. However, by our approximative point of view, we understand it as only four-valued mapping. This situation is expressed in Fig. 10 by which four arrows are left from L2 of white circle and of black diamond, respectively.
Let us study features of the orbits in the diagram of Fig. 10. We see two cyclic transition in this diagram. One is L2 with (+, +) → L7 with (+, −) → L2 with (−, −) → L7 with (−, +) → L2 with (+, +), which have already mentioned. The other is L7 with (+, +) → L2 with (+, −) → L7 with (−, −) → L2 with (−, +) → L7 with (+, +). If we once leave the former cycle at a branching point, we cannot go back to it and the orbit eventually stays the later cycle. Hence, the former one looks "unstable" and the later does "stable". Because of this, it is claimed that the effect of the indeteminate solution is not serious in this case. Figure 11 shows the diagram for "sticking" positional type with α 2 < 0. We find the self-loop at P1 and L4. We again can be summarized infinitely multi-valued mapping as four-valued one. However, we must encounter the indeterminate solutions with an infinite number of times. Hence, this diagram is rather complex than Fig. 10. Figure 12 shows the diagram for the "separating" positional type with α 2 < 0. In a similar way to Fig. 10, we see a "stable" cycle and an "unstable" one. We have obtained all diagrams for α 2 < 0 by Figs. 9, 10, 11 and 12. Note that we see all sets (P1, L1, S1, and so on, with four signs respectively) only once. Hence, we find that we have traced all orbits on the phase plane for α 2 < 0. We give some results for α 2 > 0. Figure 13 shows the approximative transition diagram which is common to all of the positional types with α 2 > 0. We see a "stable" cycle and an "unstable" one. Figure 14 shows the diagram which corresponds to the "sticking" positional type with α 2 > 0. Figure 15 shows the diagram which corresponds to the "separating" positional type with α 2 > 0. We see a "stable" cycle and an "unstable" one. We remain the "overlapping" positional type. Although we have studied it, its explanation is rather   complicated (We shall report it in forthcoming paper). If we consider Figure 9, 13, 14 and 15 and the "overlapping" type, we see all sets only once and therefore we find that we have traced all orbits on the phase plane for α 2 > 0.
In closing this section, we give Table 10 which summarizes the relationship among the sign of α 2 , type of the positional relationship, and the transition diagrams.

Concluding remarks
In [11], the ultradiscrete analogue of the hard-spring equation with parity variables was presented, and its solutions for various (but, not all) initial values were classified into some types. It was quite helpful to regard an infinite number of values as one point on the solution orbit for studying the indeterminate solution.
In this study, we focused on the same equation as [11] and we have applied further approximative method depeloped in [13], that is, reinterpreting the pUD system as the mapping which maps a set on the phase plane to the other. Moreover, the positional relationship of "cones" on the phase plane closely effects the solution behaviour. This relation is determined only by the magnitude relationship of the system parameters α 1 and α 2 , which was not pointed out in [11]. From these viewpoints, certain new perceptions have been obtained for the solutions of the pUD system as follows. Because all domains on the phase plane appeared, all initial values have been examined, which completely cover those of [11]. As a result, the condition under which an indeteminate solution does not appear has been clarified, that is, to start from a point on the phase plane included in Fig. 9. When the positional relationship of cones is "separating" or "sticking", we have captured an infinite number of branches in usual sense as a finite number of ones. This result implies that the pUD equation is still a kind of cellular automaton, because it has a finite number of "states", but its dynamics is not always "auto" because of non-uniqueness.
Furthermore, the approximative way enables us to understand the earlier work from a broader point of view. For example, in [11], the existence of exactly fourperiodic solution was pointed out. However, by our analysis, such solutions appear in several transition diagrams. Namely, all orbits of Fig. 9 and the "stable" orbit of Figs. 10, 12, 13, and 15 . Note that they were not distinguished in [11]. This study has clarified that they belong to different transition diagrams each other. We have also found that almost solution orbits (as sets) show cyclic behaviour on the phase plane with help of the transition diagrams.
Our method reported in this paper enables us to easily compare some pUD equations each other. For example, it is possible to compare the pUD analogues of several difference equations for the same differential equation, which may give a new contribution to studying nonlinear systems. It is also an interesting problem to investigate differences between integrable and non-integrable pUD systems by this approximative method.
We comment one issue in [11], that is, how to analyze the solution with "diagram of unbounded type," which seems to need an infinite number of cases for tracing the orbit. We have found that we encounter this type of solution if and only if the the positional relationship is "overlapping" and understood why such many cases appear. We shall report its details in the forthcoming paper.