11—
Non-Linear Transformation Studies On Electronic Computers:
With P. R. Stein (LADC-5688, 1963)
This paper is a continuation of the study initiated and reported in the preceding chapter on Quadratic Transformations. Interactions of polynomial transformations, particularly cubic transformations in three variables, as well as the asymptotic and ergodic properties of the sequences of iterated points are considered. The theme of the computational study of difficult problems in pure mathematics is exemplified. (Eds.)
Introduction
This paper will deal with properties of certain non-linear transformations in Euclidean spaces--mostly in two or three dimensions. In the main they will be of very special and simple algebraic form. We shall be principally interested in the iteration of such transformations and in the asymptotic and ergodic properties of the sequence of iterated points. Very little seems to be known, even on the purely topological level, about the properties of specific non-linear transformations, even when these are bounded and continuous or analytic. The transformations we study in this paper are in fact bounded and continuous, but in general many-to-one, i.e., not necessarily homeomorphisms. In one dimension such transformations are simply functions with values lying in the domain of definition; for example, if f(x) is continuous and nonnegative in the interval [0.1] and max[f(x)] < 1, then x' f(x) is a
transformation* of the type considered. Even in one dimension, however, nothing resembling a complete theory of the ergodic properties of the iterated transformation exists. On the algebraic side, we study in this paper the invariant points (fixed points), finite sets (periods)-and invariant subsets (curves) of these transformations-together with the means of obtaining them constructively. The topological properties of two (not necessarily one-dimensional) transformations S(p), T(p) are identical under a homeomorphism H: when S(p) = H[T[H-l(p)]]. When S and T are themselves homeomorphisms-and for one dimen sion-necessary and sufficient conditions for conjugacy are known.1 When S and T are one-dimensional, but not necessarily one-to-one, it is possible to give a set of necessary conditions for conjugacy; no meaningful sufficient conditions, however, are known.
For example, the set of fixed points of S has to be topologically equivalent to those of T. The same must hold for the set of periodic points, i.e., points such that the nth power of the transformation returns the point to its original position. The attractive and repellent fixed points must correspond, etc. These conditions are known from the corresponding study of homeomorphisms. For many-to-one transformations one may generalize these conditions by considering the tree of a point. For a given transformation T we define the tree of a point P as the smallest set Z of points such that: a) P belongs to Z. b) If a point Q belongs to Z, then T(Q) belong to Z. c) If Q belongs to Z, then all points of the form T-1 (Q) belong to Z.
Obviously, for two transformations to be conjugate, the trees of corresponding points must be combinatorially equivalent and, in addition, their topological interrelations must be the same.**
The present study was initiated several years ago2 with the consideration of certain homogeneous, quadratic transformations which we called binary reaction systems. A typical example is the following: x1 = x2 + 2x132 , I = 2x1X3 + 2X2X3, (1) 2 X3 - X1 ,
* Here and throughout the paper a primed variable always represents the value obtained on the next iterative step. In a more explicit notation, the above equation would read: x( n+ l) = f(x(n)). ** One-dimensional transformations are considered in more detail in Appendix I.
where we consider initial points P with coordinates xl,x2,X3 satisfying: O < Xi < 1, x + x2+x3= 1 . (2) Since xI + x2 + x3=(Xl +x 2 +x3)2
the transformation (1) maps the two-dimensional region (2) into some sub-region of itself. The choice of these transformations was motivated by certain physical and biophysical considerations. For example, the set of equations (1) could be interpreted as determining the composition of a hypothetical population whose individuals are of three types, conventionally labeled 1, 2, and 3. The xi would then represent the fraction of the total population which consists of individuals of type "i." The transformation can be thought of as a mathematical transcription of the mating rule:
type 2 and type 2 produce type 1 , type 3 and type 3 produce type 1 , type I and type 2 produce type , (3) type I and type 3 produce type 2, type 2 and type 3 produce type 2, type 1 and type 1 produce type 3 .
For any assigned initial composition, i.e., any initial vector (x1 , x2, x3 ) satisfying (2), we may then ask: What is the final (or limiting) composition of the population after infinitely many "generations," that is, after infinitely many matings according to the scheme (3)? In the present context, a mating rule can be defined as a system of three non-linear first-order difference equations of the form: x1= fi(x1,X 2,X 3), x2 = f2(xl, 2, x3), (4) X3= f3 (Xl, x2, x3)
where each fi is the sum of some subset of the six homogeneous monomials x2, x2 , x 3, 2x1X2, 21 Xl3 , 2x2x3, and each such term must belong to one and only one fi. Two transformations are called equivalent if they are conjugate under the (linear) transformation defined by a given permutation of the indices 1,2,3. (This is the only linear homeomorphism which preserves the homogeneous quadratic character of the transformation.) Under this definition of equivalence, it turns out that there are 97 inequivalent transformations of the above type. It
quickly becomes apparent that, despite their formal simplicity, these transformations are very difficult to study analytically, particularly if one is interested in their iterative properties. For example, for most initial points in the region of definition, the sequence of iterates generated by repeated application of the transformation given by equation (1) converges to a set of three points: P2 = T(pl) , p3 = T(p 2) , (5) P1 = T(p3)
Using a standard terminology to be explained in detail below, we say that the "limit set" is a "period of order three." It is clear by inspection of transformation (1) that another limit set exists; if we write pi = (XI = 1, x2 = X3 = 0), P2 = (XI = X2 =0, X3 = 1), then P2 = T(pi), pi = T(p2) (6) In addition there is the algebraic fixed point of (1): p = T(p) . (7) The general initial vector, however, always leads to (5). Certain other quadratic transformations show an even more complicated behavior.
An example is the transformation: x1= x2+ x3 + 2X1X2, x2 = X12+ 2x2x3, x3 = 2xx3 (8)
This bears a close formal relationship to (1); in fact, they differ only by the exchange of a single term. The limit sets, however, are quite different. Transformation (8) has an attractive fixed point with coordinates: 1 = , x3 . (9) 31=2, 2- 4 '3= 4 (9)
It also has a limit set of the type (6) with pi = (1, 0, 0), P2 = (0, 1,0). In this case, both limit sets are observed. It is found experimentally that the set of initial points leading to (9) is separated from those leading to the oscillatory limiting behavior (6) by a closed curve surrounding the fixed point (figure 2 of chapter 10). The analytical nature of this boundary curve remains unknown.
In view of the complicated behavior exhibited by these examples, we felt it would be useful to study these transformations numerically,
making use of the powerful computational aid afforded by electronic computing machines. From one point of view our present paper may be looked on as an introduction, through our special problems, to modern techniques in experimental mathematics with the electronic computer as an essential tool. Over the past decade these machines have been extensively employed in solution of otherwise intractable problems arising in the physical sciences. In addition to solving the particular practical problem under consideration, this work has in some cases resulted in significant theoretical advances. Correspondingly, attempts to solve difficult physical problems have led to considerable improvements in the logical and technical design of computers themselves. In contrast, the use of electronic computers in pure mathematics has been relatively rare.* This may be partly due to a certain natural conservatism; in our opinion, however, the neglect of this important new research tool by many mathematicians is due simply to lack of information. In other words, the average mathematician does not yet realize what computers can do. It is our hope that the present paper will help to demonstrate the effectiveness of high-speed computational techniques
in dealing with at least one class of difficult mathematical problems. With this end in mind, we have devoted the first section of our paper to a brief discussion of how computing machines can be used to study problems in pure mathematics. Much of this section is introductory in character, and is meant primarily for those readers who have had no firsthand experience in the use of computers. It also includes, however, a description of the numerical techniques used in this study; these may be of interest even to seasoned practitioners.
After our study of quadratic transformations in three variables,** we decided to investigate the iterative properties of other classes of polynomial transformations. As a natural generalization of the quadratics described above, we consider transformations of the form: xi -fi(x,...,Xk) (i = to k), (10) where the fi are disjoint sums of the homogeneous monomials which arise on expanding the expression: / k\ m F - xi (11) i=l
* Perhaps the greatest computational effort has been expended on problems in number theory. See refs. 2 and 5. ** We shall not discuss this work here. Full details are contained in the above references. That report contains, in addition, some fragmentary results on a few particular quadratic transformations in higher dimensional spaces.
The number of such terms, each taken with its full multinomial coefficient, is Nm ().m + k- )(12) By construction, kY,fi=F ,i=l so that if we take k xi= 1, xi > 0, (13) i=l
the (additive) normalization of the xi is preserved. We are then dealing with positive transformations in a bounded portion of the Euclidean space of k - 1 dimensions, i.e., just the hyperplane defined by (13). If m = 2, k = 3, these transformations are the 97 quadratics in three variables introduced above. The bulk of the present paper is devoted to the case m = k = 3, i.e., cubic transformations in three variables; there are 9370 independent transformations of this form. We have also examined the 34337 quadratic transformations in four variables, but our analysis of the results is not yet complete (January, 1963); for this case (m = 2, k = 4) we include only some statistical observations and a few interesting examples. These three cases-m = 2, k = 3; m 2, k = 4; m = 3, k = 3--are the only ones for which an exhaustive survey is at present feasible. For other values of m and k the number of transformations to be studied is much too large.
The determination of the exact number Tkm of inequivalent transformations for arbitrary m and k is an unsolved combinatorial problem. It can, of course, be reduced to enumerating those transformations which are invariant under one or more operations of the symmetric group on the k indices, but no convenient way of doing this is known. The problem, however, is not of much practical significance. A lower limit T*m to the number Tk of inequivalent transformations is given by: TkF = Sk, (14)
where Sj is the Stirling number of the second kind. S, is also the number of ways of putting i objects in j identical boxes, no box being left empty. This underestimates Tk by assuming in effect that each transformation has k! non-identical copies, i.e., that no transformation is invariant under any permutation (except the identity). The following table illustrates the trend:
The Tk were obtained by direct enumeration-using, of course, all known shortcuts. For m = 2, k = 4, this enumeration was actually performed on a computer. In view of the huge values of the Tk'm in the lower half of this table, it is unlikely that anyone will be interested in attempting a comprehensive numerical study of these transformations for values of m and k larger than those we have considered.
A general discussion of our results for the cubics in three variables and the quadratics in four variables is given in section II; the reader will also find there formal definitions of a few basic concepts and an explanation of the special terminology employed throughout the paper. Perhaps the most interesting result of this study is our discovery of limit sets of an extremely "pathological" appearance. The existence of such limit sets was quite unexpected,*-and is indeed rather surprising in view of the essential simplicity of the generating transformations. Sections III and IV are concerned with the effect--on the iterative properties of our transformations-of two types of structural generalization. Specifically, in section III we consider the one-parameter generalization-called by us the "At-modification"-which consists in replacing equation (10) by: i = (1- At)Xi + tfi(xl,... ,Xk) (i = I to k), 0 < At < 1 . (12a)
This generalization has the special property of leaving the fixed points of the transformation invariant, although their character-i.e., whether they are attractive or repellent-may be altered. The detailed
* Quadratic transformations in three variables apparently do not exhibit similar pathologies.
discussion of the behavior of such transformations under variations of the parameter At is limited to the cubic case.
Section IV describes the result of introducing small variations in the coefficients of the monomials which make up the various fi. Again we deal only with the cubic case, and indeed only with a few interesting examples chosen from our basic set of 9370 transformations. Let us denote the Nk monomials (e.g., x3 , 3x2x2,6x1x2x3,...) in the expansion of (11) by the symbol Mj, j = 1 to Nk. The assignment of a particular index to a particular monomial is arbitrary.
Then we have k fi= E dijMj,<i<k, (13a) j=l with dij = 1 or 0, (14a) k E i= 1. (15) i=l
The generalization then consists in relaxing the restriction (14). If this were done subject only to the condition that the dij all be nonnegative, we should be dealing with a (k - 1)Nk parameter family of positive, bounded, homogeneous polynomial transformations. At present nothing significant can be said about this class as a whole. As explained in section IV, our procedure has been to study one-parameter families of transformations which are in a certain sense "close" to some particular transformation of our original set.*
In section V we give a brief, heuristic discussion of the connection between our transformations-which are really first-order non-linear difference equations--and differential equations in the plane. Our conclusion is that the connection is not, in fact, very close, and that the techniques so far developed for treating non-linear differential equations do not seem suitable for handling the problems discussed in this paper.
* Some analogous but rather unsystematic investigations were carried out on quadratics in three variables, and are contained in the report cited in ref. 2. Subsequent to the appearance of that report we made some studies (unpublished) on quadratics with randomly chosen positive coefficients satisfying (15). For quadratics (at least in three variables), the conclusion seems to be that such randomly chosen transformations are most likely to lead under iteration to simple convergence for almost all initial points.
The final section of our main text-section VI--contains a description of a class of piece-wise linear transformations on the unit square. These transformations exhibit interesting analogies with our polynomial transformations in three variables. Relatively little work has been done on this "two-dimensional broken linear" case, but the preliminary results we report seem to indicate that a detailed study might prove worthwhile.
There are two appendices: Appendix I is largely devoted to an extended discussion of certain non-linear transformations in one dimension, on the unit interval. Some of these are special cases of our cubics in three variables; others originated independently of our principal study. It is perhaps rather surprising how little can be said theoretically even about this simple one-dimensional case. It turns out that some of the same phenomena are observed in one dimension as are found in the plane-e.g., the apparently discontinuous behavior of limit sets as a function of a monotonically varying parameter. Of course, the repeated iteration of a one-dimensional transformation is a much simpler matter than the corresponding process in several dimensions. However, as we soon discovered, great care must be taken to avoid the phenomenon of "spurious convergence." This point is discussed in some detail and a few-rather alarming-examples are given.
Appendix II contains the bulk of the photographic evidence- including the "pathology" of the limit sets-on which the discussion of sections III and IV is based. These pictures, together with others scattered throughout the main body of the text, constitute in a sense the unique contribution of this paper. In retrospect, it seems unlikely that our investigation could have been successfully carried out without the visual aid afforded by the oscilloscope and the polaroid camera. Put in the simplest terms, unless one knows precisely what one is looking for, mere lists of numbers are essentially useless. Automatic plotting devices however, such as the oscilloscope, allow one to tell at a glance what is happening. Very often the picture itself will suggest some change in the course of the investigation-for example, the variation of some hitherto neglected parameter. The indicated modification can often be effected in a few seconds and the result observed on the spot.*
Visual display is of very great value when one is in effect studying sets of points in the plane; when one passes to three dimensions automatic plotting ceases to be merely a convenience and becomes essential. * This interaction of man and computing machine has sometimes been referred to as "synergesis."3
A glance at our pictures of three-dimensional limit sets--the result of iterating certain quadratic transformations in four variables-should convince even the most skeptical reader. In our opinion, it would be virtually impossible to make sense out of a mere numerical listing of coordinates of the points plotted in these photographs.
Of the many who have helped with this work, there are three to whom we are particularly indebted: Cerda Evans, Verna Gardiner, and Dorothy Williamson. These ladies did the actual coding and supervised all the machine calculations. Without their help this paper could not have been written.
I—
The Role of the Computing Machine
1. The use of electronic computers for the solution of complicated or tedious problems, usually of practical origin, is by now familiar. Typical computer tasks are: the evaluation of integrals, the solution of large systems of linear equations, the solution of minimax problems (linear programming), the treatment of complicated boundary value or initial value problems, etc. One of the more impressive jobs that computers have done is to calculate the time history of immensely complicated physical systems (e.g., involvinig hydrodynamical motions, magnetic fields, etc.). Recently there has been considerable interest in using computers to attack problems of a less applied nature, for example those arising in combinatorial analysis4 and number theory.5 This work often takes on an experimental flavor; such experimentation has led to results of considerable interest, for example, the construction of certain types of mutually orthogonal latin squares.6 Computers can also be used to investigate formal mathematical systems,7 to reduce symbolic expressions,8 and with less success-to study games of "skill" like chess.u
The use of computing machines that we describe in the present paper differs in two respects from the examples just cited. On the one hand, our study is not essentially combinatorial in character, but falls rather in the domain of algebra and real variable function theory. On the other hand, we are not attempting to "solve" some welldefined problem; instead we investigate via repeated trials the asymptotic properties of certain non-linear transformations, usually without any advance knowledge of what we may find in a given case. Even "after the fact," so to speak, it is difficult to classify these asymptotic properties in a meaningful fashion; the broadness of the categories we employ for this classification* is merely a measure of our lack of insight into the structure of the observed limit sets.
* See section II.
Faced with this situation, one may ask the question: how does one recognize "convergence" --- -i.e., the existence of an invariant set -when one has no a priori numerical criteria to apply? We can only supply a partial answer to this question, but that answer has the advantage of simplicity, viz.: "use your eyes." The practical application of this "technique" involves, of course,* the use of automatic plotting devices.
2. Roughly speaking, computing machines are devices which perform the four elementary arithmetic operations on numbers in a certain--not necessarily simple -sequence. This sequence of operations is called the "program," and consists of a set of logical commands, both of the sequential ("do this and then do that") and of the branching ("if this holds, then do that") type. The program is composed by an investigator (the "programmer") and must therefore reflect his own limitations. Nevertheless, the machine may easily produce results quite unanticipated by the programmer, even if the program is essentially deterministic in nature.** A classic example-which happens to be relevant here-is the step-by-step application of some recurrence relation which generates a sequence whose trend the programmer cannot determine in advance. As an example, we may cite the following one-step recursion in a single variable: Yn+l = Wn(3 - 3Wn + (wn), Wn 3yn(1 - Yn) (1)
Given some initial 0 < yo < 1, we may ask: what is the result of applying the rule (1) N times, where N is some larger number, say 105? This particular transformation is discussed in detail in appendix I; here we quote three examples for the purpose of illustration.
(a) If a( = 0.99004, then for almost all yo the sequence of iterates produced by (1) converges (in < 105 steps) to a period of order 14. (b) If a = 0.99005, the corresponding limit set is a period of order 28.***
* Hand plotting is in general highly impractical, and clearly relinquishes the principal advantage of machine computation: SPEED. ** Strictly speaking, all programs used on digital computers are deterministic in nature: even when randomn numbers are employed, these are generated according to some fixed algorithm so that the sequence is in principle known. *** These results were found by using IBM "STRETCH" computer. The periods are exact to within the accuracy of that machine, i.e., 48 binary digits (-15 decimals). See further in appendix I.
(c) If a = 0.99008, no finite period is observed after N = 5 x 105 steps.
So far as we are aware, this behavior could not be predicted by current analytical or algebraic techniques. Such phenomena are easy to study on a computer, however, because of the great speed with which it can carry out the (relatively simple) operations implied by an expression such as (1) above. In fact, 200000 iterations of this transformation takes slightly less than one minute on a really fast computer.*
3. As we mentioned in the introduction, the principal content of this paper is the study of the asymptotic properties of certain nonlinear transformations of relatively simple form. This means that, if T is such a transformation, we examine the sequence T(p), T2 (p), T3 (p), ...
for various initial points p lying in the domain of T. The mathematical object of interest to us is the set (or sets) of points to which these sequences converge. In the absence of any general analytical technique for calculating these "limit sets," we must have recourse to "bruteforce" methods.
Some non-linear transformations which appear morphologically similar to those considered here can in fact be completely analyzed by elementary methods. We discovered one such case in the course of some earlier work on biological systems. It is described in our report on quadratics in three variables (see chapter 10). We restate these special results here: Let N C' = 7" CkC , 1 l<i<N, (2) k,m=l with coefficients satisfying ykm= im > 0, min (k, m)<i< max (k, m), kcmmke (3) ,k= 7y > 0, otherwise, k E y= (4) i=m
* This figure applies to STRETCH and includes all additional "diagnostic" operations such as checking for "convergence," etc.
kE m (5) (5) t=m2 We normalize the Ci by N 0 < C, < 1, all i, Ci= 1. (6) i=l
This property is clearly preserved under iteration. With the coefficients defined as above, there exists a linear invariant: N-1N-1-y(N - i)Ci = (N -i)C:. (7) i=ii=l
Given an initial vector (C(), C() , ..., C() whose coordinates satisfy (6), a is explicitly determined. It can then be shown that every initial vector satisfying (6) converges to a definite fixed point which is determined as follows. For the given value of a, there is one value of the index j such that N-j > > N-j-1. (8)
The fixed point is then explicitly given by Cj = a -(N -j-1), Cj+1 = N -j - a, all other Ci = 0. (9)
Note that the fixed point is independent of the values of the coefficients ykm.
As simple examples of coefficients satisfying (3), (4), and (5), we may mention km 1 I k-m ,Y. (10) 2k-ml i- min (k, m) and { ik-ml+l if min (k,m)< i < max (k,m), (11) 0, otherwise.
For a fuller discussion of this transformation and its possible applications, we refer the reader to the original report.
The term "brute-force" refers to the fact that, in order to determine the convergence properties of some transformation T belonging to our class, we must in general actually evaluate Tk(p) for k = 1, 2,..., N, where N is likely to be quite large, sufficiently large that
is, so that we can observe convergence* to the limit set. To make matters clear, let us consider a specific example. We choose the cubic transformation I1 3 1= X3 + 3xx3 + 3x3 x2 + 6x1 X2 X3 x2 = x + 3X2X3 + 3x3x2, (12) 3 = x2+ 3xx12 + 3X2X . We take some initial point p = (X1, x2 , x3 ) whose coordinates satisfy: x1+x2+±x3= 1, 0<xi < 1, i= 1,2,3. (13)
The program then instructs the computing machine to evaluate the right hand side of (12), thus producing a new point p' = (xI,x, x3); the coordinates of p' again of course satisfy (13). p' is then set to p, and the process is repeated. The iteration proceeds in this fashion until either some finite limit set is found** or an invariant set-presumably infinite-is "observed." The observation consists in looking at successive groups of consecutive iterates-in practice we have usually taken 900 points at a time-until no qualitative visual change is noted over a sample of several successive such groups of points. Since the transformation (12) is really two-dimensional, we may plot the successive points p in the plane. Accordingly, we define new coordinates S, a by the linear transformation*** 1+ x1 - 3 2 (14) S=a- =-. (14) 2 2
The domain of the transformation is then the 450 isosceles triangle: 1 2 0 S I
* "Convergence" must be of course understood in some approximate numerical sense. Our usual criteria are set forth in the next subsection. ** See the next subsection. *** These are the coordinates employed in our earlier work on quadratics in three variables; we have retained them more for historical reasons than for any particular advantage they may possess.
In terms of these new variables, (12) takes the form: Si -Is3 -15S23 2 3 3S' =-1 S - S2a- Sa +-a+ 6Sa -3a + l- F(S, a), 2 2 2 2 3 139 7(15)a'= 3 +Sa + Sa2- a3 - 6Saa + 3a G(S, a) 2 2 2 2
The computer is instructed to store 900 successive points p(S(n),a(n)), p(S(n+1), a(n+l)), ..., and, when the last point has been calculated, to plot all 900 points on our oscilloscope screen.* If we choose, we may then photograph the resulting pattern with a polaroid camera. Such a photograph is shown in figure 1. Here one sees 900 successive high-order iterates (n =2700 to 3600) of the initial point, S = 1/2, a = 0.17. For convenience, the triangle of reference is also shown. Fig. 1
This calculation --as well as all others which produced the photographs in this paper--was performed on the Los Alamos Laboratory's MANIAC II computer.** MANIAC II requires about 15 seconds to calculate 900 iterates of a point by repeated applications of
* The points are actually plotted in the order in which they are calculated, the whole pattern being replotted as many times as we wish. Actually, the plotting of 900 points is effectively instantaneous so far as the human eye is concerned. If we wish to see the points plotted in succession, we must introduce artificial time delays between the plotting of successive points. ** For the use of other computing machines in this work, see the next subsection.
a cubic transformation like (12) above. This figure includes the time spent in examining the successive points for simple convergence, as well as other "diagnostic" operations.* The actual numerical values of the coordinates may be printed out whenever desired by simply flipping a switch. On MANIAC II a decimal number is normally limited to eight significant figures. In the present paper, when there is occasion to quote numerical values obtained from MANIAC II print-outs, we shall generally reproduce them to seven figures without further specifying their accuracy.
Computer programs are, of course, not limited to generating sequences of numbers from an iterative formula such as (12). A considerable amount of sophistication can be incorporated into such a program so as to allow the machine to make "decisions" in the course of the calculation. It can, in fact, examine any property or any functional of the data that the programmer can describe in appropriate terms. One problem that is met with frequently in this work is to determine the points in a sequence of iterates that lie closest to some point, say within some chosen angle or set of angles. This sort of experiment is frequently of help in elucidating the local structure of a complicated limit set. Then again we may want to determine the average values of S and a, i.e., ergodic means, taken over the sequence. To achieve any sort of accuracy in such problems** we may be required to go to 50000 or even 100000 iterations. One saving feature is that several such diagnostic experiments can be carried out simultaneously. There are, however, special questions that must be dealt with by special programs. One such question arises in connection with our illustrative transformation (12). The complicated limit set shown in figure 1 is not the only one observed. This transformation has an attractive fixed point at:
S = F(S, a) = 0.6259977, a = G(S, a) = 0.1107896 ; (16) indeed, the eigenvalues of the jacobian matrix*** evaluated at this point are complex, with \A12 = 0.4366967. Consequently, there must be a neighborhood of this point in which all sequences will converge to it. The only way to find the boundary of this neighborhood is by
* For reasons of accuracy, the calculations are performed in the xi coordinates; the transformation (14) to the S, a coordinates is carried out only for plotting purposes. ** More properly, to have confidence in the results. The accuracy cannot always be satisfactorily estimated. *** The criterion for the nature of a fixed point is discussed in Section III.
trial and error. This is a time-consuming job, even for an electronic computer; if one picks a point close to the boundary of the region of convergence, several hundred-or even several thousand-iterations may be required before one can tell whether the chosen point lies inside or outside the region. Figure 2 shows the approximate boundary for the present case, drawn through 107 experimentally determined points. One of these is known to one part in 107, while the others have been determined only to I part in 104.* Fig. 2
4. General Procedure
a. Cubic transformation in three variables. Enough has been said above to make clear the necessity of using an electronic computer in such investigations. We must now say something about the systematic aspects of the study. All 9370 cubic transformations were initially studied on an IBM 7090.** First a complete list of inequivalent cubics was prepared-this was also done on the 7090-incidentally serving to check our original pencil and paper enumeration. Then by a completely automatic procedure, each transformation was taken in turn and four randomly-generated initial points were each taken as the start of an iterative sequence. For each point the iteration was continued until either convergence to a finite set of points was "observed" or 10000 iterations had been performed. By "observed" we mean that the machine sensed convergence to a fixed point or to a finite period of order
* The point S = 1/2, a = 0.2952833 lies in the region of convergence, while the point S = 1/2, a = 0.2952834 gives rise to a sequence which converges to the class IV limit set (see definition in section II). ** This computer is approximately five times as fast as MANIAC II.
<300. More precisely, the computer was programmed to test whether the following conditions were satisfied X(n')-_) < 107, i = ,2, 3 . (17) \ -xi\<10- i=1,2,3. (17)
If (17) is satisfied, a finite limit set has been reached to within the indicated accuracy. For n = n1 + I this means convergence to a fixed point ("simple convergence"). Otherwise, the limit set is a period of order n - ni. In practice, values of the xi were stored at fixed time steps nl = 300, 600,..., the test (17) being performed on each step. If "convergence" was found, the appropriate values of the xi were printed out and the next random initial point was used, etc. If no such convergence was found after 10000 steps, the values of the iterates for the last few steps were printed, and the computer proceeded as before.
When all the cubic transformations had been studied in this fashion, the "interesting" cases--i.e., those in which no convergence was observed-were examined one by one on MANIAC II, where the visual oscilloscope display could be consulted. Many cases of apparent non-convergence turned out in fact to be convergent with the iteration carried further. It should be stressed that the restriction to 10000 iterations, which we imposed in the course of the systematic, fully automatic survey of all cubic transformations, was merely one of convenience; without some such reasonable limitation, the automatic survey would have taken too long. The same remark applies to the decision that only four randomly generated initial points be taken for each case. Past experience has shown that this last restriction is not unreasonable when a complete survey of transformations is contemplated. By this we mean that the behavior of an arbitrary transformation of our class is "likely" to be defined even if iterates of only four random points are studied. To be sure, in some cases the limit set depends in a very complicated way on the initial point; for such a transformation this crude sampling technique is not adequate. In these cases, however, the four random trials are likely to produce two difference limit sets; this in itself is an indication that the transformation in question should be studied in more detail.
For the detailed examination of a given transformation, many relatively sophisticated MANIAC programs are available. We may, in effect, study any properties of the transformation that seem of interest. Typically, these may include:
1. Determination of non-attractive fixed points (see section III). 2. Checking for periodicity. 3. Exhibiting some qualitative properties of the mapping, e.g., by showing the images under the transformation of a family of lines.
4. Determining the dimensions of the limit set. 5. Verifying that low-order periods are attractive (see section III). 6. Examining the dependence of the limit set on the initial point. We cannot expatiate here on the actual procedures involved; sufficient to say that the use of visual display (i.e., the oscilloscope plot) is an essential tool in all this analysis.
b. Quadratic in four variables. All (34337) inequivalent transformations of this class were studied by the same fully automatic method as that used to study the cubics. For this purpose a faster machine than the IBM 7090 was clearly required; we were fortunate enough to have access to the IBM 7030 STRETCH computer, which is approximately 4 times as fast as the 7090 and 20 times as fast as MANIAC II. Only partial results are reported in section II, since our analysis of the STRETCH print-outs is not complete.
The detailed study of a given quadratic in four variables is more difficult than the corresponding analysis for the three variable cubics: the domain is three-dimensional, being in fact the tetrahedron defined by 4i = 1, 0 < i<1, 1 <i <4. i=l
Thus a meaningful visual display involves plotting some properly chosen projection of the three-dimensional limit set. In some cases it may require several trials before an appropriately "revealing" viewing angle is found; consequently it was not feasible to plot every potentially "interesting" limit set in this fashion, and some sort of selective procedure had to be resorted to. The method we chose was to look at three plane projections first--e.g., xl versus x2, xl versus x3 , and x2 versus X3. It turns out that one soon develops a feeling for the "interesting" case even without being able to build up an image of the actual three-dimensional configuration from the plane "slices." More serious than this purely technical difficulty is that resulting from the generally more complicated dependence of the limit set on the initial point: it turns out that in these transformations one is much more likely to miss something by restricting one's self to a few randomly generated initial points. At the present time, lacking any local or structural criteria for the prediction of asymptotic behavior, we see no way to overcome this difficulty.
II—
Limit Sets
1. Abbreviated Notation for Transformations. In order to have a convenient way of referring to a particular transformation without having to reproduce its explicit form, we introduce at this point a simple shorthand notation. As already noted in the introduction, our cubic transformations in three variables may be written in the form: 10 x' = dijMj, i=1,2,3, (1) j=1 with dij =Oor 1, alli,j, (2) and 3 dij = 1, all j, (3) i=l
where the Mj are the separate terms in the expansion of (x1+x2+x3)3 . We now choose the following conventional ordering of the Mj. Ml = x31, M2= 2 3= , M4== 3xilx2 , M5= 3xlx2 , M6= 3x2x1, M 7= 3x 2x2 , M8= 3x3x2 , (4) Mg = 3x3x2,M1o= 6x 1x 2x3 .
Any cubic transformation of our class is then completely determined by specifying which terms Mj or, equivalently, which indices j, appear in the first two lines of the schema (1). Let us call the set of indices belonging to the first line C1 and those belonging to the second line C2; C3 is of course the complement of C1 + C2 with respect to the full set {1,2,..., 10} and need not be written down. Thus, for example, the transformation:
xi = x3+ 3XX22 + 3Xlx2 + 3 x2x 2 + 3x 3x22 + 6 1X2X3 , 2 = 1+ 3 3X1, (5) 3 + 3x2x1 x3 = x2 + 32, would appear in the form: C= {3, 4, 5, 7, 9, 10} (6) C 2 =1, 8}
An analogous notation may be adopted mutatis mutandis for quadratics in four variables. Any such transformation can be written in the form: 10 x' =- dijFj, i-= 1,2,3,4, j=1 dij = 0 or 1, all i,j, (7) 4 dij = 1, all j. i=1
Our conventional assignment of indices to the Fj is as follows: F1= x,F2 = F3= x3, F 4= x, F 5= 2x1x2, F 6= 2x 1x3, F 7= 2x 1x 4, F 8=2 x2x3, F9= 2x2x4, (8) F0o= 2x3x4. Let Qk denote the set of indices belonging to the kth line of the schema (7). Then any such transformation is specified by writing down three of the four Qk thus: Q1= {2,8,9}, TQ1Q2Q: Q2 = {3, 7, 10}, (9) Q3= {5, 6}
represents the transformation xl = X2 + 2x2x3+ 2x2, X4,x2+ 2xlx 4 + 2x3, 4, ,( Z3 = 2X1x2 + 2xlX3 , 2 2 X4=xXX 4 This notation will be used extensively throughout the paper.
2. Limit Set Terminology. By a limit set Lp(T) we shall mean the set of all points of the region of definition* which are limit points, * For cubics this is the S, a triangle introduced in section I; for quadratics it is the tetrahedron 4 i = 1, Xi> 0 i=1
in the ordinary sense, of the set Tn(p), n = 1, 2,.. , for fixed p. It may happen that Lp(T) is independent of the initial point p; Lp(T) L(T) could then be called the limit set of T. In general, L(T) will only be defined for interior points p, since points on the boundary frequently* behave in a rather special way.
Thus, for example, if po is a unique fixed point of the transformation: T(po) = po, and if the iterated images Tn(p) of all interior points p converge to po, then L(T) {po}. If po,Pl,...,Pk-1 form a system S of k points such that T(pi) = pi+l, i = 1, 2,... (mod k), and if for all interior points p, limn-,, Tnk(p) is one of these points, then Lp(T) S.
It might happen that the interior points divide into a finite number of classes C1, C2, ..., Cr such that all points p belonging to the same class Lp(T) forms the same set; we should then have a finite number of limit sets L1, L2, ..., Lr. Some of these may contain a finite number of points, others may be infinite. For convenience we shall usually refer to a finite limit set containing k distinct points as a period of order k.
Although a given finite limit set belonging to some transformation T may legitimately be considered a "property" of that transformation, it is in no sense characteristic; many different transformations of our type may possess the same limit set, even for the same set of initial points. It should also be stressed that not every set of points S {pi, ..., Pk} such that T(pi) = Pi+l (i mod k) is properly a limit set. Such a set of points, each of which is a solution of the equation Tk(p) = p, must have the additional property that there exists a set of initial points whose iterated images converge to S. Finite sets S which have this property are conventionally termed attractive. Thus, we should properly refer to a finite limit set of k points as an attractive period of order k. In the sequel we shall usually omit the word attractive when the context makes it clear that this is what is meant.
There is, of course, no structure problem so far as finite limit sets are concerned; they are completely described by giving the coordinates of their constituent points. For infinite limit sets, the situation is different. On the basis of our numerical work alone, we cannot say with certainty that our transformations have such limit sets; the sets may in fact be finite (with an enormous number of points in them), but the presumption that they are infinite is very strong. For any observed infinite limit set we can at most say that it is not a period of order less than some very large k. Granting, however, that we are dealing with infinite sets, and that we may infer some of their properties by
* i.e., often enough to make it worthwhile excluding them in the definition of L(T).
examining a sufficiently large finite subset* we may attempt to classify them according to their macroscopic morphological properties.
3. Infinite Limit Sets for Cubics in Three Variables. On the basis of our empirical study of cubic transformations, we may make a rough division of infinite limit sets into four classes:
Class I. This includes all limit sets that appear to have the form of one or more closed curves. Figures 3 through 6 will serve as examples of this class. The detailed structure of these "curves" has been studied numerically in some cases, but there are as yet no theoretical arguments to the effect that these are really one-dimensional continua. Fig. 3 Fig. 4 Fig. 5 Fig. 6
* This assumption underlies all our numerical work.
To illustrate one type of numerical study that we have carried out on these limit sets we cite the case of figure 3. This shows the "infinite" limit set L(T) belonging to the transformations:
C = {2, 5, 7,9, 10 , 2 C2 = {1,3,6,8} )
In the S, a coordinates, this takes the form: S/ 3 3 15 23 3 3 1 S' 3 b- S2a - -Sa2 + a3 - 3S2 - 3a2 + qS + -a + , 2 2 2 2 2 2 2' (12) 3 3 9 1 3 3 1 - + 2a- Sa2+ a3 +3S2+3a2 - S-a+ 2 2 2 2 2 2 2
There is a (repellent) fixed point at: So = 0.6149341, ao = 0.1943821 . (13)
To six decimal places, the overall bounds on the curve are* Smax = 0.816878 at a = 0.058022, Smin = 0.411270 at a = 0.204391 (14) amax = 0.435861 at S = 0.552246, amin = 0.017750 at S = 0.728386 .
To five decimal places, the average value of coordinates is found to be* 1 1 N NS-N S S)- 0.62231, a -Na() = 0.20772. (15) i=l i=l
This set L(T) is the only infinite limit set the transformation seems to possess [the pair (S = 1, a = 0), (S = 1/2, a = 1/2) turns out to be an attractive period for this transformation.] For "most" initial points, the sequence of iterates converges to L(T). If we choose as our initial point some p C L(T), the curve will be traced out by successive images of p, though not in a continuous fashion. If, however, we look only at successive iterates of the 71 st power of T: T71 , the curve is indeed generated in a relatively continuous fashion; the successive points T(71n)(p), n = 1, 2,..., lie close to each other and trace the curve in a clockwise sense. This is illustrated in figure 7, where 246 successive values of T(71)(p) are plotted. It is striking that the nonuniform density of
* These results were obtained by carrying out N = 9600 iterations, starting from a point on the curve with coordinates: S = 0.5841326, a = 0.4125823.
points along the curve-as shown in figure 3-is reproduced by this sequence of iterates. We are thus led to the conjecture that L(T) and L(T(71 )) coincide.
It is, of course, by no means generally true that T(k) and T will have the same limit set for an arbitrary T of our type (cf. the case of periodic limit sets where k is a multiple of the period). Further experiments have convinced us that L(T) _ L(T(k)) for all k in this case. If this is so, the set is certainly infinite. That it is a continuum is also very probable.
The presumption that L(T) is one-dimensional is supported by the following experiment. We choose a point po which seems to lie, with all available precision, on a convex portion* of the curve, and obtain 100000 iterated images of it, keeping track of those iterates which lie closest to p0. We find that the two points Pl and P2 of closest approach lie in opposite quadrants with respect to Po, and that the slopes of the two line segments (pi, po) and (p2,Po) are the same to within a fraction of a percent. This suggests (1) that the limit set is a curve, and (2) that the curve probably has a continuous derivative at Po **
Limit sets consisting of several separate curves (figures 8, 9, 10) may in principle be treated in the same manner, although it is then no longer true that T(k) will have the same limit set as T for all k. For example, if L(T) consists of three separate curves, L(T( 3 )) will coincide with only one of these-which curve depends, of course, on the initial point.
Class II. This class consists of those infinite limit sets all points of which lie on a pair of boundaries of the (S, a) triangle. Alternate iterates lie on alternate sides, hence the square of the transformation will have a limit set confined to one side of the triangle. T(2)(p) is then strictly one-dimensional for all p situated on one or the other of the two sides in question. The situation is illustrated in figures 11 and 12. There seem to be only a few such one-dimensional limit sets possible within our class of cubic transformations. Correspondingly, many different cubic transformations lead to the same pair of onedimensional transformations when the set of initial points is restricted to a pair of sides of the (S, a) triangle.
* Overall convexity is rarely, if ever, a property of these limit sets. ** We do not conjecture that the derivative exists at every point, but we think it likely that the number of points where the derivative does not exist is at most a set of measure zero.
blank
For example, every transformation of the form:
with non-negative ai, bi, ci satisfying:
will lead to the pair of one-variable polynomial transformations of 6th order:
In other words, transformations T of the form (16) have the property that T(2) transforms each of the lines X3= 0 and xi = 0 into subsets of themselves (in the S, a coordinates, these lines are respectively the boundaries S+ a = I and S = a). The study of such one-dimensional transformations is much easier than that of the original plane transformations, but there are certain serious computational pitfalls connected with high-order iteration (see appendix I).
Class III. The limit sets constituting this class will be referred to as pseudo-periods. They consist of relatively dense clusters of points localized at a finite set of centers, with a few scattered points in between (figure 13). Such limit sets have not been observed for our original cubic transformations with integer coefficients; they are, however, a prominent feature of the more general transformations discussed in sections III and IV.
ClassIV. In this class we place all infinite limit sets not included in the first three classes. Viewed on the oscilloscope they appear as very complicated distributions of points with no recognizable orderly structure. Some examples are shown in figures 14 through 17. A few other examples will be discussed in detail in the following sections. For illustrative purposes, however, we include here a few remarks about figure 17.
This limit set belongs to the transformation:
blank
As is evident from the photograph, it consists of seven separated pieces; each of these is invariant under the 7th power of the transformation. Extensive experimentation indicates that the gaps are really there. There appears to be no orderly structure within the separate pieces; in figure 18 we show about 385 consecutive images T(7 )(p) (in the upper left-hand piece of the limit set) of some p lying in this subset.
Fig. 18
4. Statistical Observations.
a. A large majority of our 9370 cubic transformations in three variables-some 75 per cent--exhibited what might be called simple convergence for all initial points tried. For these the limit sets consist of a single point, i.e., a fixed point of the transformation. In many cases there are two such attractive fixed points, but we have not found a case in which both such points are interior to the (S, a) triangle.
We exclude here a few trivial cases such as the following. Consider T C1 = {1,2,3,4,5,10} , Tc,' -{C 2 = 6,9}. (
Explicitly, the second and third line read: x' = 3x2(x2 + x2 x3 ) (22 x3= 3X3 (X1 + x 2x 3) Thus 22 (23) x3 X3
so that this ratio is fixed by the initial value, and we have a continuum of fixed points. Setting x 2/x 3 _ r, we find that the fixed point is given by 1 +r- (1+r2 )/3 (24) x3 = (24) 1 + 3r + r3 with, of course, x2 = Tx3,1 = 1 - (1 + r)x3. (25)
If we consider the transformation derived from the above by interchanging the right-hand sides of (22), we shall have: 2X XI 22-1 3- (26) yielding a corresponding continuum of limit sets which are periods of order two.
b. About 16.5 percent of the transformations seem to have only finite (periodic) limits sets; not surprisingly, most of these are of order two. More than half of the latter are of a trivial nature, that is, two vertices of the triangle permute under T. Less than 20 cases were found for which the limit set was a period of order k > 3. High-order periods are, however, frequently encountered in the study of the generalized transformations discussed in sections III and IV.
c. Some 5 percent of the cases were found to have several (i.e., two, rarely three) distinct finite limit sets of the types described above. For a given transformation it would in principle be possible to determine numerically the set of initial points whose iterated images converge to a particular one of the several limit sets; lack of time has prevented us from doing this except in a few cases. We only remark that there is in general no reason to suppose that the boundary of such a set of initial points is simple.
d. The remaining 3.5 percent, some 334 transformations, possess infinite limit sets. Most of these (roughly three-quarters of them) belong to class I, that is, they look like closed curves. Perhaps 5 percent of the rest belong to class II, the 20 percent residuum being of class IV type. As mentioned above, no examples of class III (pseudo-periods) limit sets were encountered in the study of our original group of cubic transformations (i.e., those with integer coefficients 1, 3, or 6).
e. No case has been found in which a transformation has two distinct class IV limit sets, although there are cases where one of several
limit sets was of class IV type. One such has already been described in section I (page 16); a more complicated example will be mentioned in section IV below.
f. We can say very little about the rate of convergence of a sequence T(n)(p) to its Lp(T). Sometimes it may be extremely rapid (10 to 20 iterations); in other cases many thousands of iterates may be required. If Lp(T) consists of a single point, Lp(T) = {po}, this rate can, of course, be calculated (for points sufficiently close to po) by solving the approximate, linear difference equations explicitly.
This is, however, not always sufficient. If the jacobian matrix, evaluated at po, has complex roots, and \A2[= 1, the linear difference equations may generate an inlvariant ellipse. Such a case was found in one of our quadratics in three variables, and is discussed in our report on that work. In the S, a coordinates, this transformation is S' = 1 -4a + 4a 2 +2aS, a' =2aS (27) with fixed point at 1 1 S= a = - (28) 2' 4
Letting 1 1 x= S--y=a-, (29) 2' 4 the linear approximation is 1 1 x' =-y + x, yy y+x . (30) This then generates the invariant ellipse: x'2+xy' + 2y2 =x 2 +xy+ 2y2 . (31) In fact, however, for the full (nonlinear) transformation, the fixed point is attractive.
5. Limit Sets for Quadratic Transformations in Four Variables. All 34337 distinct systems of this type have been investigated on the STRETCH computer, as described in section I above.* A preliminary survey of the results indicates that only about 2 percent of
* The computing time required for the whole study was only a fraction of what one would predict on the basis of 7 seconds for 105 iterations-an average for these recurrence relations as actually coded---because a large majority of cases "converged" in a few (~50) steps.
these transformations possess infinite limit sets. The finite limit sets need no special comment; they are of the same sort as those found in cubics in three variables-except, of course, that they are not in general plane sets. A few periods of rather high-order (more than 100 points!) were found, as well as a fair number of cases with 10 to 80 points. This probably should be expected in view of the greater variety of possible algebraic structures.
We are not yet in a position to classify the infinite limit sets as we have in the case of the cubics. Perhaps the closest analogy to sets of the class I type are those which appear to be closed curves in space. These are illustrated in figures 19 through 23. They are shown in convenient projections; the "coordinate system" in the center of the picture merely indicates the orientation relative to the viewer, who is conceived of as stationed at a certain distance from the origin along the y axis.* Figure 19 shows a limit set belonging to the transformation: Q1 = {1,3,4}, TQ12Q3: Q2 = {5,6,8}, (32) Q3 = {7, 10}
Presumably what one sees is a twisted space curve. In figure 20, the limit set consists of two plane curves, one of which lies in the (x1 , x3) plane, the other lying in a plane inclined at 450 with respect to the first. The corresponding transformation is Q1 {1,2,9}, TQ1Q2Q 3: Q2 = {4, 7, 10} , (33) Q3 {3,5,6}.
The observed limit set is at least consistent with the fact that (33) evidently transforms these planes into each other (so that the points lie alternately on the separate curves). More complicated twisted curves are possible (figure 21). We have also found quite implausible looking limit sets like that shown in figure 22. As a final example, we cite the transformation:
* The "reference system" (x, y, z) is parallel to, but displaced relative to, the actual coordinate system (x1,x2,x3). The origin of the (x, y, z) system is in the (approximate) center of the picture; that of the (x1,x2,x3) system is in the lower left-hand corner.
blank
This has at least two infinite limit sets, one of which may be of class IV type (not shown); the other (figure 23) is a "curve" of unknown structure.
At the time of this writing (January, 1963) we are unable to say anything more specific about the limit sets for quadratics in four variables; to date, less than one-third of the seven hundred or so potentially "interesting" cases have been looked at on the oscilloscope.
III—
The "At-Modification"
1. We discuss here a particular one-parameter generalization of our original cubic transformations in three variables which we have called the At-modification.* It consists in replacing the usual difference equations:** 10 x'=dijMj, i = 12,3, (1) j=l by 10 x' = (1- At)xi + At dijMj , (2) j=l with O < At< 1. (3)
If At = 1, we recover the original set (1); At = 0 is excluded, since the equations then become the trivial identity transformation.
The abbreviated notation of section II is extended in an obvious manner to cover this case. Thus, if (1) is represented symbolically by TC1C2, (2) may be symbolized by Tc 1c 2(At). For a given p, the limit set will correspondingly be denoted by Lp(T), Lp(At)(T).
* This has already been mentioned in the Introduction, equation (12). The modification can, of course, be introduced for the general case [equations (10) through (13) of the Introduction.]
** The Mj are defined in section II, subsection 1, equation (4). Unless otherwise stated, equations (2) and (3) of section II are assumed to hold, as well as the condition Zi = 1, xi> 0, for all initial points p = (xl,x2 ,x3 ) .
In the S,a coordinates, (1) appears in the form S' = F(S, a), a' = G(S, a) . (4) Correspondingly, (2) reads S' = (1 - At)S + AtF(S, a), a' = ( - At)a+ AtG(S, a) . (5)
It is clear that the fixed points of (5) coincide with those of (4). As we shall see below, this fact enables us to find these fixed points by simple iteration, thus avoiding the unpleasant algebra involved in eliminating one variable from the pair of general cubics S = F(S,a), a = G(S,a) , (6)
and then solving for the roots of the resulting high-order (< 9) polynomial. One can look upon (5) at the simplest (and most naive) finite difference scheme for approximating the first-order differential system dS da d = -S+ F(S, a), d=-a +G(S, a). (7) dt dt
The analogy between (5) and (7) is not, however, very close;* consequently it is better to discuss (5) on its own merits. The effect of setting At < 1 (but > 0) on a single iteration is easy to see. Let us take a particular point S, a; the image produced by (4) will be denoted, as usual, by S', a', while we shall call the corresponding image under (5) S' mod, a' mod. Then S' rod- S' = At(S'- S), a' mod- a' = At(a'- a) . (8)
In other words, the length of the iterative step is altered, while the direction remains the same. What happens on repeated iteration is, however, not all obvious. One expects that the limit set Lp(,t)(T) will in general have smaller diameter as we decrease At, but we cannot at present predict its structure as a function of At, even relative to the (observed) structure of Lp(T). It is worthwhile illustrating this in a particular case. Consider the transformation.
* For further discussion on this point see section V.
In the S, a coordinates, this reads explicitly: 3 S2+ 3a2 S' = S3 - 6 S2a - 3Sa2 + 4a3- S2_ +3Sa-a2 +1. TA': 2 2 (A) a' =-S3 + 3Sa2 + 2a3 + S-3Sa + 2a 2 2
Since we shall refer to this transformation quite often in the sequel, we have given it the distinctive label (A). TA has one interior* fixed point (repellent), whose coordinates are: So = 0.5885696, ao = 0.1388662 . (10)
There are two infinite limit sets; these are shown in figure 24 and in figure 11 of section II. At the moment, we shall not be concerned with the limit set shown in figure 11; this is evidently of class II type and can therefore be studied in one-dimensional form.** The limit set shown in figure 24--which we shall henceforth refer to as L(TA) appears as an irregular pattern surrounding the fixed point (shown superimposed on the picture). Figure 25 again shows L(TA)-this time enlarged by a factor 3, while figure 26 shows a portion of the upper left-hand corner*** enlarged about 14 times. Figure 24 shows 900 consecutive iterates, while figure 25 shows these same 900 points plus 1800 more. For comparison, in figure 27 we plot just 50 consecutive iterates. The approximate outer dimensions of L(TA) aret
Smax at a 077251, Smin at a 204610, amax at S 491266, () amin at S
We now contrast with L(TA) the limit sets L(At)(TA) belonging to the generalized transformation TA(Xt). If we set At = 0.9931, we get a limit set entirely different from L(TA) (from the same initial point). This is shown in figure 28. It exhibits what we have called "pseudo-periodic" structure, that is, almost all the iterated images
* There is another repellent fixed point at a vertex of the triangle, namely S=a= . ** See section II. *** This is the region 0.455 <S< 0.525, 0.225 < a < 0.278. t These results are based on a calculation with N = 9600 iterations.
blank
of the initial point p are concentrated in the neighborhood of seven distinct "centers"-an example of class III limit set.*
With a very small change in At-namely, by setting At = 0.9930we find instead a period of order 7. This is shown in figure 29. As we decrease At in small steps down to At = 0.9772 (figure 30), the corresponding L(At)(TA) remains a period of order 7; the coordinates of the individual points appear to change continuously with At. For At = 0.97713, L(At)(TA) is again a pseudo-period, and this character persists down to At = 0.9770 (figure 13 of section II). Below** At = 0.9770 L(At)(TA) is a closed curve*** around the fixed point which shrinks in more or less continuous fashion as At is decreased. Figures 31, 32, and 33 illustrate, respectively, the limit sets for At = 0.97, 0.94, and 0.92. Finally, for At < 0.9180154 (see below), the limit set consists of a single point-the fixed point (10).t This peculiar behavior of L(At)(TA) as a function of At is not an isolated instance, nor is it by any means among the most extreme examples we have encountered (see section IV for a considerably more "pathological" case). Within the class of cubic transformations we have studied, it seems to be an empirical rule that the more pathological the limit set looks for At = 1, the more complicated will be its behavior as At is decreased.
2. Attractive and Repellent Fixed Points. The fixed points of a cubic transformation in the standard form (4) are those real roots of the algebraic system (6) which lie in or on the boundary of the S, a triangle.tt We are interested both in finding the values of the coordinates of these fixed points and in determining whether the points are attractive or repellent. By attractive we mean as always that, for any point p in a sufficiently small neighborhood of the fixed point po, the sequence T()(p) will converge to po. A general criterion for the attractiveness of a fixed point has been given by Ostrowski,10 viz: let lAmaxl be the largest eigenvalue in absolute value of the jacobian matrix evaluated at the fixed point. Then if )Amaxl < 1, the point
* See section II for this classification scheme. ** We have not attempted to find the critical values of At with greater precision, though this could in principle be done to, say, 7 decimal places on MANIAC II. *** This is, of course, only a conjecture. See the discussion in section II.
t In the language of functional analysis, TA(at) is a shrinking operator in this range of At values. tt Brouwer's theorem assures us that there is at least one fixed point.
blank
is attractive; if jAmaxI > 1 the point is repellent. The theorem says nothing about the case nor Amaxl = 1, nor does it yield a method for determining theoretically the appropriate neighborhood. For the two-variable transformation (4), we may give the eigenvalues of the jacobian matrix explicitly: T ±T-4J (12) 2
where To is the trace: OF OG To = ~S + 0a SSo (13) a=a and Jo is the jacobian: OF oGOG OF d OS Oa OS da S=So' (14) For the modified system (5), we find correspondingly: Amod = 1 -t + (To± ) . (15)
If the roots are complex, i.e., if To2 < 4Jo, we have I Amod 1 - \t(2 -To) + A2 (1 -To + J) (16) Defining Atlim as the value of At for which I Amod 12= 1, we obtain Atim = 2-To (17) 1 - TO+ Jo Thus, for the case of complex roots, we may make the fixed point (So, ao) attractive by choosing At such that 0 < At < Atlim . (18) Similarly, if Amax is real and negative, and l Amax [ > 1, Atlim -+2 - mx (19)
It is clear that this artifice will not work if Amax > 1. Such a situation arises in the one-dimensional case: x = x3+ 3xx2 , x + 3 l, x + x 1 (20) that is xl = x2(3- 2xi) (21)
The fixed points are xl= 0, 1/2, 1; at these points the derivative dx ldxl, has the values 0, 3/2, 0. Clearly both xl = 0 and xl = 1 are attractive fixed points; for all x(l)> 1/2, x(n - 0, while for all x( ) > 1/2, (n) - 1. The interior fixed point xl = 1/2 is repellent and cannot be made attractive by using the At-modification. The corresponding situation does not seem to occur for any of our cubic transformations in the plane.
In practice, all one has to do to obtain the numerical value of a repellent fixed point is to choose a sufficiently small At and iterate; on a computer, this calculation requires only a few seconds.*
3.Attractive Periods. The set of points constituting a period of order k are fixed points of T(k). Thus one may test whether a periodic limit set is attractive by applying Ostrowski's criterion to T(k). Let as(n)aS(-)jo(n) oS 1 a J(n) _ as(n)da(n)(22).S3 a -(So, a0 ) OS(n)OS(n) aS(n-l)aa(n-1) (23) Jn-1-Oa(n)Qa(n) _aS(n-I) Oa(n-1) S(n-1), a(n-1)
where, e.g., (S° , a° ) is a fixed point of T(k). Then, by the chain rule: J(k) = k-1 x Jk-2 ... x J (24)
Thus J(k) is easily obtained by evaluating (24) over the periodic set in question; the application of Ostrowski's criterion is then immediate.
* Early in this investigation we made the "mistake" of taking At > Atlim in a few cases, and thereby discovered the interesting limit sets Lp(/t)(T).
We have often used this technique to convince ourselves that the periods are really limit sets and not the result of spurious or accidental convergence.*
IV—
Modification of the Coefficients
1. In this section we present some results on the effect of modifying the original integer-valued coefficients of our cubic transformations in three variables. That is, we consider, as before, transformations of the form: 10 ' =dijMj, (1) j=l 3 d=ijl, (2) i=l
but we no longer require that the dij all be 1 or 0. As already remarked in the introduction, if we impose on the dij only the additional condition: 0 <dij< 1 (3)
then (1), (2), and (3) define a class of cubic transformations depending on 20 parameters, e.g., dlj, d2j (j = 1 to 10).** Since we are unable to formulate a complete theory for the finite subclass of transformations characterized by the restrictions: dij = 0 or 1, all i,j, it is clear a fortiori that we do not have a theory for the infinite class.
In this paper we limit ourselves to showing how an experimental study of some special cases can help to throw light on the properties of our original cubic transformations.
In effect, what we do is study certain transformations which are "close to" some particular transformations of our basic type. A natural
* This technique has actually been used for periods with orders k as large as 148. For very large k the method might fail owing to round-off errors or other numerical inaccuracies.
** For the definition of the cubic monomials Mj, see equation (4) of section II. The domain of this class of transformations is again the region 3 xi=- 1, 0 < Xi 1. i=l
way to define a transformation close to some given Tc,c2 would be to choose its coefficients as follows: dij-1 - eij,j C Ci , (4) dij = eij3, j Ci
where the dij must satisfy (2) and (3), and the eij are small. This class of transformations, defined with respect to some Tc,C2, is still too extensive to study, even if the various eij are restricted to a few discrete values. What we have actually done is to consider 20 such transformations, each of which depends on a single parameter e. We denote these by the symbol: T(r,s, 1< r < 10, 0 <e< , s = 0, 1. (5)
It is understood that these transformations are only defined relative to some TC 1c 2. For convenience we shall generally refer to the transformations T(r,s), defined relative to some Tc,c2 as associated transformations. The coefficients of the T(,,,), are specified as follows: For j ~ r: dij = 1, j c Ci, (6) dij = , j f Ci For j = r: dij - -e, rC 1, dir(1-s)e, r C1 , d 2r = 1-e, r C 2,d2r = se, r 02 C2
In words: T(r,s) is formed from Tc,c2 by the replacement Mr (1 -e)M wherever the term Mr occurs, and by adding eMr to one of the other two lines of the three-line schema. As an example, consider the transformation TA introduced in section III: TCA = {3, 5,8,9, 10}, C2 = {1,2,8}.
Relative to TA, T(5,1)E would read: x' = M3+ (1 - e)M5 + M7+ M9+ M 10, X2' = M1 + M2 + eMs + M8 (9) X3= M4 + M6,
while T(4 ,0) would take the form: x1=M3+ eM4 + Ms + M7 + M9+ Mo , 2 =M1+ M2+ M 8, (10) 3 = (1- C)M4+ M6,
In the S,a coordinates, the T(r,s), can be written: S' = F(S,a) + efrs(S,a), a' = G(S,a) + eg,,(S,a) ; (11)
the original Tc,c 2 is obtained from (11) by setting e = 0. For the two examples given above, we have: frs = f51-(S - a)(1 -S - a)2 e(5, 31)S : 2 (12) 951=-f51; frs = f40 = 12a2 (S - a) (13) T(4,0)e : ' (13) 940 = 0 -
It turns out that for these one-term modifications T(r,s)e we always have gr = ±frs or 0. frs can further be factored into a numerical coefficient crs and a function Mr(S,a); the Mr are of course just the original cubic monomials expressed in terms of S and a. The crs and g,s are determined as follows:
For r C C sgrs 0, Crs -1, rs -frs Crs -; ( for r C Cgrs -frs, Crs rs frs, Crs - ( for r t cl,r t C s rs O, Crs 1, s s frs Crs(
2. We have studied the modified transformations T(r,s)e for a variety of our original cubic Tc 1C 2 that happen to have infinite limit sets.
Our usual procedure has been to vary e in steps of 1/100 in the range 1/100 < e < 1/10 for a given T(,,s) relative to a given Tc,c 2, although on occasion intermediate values of e have been used. Only for the transformation TA [equation (8) above] have we looked at all 20 modified transformations. For a few other TC 1C2 we have limited ourselves to selecting certain of the associated T(r,s)e for detailed study. Since this selection has generally been made on intuitive grounds, we cannot claim that the most "interesting" modifications of the original transformations have always been considered. Nevertheless, this part of our study has proved most revealing, especially as regards the structure of class IV limit sets.
Before describing the results, we insert a few remarks on the difference between the two types of generalizations we have considered, the At-modification of section III and the associated transformations Tr,,s). The At-modification is essentially nothing but the application of a technique frequently employed in the practical solution of nonlinear equations by iterative methods; it is, in fact, one way-perhaps the simplest--of introducing a linear convergence factor. Apart from our use of this device for obtaining the coordinates of the fixed point, our principal interest is in small convergence factors (At close to unity)-too small, in fact, to produce convergence to the fixed point. In view of the fact that the At-modified transformation TC 1C 2(At) has precisely the same fixed point as the original transformation Tc1C 2, one might expect that there exists a close relationship between the corresponding limit sets Lp(At)(T) and Lp(T). In some sense this is true, as the examples given in section III show (see also below, subsection 4). We may express this more formally as follows:
We define a sequence of transformations Tc1c2(Ati)=TAt with corresponding limit sets Lp(Ati)(T) by some convenient rule: Ato = Atlim, Ati = - t (17)
The sequence Tt,, i = 1, 2, . . ., clearly converges to Tc,C2. We then formulate the following conjecture:
Given a Tc,c2and a 6 > 0, then, for all p in the triangle, there exists an N(p) such that, for i > N(p) and for all x C Lp(T), there exists a y C Lp(At )(T) satisfying I y - x < 6.
The modification of Tc1 c2 defined by the associated transformations T(r,s) differs from the At-modification in several respects. In the first place, the perturbation introduced is not linear. Furthermore,
the fixed points of T(r,s)c are in general not the same as those of TC0 C2 (fixed points on the boundary of the triangle may, of course, be common to T(,s),e and Tc,c2 for some pairs r, s).* Finally, each pair r, s must be treated separately; for fixed e, perturbations of different terms of TC1C2 may lead to quite different limit sets. Nevertheless, a conjecture analogous to that formulated for the sequence T xt would most probably turn out to be correct.
3. Limits Sets of the TransformationsT(r,s)e Associated withTA. Since we usually deal with values of e of the form: ei = 1 I, << 10, (18) we introduce a symbol to denote a set of such values: I(i,j)- {en}, i< n < j . (19) In addition, Rk will denote the closed interval of e: (r,s) Rk8 R-±k + k -Rk k(r,s) = [ (r,s) (rs)] ( (rs) < <+(rs) (20)
for which the limit sets L(r,,) of T(r,s), are periods of order k. The photographs illustrating the examples that follow will be found, suitably labelled, in subsection 2 of appendix II.
There is one significant feature common to all the T(,s)E associated with TA; for every pair r, s at least one periodic limit set-that is, a period of order k > 1- -was found in the range I(1, 10). The order of periodicity of most frequent occurrence was k = 7. Thus, for example, for (r,s) = (10,0), we found periodic limit sets with k = 7 over the range 1(3, 10), and the case (r, s) = (6, 1) behaves in the same fashion over the same range. For both series of associated transformations, the limit sets for e = 0.01 are of class IV type and closely resemble L(TA). At e = 0.02, bright spots show up in the pattern (figure A-1); this usually indicates that one is near a period, i.e., that a relatively small change in e will yield a transformation having a finite limit set. In the notation (20), this would be written: -R7 -0.02 < 1. In (r,s)
* The new fixed point S(r,s) S = S + As, a(r,s)e = a + Aa, calculated to first order in (AS, Aa), has both AS and Aa proportional to e. The ration AS/Aa in this order is therefore independent of e, though not of r, s. There are in fact, 6 possible directions of displacement, two for each of the three cases: grs = frs, grs =-frs, rs = 0, cf. equations (14), (15), and (16).
these two examples it happens that a period of order 7 is observed over the range 1(3, 10), that is: I(3, 10) cR L). This is not generally the case. Thus for the case (r, s) = (9, 0), I(2, 9) cR70, whereas L(9 ,o)0.01 and L(9 ,o)0.10 are of class IV type and are morphologically similar to L(TA). It may be recalled (section III) that an analogous behavior was observed for the At-modified transformations TA(at), namely that LAt(TA) was found to be a period of order 7 for a particular range of values of At (0.9930 < At < 0.9772), and different in character (actually, of class III type) outside the range on both sides.
Periodic limit sets of order 7 have been found for some range I(i, j) of e in 9 out of the 20 possible cases. For one of these, (r, s) = (2, 1), 1(4, 7) CR(2), while L(21)o0.io is periodic with k = 28 (figures A-2, A-3). In the transition region, i.e., for +R2, ) < e < -R(28 ) the -C) -cl 1 I ICIV(2,1)(2,1)2 limit sets are infinite. These are shown in figures A-4 and A-5 for the range 1(8,9). They look like pseudo-periods, but, when suitably enlarged, they are seen to be of class I type (figures A-6, A-7). In these pictures one clearly sees with increasing e the onset of instabilityto use an expression from mechanics-and the eventual attainment of a different stable state. The transition region at the lower end of the range also contains infinite limit sets. Figures A-8 and A-9 show (2,1)0.03, first to normal scale, then enlarged. It is manifestly a class III limit set.
For other T(r,s), periods of order k > 7 are found for certain ranges of the parameter e, viz.: k = 9, 16, 23, 30, 37, 46, 62, 148. In two cases, two periods of relatively prime order are found in different sub-ranges of I(1, 10). Thus T(5 ,l)e has two periodic limit sets, one with kl = 23 for e = 0.01 and one with k2 = 16 over 1(9, 10). Similarly, T(1,i1 ) has a periodic limit set with kl 16 for e = 0.01, and one with k2 = 9 for e = 1.10. In these cases the dependence of the limit set on e in the transition region +fRkl< e < Rk2 is more complicated than that t(r,s) (< < r,s) described above. For e-values in this region and sufficiently close to the end-points we observe the expected pseudo-periodic limit sets. For values of e not too close to either boundary the limit set may be either of class IV or of class I type. Figure A-10 shows L(5 ,1 )0.04 to normal scale; in figure A-11, a portion of the limit set is shown enlarged.
We conclude this subsection with two further examples. These illustrate a phenomenon previously mentioned in our general description of limit sets (section II), namely the coexistence of finite periods and class IV sets. Figures A-12 and A-13 show two distinct limit sets belonging to T(3 ,1 )0.01 . One is a period of order k = 23, while the other is a class IV set closely resembling L(TA). The same phenomenon
is perhaps more strikingly illustrated by the case of T(5 ,0)0.02. Here we find both a class IV limit set and a period of order k = 148 (figures A-14 and A-15). We can say virtually nothing in this case about the dependence of the limit set on the initial point. Current computing facilities and techniques are not sufficiently powerful to effect an acceptably accurate determination of the respective regions of convergence without using prohibitive amounts of computing time. We have, however, carried out a few numerical experiments, the results of which certainly confirm our first impression that the geometrical structure of these regions is immensely complicated.
4.Study of theAssociated Transformations for OtherTcIC 2. In this subsection we discuss a few additional examples to illustrate the dependence of infinite limit sets on the parameter e. The relevant photographs and tables will be found in appendix II.
For our first example, we choose the transformation: T Ci-{2,4,6, 7, 9}, TCI2=TBC2 = 5,8,10} (21)
The class IV limit set L(TB) belonging to this transformation is shown in figure B-1. As is evident, it consists of three separate pieces. Each of these is, of course, a limit set for T ). It is instructive to compare the limits sets LAt(TB) with those belonging to certain of the associated T(r,s),. In appendix II we list the results for only one case; (r, s) = (1,0). The limit sets LAt(TB) and L( 1,o), are described in table B. There are (at least) three ranges of At values for which LAt(TB) is periodic; for At close to unity the behavior of LAt(TB) as a function of At is rather wild. As At approaches Atlim = 0.854320 the (class I) limit set shrinks in a continuous manner. The behavior of L(1,0) as e is varied over 1(1, 10) is, if anything, more "pathological"; there are at least six different intervals R(l'°) for which the limit set is periodic, and each period has a different order. Note the similarity in appearance between the two class IV limit sets: LAt(TB) (At = 0.994) and L(1,o)o.oi.
The next two examples may be taken together: C: {2, 7,8,9,101 T C2 {4, 5, 6} (22) T C = {2,5,7,8,9} (23) C2={4,5,10}.
The basic class IV limit sets L(TD) and L(TE) are shown in figures D-l and E-l; their morphological resemblance is apparent. The behavior of the LAt and L(1,0)i for these two cases is set forth il the tables and photographs of appendix II. Detailed comment is perhaps superfluous at this state of our knowledge; we limit ourselves to drawing attention to the following comparisons:
1. Compare L Tt(TD) (At = 0.97) with L(l,o)o.lo(TD). 2. Compare LAt(TE) (At = 0.97) and At = 0.96) with L(1,o)o.o 9(TE) and L(1,o)o.lo(TE).
5. The original transformations TB, TD, TE are closely related from the point of view of formal structure. TD and TE differ by exchange of a single term between the defining sets C1 and C2, while each of these goes over into TB under the simultaneous interchange of two terms between C1 and C2. A comparison of the associated limit sets for TD and TE shows that the initial similarity of L(TD) to L(TE) is roughly preserved under perturbation. This suggests the possibility that some meaningful classification based on algebraic form might be devised.* Of even greater interest is the correspondence, in these examples, between the LAt and the L(1,0), for some ranges of the respective parameters. We are not at present in a position to draw any significant conclusions from the existence of this correspondence; it seems likely, however, that a closer study of these examples would yield criteria enabling one to predict such behavior.
6. There is one property of these transformations which may safely be inferred from the data, namely, that they are close to transformations having periodic limit sets (for some common set of initial points), where close is to be interpreted with reference to some appropriate parameter space---e.g., a range of e values of At values. Their limit sets are "close" to periods, not in the sense that pseudo-periods are, but rather by virtue of the fact that they contain points which lie closeperhaps arbitrarily close--to a set of algebraic solutions of Tk(p) = p. In other words, the Hausdorff distance between the set of period points and the limit set L is small. In this connection, the following piece
* The difference in behavior of Lat(TB) on the one hand and LAt(TD), LAt(TE) on the other is undoubtedly due in part to the fact that in the first case the jacobian matrix has complex eigenvalues at the fixed point, while for TD and TE the eigenvalues are real; this is probably sufficient to explain the qualitative difference of behavior of the corresponding LAt as At - Atlim for At - Atlim sufficiently small.
of evidence may be presented. Consider the transformation T(1.o)o.oi associated with TE, for which we have observed that the sequence T((,)ool(P ), n = 1,2,..., converges to a period of order k = 10 for almost all p. Let us choose a p close to the fixed point. If we then examine the sequence for n = 1, 2,..., N, where N is sufficiently large, we find that the images T() (p) of p have traced out a pattern which closely resembles the original class IV limit set L(TE) of figure E-1. This is shown in figure E-2. The bright spots are the points belonging to the periodic limit set L1 ,o)o.01 . Presumably this means that the effect of introducing a small perturbation into TE, of the form specified by T(.0o).o 01,* is to make the limit set L(TE) contract to 10 points. Alternatively, we could say that, as e - 0, the periodic limit set L(l,o)E(TE) spreads out until it becomes the class IV limit set L(TE).
This and other similar examples suggest that it might be useful to consider the periodic limit sets as fundamental, the hope being that one could develop an appropriate perturbation method, taking these periods as the unperturbed states. The effect of a small change of a parameter (in the direction of instability) is then simply to make the period non-attractive. This can in principle be studied by purely algebraic methods. Determining the structure of the resulting limit set-the perturbed state-is of course a more difficult matter.
In some cases this may amount to nothing more than the development of improved techniques for handling algebraic expressions of very high order. To clarify this statement, we offer one further example. Consider the Atmodified transformations TAt(TA), where TA is a transformation introduced in subsection 2 above. For At = 0.99300, LAt(TA) is a period of order 7. With a very small change in At- namely, for At = 0.99301-the limit set is of class III type, a pseudo-period. Rather than investigating TA(,t) let us turn our attention to the seventh power of the nmodified transformation, TA(at)(p) (At = 0.99301). If we choose our initial point p sufficiently close to one of the (repellent) fixed points of T(,t) ** we find that the first 516 iterated
* If TE is written in the form: S' =F(S, a), a' =G(S, a), then T(1,o), is S' = F(S, a) + 6(S - a)3 , a' = G(S, a). ** The actual values are not known: we have not yet developed good techniques for finding the coordinates of the points of a non-attractive period. The initial point for this example was taken as: S = 0.7034477, a = 0.1159449, chosen on the basis of some simple numerical experimentation. It is close to one of the periodic points belonging to the limit set LA(at) (At = 0.99300), viz.: S = 0.7037400, a = 0.1157123.
images of p, T(Z7t) (P), n = 1, 2,...,516, lie on an almost exact straight line in the S,a triangle. This is shown in figure A-16. The initial point p is at the lower right, and the successive images trace out the line continuously from right to left.* If we continue the iteration, we find that the later images deviate from the straight line, then oscillate in position, and finally settle down to generate another straight line with a different end-point--presumably very close to another fixed point of T(7) It is clear that if one had powerful enough algebraic tools, one could calculate this linear behavior.
7. We close this section with two remarks: 1) A study of the TAt and T(r,s)e associated with those Tc1C2 which have only class I limit sets indicates that the latter are much more stable with respect to these one-parameter modifications than are the limit sets discussed above. 2) Even these unstable limit sets appear to be stable with respect to some one-parameter perturbations of the corresponding transformations. Thus, the transformations T(3 ,0) associated with TE have limit sets visually identical with L(TE) over the whole range I(1, 10). Anomalies such as these make general pronouncements about absolute stability (or instability) impossible.
To illustrate: one rlight be tempted to explain the observed stability in this case as follows: Explicitly, T(3,0)e has the form: S' = F(S, a) +e(1- S-a)3 , a' = G(S,a) . (24)
Now the density of L(TE) is relatively large near the right-hand boundary of the triangle, S + a = 1. The perturbing term, however, vanishes on this line. Thus the transformation is on the average very little altered by the perturbation. But this "explanation" becomes less convincing when one looks at other transformations associated with TE. T(2 ,l1), for example has the form:
One would expect the same argument to apply here, but in fact the limit sets only resemble L(TE) over the two ranges 1(1, 2) and 1(6, 10). In between, we get the familiar periodic and pseudo-periodic behavior.
* The final point plotted has coordinates: S = 0.7030206, a = 0.11628136, so the slope of the line is roughly Aa/AS r -0.713. For this photograph, the scaling factor is approximately 2340.
V—
Relation to the Theory of Differential Equations
1. As we remarked in section III, the non-linear transformations discussed in this paper exhibit certain analogies with systems of differential equations. In the following we confine ourselves to discussing the plane case.
An important study in the theory of differential equations, particularly as applied to non-linear mechanics, is that of so-called autonomous systems11 '12 '13dx dy d=P(x,y), y = Q(X). (1)
The theory, initiated by Poincare, seeks to determine the properties of the solutions of (1) under very general conditions, and to deduce such properties for particular cases without actually solving the equations explicitly (i.e., obtaining the general integral). In particular, the trajectories, given parametrically as a function of t: x = x(t), = y(t) (2)
are investigated from a topological point of view. Fundamental is the classification of the singular points of the system (1), that is, the points x,y, where P(x,y) = Q(x,y) = 0. The behavior of trajectories in the neighborhood of singular points can be found by consideration of the linear approximation to (1); the real object of the theory, however, is to characterize and, where possible, predict behavior in the large. One of the most interesting phenomena connected with behavior in the large is the existence of closed trajectories, or limit cycles. The theorem of Poincare and Bendixson14 gives sufficient conditions for the existence of such. Unfortunately, the fulfillment of these conditions in particular cases is often hard to verify; to date no satisfactory theoretical method for dealing with an arbitrary given system has been found.*
2. If we write our general two-dimensional system of non-linear difference equations in the form: S— - S(n-)S(n-l) + F(S(nl) a(n-l)) ~~~~~~~At ~~~~(3) a( a(-) a(n ) + G(S(n-), a(n-)) At
* See reference 15. The practical applications are largely confined to stability theory. Also reference 13 and the literature there cited.
the analogy with (1) is evident. The fixed points of (3) correspond to the singular points of (1), and the behavior of solutions in the neighborhood of a fixed point can be investigated via the linear approximation; this procedure, in fact, yields Ostrowski's criterion (see section III). If the fixed point is attractive, the asymptotic solution in its neighborhood can of course be obtained. In the case of repellent fixed points (or if the initial point is outside the region of attraction of all attractive fixed points), the sequence of iterates sometimes converges to a limit set which appears to resemble a Poincar6 limit cycle, i.e., a closed curve. In other cases, finite limit sets (periods) are obtained; on the other hand, one may observe limit sets of quite ambiguous geometrical, not to say topological, structure. These last two alternatives have no analogues in the case of differential equations.
In fact, the analogy between (3) and (1) is more apparent than real. The significant distinction lies, perhaps, in the fact that for our difference equations there is nothing corresponding to the trajectories of (1); successive iterates do not in general lie close to each other. This fact makes it difficult to use topological arguments to determine the character of the limit set. For sufficiently small At the sequence of iterates may resemble a trajectory to some extent, but the limit as At 0 is almost certain to be a single point.*
VI—
Broken-Linear Transformations in Two-Dimensions
1. For certain special quadratic transformations in one dimension one can give an almost complete discussion of the iterative properties; this is possible because these transformations are conjugate to piecewise linear (broken-linear) mappings of the interval into itself. For example, the transformation: x' = g(x), where g(x) = 2x, 0<x< 1/2; g(x) = 2 - 2x, 1/2 <x< 1.** The iterative properties of the latter can be obtained from a study of the law of large numbers for the elementary case of Bernoulli. Stated differently, the behavior of iterates of this simple quadratic transformation turns out to depend on combinatorial rather than analytic properties of the function. With this in mind, we tried to see whether an analogous situation would obtain in two dimensions. Our non-linear, polynomial transformations of a triangle into itself might, we thought, be similar to suitably chosen broken-linear mappings of a square into itself, at least as regards their
* It may happen that some power T(n of a transformation more closely resembles a trajectory; cf. the example cited in section II. ** See further in appendix I.
asymptotic behavior. One simple generalization to two dimensions of broken-linear transformations in one variable is a mapping: x = f(x,y), = g(x,y) , (1)
where each of the functions f and g is linear in regions of the plane. In other words, the graphs of these functions consist of planes fitted together to form pyramidal surfaces. The motivation for studying such transformations is the hope that their iterative properties will turn out to depend only on the folding of the plane along straight lines or, more specifically, on the combinatorics of the overlap of the various linear regions which is generated by the mapping. The simplest nontrivial case to investigate consists in taking f(x,y) as a function defined by choosing a point in the square and making f maximum at this point, the function being linear in the four triangles into which the square is divided. g(x, y) is defined in an analogous manner.
Each of the functions f(x, y), g(x, y) is thus made to depend on three parameters. Thus for f we choose a point xl, yl in the square and erect a perpendicular of height 0 <dl< 1 at this point; this defines a surface consisting of four intersecting planes. The transformation can then be given explicitly as follows:
where the regions I to IV are specified by the bounding lines:
then region I is bounded by L 1, L 2, and x = 0, region II is bounded by L2, L3, and y = I , (region III is bounded by L 3, L 4, and region III is bounded by L3,L4, and x = 1 , region IV is bounded by L4, L1 , and y = 0 . Analogous equations hold for y' = g(x, y), with parameters x2 , Y2, d 2.
Of the several transformations of this type that we have studied numerically we mention only the following: 1 1 T Xl-= , y i , dl =0.95, T •' 3 3 0 -95(5) x2= 0.6,y2= 0.5, x = 0.5,yi= 0.9 , di= ,5 2=0.3,Y2= 0.7,d2= 0.8
and the one-parameter family: X1 iYi= Z, di= 1, x22 =2 = - z,d2 = 1.,
The limit sets are shown in figures H-l, H-2, and H-3 through H-17 of appendix II.
L(T2) (figure H-2) is, in a sense, analogous to the class I limit sets we observed for some of our cubic transformations in three variables. In contrast, L(T1) (figure H-1) represents a new phenomenon-a connected "curve" (in this case, a collection of line segments) that does not close. More interesting, however, is the behavior of the limit sets L(Tz) as z varies from 0.49 down to 0.01.* Initially L(Tz) is an open cycle like
* The case z = 1/2 has not been studied. The reason for this is technical. Straightforward iteration will always produce sequences which degenerate to zero in a finite number of steps, owing to the fact that every iteration involves multiplication by 2. In a binary machine, this operation is a "shift" to the left. A sufficiently long chain of such left shifts will always result in zero. If one wants to study this case, one must replace multiplication by 2 by some arithmetically equivalent operations, e.g., multiplication by C/(1/2C), where C is not a power of 2.
For this case (z=1/2), the problem becomes, of course, one-dimensional. The iterates remain on the line x = y, and the limit set is identical with that of the transformation x' = g(x) introduced in subsection 1 above.
L(Ti). With decreasing z, the limit set becomes more complex, until it resembles a class IV limit set (e.g., z = 0.27). With further decrease of z, the limit set appears to contract; the tails get shorter, and the points cover some sub-region of the square more and more densely. One interesting question-unanswered at the time of this writing-is: does Tz become ergodic for some range of z values? In figure H-10 we see 1000 (consecutive) points belonging to L(Tz) for z = 1/4. Figure H-11 shows these same 1000 points together with the next 2000 points. It is evident that the region containing L(Tz) is filling in. In our opinion, this is a strong indication of ergodicity.
For lower values of z, the same ergodic behavior is observed, until, at z = 0.15, the limit set splits into four disconnected pieces (figure H-15). As z is further decreased, these four pieces shrink; by the time we reach z = 0.01, the limit is nearly a single point.
2. It is clear that one can devise broken-linear transformations that are dense in the unit square; one may take, for example, product transformations with independent coefficients and use the one-dimensional result for each factor. For transformations of the type considered in this section, however, it is not easy to determine a priori what the limit set will be. It should be emphasized that there is no hope of demonstrating that our polynomial transformations in three variables are exactly conjugate to some two-dimensional broken-linear transformations. Presumably there is no such conjugacy. Nevertheless, one might hope that a somewhat weaker notion of equivalence than that of strict conjugacy could be introduced.
One suggestion along these lines is the following: define two transformations T and S to be asymptotically similar if for almost every initial point p the limit set Lp(T) is topologically equivalent to Lp,(S) for some suitably chosen p', and vice versa. Thus, for example: if for any two transformations T, S, the sets of iterates T , Sn) are dense in the (common) domain of definition for almost every initial point, then T and S are asymptotically similar in the above sense. Another special case in which two transformations, T and S, are asymptotically similar is when each transformation possesses just one attractive fixed point, the region of attraction being, in both cases, the whole space.
As we remarked above, in the case of broken-linear transformations the asymptotic behavior of iterates depends only on the combinatorial structure of the subdivisions of the fundamental regions (triangles) under repeated folding. Just how complicated this can be is shown by the behavior of the limit sets L(Tz) belonging to the one-parameter family Tz discussed above. To date we have not managed to devise any good method for handling the Boolean algebra of these iterated
intersections.
Appendix I
In this appendix we collect some general remarks about the process of iterating transformations, particularly in one dimension. We also discuss, in some detail, a few special one-dimensional transformations which we have had occasion to study.
1. One of the first simple transformations whose iterative properties were established is the following: x'= f(x) = 4x(1 - x) To obtain these properties we consider, instead of (1), the broken-linear transformation:* 2x for x< 1/2 , x: =g(x), g(x) = 2(1-x) forx > 1/2 (2) 2(1\- x) for xa> 1/2
The study of the iterates of this transformation is equivalent to investigating the iterates of a function S(x) defined as follows: if x = O.a 1a 2a 3... an..., where the ai are either 0 or 1, (3) then S(x) =.a2a3a4 .... (4)
In other words, S(x) is merely a left shift of the binary word x by one place. The iterative properties of S(x) are in turn deducible from the law of large numbers in the case of Bernoulli. In effect S(i)(x) falls into the first half of the interval if and only if ai = 0. The ergodic average E F(l)[S)(x) i=l is therefore the same as the fraction of ones among the ai for 1 <i < N. F(I) is the characteristic function of the interval [0,1/2].
The relation between (1) and (2) is that of conjugacy: there is a homeomorphism h(x) of the interval [0, 1] with itself such that g(x) = h[f[h-(x)]] . (5)
* This transformation has already been mentioned in section VI.
Thus the study of the iterates of the quadratic transformation (1) reduces to the corresponding study for the broken-linear transformation (2). In this case, h(x) can be written down explicitly:15 2 h(x) =- sin-l(V/) . (6) 7r
2. The Set of Exceptional Points. In the case of the function f(x) = 4x(1 - x), it is true, then, for almost every* initial point, that the sequence of iterated images will be everywhere dense in the interval, and what is more, the ergodic limit can be explicitly computed; it is positive for every sub-interval.
There exist, however, initial points x such that the sequence x, f(x), f[f(x)],... is not dense in the whole interval [0,1]. Obviously, all periodic points, i.e., points such that, for some n, f(n)(x) = x, are of this sort. It is interesting to notice, however, that there exist points x for which the sequence f(i)(x), i = 1,2,..., is infinite without being dense; there are, in fact, non-countably many such points. To show this we consider the equivalent problem of exhibiting such points for the function S(x) introduced above. The construction then proceeds as follows. Consider a point x = 0.ala2 ... an ... We define a set Z consisting of all those x's which have an = an+i =an+2 for all n of the form n = 3i. In other words, the set Z consists of points which have every binary digit repeated three times, the sequence being otherwise arbitrary. Consider now the transformation x' = S(x), where S(x), as defined in (4) above, is a shift of x one index to the left. We now look at the subinterval from 0.010 to 0.011. Starting with any point in Z, it is clear that no iterated image will fall in this sub-interval; no three successive binary digits of a point in Z are of the form (010). It is easy to see that Z contains non-countably many points; in particular, it contains non-periodic points, so that the set of images S(i)(x), i = 1,2, ... is infinite, but not dense in [0,1].
Presumably one can find points in Z for which the ergodic limit exists. The measure of the set Z is zero, but relative to Z, the set S of those points for which the sojourn time exists still form a majority. We may define majority either in the sense of Baire category or as follows. Take points na (mod 1), n = 1,2,..., where a is an irrational constant. Consider the set N1 of those n's for which na (mod 1) belongs to Z, and also the set N 2 of those n's for which na (mod 1) belongs to S. We then say that the points belonging to S form a majority of those points belonging to Z if the relative frequency of N2 in N1 is one.
* Almost every is to be understood in the sense of Lebesgue measure.
The behavior exhibited by the points belonging to Z is more general, in the sense of measure, for some other transformations of the interval [0,1]. It is possible to give examples of continuous functions such that, for almost every point, the iterated sequence will be nowhere dense in the interval, although the sequence does not converge to a fixed point.
3. A Remark on Conjugacy Theorem.Let g(x) be the broken-linear function of equation (2) below, i.e., 2x, o<x<g(x)-21 (2) 0 g(x)= 2(1 - x), 2 >x > 1
Let t(x) be a convex function on [0.1] which transforms the interval into itself, and such that t(O) = t(1) = O. For some p in the interval, we must have t(p) = 1; by convexity, there is only one such point. Consider the lower tree * generated by the point 1. The necessary and sufficient condition that t(x) be conjugate to g(x) is that this tree be combinatorially the same as that generated by 1 under g(x), and that the closure of this set of points be the whole interval, i.e., that the tree be dense in [0,1].
The condition is obviously necessary, since the point 1 generates a tree under g(x) which is simply the set of binary rational points. Under any homeomorphism h(x) which has to effect the conjugacy, the point 1/2 must go over into p, and our assertion follows.
To prove sufficiency, we construct h(x) in the following manner. We take h(1/2) = p by definition. We next chose h(1/4) to be the smaller of the two values of t-1 (p); h(3/4) is then by definition the larger of these two values. We then take h(3/8) to be the smaller of the two values of t- [h(3/4)], and so on. Proceeding in this fashion, we thus construct a function h(x) defined for all binary rationals. It remains to prove that we can define it for all x by passage to the limit. This, however, follows from the assumption that these points are dense in [0,1] and that their order is preserved. The function h(x) will obviously be monotonic, and, being continuous, will possess an inverse h-l(x). From our construction it then follows that h[g(x)] = t[h(x)].
* By the lower tree of p (under f(x)) we understand the smallest set of all points with the following properties: (a) The set contains the given point p. (b) If a point belongs to the set, then so do all its counter-images under f.
4.Broken-Linear Transformations. In one dimension these are functions f(x) that are continuous on [0,1] and linear in subintervals of [0,1]. We assure that the graph of the function has a finite number of vertices, i.e., that f(x) consists of a finite number of lines fitted together continuously. For these broken-linear transformations one certainly expects that the ergodic limit exists for almost every point. For example: if one considers the sequence of iterated images T(W)(p), then, for almost every initial point p, the time of sojourn should exist for all sub-intervals, i.e., N lim N fR[T(i)(p)] i=l
should exist for almost all p and all measurable sets R; here fR is the characteristic function of the set R.* The value of this limit may indeed depend on the initial point p; it is likely, however, that all the points of the interval can be divided into a finite number of classes such that, within each class, the value of the limit is the same.**
There is another finiteness property that these transformations may possess. Given n, consider all broken-linear transformations which have at most n pieces (i.e., the space divides into n regions, in each of which the transformation is linear). Then it may be conjectured that there are only a finite number of different types of such transformations, where any two transformations of the same type are asymptotically similar (in the sense defined above, section VI). In other words, according to this conjecture, the type (or class) that a given transformation belongs to does not depend on the precise numerical values of
* We should perhaps mention here a more general conjecture. Suppose T is a polynomial transformation of the sort described by equations (10) to (13) of the Introduction. We then conjecture that the sequence of iterated images T(n)(p) has the following property: Let C be any cone of directions in n-space, and let fc(p) be the characteristic function of this cone, i.e., fc(P) = 0 if p does not belong to C, fc(p) = 1 otherwise. Then, for almost every p, exists.(l6)
** See pages 71, 72 of S. Ulam, A Collection of Mathematical Problems, New York 1960. An analogous conjecture can be made concerning the actual limit sets Lp(T) for four cubics in three variables.
the coordinates of those points where the derivative is undefined (corner points), but is determined only by the combinatorial structure of the subdivision of space into linear domains. In one dimension this means that the type of a transformation is determined by the number and interrelation of the nodes in the graph of the function, and not by their precise location.
5. Numerical Accuracy. The machines we use to compute the iteration process work with a fixed number of significant digits;* in MANIAC II, for example, this number is eight. It is therefore clear that any direct, single-step iterative process carried out on this computer will exhibit a period in not more than 108 steps. Given an algorithm which is iterative and of first order [i.e., the nth step depends only on the (n - )th], the process will, with great probability, exhibit a period which is much shorter. Statistically, one can reason as follows. If we assume a random distribution of, say, the last four digits of all computed numbers, then the probability that the cycle will close long before the full theoretical run of 108 steps is extremely close to one. Indeed, after producing numbers A1, A2,...,Ak, the chance that Ak+1 will be equal to one of the preceding numbers is of the order of k/108 . If we continue the calculations up to A2k, the chance that at least two numbers in the chain will coincide is approximately
/ k\k This is practically equal to 1 1-- if k -104 e
Clearly, going to 3k, 4k, ..., the probability that the chain will be cyclic gets very close to one. The situation is quite different in an iterative process involving two or more variables, e.g., computing (Ak+l, Bk+l) from (Ak, Bk). On a probabilitic basis alone, one expects to encounter periods of length ~ 108 (the maximum possible being 1016). In practice, this means that in two dimensions one does not expect to encounter periods of accidental or false periodicity unless one generates very long sequences.
* For any particular machine one can increase the accuracy by restoring to so-called multiprecision arithmetic. This can generally be done only at the cost of considerable loss in speed.
As the argument above shows, fortuitous periodicity can be a very real danger in one-dimensional iterative calculations. Indeed, we came across a striking example in the course of studying the asymptotic properties of the transformation. y = sin ry, 0 < y 1. (7)
This classical transformation is symmetric about y = 1/2, and maps the unit interval into itself. Furthermore, both fixed points-y = 0 and y = 0.7364845-are repellent. Now a simple argument shows that for such a function there cannot be any attractive periods of any finite order, that is, the derivative dTk (y) IY=Y dy
evaluated at any periodic point T(k)(yp) = yp is always greater than one in absolute value.* Nevertheless, when we performed the iteration on MANIAC II, the limit set we observed was invariably one of two finite periods, the first of order 1578, the second of order 6168. These periods were exact to the last available binary digit! This spurious convergence is undoubtedly produced by the complicated interaction of several factors, e.g., the particular machine algorithms for multiplication and round-off, our choice of finite approximation to the function sin7y, and so forth. When we repeated the calculation on the STRETCH computer (which works with about 15 significant figures), there was no observable tendency toward convergence to such periods.
If one is interested in obtaining the asymptotic distribution of iterates under some transformation like (7), one must either provide for more significant figures or resort to ingenious devices. One such trick which appears to work well and is relatively convenient-consists essentially in computing the inverse transformation. Since the functions
* It is essential for this argument that the transformation maps the whole interval into itself. The number of distinct periods of order k can then be easily enumerated. The conclusion follows on noting that dydT(k) (y)
must be the same for all points y belonging to the same period. If, however, the maximum value of the function is less than one, T(k) (y) can have relative minima above the y axis; then, indeed, there may exist attractive periods, i.e., the line y' _ T(k)(y) = y may intersect the curve sufficiently close to an extremum so that the kth order fixed point is attractive.
under consideration are two-valued, this involves introducing a random choice at each step. Specifically, taking some initial point po, we compute the sequence o, f-1)(po), f-2(po),....
Here the symbol f-k(po) implicitly contains the prescription that we choose one of the two values of the true inverse at random; thus the sequence f-i(p)_= f[f-(i-1)(p)], i = 1,2,..,
implies a sequence of random decisions as to which counter-image to choose at each step. If the calculation is carried out in this fashion, the chance of falling into an exact period is vanishingly small until we reach a chain length in the neighborhood of the theoretical maximum. Once having obtained our inverse sequence, we can conceptually invert it, that is pretend that we started with the last point and proceeded to po by direct iteration.*'**
6. In this sub-section we present some results of a study of a particular one-dimensional transformation:*** T,: y' = W(3 - 3W + aW2 ), W =3y(1 - y) . (8) This arises in a natural manner from a certain sub-class of our generalized cubic transformations in three variables: 10 Xi=dijMj, (9) j=l
' One cannot, of course, actually reverse the calculation and expect to reproduce the sequence. If one could, there would be no need for this stochastic device.
** This procedure presupposes a good method for generating random numbers. There are several of these which are well suited for use in automatic digital computers. One of the most common-and, in fact, the method used by us-is to generate the numbers by the chain: Rn+l - RoR,(mod 2k), where k is the binary word length of the machine, and Ro is some properly chosen constant. On MANIAC II, Ro = 513. This chain closes, but its length is greater than 1012. The lowest order binary digits are themselves not random, but this makes no difference in practice.
*** This transformation was introduced by the way of illustration in section I.
3 dij = 1, all j, 0 < dij 1 . (10) i=l Namely, we choose certain of the dij as follows: d13 -= , d17 = d19 = d34 = d 36 = 1, d 31 = 0.
The rest are arbitrary, except that they must, of course satisfy (10). If we restrict ourselves to the sub-class of initial points such that x3 = 0 (i.e., the side S + a = 1 of our reference triangle), the second power of the transformation can be written in the form (8). Ta is symmetric about y = 1/2, but it does not map the whole interval [0,1] into itself; its maximum value is Ymax = T ( 4 (11)
so this is the right-hand boundary of the invariant sub-region. The lefthand boundary is then the image of y' . In the range 1> a >0.9, the fixed point varies from yo = 0.8224922 at a = 1 down to yo = 0.8193719 at a = 0.90. Over this range, the derivative at yo is negative and greater than 1 in absolute value; thus the fixed point is repellent.
We have studied this transformation as a function of the two parameters-the initial point yo and the coefficient a--on the STRETCH computer. To study the asymptotic distribution, we divide the interval into 10000 equal parts and have the machine keep track of the number of points which fall into each sub-interval over the iteration history. Such a history was usually taken to be a sequence of length n x 105, with 4 <n< 9.* Should a period (of order k< 3 x 105) be detected by the machine, the calculation is automatically terminated. If the order of the period is not too great (k< 500), the values of the periodic points are printed out, and the value of
T(k) (y) dy is calculated over the period; if dy ( o 1 ,
* As we have already remarked (section I), the calculation of 2 x 105 iterates of this transformation requires about I minute on STRETCH.
this fact is strong evidence that the observed period is actually a limit set.
Our numerical investigation of (8) has been mostly restricted to the range 0.98 <ar 1.* Even in this restricted parameter range, the observed asymptotic behavior is of bewildering complexity. For a = 1, the distribution of iterates in the interval is extremely nonuniform. There are large peaks at the end-points of the allowed subinterval, with much complicated fine-structure in between, i.e., many relative maxima as well as sizeable intervals in which the distribution is locally uniform. This general behavior persists as a is decreased down to a = 0.9902. At a = 0.9901, however, a dramatic change takes place; most of the points concentrate in a few small sub-intervals. The limit set is, in fact, a pseudo-period.** At a = 0.99009, the pseudoperiodic behavior is still evident, but the occupied sub-intervals are larger. They contract, however, for a = 0.99008, and at a = 0.990079 a period of order k = 42 is found. Actually, for this particular value of the parameter a, there appear to be two possible periodic limit sets, with k = 42 and k = 84 respectively. The dependence on the initial point is, of course, quite complicated.*** For example, as yo is varied
over the values yo = 0.11, 0.12, 0.13,..., 0.21, the period with k = 42 is found in all except three cases, namely yo = 0.11, 0.16, 0.18, which apparently lead to the periodic limit set with k = 84. Both periods are numerically well-attested; they are exact to the last binary digit and have (d T()(y) < t
* Below a = 0.98, various periodic limit sets of order k< 24 were found on MANIAC II. We have not studied this parameter range on STRETCH. ** It may actually be a period with order k = 63049. This is the result indicated by the machine. At present we have no way of verifying this.
*** Subsequent numerical work has shown that this yo dependence is spurious. The calculations were performed with multi-precision arithmetic; in some cases as many as its decimal places were retained!
t It is difficult to decide whether this behavior is real or not. Some supporting evidence for its reality is the observed behavior for values of a very close to this critical value a = 0.990079. We find that, with a = 0.9900789, only a period with k = 42 is obtained, while on the other side, a = 0.9900791, the limit set is always a period with k = 84.
As we continue down in a, the limit sets are pseudo-periodic* until we reach a = 0.99007, for which a period of order k = 112 is found. This appears to be the approximate right-hand boundary of a periodic belt; that is, in the range 0.98990 <a< 0.99007 there are only periodic limit sets for this transformation. All these have orders which are multiples of 14 (with k< 112), except for the (approximate) lefthand boundary of the a-range (a = 0.9899) where a period of order 7 exists. The dependence of these limit sets on the initial point is complicated, and will not be reproduced here.
At a = 0.98988 there is another large discontinuity in behavior; we again find an asymptotic distribution which covers the whole allowable sub-interval in non-uniform fashion. No more periodic limit sets are found as a is decreased to a = 0.9800. The only phenomenon of note is the splitting of the limit set into two parts; this occurs somewhere between a( = 0.986 and a = 0.985. The resulting gap-which contains the fixed point-continues to widen as a -a 0.980.**
The results of this investigation-which is rather incidental and subordinate to the larger study reported in the present paper-clearly show that there is a great deal to be learned about the asymptotics of iterative processes, even in one dimension. It seems that "pathological" behavior is not a property of higher-dimensional systems alone. With regard to our study of the periodic limit sets, it may be argued that what we have really done is to investigate, in a rather indirect manner, the behavior of the roots of high-order algebraic equations as a function
* In this range, the machine detected some exact periods of huge ordere.g., k = 295148 (ac = 0.990079)! There seems to be no compelling reason to take this at face value. ** Such gaps have been observed in other cases. One interesting example is the transformation:
y' = W2 (3 - 2W), W = 3 y(l - y), which is also a special case of one of our cubics in three variables. For this transformation, the limit set consists of four separate pieces, as follows: (I) 0.3455435 < y < 0.4086018, (II) 0.4296986 < y < 0.5791385, (III) 0.7562830 < y < 0.8146459, (IV) 0.8220964 < y< 0.84375.
In the present case, i.e., that of T,, the existence of the gap can easily be predicted. Let yo(a) be the fixed point. There will clearly be a gap providing T2(ymax > yo(a), since in this case only one of the two inverses of yo lies in the allowed interval. Of course, the determination of the critical value ac of a from the equation TV2(Ymax) = Y(ff) would in any event have to be carried out numerically.
of their coefficients. This is certainly true, at least in part. It is therefore of interest to observe that the iterative method seems at present to be the only effective tool for treating this purely algebraic problem.
Appendix II
1. In this appendix we collect the photographs and tables illustrating the phenomena discussed in section IV and VI. The notation used has been described in sections II, III, IV, and VI. For convenience, some of the transformations are written out explicitly.
2. Modifications of the TransformationTA. In shorthand notation, this transformation is C1= {3, 5, 7, 9, 10} , C2 = {1,2,8} 1
In the S, a coordinates, this reads: 3S2 3 2 +S' = F(S, a) S3-6S2a-3Sa2 + 4a3+ 3Sa- +1, 2 2 a' = G(S, a) -S 3 + 3Sa2 + 2a3+S2-3Sa +a2 . (2) 2 2
The (repellent) fixed point has coordinates: So = 0.5885696, ao = 0.1388662.
The generalized transformations based on TA may be written: S' = (1 -At)S + At{F(S, a) + efrs(S, a)}, a' = (1 - At)a + At{G(S, a) + egs(S, a)} ;
the original transformation is recovered by setting At = 1, e = 0.
TABLE A Figure number At e r, s Comments A-6, Class IV limit set A-2, Period: k A-2, Period: k A-2, Class I limit set A-2, Class I limit set A-2, Scaled plot of part of A- Scaling factor .(< S< 485, < a < A-2, Scaled plot of part of A- Scaling factor ~(<S< 492, <a < A-2, Pseudo-period A-2, Scaled plot of part of A- Scaling factor ~(<S< 620, <a< A-5, Class IV limit set A-5, Scaled plot of part of A- Scaling factor ~(,<S< 64, < a < A-3, Period: k A-3, Class IV limit set A-5, Period: k A-5, Class IV limit set A-successive iterates of T()p) Initia point: So 7034477, ao Scaling factor ~
Modifications ofTB. In shorthand notation, this transforma- tion is C{2,4,6,7,, C{5,8, . In the S,a coordinates, it takes the form S' F(S, a) -a-S-a- aS a ( a'G(S, a) - -S - Sa - a-S- a.
The coordinates of the fixed point are So = 0.6887703, ao = 0.1592083.
The only associated T(r,s)e discussed is the case (r,s) = (1,0); for this case, the generalized transformation, written in the form of equation (4), has frs = f0o = (S-a) (8) grs = 910= 0 .
In table B below, all scaled plots show the region: 0.30 <S< 0.65, 0.20 < a < 0.36, i.e., the upper left-hand piece of the complete limit set (shown in figure B-1 for the unmodified transformation). The scale factor is ~ 2.9.
TABLE B Figure number At e r, s Comments B-Class IV limit set (points) B-- Scaled plot of part of figure B-B-Class IV B-Class IV Compare figure B- B-- Class IV B-- Period: k (points in this piece) B-- Class IV B-1, Class IV: separate pieces (shown here) Compare figure B-B-1, Period: k (points shown here) B-1, Class IV B-ll 1, Period: k (points shown here) B-1, Class IV B-1, Period: k (points shown here) B-1, Period: k (points shown here) B-Class IV; shows "transition" from period with k (At
blank
blank
blank
blank
blank
blank
4.Modifications ofTD andTE. These transformations are given by the schemes: T C1 - {2, 7,8,9, 10}, C2 {= 4,5,6} , C1- = {2,5, 7,8,9} (10) C 2 = {4,6,10} .
In the S, a coordinates, these read explicitly: S'= S3 + S2a+ Sa2+ a3- 6Sa-6a2 + S + a T2 2 2 2 2 2 3 9 2 2 9 32 2 3 a= S3+ S2a- -Sa2 - a-3 + 3a2 + S-a, (11)
with fixed points: So = 0.6525211 (12) ao = 0.3056821; S' = 9S2 a - a3 - 3S2 - 12Sa + 3a2 +3S + 3a (13), a' =-3S2 a + 3a3 + 6Sa - 6a2
with fixed point: So = 0.6444612, (14) ao = 0.3219578 .
TABLE D Figure number At e r,sComments D-- D-Compare figure D-D-96, Limit set for At superimposed on L(TD) D-- D-1, D-1, Compare figure D- TABLE E Figure number At e r, s Comments E-- E-1, Shows convergence to periodic limit (set k from initial point close to fixed point E-- Compare with E-5, E-E-Compare with E-5, E-E-1, E-1,
blank
blank
5. Broken Linear Transformations. These are described in section VI. Figures H-1 and H-2 show, respectively, 1000 points in the limit sets of T1 and T2 . The latter are specified as follows: 1 1 Ti l= , Y 1, dl 1 0.95, (15) X2 = 0.6, Y2 = 0.5, d2 = 0.95,
with fixed point: 1 2 xo= 2' Yo= 3 (16) and T x = 0.5, yi = 0.9 , d = (1 ,7 2 2 = 0.3, Y2 = 0.7, d 2 = 0.8 , with fixed point: 80 72 43' o= 143 (18)
The remaining figures, H-3 through H-17, show the limit sets belonging to the one-parameter family T,: T Xi=Y1 = z, d = 11 x2 =2 = 1-z, d2 = 19
blank
blank
blank
The identification is given by the following table (all figures except H-11 show 1000 points):
TABLE H Figure number z Comments H-For z - 5, see the remarks in section VI, note H-H-Compare figure H-H-H-H-Like "class IV" limit set H-Like "class IV" limit set H-H- Points of H-plus the next consecutive iterates H-H-H-H-H-H-Fixed point becomes attractive below z
References
1. J. Schreier und S. Ulam, Eine Bemerkung fber die Gruppe der topologischen Abbildungen der Krieislinie auf sich selbst, Studia Math. 5 (1935).
2. Menzel, Stein and Ulam, Quadratic Transformations, Part I, Los Alamos Report LA-2305, May, 1959 (Available from the Office of Technical Service, U.S. Dept. of Commerce, Washington, D.C.)
3. See for example, S. M. Ulam, On some new possibilities ... computing machines, I.B.M. Research Report 68, 1957. 4. See, e.g., Proceedings of Symposia on Applied Mathematics, Vol. X, American Mathematical Society (1958), or the article by Marshall Hall, Jr. in Surveys of Applied Mathematics, Vol. IV, 1958.
5. In the absence of a comprehensive reference, we refer the reader to the recent issues of Mathematics of Computation (1960-1962). See, e.g., Vol. 16, No. 80, October 1962, especially the article by D. H. Lehmer, et al.
6. See the article by E. T. Parker in Proceedings of Symposia in Pure Mathematics, Vol. VI, American Mathematical Society (1962). Parker's original construction supplied the final step in the disproof of a famous conjecture by Euler.
7. See H. Gelernter et al., Proceedings of the Western Joint Computer Conference, San Francisco, May 1960, pp. 143-149.
8. M. B. Wells, Proceedings of the IFIP Congress (1962). 9. C. Shannon, Philosophical Magazine 41, March, 1950. See also Kister, Stein, Ulam, Walden and Wells, Journal of the Association for Computing Machinery, Vol. 4, Number 2, April 1957.
10. A. Ostrowski, Solution of Equations and Systems of Equations, Chapter 18, (New York 1960).
11. A general reference is N. Minorsky, Non-linear Oscillations, Princeton 1962. More detail on theoretical points can be found in 12.
12. S. Lefschetz, Differential Equations: Geometric Theory, New York 1957.
13. Nemitsky and Stepanov, Qualitative Theory of Differential Equations, Princeton 1960.
14. See reference in footnote 11.
15. This result was first published by S. Ulam and J. von Neumann, American Mathematical Society Bulletin, Vol. 53 (1947), p. 1120, Abstract 403. See also the work of O. Rechard, Duke Mathematical Journal, Vol. 23 (1956), pp. 477-488.
16. Cf. the article by S. Ulam in Modern Mathematics for the Engineer, second series, edited by E. F. Beckenbach, New York 1961, p.280 .