Abstract

A one-dimensional dynamical system of 64 particles with forces between neighbors containing nonlinear terms has been studied on the Los Alamos computer MANIAC I. The nonlinear terms considered are quadratic, cubic, and broken linear types. The results are analyzed into Fourier components and plotted as a function of time.

The results show very little, if any, tendency toward equipartition of energy among the degrees of freedom.*

*After the untimely death of Professor Fermi in November 1954, the calculations were continued in Los Alamos. The last few examples were calculated in 1955.

This report is intended to be the first one of a series dealing with the behavior of certain physical systems where the nonlinearity is introduced as a perturbation to a primarily linear problem. The behavior of the systems is to be studied for times which are long compared to the characteristic periods of the corresponding linear problems.

The problems in question do not seem to admit of analytic solutions in closed form, and heuristic work was performed numerically on a fast electronic computing machine (MANIAC I at Los Alamos).* The ergodic behavior of such systems was studied with the primary aim of establishing, experimentally, the rate of approach to the equipartition of energy among the various degrees of freedom of the system. Several problems will be considered in order of increasing complexity. This paper is devoted to the first one only.

We imagine a one-dimensional continuum with the ends kept fixed and with forces acting on the elements of this string. In addition to the usual linear term expressing the dependence of the force on the displacement of the element, this force contains higher order terms. For the purpose of numerical work this continuum is replaced by a finite number of points (at most 64 in our actual computation) so that the partial differential equation defining the motion of this string is replaced by a finite number of total differential equations. We have, therefore, a dynamical system of 64 particles with forces acting between neighbors with fixed end points. If xi denotes the displacement of the i-th point from its original position, and a denotes the cofficient of the quadratic term in the force between the neighboring mass points and f³ that of the cubic term, the equations were either ,i==(Zi+i + Xi-l - 2Xi) + a[(Xi+l - Xi)² - (Xi - xi₁ )² ] i 1, 2,... 64, or xi =(Xi+i + xi-i - 2xi) + 13 [(xi+₁ - xi)³ - (x_i- Xi_)^3] i = 1, 2,...64,

a and 3 were chosen so that at the maximum displacement the nonlinear term was small, e.g., of the order of one-tenth of the linear term. The corresponding partial differential equation obtained by letting the number of particles become infinite is the usual wave equation plus nonlinear terms of a complicated nature.

* We thank Miss Mary Tsingou for efficient coding of the problems and for running the computations on the Los Alamos MANIAC machine.

Another case studied recently was :xi = (6I(Xi+ I - Xi) - 62 (Xi - Xi--1) + c (3)

where the parameters 61, 62, c were not constant but assumed different values depending on whether or not the quantities in parentheses were less than or greater than a certain value fixed in advance. This prescription amounts to assuming the force as a broken linear function of the displacement. This broken linear function imitates to some extent a cubic dependence. We show the graphs representing the force as a function of displacement in three cases. Quadratic Cubic Broken Linear

The solution to the corresponding linear problem is a periodic vibration of the string. If the initial position of the string is, say, a single sine wave, the string will oscillate in this mode indefinitely. Starting with the string in a simple configuration, for example in the first mode (or in other problems, starting with a combination of a few low modes), the purpose of our computations was to see how, due to nonlinear forces perturbing the periodic linear solution, the string would assume more and more complicated shapes, and, for t tending to infinity, would get into states where all the Fourier modes acquire increasing importance. In order to see this, the shape of the string, that is to say, x as a function of i and the kinetic energy as a function i were analyzed periodically in Fourier series. Since the problem can be considered one of dynamics, this analysis amounts to a Lagrangian change of variables: instead of the original xi and xi, i = 1, 2, ... 64, we may introduce ak and ak, k = 1, 2,... 64, where ak = xi sin 64 (4) The sum of kinetic and potential energies in the problem with a quadratic force is

Ekin Epot 1 .2(xi+l - Xi)² + (xi-i-1)²i + Ei- i + (5a) iWkink222 7 k Eak +E -pt = + 2a sin (5b) a, ₂k+n 2128 if we neglect the contributions to potential energy from the quadratic or higher terms in the force. This amounts in our case to at most a few per cent.

The calculation of the motion was performed in the x variables, and every few hundred cycles the quantities referring to the a variables were computed by the above formulas. It should be noted here that the calculation of the motion could be performed directly in ak and ak. The formulas, however, become unwieldy and the computation, even on an electronic computer, would take a long time. The computation in the ak variables could have been more instructive for the purpose of observing directly the interaction between the ak's. It is proposed to do a few such calculations in the near future to observe more directly the properties of the equations for ak.

Let us say here that the results of our computations show features which were, from the beginning, surprising to us. Instead of a gradual, continuous flow of energy from the first mode to the higher modes, all of the problems show an entirely different behavior. Starting in one problem with a quadratic force and a pure sine wave as the initial position of the string, we indeed observe initially a gradual increase of energy in the higher modes as predicted (e.g., by Rayleigh in an infinitesimal analysis). Mode 2 starts increasing first, followed by mode 3, and so on. Later on, however, this gradual sharing of energy among successive modes ceases. Instead, it is one or the other mode that predominates. For example, mode 2 decides, as it were, to increase rather rapidly at the cost of all other modes and becomes predominant. At one time, it has more energy than all the others put together! Then mode 3 undertakes this role. It is only the first few modes which exchange energy among themselves and they do this in a rather regular fashion. Finally, at a later time mode I comes back to within one per cent of its initial value so that the system seems to be almost periodic. All our problems have at least this one feature in common. Instead of gradual increase of all the higher modes, the energy is exchanged, essentially, among only a certain few. It is, therefore, very hard to observe the rate of "thermalization" or mixing in our problem, and this was the initial purpose of the calculation.

If one should look at the problem from the point of view of statistical mechanics, the situation could be described as follows: the phase

space of a point representing our entire system has a great number of dimensions. Only a very small part of its volume is represented by the regions where only one or a few out of all possible Fourier modes have divided among themselves almost all the available energy. If our system with nonlinear forces acting between the neighboring points should serve as a good example of a transformation of the phase space which is ergodic or metrically transitive, then the trajectory of almost every point should be everywhere dense in the whole phase space. With overwhelming probability this should also be true of the point which at time t = 0 represents our initial configuration, and this point should spend most of its time in regions corresponding to the equipartition of energy among various degrees of freedom. As will be seen from the results this seems hardly the case. We have plotted (Figs. 1-9) the ergodic sojourn times in certain subsets of our phase space. These may show a tendency to approach limits as guaranteed by the ergodic theorem. These limits, however, do not seem to correspond to equipartition even in the time average. Certainly, there seems to be very little, if any, tendency towards equipartition of energy among all degrees of freedom at a given time. In other words, the systems certainly do not show mixing.*

The general features of our computation are these: in each problem, the system was started from rest to time t = 0. The derivatives in time, of course, were replaced for the purpose of numerical work by difference expressions. The length of time cycle used varied somewhat from problem to problem. What corresponded in the linear problem to a full period of the motion was divided into a large number of time cycles (up to 500) in the computation. Each problem ran through many "would-be-periods" of the linear problem, so the number of time cycles in each computation ran to many thousands. That is to say, the number of swings of the string was of the order of several hundred, if by a swing we understand the period of the initial configuration in the corresponding linear problem. The distribution of energy in the Fourier modes was noted after every few hundred of the computation cycles. The accuracy of the numerical work was checked by the constancy of the quantity representing the total energy. In some cases, for checking purposes, the corresponding linear problems were run and these behaved correctly within one per cent or so, even after 10,000 or more cycles.

It is not easy to summarize the results of the various special cases. One feature which they have in common is familiar from certain problems in mechanics of systems with a few degrees of freedom. In the

* One should distinguish between metric transitivity or ergodic behavior and the stronger property of mixing.

compound pendulum problem one has a transformation of energy from one degree of freedom to another and back again, and not a continually increasing sharing of energy between the two. What is perhaps surprising in our problem is that this kind of behavior still appears in systems with, say, 16 or more degrees of freedom.

What is suggested by these special results is that in certain problems which are approximately linear, the existence of quasi-states may be conjectured.

In a linear problem the tendency of the system to approach a fixed "state" amounts, mathematically, to convergence of iterates of a transformation in accordance with an algebraic theorem due to Frobenius and Perron. This theorem may be stated roughly in the following way. Let A be a matrix with positive elements. Consider the linear transformation of the n-dimensional space defined by this matrix. One can assert that if x is any vector with all of its components positive, and if A is applied repeatedly to this vector, the directions of the vectors xi, A(x), ..., Ai(x), ..., will approach that of a fixed vector xo in such a way that A(Xo) = A(2o). This eigenvector is unique among all vectors with all their components non-negative. If we consider a linear problem and apply this theorem, we shall expect the system to approach a steady state described by the invariant vector. Such behavior is in a sense diametrically opposite to an ergodic motion and is due to a very special character, linearity of the transformations of the phase space. The results of our calculation on the nonlinear vibrating string suggest that in the case of transformations which are approximately linear, differing from linear ones by terms which are very simple in the algebraic sense (quadratic or cubic in our case), something analogous to the convergence to eigenstates may obtain.

One could perhaps conjecture a corresponding theorem. Let Q be a transformation of a n-dimensional space which is nonlinear but is still rather simple algebraically (let us say, quadratic in all the coordinates). Consider any vector x and the iterates of the transformation Q acting on the vector x. In general, there will be no question of convergence of these vectors Qn(x) to a fixed direction.

But a weaker statement is perhaps true. The directions of the vectors Qn(x) sweep out certain cones C, or solid angles in space in such a fashion that the time averages, i.e., the time spent by Qn(x) in C, exist for n -— oo. These time averages may depend on the initial 2 but are able to assume only a finite number of different values, given C,. In other words, the space of all direction divides into a finite number of regions Ri, i = 1,... k, such that for all vectors x taken from any one of these regions the percentage of time spent by images of x under the Q" is the same in any C,.

The graphs which follow show the behavior of the energy residing in various modes as a function of time; for example, in Fig. I the energy content of each of the first 5 modes is plotted. The abscissa is time measured in computational cycles, bt, although figure captions give 6t² since this is the term involved directly in the computation of the acceleration of each point. In all problems the mass of each point is assumed to be unity; the amplitude of the displacement of each point is normalized to a maximum of 1. N denotes the number of points and therefore the number of modes present in the calculation, a denotes the coefficient of the quadratic term, and / that of the cubic term in the force between neighboring mass points.

We repeat that in all our problems we started the calculation from the string at rest at t = 0. The ends of the string are kept fixed.

t In Thousands Of Cycles

Fig. 1. The quantity plotted is the energy (kinetic plus potential in each of the first five modes). The units for energy are arbitrary. N = 32; a = 1/4; 8t² = 1/8. The initial form of the string was a single sine wave. The higher modes never exceeded in energy 20 of our units. About 30,000 computation cycles were calculated.

'ig. 2. Same conditions as Fig. 1, but the quadratic term in the force was stronger. = 1. About 14,000 cycles were computed.

Fig. 3. Same conditions as in Fig. 1, but the initial configuration of the string was a "saw-tooth" triangular-shaped wave. Already at t = 0, therefore, energy was present in some modes other than 1. However, modes 5 and higher never exceeded 40 of our units.

Fig. 4. The initial configuration assumed was a single sine wave; the force had a cubic term with 3 = 8 and 6t² = 1/8. Since a cubic force acts symmetrically (in contrast to a quadratic force), the string will forever keep its symmetry and the effective number of particles for the computation N = 16. The even modes will have energy 0.

Fig. 5.N = 32; 8t² = 1/64; / = 1/16. The initial configuration was a combination of 2 modes. The initial energy was chosen to be 2/3 in mode 5 and 1/3 in mode 7.

Fig. 6. 5t² = 2^-6 . The force was taken as a broken linear function of displacement. The amplitude at which the slope changes was taken as 2-5 + 2^-7 of the maximum amplitude. After this cut-off value, the force was assumed still linear but the slope increased by 25 per cent. The effective N = 16.

Fig. 7. 6t² = 2-. Force is again broken linear function with the same cut-off, but the slope after that increased by 50 percent instead of the 25 percent charge as in problem 6. The effective N = 16.

Position Of The Mass Point Fig. 8. This drawing shows not the energy but the actual shapes, i.e., the displacement of the string at various times (in cycles) indicated on each curve. The problem is that of Fig. 1.

t In Thousands Of Cycles Fig. 9. This graph refers to the problem of Fig. 6. The curves, numbered 1, 2, 3, 4, show the time averages of the kinetic energy contained in the first 4 modes as a function of time. In other words, the quantity is i E Z1 i=l v is the cycle number, k = 1, 3, 5, 7.