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.