- Email: [email protected]
A method of point mapping under cell reference for global analysis of nonlinear dynamical systems
16 may 1994 PHYSICS PhysicsLettersA
ELSEVIER
LETTERS
A
188 (1994) 137-145
A method of point mapping under cell reference for global analysis of nonlinear dynamical systems Jun Jiang, Jian-xue Xu Institute of EngineeringMechanics,Xi’an Jiaotong University,Xi’an, China Received 29 October 1993; revised manuscript received 18 January 1994; accepted for publication 9 March 1994 Communicated by A.P. Fordy
In this paper a criterion for preserving the system characteristics for numerical global analysis is introduced. By setting a cell coordinate in state space as a reference, a method of point mapping with this criterion can be implemented efftciently, which preserves system accuracy while reducing a lot the computation.
1. Introduction In the analysis of many finite dimensional nonlinear dynamical systems, one often casts them in the form x(n+l)=G(x(n)),
(1)
where x(n) is an N-dimensional state vector, the vector mapping G maps x(n) to its image point x( n + 1) after a time duration r and x(n) means x( nr). As different initial conditions can lead to different responses for point mapping ( 1) , an individual trajectory must be allowed to settle on a steady state response for every initial point within a chosen region of state space in order to fully characterize the system’s response. Clearly, this will quickly become extremely costly and time consuming. In an effort to relieve this difficulty, Hsu [ 1 ] introduced the cell mapping method, which performs global analysis on a discrete cell state space constructed by overlaying an array of cells on a given region of state space. Simple cell mapping (SCM) [ 21 and generalized cell mapping [ 3,4] are devised under this framework. In a cell mapping the finite cell size places an immediate limit on both the accuracy of the system response and the resolution of the fine structure that is possible. An approach that attempts to overcome the drawback of cell mapping was proposed by Tongue [ 5 ] under the name of interpolated cell mapping (ICM). In this technique interpolated trajectories are created on the basis of the finite grid points in a chosen region of state space and their image points under one step of the point mapping. In fact, ICM defines an approximated mapping of G and the accuracy is restricted by the distance between the grid points. The point mapping method can achieve highly accurate results often referred to as exact solutions, but it is not economic in its implementation. What can we do to improve the point mapping method: This question provides basically the motivation for the study of the present paper. Elsevier Science B.V. SSDZO375-9601(94)00204-3
138
J. Jiang, J. Xu /Physicsbt&rsA
188 (1994) 137-145
2. Criterion of the persistence of system characteristics
Theoretically, a trajectory in state space cannot intersect another trajectory except at a steady state orbit. From this point of view, in practical computation of a global analysis by the point mapping method, each of a finite number of initial points is individually used as a start of a trajectory. Once the steady state orbit that the particular initial point leads to is determined, the global character of that point in state space is marked. But no intersection of trajectories except at a steady state orbit does not mean that two closely located trajectories must have different dynamical characters. In fact, we are more interested in the final dynamical characters of a particular point in state space than in the details of its evolution. In a practical analysis, we can suppose that two trajectories will reside at the same steady state motion if they approach each other within a distance much smaller than the distance of their initial positions. If all the initial points keep definite distances away from each other globally, just as is achieved by the point mapping method, a distance much smaller than the minimal distance between the initial points can be set to do the same work as above in the case of two trajectories. From the above point of view, we propose a criterion of persistence of system characteristics in the global numerical analysis of the dynamics. By a criterion of persistence of system characteristics or a criterion of persistence we mean the following: If a trajectory comes close to an existing trajectory, whose dynamical property is known, at a distance less than the distance set in the beginning by a numerical analyst, then this trajectory will be regarded to have the same dynamical character as that of the existing trajectory. The distance we set is calledproposeddistance and denoted as di,. If an exact result can be achieved under the use of a proposed distance in this criterion, this proposed distance will be denoted as character preserving proposed distance. Remarks. ( 1) Although two arbitrarily closely located points on the fractal basin boundaries may have quite distinct global properties, the authors think, from the point of view of global numerical analysis, that it is not generic that two trajectories, belonging to two distinct basins of attraction and with their initial conditions being a finite distance away from each other, could approach each other infinitesimally close, or we think that two trajectories of this kind can at best approach each other to some extent, namely, they will always be away from each other a finite distance (see the examples in Section 4 below). (2) In the numerically calculated point mapping method, there really exists a criterion of persistence. For instance, when two trajectories come very close within a distance less than the round-off error of the computer, they will be taken as the same even if their “real” dynamics may be quite different. Analogous, the criterion, performed by the inherent features of the computer, also has a so-called proposed distance, which is the roundoff error of the computer and is unchangeable. We shall call it the passive criterion of persistence. (3) Obviously, the character preserving the proposed distance is dependent on the discretization or the distances between the initial points and the system studied. In almost all of our computations in the global numerical analysis, there might be an infinite number of character preserving proposed distances in determination of the dynamics of a given system with a definite discretization. The smallest character preserving proposed distance is, of course, the round-off error of the computer, as used passively in the numerically calculated point mapping method. But it is impossible to find the largest character preserving proposed distance, which will bring about the relatively most efficient computation.
Using such a criterion of presistence, we have more choice regarding the system accuracy and computation work. As shown in the applications, exact results can be obtained and computation is greatly reduced if a character preserving proposed distance is adopted in the analysis. Unfortunately, the point mapping method currently used is not improved if a new criterion, like the one of system characterstic persistence, is added beside the already existing criterion of periodicity. For instance, each newly generated point on an extended trajectory must be compared by the criterion of persistence with all points on the previously calculated trajectories in order
J. Jiang, J. Xu /Physics LettersA 188 (1994) 137-145
139
to decide whether it might be thought to have the same character as one of the previously calculated trajectories. This increases the computation time as well as the memory required a great deal. To overcome these difficulties, a new method of point mapping under cell reference is proposed in the next section.
3. Method of point mapping under cell reference To make the point mapping with the criterion of persistence of system characteristics implementable, we try to build a reference which is much easier to be dealt with and can provide guidance on whether and with which point it is necessary to judge the newly generated point of processing trajectory by the criterion. It is fortunate that a cell coordinate space, developed by Hsu [ 1,2], to cast a continuous state space into a discrete cell space for cell-to-cell mapping is readily available and can act as the reference coordinates for our purpose. Let Xi, i= 1, 2, .... N, be the state variables of the state space. To construct a cell coordinate space, divide the coordinate axis of a state variable xi into a countable infinity of intervals of uniform size hi. The interval zi along the xl-axis is defined to be one which contains all xi satisfying (zi-&)hi
(zi+f)hi
a
(2)
As defined in Ref. [ 11, Zi is an integer. A N-tuple Zi, i= 1, 2, .... N, is then called a cell vector of the state space denoted by Z. In this way a cell coordinate system ZN is created. For simplicity, condition (2) to determine in which cell a given point is located is rewritten in a form of a mapping: z=C(x),
zcZN, xeRN,
(3)
where C: RN+ZN. This mapping means that for a point in state space a cell in the cell coordinate system can be found corresponding to it. Since a grid of initial points is only set in a chosen region in the state space, the cells in the corresponding part of the cell coordinates are finite, say N,. Moreover, it is not difficult to identify these N, cells sequentially in a one-dimensional way and, therefore, label them by positive integer 1,2, .... N, or express it in a mapping c=N,(z),
ZEZN,
(4)
where NU:ZN-+S, S= {1,2, .... N,}, c is an integer belonging to S. By now a point x in state space has been related to a positive integer belonging to S. Below, we always refer to a point or a trajectory in a Euclidean state space and refer to a cell or a sequence as the one from the scope of reference cell coordinates on the same state space. Before this cell coordinate system can be used as a reference, we must assign each of the cells in S with an identijkation parameter Id(c) and an identification point Y(c, N). We call the vector R(c)={Id(c), Y(c, N)} an identification vector of cell c. In implementation, Id(c) is given by a one-dimensional integer array of size N, and Y( c, N) by a two-dimensional real-valued array in size of N, x N. Similar to the group number in SCM [ I], each of the cells in S has four possible states identified by the value of Id(c) . To be in the first state the cells have not been passed by any trajectory and are characterized by Id(c) = 0, called virgin cells. The cells in the second state are denoted as cells underprocessing, which are only passed by a trajectory currently being processed and identified by Id(c) = - n,, where n, is a positive integer larger than the maximal number of attractors thought to be captured in the analysis. Cells in the third state are those which have been passed by some trajectories whose global properties have been determined and are characterized by Id(c) = n, 0 < n < n,, and this number, called a global characteristic symbolic number, indicates that a trajectory (among others passed through this cell) will approach the nth attractor in this analysis. These kind of cells are called processed cells. Cells in a fourth state are processed cells before they are passed by a currently evolving trajectory to become cells called processed cells under processing characterized by Id(c) = - n. Y( c, N) is just a
J. Jiang, J. Xu /Physics LettersA 188 (1994) 137-145
140
representation point of the corresponding state of the cell. For instance, a virgin cell will have no identification point, Y( c, N) = 0.0; a cell under processing chooses the most recently passed point of the current calculating trajectory as its Y(c, N); a processed cell or a processed cell under processing takes a point within it, whose character is determined, as its Y( c, N). By now we have finished the work of establishing a reference, so a point mapping with this criterion of persistence can be carried out efficiently in the following way. For a given initial point, x0, one first determines in which cell it is located through mapping (3) and (4), say co, the corresponding identification vector R ( co) = {Id ( co), Y( co, N) } is readily available, then updates cell co to a corresponding cell under processing by reassigning Id ( co) and Y( co, N) as defined above, finally brings x0 into ( 1) to get its image point G’ (x0). Repeating the same procedure, one gets a trajectory and a sequence of cells with corresponding identification vectors. This procedure is illustrated by the following diagram, -
x0
G’bo)
d
1( co
I
R;coI~Nco)
Cl
&d~Ncd
G2(xo) d
-
I
c2 d R(c2)%~2)
I
...-
G”(x,) w!
-
GI
... (5)
I &~~Wrn~
U in the above diagram stands for the update procedure to make a cell, determined by a newly generated point G’(x,) =O. To become a cell under processing in case x does not satisfy one of the following conditions to end the current processing trajectory. At each step in generating (5 ) there are four possibilities (the initial step which is a special case of the other step with only two possible cases): ( 1) The newly generated point G’(x,) is in the cell Ciwhich is such that Id( Cl)= 0 indicating that the cell is a virgin cell. In this case we continue forward to locate the next point G’+ ’ (x0) in the processing trajectory. However, before doing that, we first set Id(ci) = - n, and Y( Ci,N) = G’ ( x0) in order to indicate that cell Ciis no longer a virgin cell but a cell under processing. (2) The newly generated point Gi(xo) =x is found to be in a cell which has a positive integer as its identification parameter. This means that this cell has been passed by the previous processing trajectories, whose global characters have been determined. Thus, we will get Id( ci) = n and Y( Ci,N) =y. The global property of point y on one of the previous processing trajectories is that it lies in the basin of attraction of the nth attractor in this analysis. In this case we will judge, by the criterion of persistence, if x approaches y close enough so that its global character as well as that of complete processing trajectory can be determined. If x and y satisfy Ilx-Yll
(6)
the current processing trajectory is terminated at this point. Thus, the global character of all the points in the presently processing trajectory is easily determined. Obviously, all points of this trajectory will be in the same basin of attraction as that of pointy. Now, the points of this trajectory can be stored on a disk together with their global characteristic symbolic number n. Before we go back to pick the next initial point to begin a new trajectory, we must deal with all the identification vectors accompanying the just finished trajectory in the following manner: reassign Id ( cj) = n for the identification vectors with Id (Cj) = - IZ,, j= 0, 1, 2, .... i, and Id (cj) = n’ for the identification vectors with Id (c,) = -n’, n’# n and keep Y(cj, N) unchanged. If criterion (6) is not met, reassign Id (ci) = - n and keep Y( ci, N) unchanged, then go forward to determine the next step image point of x0, Gi+‘(xo). (3) The newly generated point Gi(xo) =x’ is just in a cell under processing with Id( Ci)= - nnr and Y( ci, N) =y’. This indicates that a current processing trajectory exhibits periodicity on the scale of a cell. To confirm if periodic motion is contained in the present trajectory, we check by the criterion of periodicity just like that in the point mapping method, i.e. that 11x’-Y’ll
=GE .
(7)
J. Jiang, J. Xu /PhysicsLettersA
ISg(1994) 137-145
141
If criterion (7) is met, a new periodic solution with a new sequential number, say n, is found and the current trajectory is terminated. All points of the trajectory are assigned by this new global characteristic symbolic number n, which is larger than those of previous trajectories. All points on this trajectory will be written out. The identification vectors are processed just like that in case (2), when a trajectory is terminated. Otherwise, set = - n,,, unchanged to continue on the next step mapping to get G’+ ’ (x0). Y(c,,N)=~‘andkeepId(c~) (4) The cell, Cl, corresponding to a newly generated point Gi(xO) =x” is a processed cell under processing with Id(ci) = -n and Y(ci, N) =y”. This indicates that the presently processing trajectory comes at least two times into this processed cell. In this case we have to judge with criteria of both ( 7 ) and ( 6 ) . Obviously, criterion (6) can be performed directly like that in case (2) but uses x” instead of x in criterion (6). To carry criterion (7) out, one would have to search back in the current trajectory to find a point which most recently passed a cell, Ci,other than x”. After finding this point replace y’ in criterion (7) by this point and do like in case ( 3 ) . Once the global characteristic symbolic numbers have been assigned to each point of the current processing trajectory an each identification vector in this procedure has been changed to a processed one, the work on this processing trajectory is finished and we go back to find a new initial point and start a new processing trajectory. After all the initial points we set are dealt with, the global dynamical behavior of a given nonlinear dynamical system is known. Remark. ( 1) In this method the cell coordinate system only serves as a reference in the analysis and the size of the cells
or the number of cells used is, in fact, not decided by the number of initial points used. This means that a cell might contain as many initial points as possible. The precision of the analysis of the present method is predominantly controlled by the proposed distance d, in the computation. (2) Usually, a maximal length of a trajectory in practical computations is set by people. If a trajectory like the one in (5) reaches this limit before a periodic attractor or a previously processed trajectory is found to go on, a nonperiodic motion is thought to be detected and we perform just like in case (3), where a periodic attractor was found. (3 ) In the computation, a trajectory may map outside the chosen region, where no identification vector is defined on the cell in that part. In this case we apply a mapping G,,,( ) =G“( ), k is a positive integer to the point outside the chosen region. If this trajectory is led back to the chosen region by this mapping, the same analysis is performed as introduced above. Otherwise, this trajectory will be thought undeterminable and is terminated with no global characteristic symbolic number assigned.
4. Examples of applications
In order to illustrate the point mapping under the cell reference method, a series of computations has been performed with different proposed distances to apply it to the system described by the forced van der Pol equation, Z+5.10(xZ -I)~+x=ll.Ocos(3.1416t).
(8)
The region of state space to be examined is defined by -3.Odx,<2.0,
-6.O~Rg4.0.
This particular system supports a steady-state period-three motion as well as a steady-state period-five motion. Moreover, the system exhibits fmctal basins of attraction (see Ref. [ 71). We choose this system rather than a system which can support nonperiodic motion so as to make a reasonable comparison between the computation work of the point mapping method and the present method. If a system supporting chaotic motion is used, the point mapping method will extend all trajectories leading to chaos to their maximal length before they are ter-
142
J. Jiang, J. Xu /Physics Letters A I88 (1994) 137-145
minated. This maximal length of a trajectory in numerical analysis is set differently by different people. All computations were performed on a LINK 386/DX33 microcomputer using an integration step length of 0.04. We overlay a grid of 100x 100 points on the chosen region defined above. Thus, the minimal distances between two initial points in the x-direction and i-direction are d,= 0.06 and d*=O. 1, respectively. They will be denoted as grid scale below. In our method cell coordinates serving as a reference must also be built in the state space. To make the dimension of the cells be also in 100x 100 within the chosen region as the grid of initial points, a discretization of L\x=O.O6 and A.&O. 1 is used. This implies that each cell contains only one initial
(4
4.0
2.0 0.0 X -2.0
-4.0
-0.0 -< X
2.0 .
.
%
0.0 ii -2.0
-4.0
. -6.0
1.0
-2.0
-1.0
0.0 X
1.0
2.0
3:0
Fig. 1. Basins of attraction represented by the states of 100 x 100 initial points and attractors:Point mapping method; (m) period-three motion, (m) period-five motion.
J. Jiang, J. Xu /Physics L&ten A 188 (1994) 137-145
143
point. But as stated in the remarks of Section 3, the cell array is not necessary of the same dimension as that of initial conditions. Fig. 1 represents an “exact” solution in the sense that each trajectory for a particular initial point was continuously mapped until the response goes on either of the steady-state motion. No criterion of persistence is used. Figs. 2-4 illustrate the present method for different proposed distances used in the criterion of persistence. The CPU time required and the final states of all 100 x 100 initial points were recorded in all cases for comparison in details, see Table 1. Fig. la shows that points in the state space are divided into two main types: Those that approach a periodthree orbit (black) and those that lead to period-five motion (gray). The actual Poincare maps of the attractors are illustrated in Fig. lb. Fig. 2 shows the results of using the present method with a proposed distance, whose components are G&= 0.006 and d&=0.01, that is, the proposed distance in the criterion of persistence is one-tenth of the grid scale. The match between this plot with the exact solution, Fig. la, is good but with discernible differences in some parts of the boundaries of the period-three and period-five motion. We will make a comparison in more detail. Table 1 lists the details of the relative accuracy of basins of attraction over exact solution as well as the relative CPU time used in these computations. It is shown quantitatively that the present method in this proposed distance achieves an accuracy of 83.6W in the basin of attraction of period-three motion and an accuracy of 89.6W to that of period-five motion using less than $ of the original CPU time. Figs. 3 and 4 display the results of applying the present method with proposed distances of one-hundredth of the grid scale and one-thousandth of it, respectively. The matches between the two plots and the exact one in Fig. 1 are extremely good or, to be more precisely, indistinguishable. From Table 1 one can find that the result in Fig. 4 is the exact one, which is obtained using only one-fourth of the CPU time used by the point mapping method. If one thinks that the accuracy of the result in Fig. 3 is enough for one’s purpose, then a further reduction of the computation can be made as shown in Table 1. In all the above computations, the only character preserving the proposed distance is the one-thousandth of the grid scale and proposed distances less than it will be, undoubtfully, the character preserving ones. The attractors determined by the present method are the exact
4.0
2.a 0.0 X -2.0
-4.0
-6.0
Fig. 2. Basins of attraction represented by the states of 100x 100 initial points: Present method with the proposed distance being one tenth of grid scale; (m) period-three motion; ( n ) period-five motion.
144
J. Jiang, J. Xu /Physics Letters A 188 (1994) 137-145
4.0
2.0
0.0 ir -2.0
-4.0
-6.0
X
X
Fig. 3. Basins of attraction represented by the states of 100x 100 initial points: Present method with the proposed distance being one hundredth of the grid scale; (m) for the period-three motion; (I) for the period-rive motion. Fig. 4. Basins of attraction represented by the states of 100x 100 initial points: Present method with the proposed distance being one thousandth of the grid scale; (m) period-three motion; (m) period-five motion. Table 1 Comparison of relative accuracy and CPU used point mapping method basin of attraction of P-3 basin of attraction of P-5 CPU time used
100% 100% 100
Point mapping under cell reference l/10 of G.S.
l/100 ofG.S
l/1000 of G.S.
83.6% 89.6% 7.50
99.1% 99.15% 16.7
100% 100% 23.2
Note: G.S. means the grid scale defined above.
ones as in Fig. lb for all cases since the criterion of periodicity used is the same as that in the point mapping method. The usefulness of the point mapping under reference for global analysis is well demonstrated by the above examples and by comparison with the exact solutions. In fact, the present method is a flexible method for global analysis with controllable accuracy and computational work.
5. Conclusions After analyzing the inefficiency of the originally introduced point mapping method in global analysis, a criterion of persistence of the system characteristics is first proposed in this paper. The necessity to explicitly propose the criterion of persistence lies in the following two aspects: (a) To make people pay attention to the reliability
J. Jiang, J. Xu /Physics LettersA 188 (1994) 137-145
145
of the results of the point mapping method, especially in the case to determine the fractal dimension of the basins of attraction, where characters of two closely located points (the distance between them may be of the order of lo-“) are needed to be identified. (b) To make it possible to find an efficient way to get exactly the same results as those of the point mapping method, as the work done in the present paper. To make the point mapping method together with this criterion implementable, a cell coordinate system developed by Hsu is adopted to establish another reference coordinate system in the state space to assist the global analysis. It is shown by direct applications that the present method is a flexible method in global analysis, which can control the achievement of the system accuracy and the expense of the computation. Moreover, the present method is much more efficient than the point mapping method formerly used in global analysis after the comparison of the CPU time used to get exactly the same results.
Acknowledgement This work was supported by the National Natural Science Foundation of China and the Provincial Natural Science Foundation of Shaanxi, China.
References [ 1 ] C.S. [2] C.S. [3] C.S. [4] C.S. [S] B.H. [6] B.H. [7] S.W.
Hsu, ASME J. Appl. Mech. 47 ( 1980) 931. Hsu and R.S. Guttalu, ASME J. Appl. Mech. 47 ( 1980) 940. Hsu, R.S. Guttalu and W.H. Zhu, ASME J. Appl. Mech. 49 (1982) 885. Hsu,ASME J. Appl. Mech. 49 (1982) 895. Tongue and K. Gu, ASME J. Appl. Mech. 55 (1988) 461. Tongue, Physica D 28 (1987) 401. McDonald, C. Grebogi, E. Ott and J.A. York, PhysicaD 17 (1985) 125.
ELSEVIER
LETTERS
A
188 (1994) 137-145
A method of point mapping under cell reference for global analysis of nonlinear dynamical systems Jun Jiang, Jian-xue Xu Institute of EngineeringMechanics,Xi’an Jiaotong University,Xi’an, China Received 29 October 1993; revised manuscript received 18 January 1994; accepted for publication 9 March 1994 Communicated by A.P. Fordy
In this paper a criterion for preserving the system characteristics for numerical global analysis is introduced. By setting a cell coordinate in state space as a reference, a method of point mapping with this criterion can be implemented efftciently, which preserves system accuracy while reducing a lot the computation.
1. Introduction In the analysis of many finite dimensional nonlinear dynamical systems, one often casts them in the form x(n+l)=G(x(n)),
(1)
where x(n) is an N-dimensional state vector, the vector mapping G maps x(n) to its image point x( n + 1) after a time duration r and x(n) means x( nr). As different initial conditions can lead to different responses for point mapping ( 1) , an individual trajectory must be allowed to settle on a steady state response for every initial point within a chosen region of state space in order to fully characterize the system’s response. Clearly, this will quickly become extremely costly and time consuming. In an effort to relieve this difficulty, Hsu [ 1 ] introduced the cell mapping method, which performs global analysis on a discrete cell state space constructed by overlaying an array of cells on a given region of state space. Simple cell mapping (SCM) [ 21 and generalized cell mapping [ 3,4] are devised under this framework. In a cell mapping the finite cell size places an immediate limit on both the accuracy of the system response and the resolution of the fine structure that is possible. An approach that attempts to overcome the drawback of cell mapping was proposed by Tongue [ 5 ] under the name of interpolated cell mapping (ICM). In this technique interpolated trajectories are created on the basis of the finite grid points in a chosen region of state space and their image points under one step of the point mapping. In fact, ICM defines an approximated mapping of G and the accuracy is restricted by the distance between the grid points. The point mapping method can achieve highly accurate results often referred to as exact solutions, but it is not economic in its implementation. What can we do to improve the point mapping method: This question provides basically the motivation for the study of the present paper. Elsevier Science B.V. SSDZO375-9601(94)00204-3
138
J. Jiang, J. Xu /Physicsbt&rsA
188 (1994) 137-145
2. Criterion of the persistence of system characteristics
Theoretically, a trajectory in state space cannot intersect another trajectory except at a steady state orbit. From this point of view, in practical computation of a global analysis by the point mapping method, each of a finite number of initial points is individually used as a start of a trajectory. Once the steady state orbit that the particular initial point leads to is determined, the global character of that point in state space is marked. But no intersection of trajectories except at a steady state orbit does not mean that two closely located trajectories must have different dynamical characters. In fact, we are more interested in the final dynamical characters of a particular point in state space than in the details of its evolution. In a practical analysis, we can suppose that two trajectories will reside at the same steady state motion if they approach each other within a distance much smaller than the distance of their initial positions. If all the initial points keep definite distances away from each other globally, just as is achieved by the point mapping method, a distance much smaller than the minimal distance between the initial points can be set to do the same work as above in the case of two trajectories. From the above point of view, we propose a criterion of persistence of system characteristics in the global numerical analysis of the dynamics. By a criterion of persistence of system characteristics or a criterion of persistence we mean the following: If a trajectory comes close to an existing trajectory, whose dynamical property is known, at a distance less than the distance set in the beginning by a numerical analyst, then this trajectory will be regarded to have the same dynamical character as that of the existing trajectory. The distance we set is calledproposeddistance and denoted as di,. If an exact result can be achieved under the use of a proposed distance in this criterion, this proposed distance will be denoted as character preserving proposed distance. Remarks. ( 1) Although two arbitrarily closely located points on the fractal basin boundaries may have quite distinct global properties, the authors think, from the point of view of global numerical analysis, that it is not generic that two trajectories, belonging to two distinct basins of attraction and with their initial conditions being a finite distance away from each other, could approach each other infinitesimally close, or we think that two trajectories of this kind can at best approach each other to some extent, namely, they will always be away from each other a finite distance (see the examples in Section 4 below). (2) In the numerically calculated point mapping method, there really exists a criterion of persistence. For instance, when two trajectories come very close within a distance less than the round-off error of the computer, they will be taken as the same even if their “real” dynamics may be quite different. Analogous, the criterion, performed by the inherent features of the computer, also has a so-called proposed distance, which is the roundoff error of the computer and is unchangeable. We shall call it the passive criterion of persistence. (3) Obviously, the character preserving the proposed distance is dependent on the discretization or the distances between the initial points and the system studied. In almost all of our computations in the global numerical analysis, there might be an infinite number of character preserving proposed distances in determination of the dynamics of a given system with a definite discretization. The smallest character preserving proposed distance is, of course, the round-off error of the computer, as used passively in the numerically calculated point mapping method. But it is impossible to find the largest character preserving proposed distance, which will bring about the relatively most efficient computation.
Using such a criterion of presistence, we have more choice regarding the system accuracy and computation work. As shown in the applications, exact results can be obtained and computation is greatly reduced if a character preserving proposed distance is adopted in the analysis. Unfortunately, the point mapping method currently used is not improved if a new criterion, like the one of system characterstic persistence, is added beside the already existing criterion of periodicity. For instance, each newly generated point on an extended trajectory must be compared by the criterion of persistence with all points on the previously calculated trajectories in order
J. Jiang, J. Xu /Physics LettersA 188 (1994) 137-145
139
to decide whether it might be thought to have the same character as one of the previously calculated trajectories. This increases the computation time as well as the memory required a great deal. To overcome these difficulties, a new method of point mapping under cell reference is proposed in the next section.
3. Method of point mapping under cell reference To make the point mapping with the criterion of persistence of system characteristics implementable, we try to build a reference which is much easier to be dealt with and can provide guidance on whether and with which point it is necessary to judge the newly generated point of processing trajectory by the criterion. It is fortunate that a cell coordinate space, developed by Hsu [ 1,2], to cast a continuous state space into a discrete cell space for cell-to-cell mapping is readily available and can act as the reference coordinates for our purpose. Let Xi, i= 1, 2, .... N, be the state variables of the state space. To construct a cell coordinate space, divide the coordinate axis of a state variable xi into a countable infinity of intervals of uniform size hi. The interval zi along the xl-axis is defined to be one which contains all xi satisfying (zi-&)hi
(zi+f)hi
a
(2)
As defined in Ref. [ 11, Zi is an integer. A N-tuple Zi, i= 1, 2, .... N, is then called a cell vector of the state space denoted by Z. In this way a cell coordinate system ZN is created. For simplicity, condition (2) to determine in which cell a given point is located is rewritten in a form of a mapping: z=C(x),
zcZN, xeRN,
(3)
where C: RN+ZN. This mapping means that for a point in state space a cell in the cell coordinate system can be found corresponding to it. Since a grid of initial points is only set in a chosen region in the state space, the cells in the corresponding part of the cell coordinates are finite, say N,. Moreover, it is not difficult to identify these N, cells sequentially in a one-dimensional way and, therefore, label them by positive integer 1,2, .... N, or express it in a mapping c=N,(z),
ZEZN,
(4)
where NU:ZN-+S, S= {1,2, .... N,}, c is an integer belonging to S. By now a point x in state space has been related to a positive integer belonging to S. Below, we always refer to a point or a trajectory in a Euclidean state space and refer to a cell or a sequence as the one from the scope of reference cell coordinates on the same state space. Before this cell coordinate system can be used as a reference, we must assign each of the cells in S with an identijkation parameter Id(c) and an identification point Y(c, N). We call the vector R(c)={Id(c), Y(c, N)} an identification vector of cell c. In implementation, Id(c) is given by a one-dimensional integer array of size N, and Y( c, N) by a two-dimensional real-valued array in size of N, x N. Similar to the group number in SCM [ I], each of the cells in S has four possible states identified by the value of Id(c) . To be in the first state the cells have not been passed by any trajectory and are characterized by Id(c) = 0, called virgin cells. The cells in the second state are denoted as cells underprocessing, which are only passed by a trajectory currently being processed and identified by Id(c) = - n,, where n, is a positive integer larger than the maximal number of attractors thought to be captured in the analysis. Cells in the third state are those which have been passed by some trajectories whose global properties have been determined and are characterized by Id(c) = n, 0 < n < n,, and this number, called a global characteristic symbolic number, indicates that a trajectory (among others passed through this cell) will approach the nth attractor in this analysis. These kind of cells are called processed cells. Cells in a fourth state are processed cells before they are passed by a currently evolving trajectory to become cells called processed cells under processing characterized by Id(c) = - n. Y( c, N) is just a
J. Jiang, J. Xu /Physics LettersA 188 (1994) 137-145
140
representation point of the corresponding state of the cell. For instance, a virgin cell will have no identification point, Y( c, N) = 0.0; a cell under processing chooses the most recently passed point of the current calculating trajectory as its Y(c, N); a processed cell or a processed cell under processing takes a point within it, whose character is determined, as its Y( c, N). By now we have finished the work of establishing a reference, so a point mapping with this criterion of persistence can be carried out efficiently in the following way. For a given initial point, x0, one first determines in which cell it is located through mapping (3) and (4), say co, the corresponding identification vector R ( co) = {Id ( co), Y( co, N) } is readily available, then updates cell co to a corresponding cell under processing by reassigning Id ( co) and Y( co, N) as defined above, finally brings x0 into ( 1) to get its image point G’ (x0). Repeating the same procedure, one gets a trajectory and a sequence of cells with corresponding identification vectors. This procedure is illustrated by the following diagram, -
x0
G’bo)
d
1( co
I
R;coI~Nco)
Cl
&d~Ncd
G2(xo) d
-
I
c2 d R(c2)%~2)
I
...-
G”(x,) w!
-
GI
... (5)
I &~~Wrn~
U in the above diagram stands for the update procedure to make a cell, determined by a newly generated point G’(x,) =O. To become a cell under processing in case x does not satisfy one of the following conditions to end the current processing trajectory. At each step in generating (5 ) there are four possibilities (the initial step which is a special case of the other step with only two possible cases): ( 1) The newly generated point G’(x,) is in the cell Ciwhich is such that Id( Cl)= 0 indicating that the cell is a virgin cell. In this case we continue forward to locate the next point G’+ ’ (x0) in the processing trajectory. However, before doing that, we first set Id(ci) = - n, and Y( Ci,N) = G’ ( x0) in order to indicate that cell Ciis no longer a virgin cell but a cell under processing. (2) The newly generated point Gi(xo) =x is found to be in a cell which has a positive integer as its identification parameter. This means that this cell has been passed by the previous processing trajectories, whose global characters have been determined. Thus, we will get Id( ci) = n and Y( Ci,N) =y. The global property of point y on one of the previous processing trajectories is that it lies in the basin of attraction of the nth attractor in this analysis. In this case we will judge, by the criterion of persistence, if x approaches y close enough so that its global character as well as that of complete processing trajectory can be determined. If x and y satisfy Ilx-Yll
(6)
the current processing trajectory is terminated at this point. Thus, the global character of all the points in the presently processing trajectory is easily determined. Obviously, all points of this trajectory will be in the same basin of attraction as that of pointy. Now, the points of this trajectory can be stored on a disk together with their global characteristic symbolic number n. Before we go back to pick the next initial point to begin a new trajectory, we must deal with all the identification vectors accompanying the just finished trajectory in the following manner: reassign Id ( cj) = n for the identification vectors with Id (Cj) = - IZ,, j= 0, 1, 2, .... i, and Id (cj) = n’ for the identification vectors with Id (c,) = -n’, n’# n and keep Y(cj, N) unchanged. If criterion (6) is not met, reassign Id (ci) = - n and keep Y( ci, N) unchanged, then go forward to determine the next step image point of x0, Gi+‘(xo). (3) The newly generated point Gi(xo) =x’ is just in a cell under processing with Id( Ci)= - nnr and Y( ci, N) =y’. This indicates that a current processing trajectory exhibits periodicity on the scale of a cell. To confirm if periodic motion is contained in the present trajectory, we check by the criterion of periodicity just like that in the point mapping method, i.e. that 11x’-Y’ll
=GE .
(7)
J. Jiang, J. Xu /PhysicsLettersA
ISg(1994) 137-145
141
If criterion (7) is met, a new periodic solution with a new sequential number, say n, is found and the current trajectory is terminated. All points of the trajectory are assigned by this new global characteristic symbolic number n, which is larger than those of previous trajectories. All points on this trajectory will be written out. The identification vectors are processed just like that in case (2), when a trajectory is terminated. Otherwise, set = - n,,, unchanged to continue on the next step mapping to get G’+ ’ (x0). Y(c,,N)=~‘andkeepId(c~) (4) The cell, Cl, corresponding to a newly generated point Gi(xO) =x” is a processed cell under processing with Id(ci) = -n and Y(ci, N) =y”. This indicates that the presently processing trajectory comes at least two times into this processed cell. In this case we have to judge with criteria of both ( 7 ) and ( 6 ) . Obviously, criterion (6) can be performed directly like that in case (2) but uses x” instead of x in criterion (6). To carry criterion (7) out, one would have to search back in the current trajectory to find a point which most recently passed a cell, Ci,other than x”. After finding this point replace y’ in criterion (7) by this point and do like in case ( 3 ) . Once the global characteristic symbolic numbers have been assigned to each point of the current processing trajectory an each identification vector in this procedure has been changed to a processed one, the work on this processing trajectory is finished and we go back to find a new initial point and start a new processing trajectory. After all the initial points we set are dealt with, the global dynamical behavior of a given nonlinear dynamical system is known. Remark. ( 1) In this method the cell coordinate system only serves as a reference in the analysis and the size of the cells
or the number of cells used is, in fact, not decided by the number of initial points used. This means that a cell might contain as many initial points as possible. The precision of the analysis of the present method is predominantly controlled by the proposed distance d, in the computation. (2) Usually, a maximal length of a trajectory in practical computations is set by people. If a trajectory like the one in (5) reaches this limit before a periodic attractor or a previously processed trajectory is found to go on, a nonperiodic motion is thought to be detected and we perform just like in case (3), where a periodic attractor was found. (3 ) In the computation, a trajectory may map outside the chosen region, where no identification vector is defined on the cell in that part. In this case we apply a mapping G,,,( ) =G“( ), k is a positive integer to the point outside the chosen region. If this trajectory is led back to the chosen region by this mapping, the same analysis is performed as introduced above. Otherwise, this trajectory will be thought undeterminable and is terminated with no global characteristic symbolic number assigned.
4. Examples of applications
In order to illustrate the point mapping under the cell reference method, a series of computations has been performed with different proposed distances to apply it to the system described by the forced van der Pol equation, Z+5.10(xZ -I)~+x=ll.Ocos(3.1416t).
(8)
The region of state space to be examined is defined by -3.Odx,<2.0,
-6.O~Rg4.0.
This particular system supports a steady-state period-three motion as well as a steady-state period-five motion. Moreover, the system exhibits fmctal basins of attraction (see Ref. [ 71). We choose this system rather than a system which can support nonperiodic motion so as to make a reasonable comparison between the computation work of the point mapping method and the present method. If a system supporting chaotic motion is used, the point mapping method will extend all trajectories leading to chaos to their maximal length before they are ter-
142
J. Jiang, J. Xu /Physics Letters A I88 (1994) 137-145
minated. This maximal length of a trajectory in numerical analysis is set differently by different people. All computations were performed on a LINK 386/DX33 microcomputer using an integration step length of 0.04. We overlay a grid of 100x 100 points on the chosen region defined above. Thus, the minimal distances between two initial points in the x-direction and i-direction are d,= 0.06 and d*=O. 1, respectively. They will be denoted as grid scale below. In our method cell coordinates serving as a reference must also be built in the state space. To make the dimension of the cells be also in 100x 100 within the chosen region as the grid of initial points, a discretization of L\x=O.O6 and A.&O. 1 is used. This implies that each cell contains only one initial
(4
4.0
2.0 0.0 X -2.0
-4.0
-0.0 -< X
2.0 .
.
%
0.0 ii -2.0
-4.0
. -6.0
1.0
-2.0
-1.0
0.0 X
1.0
2.0
3:0
Fig. 1. Basins of attraction represented by the states of 100 x 100 initial points and attractors:Point mapping method; (m) period-three motion, (m) period-five motion.
J. Jiang, J. Xu /Physics L&ten A 188 (1994) 137-145
143
point. But as stated in the remarks of Section 3, the cell array is not necessary of the same dimension as that of initial conditions. Fig. 1 represents an “exact” solution in the sense that each trajectory for a particular initial point was continuously mapped until the response goes on either of the steady-state motion. No criterion of persistence is used. Figs. 2-4 illustrate the present method for different proposed distances used in the criterion of persistence. The CPU time required and the final states of all 100 x 100 initial points were recorded in all cases for comparison in details, see Table 1. Fig. la shows that points in the state space are divided into two main types: Those that approach a periodthree orbit (black) and those that lead to period-five motion (gray). The actual Poincare maps of the attractors are illustrated in Fig. lb. Fig. 2 shows the results of using the present method with a proposed distance, whose components are G&= 0.006 and d&=0.01, that is, the proposed distance in the criterion of persistence is one-tenth of the grid scale. The match between this plot with the exact solution, Fig. la, is good but with discernible differences in some parts of the boundaries of the period-three and period-five motion. We will make a comparison in more detail. Table 1 lists the details of the relative accuracy of basins of attraction over exact solution as well as the relative CPU time used in these computations. It is shown quantitatively that the present method in this proposed distance achieves an accuracy of 83.6W in the basin of attraction of period-three motion and an accuracy of 89.6W to that of period-five motion using less than $ of the original CPU time. Figs. 3 and 4 display the results of applying the present method with proposed distances of one-hundredth of the grid scale and one-thousandth of it, respectively. The matches between the two plots and the exact one in Fig. 1 are extremely good or, to be more precisely, indistinguishable. From Table 1 one can find that the result in Fig. 4 is the exact one, which is obtained using only one-fourth of the CPU time used by the point mapping method. If one thinks that the accuracy of the result in Fig. 3 is enough for one’s purpose, then a further reduction of the computation can be made as shown in Table 1. In all the above computations, the only character preserving the proposed distance is the one-thousandth of the grid scale and proposed distances less than it will be, undoubtfully, the character preserving ones. The attractors determined by the present method are the exact
4.0
2.a 0.0 X -2.0
-4.0
-6.0
Fig. 2. Basins of attraction represented by the states of 100x 100 initial points: Present method with the proposed distance being one tenth of grid scale; (m) period-three motion; ( n ) period-five motion.
144
J. Jiang, J. Xu /Physics Letters A 188 (1994) 137-145
4.0
2.0
0.0 ir -2.0
-4.0
-6.0
X
X
Fig. 3. Basins of attraction represented by the states of 100x 100 initial points: Present method with the proposed distance being one hundredth of the grid scale; (m) for the period-three motion; (I) for the period-rive motion. Fig. 4. Basins of attraction represented by the states of 100x 100 initial points: Present method with the proposed distance being one thousandth of the grid scale; (m) period-three motion; (m) period-five motion. Table 1 Comparison of relative accuracy and CPU used point mapping method basin of attraction of P-3 basin of attraction of P-5 CPU time used
100% 100% 100
Point mapping under cell reference l/10 of G.S.
l/100 ofG.S
l/1000 of G.S.
83.6% 89.6% 7.50
99.1% 99.15% 16.7
100% 100% 23.2
Note: G.S. means the grid scale defined above.
ones as in Fig. lb for all cases since the criterion of periodicity used is the same as that in the point mapping method. The usefulness of the point mapping under reference for global analysis is well demonstrated by the above examples and by comparison with the exact solutions. In fact, the present method is a flexible method for global analysis with controllable accuracy and computational work.
5. Conclusions After analyzing the inefficiency of the originally introduced point mapping method in global analysis, a criterion of persistence of the system characteristics is first proposed in this paper. The necessity to explicitly propose the criterion of persistence lies in the following two aspects: (a) To make people pay attention to the reliability
J. Jiang, J. Xu /Physics LettersA 188 (1994) 137-145
145
of the results of the point mapping method, especially in the case to determine the fractal dimension of the basins of attraction, where characters of two closely located points (the distance between them may be of the order of lo-“) are needed to be identified. (b) To make it possible to find an efficient way to get exactly the same results as those of the point mapping method, as the work done in the present paper. To make the point mapping method together with this criterion implementable, a cell coordinate system developed by Hsu is adopted to establish another reference coordinate system in the state space to assist the global analysis. It is shown by direct applications that the present method is a flexible method in global analysis, which can control the achievement of the system accuracy and the expense of the computation. Moreover, the present method is much more efficient than the point mapping method formerly used in global analysis after the comparison of the CPU time used to get exactly the same results.
Acknowledgement This work was supported by the National Natural Science Foundation of China and the Provincial Natural Science Foundation of Shaanxi, China.
References [ 1 ] C.S. [2] C.S. [3] C.S. [4] C.S. [S] B.H. [6] B.H. [7] S.W.
Hsu, ASME J. Appl. Mech. 47 ( 1980) 931. Hsu and R.S. Guttalu, ASME J. Appl. Mech. 47 ( 1980) 940. Hsu, R.S. Guttalu and W.H. Zhu, ASME J. Appl. Mech. 49 (1982) 885. Hsu,ASME J. Appl. Mech. 49 (1982) 895. Tongue and K. Gu, ASME J. Appl. Mech. 55 (1988) 461. Tongue, Physica D 28 (1987) 401. McDonald, C. Grebogi, E. Ott and J.A. York, PhysicaD 17 (1985) 125.