10—
Quadratic Transformations Part I:
With P. R. Stein and M. T. Menzel (LA-2305, March 1959)

This report is the original study of properties of iterations of non-linear transformations. It has given rise to a very large body of work, and by now, an extensive literature is still appearing on the subject at an ever increasing rate. (Author's note).

Abstract

This report deals with the properties of a restricted class of homogeneous quadratic transformations, with interesting physical and biological analogues, which we have called Binary Reaction Systems. All possible transformations of this class in 3 variables have been studied numerically on a computing machine, and the limiting behavior of random initial vectors under iteration of each of these transformations is tabulated. Some examples of 4-variable Binary Reactions Systems are also studied, and a few generalizations of the notion of Binary Reaction System are investigated for particular cases. Some remarks and results concerning the behavior in the large are presented, and examples of the mode of approach to the limit are given. Several of the more interesting phenomena are illustrated graphically.

The appendix deals with a different class of homogeneous quadratic transformations (of arbitrary dimension) which arise naturally from the study of a

simple evolutionary mode. For this class of transformations, the limiting behavior of arbitrary vectors under iteration can be given explicitly.

Introduction

This report summarizes and discusses some recent studies of the properties of quadratic transformations in several variables under iteration. This report is of an interim nature, and consists mainly in a presentation of "experimental" (i.e., numerical) results. A general theory of quadratic transformations (in contrast to the linear case) is essentially non-existent, and from the theoretical point of view, the work summarized below does little to improve the situation. It can only be hoped that as more facts become known, some outlines of a theory-at least a classification or "descriptive theory"-will emerge.

The motivation for the considerations which follow lies in the combinatorial problems suggested by genetic or biological systems. One has to deal with large populations of individuals (or particles) present in a given generation. Those may combine in pairs and produce, in the next generation, new particles. Suppose the original particles are of N different types. Given a rule for the type i (i = 1,... N) produced by mating of individuals of type j and k, the proportion or fraction xi of a given type in the next generation will be a quadratic function of the two fractions zj and xk.

More generally, one could consider a system (gas) of physical particles with N possible characteristics which collide in pairs and produce through the collision, say, a pair of particles with new characteristics. There could be many different values of the momenta, but possibly the "type" of the particle resulting from the collision could be different from the original ones.

Our present considerations concern the averages or expected values of the fraction xi in the next generation. The role of fluctuations or deviations from the fractions given by the quadratic formulae will be studied in a subsequent report.

As suggested in the title, this report is the first of a series (one might add, of unspecified length). By the time these remarks appear in print, much of the work will probably have been generalized and the results extended. It is hoped that at least some of the tentative conclusions presented below will stand.

I—
Homogeneous Quadratic Transformations

We begin by defining a homogeneous quadratic transformation as a set of N coupled non-linear first-order difference equations of the form: N Xi=Xi zke i=l,...N (1) k,Q=l where the -y are some real numerical coefficients with the property: Ik = ijk (2) In the present work we restrict ourselves to the case of non-negative coefficients: 7" > 0, all k, , i (3) Then if we choose the Xj as non-negative real numbers, the xj will also have this property. In the following we shall always restrict our Xj in this manner.

In order that the system (1) be of dimension N, we insist that for a given index i, not all ye can vanish, i.e., N k > 0 (4) k,Q=l

Systems of the form (1) (not necessarily with the restrictions of reality or non-negativeness) can be considered from two points of view:

(a) As the transformation of all vectors ( ) into vectors XN X IN or, more generally, as a mapping of some specified region X into a region X'. (b) As a set of difference equations which determine the value of X1I a vector at "time" n from the value of some initial vector X(^°) XN (XN)

If we take the first point of view we are, in effect, studying a single iteration, the "mapping" in explicit form. The literature does contain some sporadic work on this problem for low values of N.⁽⁾ On the other hand, the second point of view does not seem to have been considered except for the case of one dimension.* The general "solution" to the problem posed by (b) is the explicit construction of the vector x(n) in terms of the initial vector x(⁰ ) and the time-variable n. Such a solution can only, in general, be presented as an iterative procedure. The most one can hope to do is to predict the limiting behavior of the system as n -— oc. Luckily, for practical applications this is usually the only thing of interest. However, even this problem has hardly been touched upon.

In the limit n - oo a variety of behaviors is possible; the vectors may, for example, converge under iteration to a limiting vector x, they may oscillate between a finite set of limit vectors xi, or they may exhibit a more or less chaotic behavior, i.e., do neither, but have an ergodic behavior, i.e., the limit in time of average 1/N Z_4 Xi will exist. (The last is familiar from Kronecker-Weyl's theorem⁴ of irrational rotations; see also Ref. 5.) The type of behavior, as well as the numerical value of the limit (if it exists) may or may not depend on the initial vector x(O) All the examples (with one exception--see Section XI) studied in this report turn out to be of the first two types; i.e., they either converge to a limit vector or converge to an oscillation between a finite set of vectors in a definite order, i.e., with a definite period.

II—
Normalization

Formally, the system (1) contains N² (N+1)/2 parameters, the yke. It proves convenient to reduce this number (somewhat arbitrarily) by postulating the condition: N S 7Z = 1, all k, e (5) i=1

* To be sure, there is an extensive literature on coupled non-linear differential equations of 1st order in 2 variables.² For N variables, a special class of differential equations is treated in V. Volterra.³

This has the great practical advantage that if we now normalize the initial vector by the conditions: N Exi - S (6) i=l then if we divide the right-hand side of each equation of (1) by the constant S, we have: N x = S (7) i=l

In biological terminology (see below) this means that we restrict ourselves to the case of a "constant population." Of course, there is no loss of generality in taking S = 1. In the sequel we therefore always impose condition (5), and furthermore, take: 0 < x(°) < 1, for all i with o) = 1 (8) i=l

This, of course, reduces the number of variables to N - 1. (For an example of a different procedure when (5) is not postulated, see the last section of the Appendix.) The xi are now restricted to lie on the positive portion of the hyperplane: N xi= 1 i=l

Even with all these restrictions, it has so far proved impossible to give any general theory of the limiting behavior of the systems (1). However, one sub-class, defined by certain reasonable restrictions on the coefficient -k, has been completely studied for all N. The results for this case are summarized in the Appendix. The rest of the report is concerned with a discussion of a different sub-class which seems to be of considerable interest, but for which no general theory as yet exists.

III—
Binary Reaction Systems

1. Equations (1) have a natural interpretation in terms of biological or genetic language. Consider a large population which consists of N different "types" of individuals.* Let Xj represent the fraction of male individuals of type j, xz be the fraction of the population of type e (we assume that there are equal numbers of males and females of every given type, hence we may represent them both by the same letter x). Then the system (1), defined by the coefficients y7e, determines the composition of the next generation, which results from a random pairing (once for each individual) of the population at the present generation. If we assume that the members of the old generation do not survive into the next, then for very large populations the expected value of the fractions of the individuals of type j will be given by the system (1). By virtue of our restrictions: N kl = I, all k, (5) i=l N yxi = 1, <xi < 1 (8) i=l

the size of the population is constant. The system (1) can therefore be looked upon simply as defining a "mating rule," i.e., determining the characteristics of the offspring from the characteristics of parents. The problem, of course, is to determine the composition of the limiting population.

2. As mentioned in II, the set of equations (1) is, in general, very difficult to study. Consequently, we decided to restrict ourselves initially to a sub-class of systems defined by the following additional restrictions on the coefficients kE'

For each pair k, , there is exactly one value of i for which /k = 1. For all other values of i the coefficient is zero; i.e., Vi = 5ii, (9) where each pair (k, c) determines one value of i _l.

* These "types" are not, of course, meant to correspond to the genotypes or phenotypes of Mendelian genetics!

This means that every term in the product (xl +₂ ...+x N)² will appear in exactly one row of the set (cross-terms appearing with the factor 2). For example, X₁ = X1 + x₂+ x₄+ 2X₁X₄2 ₄XX + 2X ₃X ₄ X₂= 2x₁X₃+ 2X₂ x3 33 = 2x₁X₂x₄₃

Since there are N(N + 1)/2 terms in the square of the sum, the number of different possible systems of this sort is clearly equal the number of ways of placing N(N + 1)/2 different objects in N boxes, no box being empty. Setting P- N(N + 1)/2, this number is easily shown to be: TN = N^P ( (N-1)P + (N-2 -. .. (10) () 2 For example, T3 = 540 T4 = 818,520.

Naturally, many of these are equivalent to each other under permutation of the indices 1, 2,... N. A lower limit to the number SN of different systems, inequivalent by permutation, is: T*TN (11) N N!

In fact, the number S _N of inequivalent systems (in this sense) will be somewhat higher than this because of the fact that some of the systems are formally invariant under certain of the N! permutations. Thus, for N = 3, T* = 90, but by actual enumeration it is found that there are 97 inequivalent systems, i.e., SN = 97. For N = 4, Tk = 34,105; at the present writing the actual number SN of inequivalent systems has not been determined by us.

We have called systems restricted by condition (9) "Binary Reaction Systems. " The reason for this name is that such a system associates with each pair (Xk,x() a unique result, say x'. Symbolically: k ® (-Dj (12)

This seems to us to be a natural and simple definition for binary reactions in which "particles" of types k and e produce by "collision"

particles of type j. (The "genetic" case is inherently more complicated; see, for instance, the Appendix.)

More general and natural, though less simple, would be a rule of the form: k e f (j,m) (13) i.e., a pair of particles produces a pair of not necessarily similar particles. We have studied a few generalizations of the simple scheme (12), though not yet in comparable detail (see below, Section VIII.2).

The reaction rule may be presented in tabular form. Consider, e.g., the system: The table for this would be:

Considered as an algebraic system with a law of composition (multiplication) given by the table, this scheme is seen to be commutative

(xixj = xjxi) but non-associative, e.g., (XlX₂)x ₃= X3X3 = X2 X ₁(X ₂X ₃) = XlXl= X1

Binary reaction systems, as defined above, are always commutative (since each product occurs in only one of the set of equations) but are not in general associative. Indeed, for N = 3 there are just five associative schemes of this sort:

This last table corresponds to the finite group in 3 variables. Although this classification is suggestive, it does not appear that these associative systems are distinguished from the non-associative ones as regards their convergence properties. This is at least the case for N = 3, and the presumption is that no particular significance will attach to the associative property for higher N either. However, if the "reaction table" possesses the properties of a group table, certain special properties are easy to establish; e.g.:

The fixed point of the transformation has coordinates xl = x2 ... XN = 1/N and is attractive, i.e., the iterates of any vector in its neighborhood converge to it.

IV—
Procedure and Results

At this point it is necessary to describe our experimental procedure and results in some detail, since they will be referred to frequently in what follows.

As stated above, for the case of 3 variables, X₁,x₂,x₃, there are 97 binary reaction systems inequivalent to each other under permutation of the labels 1, 2, 3. Each one of these has been studied numerically (on an IBM 704) by having the machine select randomly three initial vectors with coordinates x(), ) x (satisfying x() = 1, 0 < x <1) 1 ,22 ,23 _i _ 2__ stsyn and letting the computer iterate the transformation in question "as

long as necessary," i.e., until some definite limiting behavior was observed.* In all but two cases such limiting behavior became evident without further analysis (one "ambiguous" case is discussed in Section IX.1; the other is mentioned in Section XI). All other transformations eventually either:

(a) Reached a stable distribution or a "fixed point" of the transformation, or (b) Oscillated between two or three fixed sets of values.

Only in (relatively) few cases did the behavior depend on the choice of the initial vector. However, the rate of convergence (used in the generalized sense to refer to both fixed points and "fixed" oscillations) often varied considerably with this choice.

The main results are contained in Table II. Here each system is written down in symbolic form, the results of iteration being given below. Each system has attached to it a conventional symbol, e.g., I.5.p, II.l.d, etc. The Roman numerals I, II, III refer to a distribution of quadratic terms on the right-hand side corresponding respectively to the three partitions of 6 into exactly 3 parts, viz.: (3,2,1), (4,1,1), (2,2,2). The other symbols refer to distributions within these main divisions, and are purely conventional. (They correspond to a particular order in which we have examined these cases on computing machines.)

A few examples will serve to illustrate the notation. (a) Consider the system: x1 = 2x _lx2 + 2x _lx ₃+ 2x ₂x ₃ x₂ = x2 + x₃ Conventional name: I.6.b , =_ 2 X3 ² In Table II this appears symbolically as: 2(12) + 2(13) + 2(23) b (11) + (33) b.d.p. ) not degenerate (22) i.f.p. given by(x = x₂ ): x= .56311573 2x⁴ + 2x³ _ ² - 3x + I = 0 x₂ = .32878482 i.f.p. m -* i.f.p. X3 = .10809945

* Although the sum xl+x ₂+x ₃= 1 is formally conserved, it is necessary to normalize at each step to avoid loss of accuracy by round-off errors in the last digit.

The notation b.d.p. means that if any of the initial x's = 1, the limiting configuration will be an oscillation between the two states: (x₁ = 0, x2 = 1, x₃ = 0) and (xl= 0, x2 = 0, x₃ = 1). This we call a "boundary double point" (b.d.p.). Since it is evident from the structure of the system which variables will assume these values, it is not necessary to specify the b.d.p. more completely. In some systems the b.d.p. will only be reached if either of some two rather than any of all 3 variables is initially equal to 1. These cases are always immediately obvious from the structure of the system.

The words "not degenerate" mean that if initially we have some x0) = 0, it will not automatically remain so for all time, i.e., that xi is not a factor of the r.h.s. of the ith equation. The notation i.f.p. (x = x₂) followed by an equation means that there exists an interior fixed point (i.f.p.), that is, a fixed point with no xi = 0, and that the value of one of the variables (in this case x₂ ) is given by the relevant root of the equation. This equation is simply gotten by suppressing the primes on the l.h.s. and eliminating two of the variables. By relevant root we mean a real root between 0 and I which satisfies the set (sometimes extraneous roots are introduced in the elimination process; these do not satisfy the original set). To the right of this is given the resulting fixed point obtained from this equation. The notation m -+ i.f.p. means that the transformations as carried out on the machine actually converged to this value (to 8 decimal places) for three random initial vectors.

(b) As a second example, consider: 1 2 + X₁=- X² + 2xXl₂+2x₁x ₃ x2 = x² + x² Conventional name: I.4.a X.= 2x₂x₃ In Table II this appears as: (11) + 2(12) + 2(13) 2 n.f.p.'s (22) + (33) b.f.p. 2(23) doubly degenerate no i.f.p. m -* n.f.p. (xl = 1)

n.f.p. means "nodal fixed point" and refers to the fact that xi = 1, for at least one value of i, is a fixed point. b.f.p. or "boundary fixed point" means that there exists a fixed point for which one x = 0. In view of the explanation of "not degenerate," the term "doubly degenerate" is self-explanatory. "No i.f.p." means that there is no interior fixed point, and m - n.f.p. (xl = 1) means that the system converged to xl = 1 for 3 randomly chosen initial vectors.

(c) Consider finally: = + 2 X1 = X1 + X3 X₂ = 2X1X3+ 2X2X3 Conventional name: III.1.f X₃ = 2X1X2+ X²

In Table II (leaving out some explanatory material with regard to the b.d.p.): ( ( ) n.f.p. b.d.p. ( not degenerate i.f.p. given by (xX: . - - x) .i.f.p. m -* i.d.p. X. xl .xl X.X .i.d.p. X.X.

In this case the transformation did not converge to the i.f.p., but rather ended up oscillating between two sets of values, that is, it achieved an "interior double point" (i.d.p.).

We hope that with these examples in mind, it will be possible to interpret Table II. In a few instances, some remarks are appended, but these are self-explanatory.

V—
Convergence Behavior

To date we know of no criterion which enables one to predict combinatorially, i.e., from an inspection of the reaction table, whether or not the limiting behavior of a given system will be true convergence (fixed point attained), oscillation between a finite set of limit vectors (periodic point), or neither. For general N, such a criterion will certainly not be simple; it would, for instance, have to take into account the various boundary fixed points and boundary periodic points which are, of necessity, present in any binary reaction system. Of course, every binary reaction system has either boundary fixed points or boundary periodic points (or both). (from Brouwer's fixed point theorem it follows that at least one fixed point must exist.) In other words, the behavior of boundary points under iteration will always have to be treated specially.

An example of this type of complication is provided by system I.1.c. Transcribed from Table II, this reads: X1 = X1+X2 + X3 X2= 2x₁x ₃+ 2x₂ 3 3₃ = 2X₁ X₂

First of all, it is clear that x(0) = 1 will lead to the n.f.p. x₁ = 1. Furthermore, it is clear that a b.d.p. exists, namely: x = 1/2 x₁ = 1/2 x2 = 0 x2 = 1/2 3 = 1/2 x3 = 0

Experimentally, for three different randomly chosen interior vectors, the systems converged to the i.f.p. given in the table.

Other special cases have a behavior more difficult to discover; e.g., the system III.l.f, in addition to the i.f.p., has the periodic solution: l = .31944846 x = .56519772 X2 = 0 x2 = .43480228 x3 = .68055154 X3 = 0 which is attained if x?f⁾ or 3x0) = 0.

For our randomly selected initial vectors, however, the system attains the i.d.p. given in the table.

One might think that things get simpler if we consider only interior initial points. Several examples, however, show that even here nothing universally true can be asserted. For example, the system I.2.m attained its i.f.p. with two different initial vectors, but went to the b.d.p. from a third initial vector.

The situation clearly gets more complicated in higher dimensions, where the classification of special boundary solutions in general depends on a complete knowledge of the behavior of lower-dimensional systems.

One may consider, in order to determine whether the fixed point of the transformation is "attractive" or "repellent," the value of the Jacobian of the transformation or the fixed point. If, for example, the absolute value of the Jacobian is > 1, then, in general, iterates of points in every neighborhood of the fixed point will diverge from it.

A summary of the convergence behavior in all 97 three-variable systems for random initial vectors is given in Table I. In 23 systems there is convergence to xi = 1, for some i. Twelve converged to a b.f.p., 15 to a b.d.p., 4 to an i.d.p., 4 to a b.t.p. (boundary triple point), one to an interior triple point, while 6 showed varying behavior depending on the initial conditions. (This class may turn out to be larger with more initial points sampled.) One does not converge at all. All the rest converged to an i.f.p.

The two systems II.1.d and II.1.f showed a continuum of i.f.p.'s and i.d.p.'s, respectively, but this behavior is easy to understand, and not specially significant (see remarks in Table II).

Except in 2 cases, systems I.2.j and III.2.3.a, convergence was numerically evident (though occasionally extremely slow). 1.2.j is particularly interesting, and is discussed in detail in Section IX.1. III.2.3.a is briefly discussed in Section XI.

VI—
The Nature of the Interior Fixed Points

1. Leaving aside for the moment the question of convergence, it is of interest to inquire into the nature of the various i.f.p. Since for a given N, there are a finite number of different systems, there are only a finite number of i.f.p. These are, of course, defined by the set of algebraic equations obtained on suppressing the primes on the left-hand side of the systems in question. As mentioned above, from Brouwer's theorem it follows that there exists at least one fixed point,

but it need not lie in the interior. Frequently these systems have no solutions such that 0 < xi< 1, all xi, which means that no i.f.p. exists. Consider, for example, system I.1.a. The set of equations defining the fixed point is: 2 2 2 X1 = C² + X2 +-X² X2 = 2x₁X₂+2x₁x₃x3 = 2X₂ X₃

Since, by definition of an i.f.p., x₃4 0, we must have x₂= 1/2; the second equation then implies: 1/2 = 2x(1 - xl), i.e., xl = 1/2, so that xl+ x₂ = 1, implying x₃= 0. Thus no i.f.p. exists.

In general, to find the i.f.p. we must eliminate two of the variables. The resulting equation is then of ₄th order in the remaining variable, say xi, although it may have factors corresponding to x 1 = 1, xi = 0,

or perhaps to extraneous roots like xi = -1. (In Table II the equation listed is always in "reduced" form, with these factors removed.) Occasionally the equation may have two real roots in the interval 0 to 1. For N = 3, in all such cases one of the roots proved to be spurious, i.e., it did not satisfy the original system. In fact, excepting the case II.l.d, mentioned above, which had a continuum of i.f.p., no system had more than one i.f.p. Although it is doubtless possible to give a complete theory of these equations for N = 3, a similar treatment for general N seems beyond reach. Here the elimination process can yield an (unreduced) equation of order ₂N-1.

2. Bounds for the i.f.p. Consider an i.f.p. satisfying: I >X1>x 2 >X ₃>...> > 0 (16)

Clearly, we lose no generality by specifying this ordering, since we can always carry out a permutation on the system so that (16) holds. For a given N, the "largest" i.f.p. will be defined as that i.f.p. for which xl< 1 has the largest numerical value as we range over all possible systems. (Since the number of systems is finite for finite N, there will always exist a largest i.f.p.) The question then arises: Given N, which system has the largest i.f.p., and what is the corresponding value of

xl? In view of the astronomical number of inequivalent systems (for even moderate values of N) it is of some interest that a partial answer can be given to this question.

For N = 3, our complete study reveals that the system possessing the largest i.f.p.-hereafter called the "maximal system"-is II.3.d, for which the defining equations are (we interchange x₂ and x₃ for convenience): X1 = X3 + 2x₁x₂+ 2x₁x₃+ 2x₂x₃X2 = X (17) X3 = X² The (unreduced) equation is clearly (x = xi): x+x² +x⁴ = 1 (18) which yields: x =.56984029 (19)

The natural generalization of this system to N dimensions is: N\2 N-1 XI=XN + ...= Xi - X_ii=l i=l X2 = X12 X3 = X22 (20) 2 XN-XN-1 The root x = xl is then given as a root of the equation: p=N-1 fN(x)-- E x² =1 (21) N=O This root converges very rapidly as N - oo; for example, N = 4: x = .566160865... N = 5: x = .566123797... (22) N = o: x = .566123792...

It is tempting to consider this last number as an N-independent bound for all binary reaction systems. Unfortunately, this is false, as will be shown below.

Consider a system with xl > 1/2 and satisfying the ordering (16). Let us assume the "skeleton": x1 = x1+ 2x₁x₂+ ... X2 = 2x _lx ₃ + ... X ₃= 2X ₁X ₄+ ... (23)XN_₁= 2X_N+ ...XN

Clearly: 1-Xl = X2+X3+X4+ ...+XN<X2 [1+₂+₄₂+-8X₃(₂2 (24)l 2xL 2 (2x,1) 1- xi But X2< 2 (25) i - Xl 1 - xl or >2>p ² (₂x)P (26) If we equate these bounds and set y = 2x₁ , we obtain the equation: yN-1 - ₂yN-2 + I = 0 (27) Calling the root of this equation yN, it is evident from (26) that: Xz< Y (28)

Clearly, yN - 2 as N - oo, i.e., the bound is N-dependent. One might suspect at first that this bound is a very weak one, and that the actual maximal system has a much lower i.f.p. However, that the bound is the best possible is proven by exhibiting a system for which yN/2 = x1 is actually obtained. In fact, such a system is:

X= X1+ 2x1x2 X2 = 2l₁X₃X3 = 2x1x4 (29) XXN_₁ = 2xlx _NN N = 1 2 x2 x j=(x .....X)² =(1X)²xNx = E i + E2 E^ = Xj=(X2+• *.1) = (-Xl) i=2 i<j=2

It is easily verified that this system has (27) as its i.f.p. equation. For N = 4 x _l= = .809016995 ( 4 (30) N = 5 x₁ = .919643378

Experimentally (N = 4, 5, 6), this i.f.p. is not attained on iteration starting from a general point, but these converge to the n.f.p. XN = 1. This is to be contrasted with the behavior of the system (20), which actually attained its i.f.p. (N = 3, 4, 5). Indeed, for system (29) it turns out that the absolute value of the Jacobian at the i.f.p. is (y = 2x₁): IJI = yN-2(² - ^y)(³¹⁾ > 1 which makes it reasonable that this i.f.p. is not attractive.

On the other hand, it is clear that for a "skeleton" of the form: Xl = ... X2 = X²+ ... X3 X2 + ... (31) x_N-XN² -_1 + -' we have:

2 x X2 >x2 X₃> X₂> xl (32) 2 >x 2>x2^N-1N N-1 - 1 N-1 whence: 1 = xl+x2 ...+x>Zx = fN(Xl) p=O

Therefore, for such a skeleton, the root of fN(x) = 1 does indeed provide an upper bound (attained for the system (20)).

At this writing it has not yet been shown that the system (29) is actually maximal. However, a weak upper bound can be obtained for all systems such that x₁ does not contain the term x2 (skeletons of the form (31) are a sub-class of these).

Namely, in this case xl<2xl (1 -xl ) + (1 -x)² = 1--2 (33) Therefore, clearly x2 + Xl 1 or xI < .61803399 (34)

However, we can do much better for this case. In fact, we can show that for xl > 1/2, we must have, under the ordering (16): Xk > z-21 (35)

which then establishes the bound (21) by the previous argument.

VII—
Periodic Limits

For a large number of 3-dimensional systems, randomly chosen initial vectors iterated to a periodic limit, i.e., the limiting behavior was an oscillation of period 2 or 3 between fixed points. Twentyfour systems exhibited this behavior for three initial vectors, while six

others achieved a similar limiting configuration for at least one choice of initial vector (with no coordinates lying on the boundary). In most of these cases the final state was of the form:

i.e., a boundary double point. A few cases of a boundary triple point were also observed (cf. Table I). Such final states we call "trivial," for the reason that the algebraic structure of the transformation alone indicates that such a final state is at least possible. We may contrast this "trivial" type of oscillatory final state with those for which the oscillation takes place between two or three interior points. The latter we call an "interior double point" (i.d.p.) or "interior triple point" (i.t.p.).

For N = 3 we found just four examples of an i.d.p. and one of an i.t.p. There was also one case of a "non-trivial" b.d.p. (system I.2.e) for which the final state was oscillatory with period 2, but between two "non-trivial" boundary points, viz.:

The existence of an interior double (or triple) point means that the second (or third) power of the transformation possesses these limit values as fixed points. The algebraic difficulty of finding such points is in general prohibitive. For example, in the unique case of the interior triple point (system 1.3.g), if we let:

then one coordinate of the triple point is determined by the set of equations:

It may be verified that the successive values of xl given in Table II indeed satisfy this set of equations.

Although there are no oscillatory limiting configurations with periods greater than 3 for N = 3, one can, of course, find such by going to higher N. Indeed, we discovered, by chance, a particularly interesting case, viz.:

This can be considered as one particular generalization of the 3variable system I.5.j. This generalization-which we applied to several of our original systems (Table III)-consists in setting x4 = x3, replacing x3 in the original system by x4, and putting the new crossterms 2x ₄(il + x ₂+ X ₃) in the top line. When we generalized the triple periodic case, I.3.g, in this manner, the resulting limiting behavior (for 3 randomly chosen initial vectors) was still periodic with period 3, but the configuration was of the "trivial" sort, i.e., the b.t.p. (1,0, 0, 0), (0,0, 1,0), (0, 0,0, 1). However, in the case of I.5.j, whose resulting generalization is given above, the limiting configuration was oscillatory with period 12. (The values of the coordinates are given in Table III.) A further generalization to 5 variables, following the same prescription, yields the system:

In this case, 3 random initial vectors achieved an oscillatory limiting configuration of period 6. (See Table III for the numerical values.) In our opinion, it is not likely that the behavior observed in these two cases could be predicted by means of any simple criteria.

VIII—
Form Stability

1. In a few cases we investigated the effect of making slight changes in the form of the equations themselves. One way of doing this, which makes the change of form depend on a single parameter, is as follows: Multiply each term by the factor I - e, and add e/2 times the term in each of the two other rows.

For example, the system 1.5.0 is:

This was now modified to:

In this case it is easy to carry out the elimination of x₂ , x₃ to obtain the equation for the i.f.p. as a function of e. There results: blx⁴+ b2^x3 + b₃x² + b₄x_I+ b ₅ = 0 (38) where the bi are given in terms of the parameter: 2- 3e a = 2 (39) as follows: b = a⁴ b2 = 4a³ b3 2a² - ³ (2 + a)3 (40) 4a ²b4 =a(1 + 2a) - 3 (2 + a) o~ o ~a2 b₅ = (2a + 1)+ (2 + a)²

For e = 0 (a = 1) we get (dividing out xl) the original i.f.p. equation: x+ + 4x² - 1 = 0 (41) At e = 2/3, a = 0, all the coefficients bi vanish. This corresponds to the set: 12X₁= x ₂= x3 = (X2 + X2 + ) + 3 (X1l2 + x1X3 + 3X23) (42)

which reaches the i.f.p. (1/3, 1,3, 1/3) on a single iteration starting from any initial vector. In fact, as e -- 2/3, all systems to which this generalization is applied tend toward this simple case (since at e = 2/3, 1- e = e/2).

For the system 1.5.0, generalized in this manner, we investigated the convergence for several values of e, viz.: e = .001, .01, .05, .1. In each case, randomly chosen initial vectors iterated to the i.f.p. predicted by equation (38). In other words, there is a sort of continuity in the convergence behavior as the form of the equation is changed in this simple manner.

Slightly more interesting is the result of applying this one-parameter generalization to the system 1.3.g. As mentioned above, random initial vectors iterated according to this transformation reached an oscillatory final state with period 3. Vectors under the generalized transformation behaved in the same fashion, but as e was increased the final-state oscillations decreased in amplitude, until at e _ .045 general initial vectors appeared to converge to the i.f.p. predicted by the corresponding i.f.p. equation. Presumably, the final state is still oscillatory (with period 3), but the oscillations are too small to observe with 8-decimalplace accuracy.*

The conclusion to be drawn from these experiments is that binary reaction systems are "stable" under small perturbations of formal structure.

*For e > 0, some initial points converged to the i.f.p. In other words, the i.f.p. (which is a function of e) appears to be attractive for e > 0. More detailed analysis would probably show that what is happening is that the area of the region of the triangle for which convergence (of a point in this region) to the i.f.p., is increasing with increasing e. Correspondingly, one is less likely to pick an initial point outside this region, i.e., a point which will iterate to the i.t.p. For further discussion of behavior in the large, see Section XI.

2. When more radical changes of form are made, we do get correspondingly greater changes in behavior. For instance, we took the above-mentioned system I.3.g and kept only the skeleton: s/ 12 12 xI = 22 + X23 + 2x1x2 + ' = 2x₁x₃+ 2x2x3+ . .. (43) 1 2x₃ X₁ -t +.

To this we added the missing terms (1/2) x2, (1/2) x² , (1/2) x3 in all possible ways (i.e., 27 ways, of which one corresponds to the original system). For 25 of the resulting 26 new transformations, random initial vectors iterated to an i.f.p. (The results are summarized in Table IV.) One system, however, gave an oscillatory final state with period 3, a behavior analogous to that of the original system, viz.:

The configuration of the final state (starting from 3 randomly chosen initial vectors) was:

It is clear that a "rule" which allows x₃@z3 - xl or x2 (with equal probability) is a rather unnatural one. A somewhat more logical modification is to assume a skeleton in which the cross-terms appear with coefficient unity, and to add the missing terms x1x2,xx13,x2x3, in all possible ways. When this was done for I.3.g, all the resulting transformations had an i.f.p., and in every case random initial vectors iterated to the i.f.p. (see Table V).

This change of rule has a natural interpretation; effectively it allows non-commutativity, since e.g., xixj may give two results, which

we can interpret as the respective results of xixj and xjxi. Formally, this is an interesting generalization of the concept of a binary reaction system (as originally defined), but the convergence behavior does not appear to be startlingly different.

This same "non-commutative" generalization was tried on the system III.2.1.b; the results offered no surprises.

IX—
A Specific Convergence Problem

1—
The Exceptional Case I.2.j

As mentioned above, there are only two systems among all the 97 which exhibit an ambiguous convergence behavior. The most interesting of these has the form: x' = xl+ x3 + 2x₁x₃ x'₂ = x + 2xlx₂ I.2.j xI3 = 2J:2:J:3

The system is degenerate, and lias two n.f.p.'s; in the following we shall ignore the behavior of tile boundary points (since these present no problem) and only discuss the behavior under iteration of interior points.

By inspection, it is evident that the system possesses the i.f.p: 1Xl = 3 = 1 (45) X2 = 2

For three randomly chosen (interior) initial points, little convergence was evident even after some 85,()()( iterations. In order to see better what was happening, wc intro(duce(ld iew coordinates: 1 + .1 - : ₃ 2 (46) ( = -r2

These are effectively Cartesian coordinates in the plane of the triangle formed by the three vertices (1,0,0) (0, 1,0) (0,0, 1).

More specifically, in the original form of the transformation, the xi's are constrained to move on the positive portion of plane x₁ + x₂ + X₃ = 1:

For algebraic convenience, we distort the triangle into a 45^° triangle with base unity. The coordinates of a point in this triangle are then (S, a), as shown in the sketch below. Here

In terms of these new coordinates, the transformation takes the form: S' = 1 - 4a + 4²+ 2aS (47) a' = 2aS and the i.f.p. is: 1 2 (48) 1 4

Note that the Jacobian of (47) is exactly 1 at the f.p. (48). For future reference, we write down the inverse of (47): - ' - (49) 1-S'-a' a= -' 2

If we now make the further transformation: x =S-1 f~~~~~~~~2 ~~(50) y =c-awe get: x x' -y + + 4² + 2xy (51) x y =y + + 2xy with i.f.p. x = y = 0.

Note that if we consider only linear terms: x 2' =-y + (52) , x y =y+ ~~2~~~~~ ~215

then an invariant ellipse exists: x12 + x^ty+ + 2y^t2 = + ^xy + ^2y2 (53)

Figure 1 shows a plot of the observed iterates in the x, y plane. The three curves show the behavior of successive iterates, the initial point of the sequence being taken respectively at n = 1550, n = 10,101, n = 75, 001. Each curve is roughly an ellipse of the form: +xy + 2y² = C (54)

Numerically, at least, C appears to -— 0 as n -—oo, or, to put it another way, the axes of the ellipse of reference are shrinking.

In order to convince ourselves that this apparent convergence was not simply the result of systematic round-off, we used the inverse transformation (48), numerically retracing our steps from n = 10, 101 to n = 1,550. All coordinates agreed to 6 decimal figures, which pretty well precludes the possibility that the observed convergence is a numerical accident.

It should be noted, by the way, that no matter how close one is to the fixed point, the quadratic terms in (51) cannot be ignored, for according to (53), the linear terms by themselves will generate a sequence of iterates which will all lie on the ellipse of reference.

For the purpose of discussion, it is convenient to transform equation (51) so that the reference curve is a circle. The linear part of the transformation can be written: ()A(x) A- (1 ) (55) y (55) 2 If we transform this with the matrix: S=( ) (56) 0 2 we find: R = S-1 AS = 4 (57) 216V-- 3 4 4

which corresponds to a rotation through an angle 0 = cos^- . The transform through S of the complete mapping (50) has, then, the form: , 3 V7 10 2 =-x + -y+ y -6xy ( 4^~~~ 4 ~~~~(58) -v 3 4 4 y' = -T-x + y+2y² + 2Vfxy

2—
The Asymptotic Behavior of the Angle of the Radius Vector under Iteration

The transformations with which we are concerned do not preserve the ordinary measure (Lebesgue measure) of the space which is mapped into itself. Moreover, they are not one-one. In some cases they shrink a neighborhood of the fixed point into a proper part of itself and the limiting image of such a region may consist of this point alone. Obviously, an invariant measure, if it is to be constructed, would have to be of Lebesgue-Stieltjes type and assign positive values, in some cases, to sets consisting of single points.

One may be interested in the behavior of the angles (with a fixed direction issuing from a fixed point) of the vectors Ti(x), where x = (x₁ , x₂ x₃ ); i = 1, 2,... more generally in the behavior of the points: T(Si-) (59) ITi(()I

on the unit sphere. Something can be said about it even in cases where T does not transform any bounded region into itself. Thus, for example, if T is an arbitrary linear transformation of the n-dimensional Euclidean space into itself, T(0) = 0. Then the ergodic limit of the average of Si exists. In other words, if C is an arbitrary "cone" of directions in space, i.e., a sub-set of the unit sphere fc(s), its characteristic function being:

fc(s) = 0 if .s C; fc(s) = 1 if s c C then for almost all s, the limit: N lim N Zfc(T'(s)) (60) i=1

exists. This follows for a general linear transformation T, from the well-known theorems giving, in fact, more precise information-in the two extreme cases: if T is an orthogonal transformation we have the Kronecker-Weyl theorem on equipartition; if T, considered as a matrix, has all coefficients positive, the Perron-Frobenius theorem asserts convergence to a unique direction. In the general case one obtains, by considering a decomposition of the space into sub-regions where one or the other behavior dominates, at least the existence of the ergodic limit.

Presumably, the theorem still holds true if T is a general homogeneous quadratic transformation of the n-dimensional space.

In our very special quadratic transformations of the plane, more can be said: The transformation of the previous section possesses the Knonecker-Weyl property: The angle described by the iterates of almost every point covers the circumference of the unit circle densely and uniformly.

We hope to show that if the linear part of a quadratic transformation Q, which has the origin as its fixed point, consists of a rotation through an irrational angle, then the iterates of Q converge to the origin provided one starts with points in a circle of sufficiently small radius.

A detailed discussion of these matters will be given in a subsequent report.

X—
Further Generalizations

In view of the impracticability of studying all possible Binary Reaction Systems for any N > 3, we thought it worthwhile to generalize a few of our 3-variable systems to higher dimension by arbitrary but fixed rules. One such generalization is mentioned in SectionVII (see I.5.j - ext - I and subsequent discussion). Another essentially different way to generate "interesting" systems is to construct for any given 3-variable transformation the corresponding "super-system." This is constructively defined as follows: We introduce nine variables, yl, y ,... yg, according to the prescription:

and substitute these in the original transformation. We then have three transformations for the three triads of variables yi+ Y2 + Y3, y4 + y5 + Y6, Y7 + Y8 + y9- Consider, e.g., the system I.2.e. The last line of the transformation now reads: y' + y8 + y' = 2 (yl + Y2 + Y3) (Y4 + Y5 +Y6)

In order to convert this into 3 separate expressions for y', y8, y', we could, for example, formally identify the variables modulo 3, i.e., YI ~ Y4 ~Y7 Y2 ~ Y5 ~ Y8 (62) Y3Y~ 6 ~ Y9

Correspondingly, on the right-hand side, we could make the identification: ² y1l4 - y² y2y5 ~ Y2 2y2y6 +2y3Y5 - 2y2y3 (63) 2yy16 + 2y3y4 ~2ylY3 2yly5 +2y2Y4~ ² Y1Y2

We can now write expressions for y',y8, y' so that, with these formal identifications, the resulting sub-system will have the same form as the original 3-variable system, i.e., Y7~ y = yl +Y2 +² Y2Y3 - 2Y14 + ² Y2y5 + 2Y2Y6 + 2y3Y5 Ys ~ Y2 = Y3 + 2yy 33 - 2Y3Y6 + 2y1Y6 +2y3y4 (64) ys y13 = 2y1y2~ 2yly5 + 2y2Y4

In this way we (arbitrarily) obtain equations for y', y8, y' in terms of the yi's. When this is done for each triad, a 9-dimensional B.R.S. results.

In the present case, using the symbolic notation of Table II: (11) + (44) + 2(47) + (22) + (55) + 2(58) + 2(23) + 2(56) + 2(59) + 2(68) (33) + (66) + 2(69) + 2(13) + 2(46) + 2(49) + 2(67) 2(12) + 2(45) + 2(48) + 2(57) (77) + 2(17) + (88) + 2(28) + 2(29) + 2(38) + 2(89) (99) + 2(39) + 2(79) + 2(19) + 2(37) (I.2.e - Super) 2(78) + 2(18) + 2(27) 2(14) + 2(25) + 2(26) + 2(35) 2(36) + 2(16) + 2(34) 2(15) + 2(24)

As a second example, we quote the result of treating system I.3.g in the same manner: (55) + 2(25) + (88) + (66) + (99) + 2(36) + 2(45) + 2(78) + 2(15) + 2(24) 2(46) + 2(79) + 2(16) + 2(34) + 2(56) + 2(89) + 2(26) + 2(35) (44) + (77) + 2(14) 2(28) + 2(58) + 2(39) + 2(69) + 2(18) + 2(48) + 2(27) + 2(57) 2(19) + 2(37) + 2(49) + 2(67) + 2(29) + 2(38) + 2(59) + 2(68) 2(17) + 2(47) (I.3.g - Super) (22) + (33) + 2(12) 2(13) + 2(23) (11)

When these transformations were iterated for randomly chosen initial vectors, the sums: Xl = Yl + Y2 + Y3, x₂= Y4 + ₅ + Y6, 33 = Y7 + Y8 + Y₉

reached, of course, as they should, the same limiting configuration as was observed in the original 3-variable system. However, the actual values of the individual yi varied with the initial configuration. The results are given in Table VI.

XI—
Properties in the Large

1. Although it is, in general, possible to discuss the behavior of points under iteration in the neighborhood of a fixed point, the iteration behavior of such points over the whole domain (positive portion of the hyperplane) can, at present, be treated only experimentally. In what follows, we shall take the variables to be S, a, and the domain to be the corresponding 45^° triangle with unit base (see Section IX.1 for a definition of this coordinate system).

As stated above, we have not found any general criteria for determining which of several possible limiting behaviors will be realized for a given 3-variable system, starting with a general point in the triangle. On the basis of a rather small sample (~ 3 random initial points for each system), it appears that, excluding boundary points, the limiting behavior is independent of the initial point for the large majority of systems. As shown in Table I, however, there are (at least) 6 systems in which this limiting behavior depends on the initial point (we exclude from consideration the "pathological" systems II.1.d and II.l.f; see the discussions under these entries in Table II). Two of these, I.2.m and I.5.h, have been examined in greater detail. What was done was to look for boundaries which separate regions of different limiting behavior. This was accomplished by programming the computing machine to "search" the whole triangle in a systematic manner. On the first pass a crude net was used (AS = Aa = .05). Then, when the boundaries had been approximately located, a more refined interval was employed in the appropriate neighborhoods. For each trial point in the triangle, sufficient iterations had to be performed to identify the limiting behavior. Despite the apparent magnitude of the task (several hundred trial points had to be followed for some 70 iterations each), a complete search (first crude, then appropriately refined) takes only about 15 minutes of computing time per system.* The results for the two

* On the average, the machine will perform 50 iterations/sec. for a 3-variable Binary Reaction System. The systematic search of the triangle is somewhat slower for various reasons connected with input-output requirements.

systems studied are shown in Figures 2 and 3. In these, each calculated boundary point is determined to within an absolute error < .0025 i.e., to within 1/4% of the length of the base of the triangle. In these figures, initial points lying in the region marked "OSC." will iterate to the appropriate boundary oscillation, while all points lying outside these regions will converge under iteration to the i.f.p. The boundaries appear to be complicated. We have not attempted to study what happens to points actually lying on the boundary curves; for this an analytical treatment is necessary.

In a few simple cases it is possible to give an analytical treatment of such boundary regions. As an example, we cite the system III.2.2.a. In the S,a coordinates (we use these to conform with our treatment of the other systems; the argument can be carried out in the original coordinates with equal ease) the transformation takes the form: S' = 2S - S2 + 3,² - 4aS (65) a' = 2a(1- S)

There are 3 n.f.p.'s, namely, (S = 0, a = 0), (S = 1, a = 0), and (S = 1/2, a = 1/2), while the i.f.p. is (S = 1/2, a = 1/6). The boundaries a = 0 and S = a are clearly transformed into themselves. In addition, there exists an invariant line: S = 3a .

All points lying on this line can easily be shown to iterate to the i.f.p.with the exception of S = a = 0, which is a non-attractive f.p. This line is shown on the diagram below. The curve S =S* is the locus of all points such that S'= S; its equation is: S*= - 4a+ (1 4a)² 12 (66) 2

It is easily shown that a point lying in the region below the line S = 3a, i.e., such that S > 3a, remains in this region under iteration; similarly for points lying above the line. Further, all points lying to the left of the curve S = S* remain to the left under iteration. It can be further shown that points lying to the right of S = S* and not situated at the corners of the triangle or on the line S = 3a will eventually cross

the curve. Furthermore, in the neighborhood of the i.f.p., the linear approximation to the transformation is (x S - 1/2, y a - 1/6): x x (67) Y =- +Y

from which it can be deduced that the i.f.p. is attractive only along the line x = 3y. Finally, using the fact that for S < 1/2, we always have a' - a > 0 and that correspondingly for S > 1/2, a' - a < 0, it can be shown that all points lying above the line S = 3a will iterate to S = a = 1/2, while all points below the line will iterate to S = 1, a = 0.

The reason why this system can be treated so simply is, of course, that the boundary curves are explicitly known. When this is not the case, it may be helpful to have before one a picture of the mapping (single iteration) of the entire triangle, as well as the curves along which a and S are stationary. A few interesting examples are given in Figures 4, 5, and 6.

2. There is one case which does not appear to converge either to a fixed point or to a finite oscillation. This is the system III.2.3.a: x₁ = X²+ 2x1x ₂ x₂ = 2 + 2x₂x₃ (68) x -+ - 2x ₁x ₃

or, in terms of the coordinates, S, a: S' = 2a+S² - 3a² (69) a' = 2(1 - S)

The situation is illustrated in Figure 7. The i.f.p. is non-attractive and the 3 corners of the triangle are attractive only along the boundaries in a clockwise direction. Under iteration, points will spiral out, approaching arbitrarily close to the boundaries, but, e.g., as the transformation a' = 2a(l - S) shows (for the bottom boundary), no point inside the triangle can ever reach the boundary. Thus a general point will continue to spiral indefinitely.

Numerically, a spurious convergence was observed owing to the fact that the several random initial points chosen rapidly iterated to within a distance less than 10^-8 from one or another boundary line.

If one transforms the triangle into the unit circle in an appropriate manner, the situation can be viewed as follows: The center is a nonattractive fixed point, and all points lying in the circle spiral outwards towards the circumference. On the circumference itself, there are 3 fixed points located, say at 0 = 0, 27r/3, 47r/3, which in turn define 3 arcs. Any point lying on one of these arcs will move under iteration in a clockwise direction, ultimately converging to the fixed point which constitutes the right-hand boundary of the arc in question. Interior points, however, can never reach the boundary. In general, the sequence of iterates of any interior point (excluding the center) does not converge.

3. We have not studied in detail the rate of approach to the limiting configuration. In the neighborhood of a fixed point, this rate is usually easy to obtain, but for oscillating configurations the algebra is more difficult. For most of the cases studied, convergence (to 8 decimal places) was either attained within 100 iterations (some were much

faster) or else not for many thousands of iterations. We shall not discuss this further except to remark that the path of approach to a fixed point may depend very critically on the initial point. As an example, we refer the reader to Figure 8. For this transformation, two of the corners of the triangle are non-attractive fixed points; initial points in their neighborhoods iterate smoothly to the i.f.p. In contrast, points in the neighborhood of the origin iterate to the i.f.p. in an oscillatory manner; the existence of a limiting line through the i.f.p. is clearly evident.

XII—
Connection with Ordinary Differential Equations

No doubt it will have occurred to the reader that Binary Reaction Systems and their extensions have an obvious connection with systems of ordinary differential equations. Consider, for instance, the system of differential equations:

where the fi (xl ... XN) are homogeneous quadratic functions of the variables. A straightforward finite difference approximation to the system is: xn+ ₁ (1 -t)x(⁾ + At ft (71) xN⁺ l) - (1 - At)X⁾ + At f ⁾

If we now restrict the fi to be disjoint partial sums of the terms in (Xl+...+ XN) ² , such that ilfi = (xi +... +N)² , i.e., just the terms that occur in our binary reaction transformation, then the above set of difference equations goes over into a Binary Reaction System for At = 1. With this restriction on the fi, the system of difference equations has the property xn = 1 for all n if Ex() = 1, independent of the time-step At. It should be observed that for At > 0, the system has the same fixed points as the corresponding Binary Reaction System.

As a very simple example, consider the system: 2 2 1 with i.f.p. xl = 2= -, 2 (72) X₂ = 2xlX₂ and b.f.p. xI = 1, x₂ = 0.

Any interior point (xif 1, 0) will iterate to the i.f.p. The corresponding differential equation is (eliminating x₂ ): dx=2x - 3xl + 1 (73) dt 1 of which the solution is: 1 /2-1-x 1/2 -Jx (0 l(74) 1/2-^ 1/2-2' x

Clearly, as t - oc, xl -, 1/2, x(°⁾# 1, 1/2. Thus, in this case, the asymptotic behavior of the differential equation is the same as that of the corresponding Binary Reaction System. This, however, is not generally so. For example, the transformation I.5.d (see Table II) possesses a boundary double point which is attained for a certain set of initial points. The corresponding differential equation system, when integrated according to the finite difference scheme (71) with At = 26, converged to the n.f.p. xl = 1 regardless of the initial point.

In such cases the Binary Reaction System, viewed as a finite difference approximation to the corresponding system of differential equations, is clearly an "unstable" scheme. We hope to discuss the point further in a subsequent report.

Explanation of Graphs

1. Figure 1 is a plot of successive iterates of a point under the transformation equation (53) of the text. The outer curve shows the iterates from cycle n = 1550 (point ii) to cycle n = 1570 (point fl); the next curve goes from n = 10,101 (point i₂) to n = 10,114 (point f2), while the innermost curve goes from n = 75,001 (point i ₃) to n = 75,011 (point f3).

2. Figures 2 and 3 show the (experimentally determined) boundaries between two types of limiting behaviors for two different Binary Reaction Systems (see discussion in Section XI). In terms of the S, a coordinates, these transformations are: 1.2.m (Figure 2): 3S² - a²S' = 1 + + 3aS - 2S S² _ 3a² a -3aS + 2a 2 I.5.h (Figure 3): 1 + a² - 3S² S= l+~2 -a+S+3aS I + S² - 3,a²(a = + a -S - aS

In each case, points in the inner region (marked "CONV.") will iterate to the fixed point, while points outside this region will reach an oscillatory final state (in I.5.h, all three outer regions are oscillatory, though only one is so marked).

3. Figures 4, 5, and 6. These figures give the S, a mapping explicitly for three different Binary Reaction Systems: I.2.j (Figure 4): S' = 1 - 4a + 4a² + 2aS (cf. equation (47) of the text) a' = 2aS I.5.h (Figure 5): See 2. above.

III.2.1.a (Figure 6): S' = S + 3a + 6aS a' = S- a + 3a² - S2

In Figures 4 and 5, the numerically labelled lines are the transforms of lines of constant a; e.g., in I.5.h the curve labelled .25 is the transform of the horizontal line a = .25. The curves labelled a' = a and S'= S are lines of constant a and constant S, respectively.

Figure 6 is plotted in a 60^° , 30^° reference triangle, i.e., in terms of the variables S and t = vJa. The labelled curves are the transforms of lines of constant t, and the arrows indicate the order of the transformed points as the relevant lines of constant t are traversed from left to right (direction of increasing S). Note the change of direction for t > v3/6. (The vertical line with two arrows labelled 1/6 is the transform of a = 1/6 or t > \/3/6. This line is doubly covered, first upwards, then downwards, as S increases.)

4. Figure 7 illustrates the "non-convergent" case. See Section XI, equation (69), and the accompanying text.

5. Figure 8 illustrates different modes of convergence to a fixed point. See part 3 of Section XI for discussion.

blank

blank

blank

blank

Table I—
Summary of Convergence Behavior of Three-Variable Binary Reaction Systems

Note: In addition, we have the following I.3.g i.t.p. III.2.3.a non-convergent

Introduction to Table II

This table summarizes the properties of all 97 Binary Reaction Systems in 3 variables. The notation is explained in Section IV of the text. For a few systems, the behavior of arbitrary vectors under iteration can be predicted theoretically owing to the fact that the system reduces to a difference equation in a single variable (e.g., I.l.b, 1.2.q, etc.).* We have not though it worthwhile to note these instances explicitly; anyone using the table will immediately discover them for himself.

For the 3-variable case, the coordinates of all fixed points could be written explicitly in terms of radicals, since the f.p. equation is at most of ₄ th degree. There are many interesting relationships between the roots of the various f.p. equations. We have not investigated these relationships systematically, although it could easily be done using standard tools (cf. L. E. Dickson).⁶

Table II—
Three-Variable Binary Reaction Systems

I.l.a (11) + (22) + (33) n.f.p. 2(12) + 2(13) b.f.p. 2(23) J degenerate no i.f.p. m -* b.f.p. I.l.b (11) + (22) + (33) n.f.p. 2(12) + 2(23) 2 b.f.p.'s 2(13) J doubly degenerate no i.f.p. m - b.f.p. (x3 = 0)

* In some other cases, one may observe that one of the variables will obviously iterate to zero, in which case the limiting behavior is also evident.

Table II (cont.) I.l.c ( n.f.p. b.f.p. j not degenerate i.f.p. given by (xX - - .or (y l) x .i.f.p. yy - .m -, i.f.p. I.a ( ( n.f.p. ( ne degenerate i.f.p. given by (x x v- xx - .i v/-(cf e)m- i.f.p. x-.(cf I) I.b ( ( n.f.p. ( ne degenerate no i.f.p. m -n.f.p.

Table II (cont.) I.c ( ( n.f.p. ( nfp not degenerate no i.f.p. m -, n.f.p. I.d ( ( n.f.p. ( b.f.p. ) degenerate no i.f.p. m - b.f.p. I.e ( ( n.f.p. ( b.d.p. (see below) J not degenerate no i.f.p. one coordinate of b.d.p. (x x given by - - x .56519772, xl . x .X0, X.b.d.p. m - b.d.p. x.43480228, x J

Table II (cont.) I.f ( ( n.f.p. ( b.f.p. J degenerate no i.f.p. m -— b.f.p. g ( ( ) n.f.p.'s ( b.f.p. J doubly degenerate no i.f.p. m -b.f.p. h ( ( n.f.p.'s ( degenerate degenerate no i.f.p. m -t n.f.p. (x i ( ( nn.f.p.'s Xl ] ( x ni.f.p. degenerate m , n.f.p. (xl

Table II (cont.) I.j ( ( nf's ( degenerate i.f.p. J m — i.f.p. (See Section IX. of text for discussion) I.k ( ( p's n.f.p.'s ( ( ot degenerate i.f.p. given by (x xi) x . - i.f.p. £.m -* i.f.p. x. ( ( n.f.p.'s ( b.f.p. J doubly degenerate no i.f.p. m -* n.f.p.(x I.m ( ( x - b.d.p. • . ( b i i.f.p. degenerate (cf I.b) I ^ ^ - (cf b) x) - .(cf. III.l.d)

Table II (cont.) I.m (cont.) m -* i.f.p. with initial vectors ( x) .x) and ) .x x .m - b.d.p. with initial vector (See Section XI for discussion) x(⁾ .x ⁾ .x . I.n ( ( bd b.d.p. ( degenerate degenerate The only f.p. is the b.f.p., which is the i.f.p. of the 2-variable system System ( .| ( 3-e .(X This f.p. is easily shown not to be attractive. And, experimentally m -* b.d.p. I.o ( ( b b.d.p. ( ( not degenerate i.f.p. given by (x x x.- - X.i.f.p. m -i.f.p. x.

Table II (cont.) I.p ( ( b.f.p. (cf. I.n) ( b.d.p. J degenerate i.f.p. given by (x xi) l.x - xi.f.p. /3- (cf. h) - (cf. II.b) m -* i.f.p. with initial vectors x(⁾ x() x x() and x?⁾ and x) .x -) x) x(⁾m -* b.d.p. with initial vector x() .x() .x . I.q ( ( b.d.p. ( b not degenerate i.f.p. given by (x xl) .x- V5-i.f.p. -Ox- - .i m - b.d.p. X.

Table II (cont.) r ( ( b.f.p. (cf. I.n) ( b.d.p. J degenerate no i.f.p. m -* b.d.p. I.a ( ( n.f.p.'s degenerate degenerate ( no i.f.p. m -* n.f.p. (x I.b ( ( ) n.f.p.'s b.f.p. ( J doubly degenerate no i.f.p. m - b.f.p. I.c ( ( n.f.p.'s b.f.p. ( J not degenerate no i.f.p. m -- n.f.p. (x xl)

Table II (cont.) I.d ( ( n.f.p. n.f.p not degenerate ( no i.f.p. m - n.f.p. I.e ( ( xl ) In.f.p. den e i.f.p. ( x i.f.p. exists m -l i.f.p. I.f ( ( n.f.p. nfp not degenerate ( i.f.p. given by (x x) .- x.i.f.p. m - i.f.p. x .J I.g ( ( }) . b.d.p. .i.f.p. not degenerate ( .

Table II (cont.) I.g (cont.) i.f.p. given by (x x - - x m -- interior triple point xnx) .73924369, x() .08071790, n(n x) .01294586, x( .37280086, x) .i.t.p. - x ) .24781045, x(n) .54648124, .. I.h ( ( ) b.f.p. (cf. I.n) b.d.p. ( J degenerate i.f.p. given by (x xl) xi . m i.f.p.) mi.f.p. 1- X . I.i ( ( b.d.p. bdp not degenerate ( 3- v/. i.f.p. given by (x xl)x .x- - .i.f.p. m - b.d.p. - .(cf. f)

Table II (cont.) I.a ( ) n.f.p.'s ( ( b.f.p. doubly degenerate no i.f.p. m -- n.f.p. (xl I.b ( n.f.p.'s ( ( degen degenerate i.f.p. given by (x x xl - V - .i.f.p -Vif..2- (cf. I.m) m -, i.f.p. X . I.c ( n.f.p.'s ( ( not degenerate no i.f.p. m -- n.f.p. (xl

Table II (cont.) I.d ( b.f.p. (cf. I.n) ( ( b.d.p. ) degenerate no i.f.p. m -. b.d.p. I.e ( b.f.p. (cf. I.n) ( ( b.d.p. J degenerate i.f.p. given by (x xx - - i.f.p. m - b.d.p. (v2-(cf a) X. I.f ( b b.d.p. ( ( ( (ot degenerate i.f.p. given by (x x .- x - .i.f.p. m -* i.f.p. X.J

Table II (cont.) I.g ( n.f.p. ( ( b.f.p. J degenerate no i.f.p. m - b.f.p. I.h ( n.f.p. ( ( b.f.p. J degenerate no i.f.p. m -- b.f.p. I.i ( n n.f.p. ( ( not degenerate i.f.p. given by (x x x.- - x. m - i.f.p. x .J a ( nfps n.f.p.'s ( triply degenerate no i.f.p. m - n.f.p. (xl

Table II (cont.) I.b ( n.f.p. b.f.p. (cf. I.n) ( b p. b.d.p. ( degenerate no i.f.p. m - n.f.p. I.c ( nfp's ( deen degenerate ( no i.f.p. m -- n.f.p. (Xa I.d ( n.f.p. ( b.d.p. (see below) ( not degenerate no i.f.p. In addition to the "trivial" b.d.p. (0, 1,, (0,0, , there is a b.d.p. which is a f.p. of the transformation X X-² , Viz. Xz X 3- xV5- - b.d.p. x This is not an attractive f.p. If xO) 0, x < (- v/2, the ini- tial vector - xi If a⁾ 0, x ⁾ > (- v/)/2, the initial vector the trivial b.d.p. This behavior has been verified numerically.

Table II (cont.) I.e ( n 's n.f.p.'s ( ( doubly degenerate ( no i.f.p. m - n.f.p. (x I.f ( n.f.p. ( b.d.p. ( J not degenerate i.f.p. given by (x x X - - Z.i.f.p. - (cf. i) m - . i.f.p. . I.g ( n.f.p. b.f.p. (cf. n) ( bd.p. b.d.p. ( degenerate no i.f.p. m - b.d.p.

Table II (cont.) I.h ( b.t.p. ( btp not degenerate ( i.f.p. given by (x xi)x .i.f.p. X . m -- i.f.p. with initial vectors x) .x)and x) .XJ ( m -- b.t.p. with initial vector x⁾ .x) .(See Section XI for discussion.) . I.i ( n.f.p. b.f.p. (cf. I.n) ( b.d.p. ( degenerate no i.f.p. m -* b.d.p.

Table II (cont.) I.j ( b.d.p. ( bdp not degenerate ( i.f.p. given by (x x) x. xxx - .i.f.p. m -* b.t.p. x .j k ( n.f.p. b.f.p. (cf. I.n) ( b.d.p. ( degenerate no i.f.p. m -- b.d.p. ( x . ( b.t.p. .cf. .n for not degenerate x. fo value of ( ) .(cf. of .r) (cf. I.r) i.f.p. given by (x x x³ - x ² - m -* i.d.p. (n (n) (n-i) x_I(x - .67046846, xi .(n(n) (n) x) .31224737, x .i.d.p. (n (n)(n) x .01728417, x .

Table II (cont.) I.m ( n.f.p. ( b.d.p. b.f.p. (cf. n) ( degenerate no i.f.p. m -* b.d.p. I.n ( ) b.t.p. ( t (not degenerate ( i.f.p. given by (x x x . x - xi.f.p. m -* i.f.p. x .(cf. for value of x o ( ) n.f.p. ( b.d.p. ( not degenerate i.f.p. given by (x x x .x _ .i.f.p. m - i.f.p. X I.p ( b.t.p. ( not degenerate (

Table II (cont.) p. (cont.) i.f.p. given by (x x x x- x- .i.f.p. m i.d.p. x. m - i.d.p. (n ) .42717428, x .xnx) .53818709, x .i.d.p. (n x) .04463863, x( . I.q ( n.f.p. ( b.f.p. (cf. n) ( b.d.p. ( degenerate no i.f.p. m -- n.f.p. with initial vectors x() x) .x) and x) .( x) .m -* b.d.p. with initial vector () . .x) . I.r ( b.t.p. ( t ( f not degenerate ( i.f.p. given by (x x xi X- x- I X .\ i.f.p. m - b.t.p. x.J(cf.

Table II (cont.) a n.f.p.'s ( ( b.f.p. ( J degenerate no i.f.p. m - b.f.p. I.b bd b.d.p. ( ( not degenerate ( i.f.p. given by (x x xl ._-x - X.i.f.p. m -* i.f.p. x . I.c n n.f.p. ( ( not degenerate ( i.f.p. given by (x X X.- i.f.p. m -÷ i.f.p. x . II..a ( ( ( } n.f.p. b.f.p. J doubly degenerate no i.f.p. m - b.f.p.

Table II (cont.) II.l.b ( ( ( n.f.p. ~~~ }degenerate no i.f.p. m -- n.f.p. II.d ( ( ( n.f.p. b.f.p.'s J doubly degenerate X There is a continuum of i.f.p.'s, since -- a. The i.f.p. corresponding to this a is Xx l a Xa) i.f.p. -a) Note that if we write xxY2, xl yi, then the system takes the form tYY which, for a general point, Y - Yl Ym -* appropriate i.f.p. II.l.f ( ( ( ) n.f.p. ~ not degenerate There is a continuum of i.d.p.'s owing to the fact that — Thus, if X a, the i.d.p. is given by x' x(

Table II (cont.) II.f (cont.) x(n xn) n _ n) (n a x)' - - 1-a) i.d.p. )' ) _(n ( xn) a (n) -(cf II.d) This system reduces to the same 2-variable system as does II.c m -* appropriate i.d.p. II.a ( ( } n.f.p. ( ) } degenerate no i.f.p. m - n.f.p. II.b ( ( } n.f.p. ( degenerate i.f.p. given by (x xx x x - .i.f.p. i.f.p. - (cf. I.p.) - II.c ( ( n n.f.p. ( not degenerate no i.f.p. m -* n.f.p.

Table II (cont.) II.d ( ( p's n.f.p.'s ( ~~~( tdoubly degenerate no i.f.p. m n.f.p. (x II.e ( ( n.f.p.'s ( b.f.p. J doubly degenerate no i.f.p. m -* b.f.p. II.f ( ( n.f.p. ( ) doubly degenerate no i.f.p. m - n.f.p. (x II.g ( ( b.f.p. (cf. I.n) ( b.d.p. Jdegenerate no i.f.p. m -* b.d.p.

Table II (cont.) II.h ( ( b.f.p. (cf. I.n) ( b.d.p. J degenerate xl } i.f.p. exists M --b.d.p. ^~~xi.f.p. m -* b.d.p. x II.i ( ( b.d.p. not degenerate xl ) i.f.p. exists x m - ifp 4- i.f.p. X(cf. II.h) II.a ( n.f.p.'s doubly degenerate ( no i.f.p. m -- n.f.p. (xl II.b ( n.f.p. ( b.d.p. ( J not degenerate no i.f.p. m - n.f.p.

Table II (cont.) II.c ( n.f.p. b.f.p. (cf. n) ( ^J bdb.d.p. ( degenerate no i.f.p. m - b.d.p. II.d ( b.t.p. ( not degenerate ( J i.f.p. given by (x xl) - -or (unreduced) x - -(cf. Eq. ( of text) or (y Xx .i.f.p. y - xm —- i.f.p. III.a ( ( ) n.f.p.'s b.f.p. ( J doubly degenerate no i.f.p. m -b.f.p.

Table II (cont.) III.l.b ( ) n.f.p.'s b.f.p. ( J degenerate no i.f.p. m - n.f.p. (x III. .c ( ( ) n.f.p. b.f.p. ( Jdegenerate xi.f.p. exists i m- i.f.p. X .i.f.p. m --- i.f.p. x- . III.l.d ( ( n.f.p. ) b.f.p. ( Jdegenerate 2- v/ - .i.f.p. exists i.f.p. m i.f.p. JV(cf. I.m) .

Table II (cont.) III.l.e ( ( ) n.f.p. ot degenerate ( i.f.p. given by(x xl)x . - 6- ^- xi.f.p. m - i.f.p. x .J III.l.f ( ( ) n.f.p. b.d.p. (see below) ( j not degenerate x. i.f.p. given by (x xx.i.f.p. I - - x)X- /.A non-trivial b.d.p. exists, which is a f.p. of the transformation x' - x [- - x] This f.p. is given by the root of (xX - - which leads to x x) 31944846, x .(n (n) (n l) xx , x.(n) .68055154, x(n £--£ Actually m - i.d.p. x(n x .17899745, x .(n, (n) . x_(n xn 44248, .i.d.p. x( x .61156007, xz(n .

Table II (cont.) III.g ( ( ) b.f.p. ( cf. n) b.d.p. ( degenerate - i .i.f.p. exists X i.f.p. m - i.f.p. - (cf. I.a) X .J III.h ( ( b b.d.p. not degenerate ( i.f.p. given by (x x) x .- - .i.f.p. m -, i.f.p. X- . III.i ( ( b b.d.p. not degenerate ( i.f.p. given by (x - - x.i.f.p. m -, i.f.p. x.J III.a ( n.f.p.'s ( not degenerate ( ) i.f.p. xiXXm - P i.f.p.

Table II (cont.) III.b ( ) n.f.p. ( b.d.p. ( J not degenerate i.f.p. xi Xxm -P i.f.p. III.e ( b.t.p. ( b.t.p. not degenerate ( i.f.p. xxxm - i.f.p. III.a ( n.f.p.'s ( doubly degenerate ( i.f.p. XXX m - i.f.p. (xl with initial vectors xl .xl ..and x...m - n.f.p. (x with initial vector x.x.(See Section XI for discussion) X .

Table II (cont.) III.b ( n.f.p. ( b.d.p. ( J not degenerate i.f.p. x xxm -- b.d.p. III.c ( ) n.f.p. ( sIb.f.p. (cf. I.n) ( bdp b.d.p. ( degenerate i.f.p. Xlxx m - i.f.p. III.e ( ) b.t.p. ( b.t.p. not degenerate ( i.f.p. lXX m -— i.f.p. III.a ( n.f.p.'s ( > triply degenerate ( i.f.p. xX Xm - no convergence (see Section XI for discussion)

Table II (cont.) III.b ( ) n.f-p b.f.p. (cf. I.n) ( b.d.p b.d.p. ( ) degenerate i.f.p. xlx x m -- n.f.p. with initial vectors xl.xl .X.and X..x.m - b.d.p. with initial vector xl .x .x. III.e ( ) b.t.p. ( b not degenerate ( i.f.p. XlX xm - b.t.p. I.f ( b.t.p. ( b.t.p. not degenerate ( i.f.p. Xxxm - b.t.p.

Introduction to Table III

This table lists a few 4-variable generalizations (and two 5-variable ones) of selected 3-variable systems. The method of generalization is explained in the text (Section VII, 1.5.j - ext - I and subsequent discussion). The basic notation is that of Table II, but degeneracy, b.f.p.'s, n.f.p.'s, etc., are not noted.

Table III—
Examples of Binary Reaction Systems for N>3

I.l.c - ext - ( ( ( ( f( xi . i.f.p. given by (x xi)x) - X. m - i.f.p. X . I.a - ext - ( ( ( ( x i.f.p. given by (x x .x- - x - x . m - .i.f.p. X.

Table III (cont.) k - ext - ( ( ( ( i.f.p. given by (x xl) X. (- - X .m -i.f.p. . I.p - ext - ( ( ( v/ xl ()x i i.f.p. m - b.d.p. -/-xs ^ (1,0,0,, (0,1,0, X . I.g - ext - ( ( ( ( i.f.p. gienby (x( xxl . i.f.p. given by ( x ) .- . m - b.t.p. (1,00,,, (0,0,1,, (0,0,0, X.

Table III (cont.) I.i - ext - ( ( ( )x .( .. i.f.p. given by (x xl) x .i.f x - .m -- i.f.p. with initial vectors: x(⁾ .x(⁾ x⁾ .and x⁾ .x⁾ - .x^® .x) .x().m -b.t.p. with initial vector: x) .( . .(.I.e - ext - ( ( ( (x i.f.p. given by (x x .xxx - X.nm - b.d.p. Z.J

Table III (cont.) m -* interior periodic points of period 12: Cycle

Table III (cont.) m -, interior periodic points of period 6: Cycle i.f.p. given by

Table III (cont.) b - ext - ( ( ~~( ( ~xl . ~(~)~ x ( .i f.p' i.f.p. given by (x x X. - - ² - m - i.f.p. This i.f.p. is very close to that of II.d - ext - (q.v.), as might be expected. I.c - ext - ( ( ( ( i.f.p. given by (x x) . xxxx X .i.f Xm - •i.f.p. X.(cf. I.b - ext - II.d - ext - ( ( ( ( xl .} i.f.p. given by (x x X .X X - ..f.. x .m - i '.f.p. X.(cf. Eq. ( of the text)

Table III (cont.) II.d - ext - ( ( ( ( ( x .i.f.p. given by (x xi) .x xxx x- x.i.f.p. (cf. Eq. ( of the text) .m i.f.p. X. III. .c - ext - ( ( ( ( x, .. f • ^ ( ^ X .i.f.p. given by (x X i.f.p. x.- - (-^x - . m -* i.f.p. III.f- ext - ( ( ( i( . i.f.p. given by (x XX .- x - xx -X ( m i.f.p. X.

Explanation of Tables IV and V

Each of these tables consists in a tabulation of the convergence behavior of random initial points under 26 transformations which are simple modifications of I.3.g. These 26 systems are generated as follows (see Section VIII.2 of text):

In Table IV, we retain the "skeleton" x1 =1/2x2 + 1/2x+ 2xx₂ + ..₂=2xlx₃+ 2x₂X₃+ ... x_' =1/2x2 + ...

To this we add in all possible ways the missing terms 1/2l² , 1/2x2, 1/2x3. Twenty-six new systems result (the 27th is identical with 1.3.g). The notation is the same as in the previous tables, but we give only the numerical results, omitting comments and i.f.p. equations.

Table V is similar, except that the skeleton is x] =x2 + x3 + x3x2 +.. 3 =3x ₂ + 2₃ + ₁2+ . x₂=1XX3+ X₂X₃+ ... x3 ~l +... and the missing terms are x₁x ₂, xlx ₃, x ₂x ₃.

Table IV—
Modifications of System I.3.g.

Case ( ( . x.i.f.p. x. m -— i.f.p. Case x . .i.f.p. m - i.f.p. .

Table IV (cont.) Case x . x .\ i.f.p. ( X .m -, i.f.p. Case ( x. x.i.f.p. X.m - i.f.p. Case ( . x.i.f.p. X.m -* i.f.p. Case ( ( x . x.i.f.p. X .m -, i.f.p. Case ( x . x.i.f.p. X . m -, i.f.p.

Table IV (cont.) Case ( xi . X., i.f.p. X .m -. i.f.p. Case xi . x .i.f.p. X. m - i.f.p. Case xl . x .i.f.p. X .m - i.f.p. Case xi . X.i.f.p. ( m , i.f.p. Case . .i.f.p. X . m -P i.f.p.

Table IV (cont.) Case ( xi . .i.f.p. .m -. i.f.p. Case x_I . x .i.f.p. x.m -i.f.p. Case xl . x .i.f.p. ( x.Jm -* i.f.p. Case xi . .i.f.p. ( x.m -* i.f.p. Case xl . x.i.f.p. x . m -, i.f.p.

Table IV (cont.) Case ( x x.i.f.p. ( .m -> i.f.p. Case ( x . x .i.f.p. X.m -- i.f.p. (Very slow convergence; IJacobiani .at i.f.p.) Case x .x - -.i.f.p. m i.f.p. X. Case xl . x.) i.f.p. . m -- i.f.p. Case . x.i.f.p. X .m - i.f.p.

Table IV (cont.) Case ( xl x .\ i.f.p. x.m -l i.f.p. Case ( . x .) i.f.p. ( X.i.f.p. given by (x xl) - - m -, i.t.p. (n () .87777286, x(n) .00766150, x(n .(n x).00011739, x (nl) .22185332, x - .i.t.p. Xn x) .12210975, X(.77048518, ( Case ( xi . x .) i.f.p. .m -i.f.p. Case ( x . .) i.f.p. ( X. rn - i.f.p.

Table V—
Modifications of System I.3.g.

Case ( ( ( ( .( ( x.i.f.p. ( X.m - i.f.p. Case ( ( ( . ( x.i.f.p. ( x.m - i.f.p. Case ( ( ( x_I .( ( x.i.f.p. ( ( ( ( X.m -, i.f.p. Case ( ( ( x ( xi.f.p. ( Xm -- i.f.p.

Table V (cont.) Case ( ( ( x . ( X.i.f.p. ( X .m -i.f.p. (slow convergence) Case ( ( ( ( ( xi ( ( ( x.i.f. ( X.m - i.f.p. (cf. Case Case ( ( ( xi .( ( x.i.f.p. ( ( X . nm - i.f.p. Case ( ( ( xi .( ( x .i.f.p. ( ( x . m - i.f.p. Case ( ( ( ( ( x ( ( x. ( ( X .m i.f.p. (cf. Case

Table V (cont.) Case ( ( ( xi ( ( x.i.f.p. ( ( X .m -- i.f.p. Case ( ( ( ( ( x .i.f.p. ( ( X .m -, i.f.p. Case ( ( ( x .i.f.p. ( ( X .m - i.f.p. Case ( ( ( ( . ( ( X.i.f.p. ( X.m -— i.f.p. Case ( ( ( ( x( ( x .i.f.p. ( X .m -i.f.p.

Table V (cont.) Case ( ( ( XXX i.f.p. ( ( ( ( m -* i.f.p. Case ( ( ( xl xx i.f.p. ( ( ( ( m -* i.f.p. Case ( ( ( xi X i.f.p. ( ( ( ( ( ( m -* i.f.p. Case ( ( ( ( xi .i.f.p. ( ( x. ( ( ( (cf. Case mn -- i.f.p. Case ( ( ( ( .( ~ ~(~.i.f.p. ( ( x.( ( ( .m -* i.f.p.

Table V (cont.) Case ( ( xl X .} . ( ( .( ( ( (cf. Case m -, i.f.p. Case ( ( x . ( x.i.f.p. ( ( X.J m - i.f.p. Case ( ( ( XX .( ( (cf. Case m -> i.f.p. Case ( ( ( ( x ( ( ( x.i.f.p. ( ( X.m -* i.f.p. Case ( ( ( ( x . ( XX .J ( ( (cf. Case m - i.f.p.

Table V (cont.) Case ( ( ( ( x .( ( ( x .i.f.p. ( ( . m -* i.f.p. Case ( ( ( ( . ( X .i.f.p. ( ( X.m -, i.f.p.

Table VI—
Super-Systems

e - Super a. Initial Configuration yi .y .y.Y .s .Y .y .Y .yg .b. Initial Configuration yl ..y .Y.Y .Y.Y ..yg .Both of these gave the final periodic configuration (Yi ( y(n) .yn y(n .) .yn) y) - .y V y(n) (n)y(n)(n n .y(n (n) (n) (n) Y Y?n) ) (n _ (n- .Y .y For the initial configuration c. yi .Y ..Y .Y .Y .y .Y .yg .07284413, the final configuration was ( y(): y(n .y- .(n) y) .y .y(n) (n) (n) (n)y(n+- 31944846 y(n)= 0 y(n) =.24574926 (n+l) (n+)0 (n+l) =Y = 2492 6 = . 290536 9 -0(n)= .24574926 )= 40 = 18905302 Y3 .18905302

Table VI (cont.) I.g- Super Initial Configurations: a. yi ...Y..Y.Y .Y.yg .b. yl .y..Y.y.Y.Y .Y .y .These both gave the final configuration (y n y(n)): y ) .(n) .y.y) .y(n) .) - .) .(n) () .y y_ (n) yn .(n .._~n) ( _ ) ~ (n)yl .y .y .(n) y(n (nya^ .Y) .Y.y(n y . .(n (n (n Y.Yn .Y .n .Y. .Initial Configuration: C. yi ..y .Y.Y.Y.y .Y .y .This gave the final configuration (yn yi: y ) .y .y() .y(n) .) .( .y) .y() .(n) _n _(n)(n y() .(n) y ._n- _(n .yn) .Y) ..( _) _(n (n Y- .Y ..(n (n (ny( ...(n y (n .Y...(n .(n (n _ .Y.YS~~~~~~n_s .

Appendix

1. In this appendix we summarize the results for a class of homogeneous quadratic transformations quite distinct from the class we have called Binary Reaction Systems. This class arises in a natural way in the study of a certain crude model of the "evolutionary process" which will be described below. It also has some mathematical interest, owing to the fact that the limiting behavior of all systems belonging to this class can be explicitly predicted.

2.The Evolutionary Model

Consider a large population in which each distinct "type" of individual is labelled by an index pair (i,j), i,j = 1,2,... ,N. Let the fraction of the male population which is of type (i,j) be denoted by xij. We shall assume that xij = xji; also, we take the number of females of type (i,j) to be equal to the number of males of this type-hence there is no need to denote the fraction of females by a separate letter. We now impose a mating rule (random mating is assumed) which states, in effect, that if individuals of type (k, ) mate with individuals of type (m,n), the progeny will be of all types (i,j) such that min(k, m)<i <max(k, m) min(e, n)<j< max(f, n)

Loosely speaking, we may call each index of a pair a "characteristic." A given mating will produce all possible children such that (1) is satisfied, the distribution of the two indices determining the progeny being the product of two identical distributions. The number of children will, of course, be proportional to the number of parents of each type. Mathematically: x(n+1) _ ,kmnxn(n))(n) (2) ^ij~/--li lj X k£ -"mn k,£,m,n

The sum in (2) is to be carried out under the restriction (1). We specify the system further by postulating: m= m mk > 0, min(k, m)<i < max(k, m) z - (3) = 0 otherwise; k i/kN 1 (4) i=m

ikrn m + k-EZkm=rn+k 2(5) i=m In addition, we normalize by taking: N Zx)- = 1,<x⁾< 1, alli,j (6) i,j=1 It is then evident that we have: N E ij i,j for all n.

The class of systems defined above is the one we shall actually discuss. However, it may be of interest to at least mention the actual evolutionary model in connection with which equation (2) arises. Evolution is assumed to take place by mutation. A type (i,j) can give rise to two new types (i+1,j) and (i,j+1) with some small probability. When we include this (linear) effect, we get a series of equations: Xij = -exij + e/2(xi-l,j + xi,j-) + Z 7kmnXktxmn (7)

Here E is taken as some small number. We actually performed many numerical experiments on systems of the form (7) with special values of the <ym satisfying (3), (4), (5). Two particularly convenient choices are: -= 2-k i-min(j,k)) (8) and k = Jj-kJ+ , min(j, k) <i< max(j, k) = 0 otherwise

Note that with our definition (3) we may take the sum in (2) as unrestricted, i.e., over-all k, , m,n = 1/2,..., N. We shall not discuss the behavior of the "mutating" system (7) in this report, but rather restrict ourselves to the pure "mating" system (2).*

It may be objected that our mating rules have nothing to do with Mendel's Laws. This is intentionally the case. The Mendelian case has been treated in great generality in a series of papers by Hilda Geiringer;⁷ , See also C. C. Li:⁸

3. If we sum over one of the indices in (2), we obtain the system: N C (^n+l) = E "kmFC (n)(n) i1,... N (10) k,m=l N C_i-xij (11) j=l and, of course, N E C(n)= 1 all n. (12) i=l

By virtue of condition (5) on the km, the system (10) possesses a linear invariant (distinct from E Ci = 1); in fact: E jC(n+ - EC()C(n)E jm = E (n)C(n)(e+ ) jm 29 j e,m j £ ,m = E7 cj) (13)

The consequences of this property are very interesting. It turns out that the existence of this linear invariant enables one to predict explicitly the limiting behavior of any initial vector (C(O) C°O),...CO)) when iterated according to (10). Using the fact that EN Ci = 1 is also an invariant, we may define an invariant: N-1 a--, (N - i)Ci (14) i=l

a is, of course, explicitly determined by the initial vector. It can then be proved that every initial vector will converge to a definite fixed point which is determined as follows: For the given value of a, there is one index j such that: N-j>a> N-j-1 (15) The f.p. is then explicitly given by: Cj = a- (N - j - 1) Cj+i = N-j-a o (16) all other Ci = 0

The f.p. is independent of the actual values of the coefficients "km providing they satisfy (3), (4), and (5).* These results can easily be referred back to the original variables xij. Defining a quantity a from (16) by: Cj - l+ (17) we find that the corresponding symmetric tensor (x(°) will converge to the final state: 1 3 (1 + a)²Xj+l,j = j,i+l- (1 + a)² (18) a²Xjj+I,j= (1 + )²

The results can also be extended to the case of M "characteristics" il, i2,- ·-iM The fraction of the population of type (il,. . .,iM) is then denoted by xii...iM. This is taken to be symmetric in all M indices; hence, there are: (N+M-1) M

types of individuals in the population. The corresponding mating rule is: X(n+1) = E \^,1S,354 S2M-1S2M il...iM...~S1S2 S₃S₄...• ' ' 'l 'Yi₂M' ' ''M S₁,S₂S ₃,S ₄(n () (19) S1S3...S₂M-1S2S₄...S2M

If we sum over M- 1 indices, we can again define:

*Of course the rate of convergence will, in general, depend on the actual values of the "km. In one case it is actually possible to solve explicitly for the iterates as functions of n. This is the case N = 2 where we find: (n)ii - n [O) 2(2 "X~l^)--2- Ix?) + or (2 - 1)] X(0) (0) 11 — lx₁₂ This case can be viewed as a generalization of the Mendelian law for a single gene. The population has the same limiting configuration in both cases, but for Mendel's rules equilibrium is reached in a single step "(Hardy's Law").

N Ci il...iM (20) i2,...,iM=l and we again get the system (10). In terms of the tensor xi, . . iM, the final state is: •• i = (1 + a)Mxj,j+ ...j = (1 + )^M = (1 + a)M (21) Xj+l,j+l...I+ = (1 + c)M Here j and a are determined as in (16) and (17).

4. The explicit nature of these results is due to the existence of the linear invariant a, which in turn is a consequence of the "meanpreserving" property (5) assumed for the coefficients yijk. However, the actual convergence properties are probably more closely connected with the "index-limiting" condition (3). (Our Binary Reaction systems do not, in general, have this property; the existence of oscillating final states may well be a consequence of allowing such mating rules as jEj - k 5 j.) Some cases have been investigated in which the conditions on the coefficients were relaxed. In particular, we have considered systems for which we no longer require: jk > 0 all j, k. As an example, we chose the case: N yi = bj,k + 6j+l,k + 6j-l,k (22) i=l

This corresponds to a sort of "selective mating scheme" in which only individuals with "nearby"- indices can mate. Of course we then no longer have E Ci = 1. In order to secure this property, we must "renormalize" at every step, i.e., set: C'i *

If this is formally written out, the system is no longer quadratic. There are then many types of fixed points possible. So far it has not been possible to give an analytical treatment of such systems.*

References

1. See, e.g., R. Sauer, Math. Ann.106, 722, (1932), 0. Baier, Math. Ann. 112, 630.

2. Coddington and Levinson, Theory of Ordinary Differential Equations, Chapters 15, 16, McGraw-Hill, New York (1955).

3. V. Volterra, Leqons sur la The6rie Mathematique de la Lutte Pour La Vie, Gauthier-Villars, Paris (1939).

4. H. Weyl, Math. Ann. 77, 313 (1916).

5. S. Ulam and J. von Neumann, Bull. Am. Math. Soc. 53, 1120 (1947).

6. L. E. Dickson, Modern Algebraic Theories, Chapter II, Sanborn, Chicago, (1930).

7. Hilda Geiriger, Ann. Math. Statistics, 15, 25-57 (1944), Genetics,33, 548-564, (1948). 8. C. C. Li, Population Genetics, Univ. of Chicago Press, Chicago (1955).

*Note, however, that if we simply take 7^j k = ,we find that every initial vector is a fixed point.

Enrico Fermi

Johnny von Neumann with 11 year old Claire Ulam

Bob Richtmyer

David Hawkins

From left to right, Ulam, Mycielski, Bednarek

Bill Beyer

"Computers can also be used to investigate... and with less success, to study games of 'skill' like chess." (See page 302.) P. R. Stein playing the first game of chess he and others had programmed with Ulam in 1957 against the MANIAC.