by Colin Pask
Figure 12.5. (a) Example of the FPU problem for eight particles. (b) The FPU system for two particles and the two modes that are found in the linear system. In the first mode, the particles oscillate together; in the second mode, they oscillate with the same frequency and amplitude but in opposite directions. Figure created by Annabelle Boag.
The equations to be solved are Newton's second law of motion (mass times acceleration equals force) for each of the particles. Let the position of the ith particle along the chain be xi and assume that masses and linear spring constants are both scaled to be 1. Then the FPU report gives the equations to be used for two choices of nonlinearity as
There is an equation like this for each particle, so for the largest number of particles considered (which was 64) we must put i = 1, 2, 3,…64. Two possible forms for the nonlinear force components are used, and the parameters α and β measure their strengths. Taking α and β as zero gives the linear problem for which the modes can be found and they are given in the FPU report. They are extensions of those given in figure 12.5 (b) but obviously more involved. They may be thought of as much like the modes in figure 12.1, with the first mode having all particles moving together. The calculation is carried out using the linear modes as a basis and discretizing the differential equations so the computer can time-step along the solutions.
The report, “Studies of Nonlinear Problems,” now describes the calculation to be made:
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.6
So the idea was to use one mode as a starting point and see how the system evolved over time. I can do no better than to quote from the report:
Let us say here that the results of our computations show features which were, from the beginning, surprising to us. [Even a Nobel Prize winner famous for his intuition was surprised!] Instead of a gradual continuous flow of energy from the first mode to the higher modes, all 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 [mode one] we indeed observe 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 the 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 that exchange energy among themselves and they do this in a rather regular fashion. Finally, at a later time mode 1 comes back to within one percent of its initial value so that the system seems to be almost periodic. All our problems have at least this one feature in common…. 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.
An example of FPU results for 32 particles is shown in figure 12.6. The system starts in mode 1 and after about 29,000 cycles (measured by the linear mode frequency) it has pretty much returned to that modal state. No wonder they were surprised! Incidentally, as part of this new approach or “mathematical experiment” there is a move toward visual presentation of results and now elaborate graphics is a feature of many papers. A whole set of FPU calculations with varying particle numbers and force parameters show similar results. In terms of the aim of exhibiting the way in which nonlinearity leads to mixing of modes and certain statistical properties for the system, the numerical experiment was a failure and the authors write that:
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.
Figure 12.6. Modal energies for the first five modes in an FPU calculation using thirty-two particles and quadratic nonlinear forces. The curves are labeled by the mode numbers. From E. Fermi, J. Pasta and S. Ulam, “Studies of Nonlinear Problems,” Los Alamos Document LA-1940 (May 1955).
While the FPU results themselves were surprising, it will come as no surprise to learn that they have led to an enormous amount of work on many diverse aspects of the problem and the system under investigation. The computer power available to modern researchers far exceeds that available for the FPU calculation and thus more extensive results have been generated. Exploring what happens when the various parameters are given different and sometimes extreme values has shown that the oscillations can take on a variety of behaviors including that chaotic form to be discussed in the next section.
The paper by Pettini and colleagues (see bibliography) explains how, as the energy of the initial state is increased, there is a transition first from a regular type of motion to a partially chaotic regime, and then to a fully chaotic situation where full mode mixing is indeed present. The continuum limit (where the number of particles approaches infinity and the space between them becomes zero) leads to a whole new area of science as briefly outlined below in section 12.4.2. This was a calculation that introduced a new approach for studying physical systems, produced intriguing results, and spawned a whole area of science; calculation 49, FPU and oscillation surprises surely must find a place in any list of important calculations.
(Any literature search will reveal a great many papers and books about the FPU calculations, their impact, and the new research areas they have generated. I recommend the paper by Porter, Zabusky, Hu, and Campbell for a well-illustrated general introduction. The topic of “computational synergetics” is authoritatively introduced and reviewed by Norman Zabusky, who has played a major part in this area of science. The Status Report book edited by G. Gallavotti contains eight up-to-date reviews and also reprints the original FPU report. The review, “The Fermi-Pasta-Ulam Problem: Fifty Years of Progress,” by Berman and Izrailev appears in volume 50 of the journal CHAOS, which contains a whole series of reviews and progress reports marking their fiftieth anniversary. A paper by Akhmediev shows how the FPU problem is relevant to observations of light waves in optical fibers. Kevrekidis gives an extensive survey of the theory of nonlinear lattice waves, an area of research spawned by the FPU work.)
12.4.2 Continuous Systems
Historically, the initial approach to the stretched-string problem considered a finite number of masses joined by elastic elements (see Kibble for a simple introduction). The limit of an infinite number of masses leads to the continuous stretched-string model as described by the differential equation (12.1). Thus it was natural to ask about the continuous form of the FPU problem. For the quadratic nonlinearity case (non-zero α in equation (12.11)), and with suitably defined displacement u, time variable, and constant parameter δ, the relevant equation is
(See the Zabusky paper for further details.) This is known as the KdV equation since it was first written down by D. J. Korteweg and G. de Vries in 1895. They were investigating the propagation of shallow water waves, in which case u refers to height of water in the wave.
The KdV equation has a remarkable solution comprising a wave of given height and shape which propagates unchanged. (This contrasts with the usual linear waves built from many Fourier components and so dispersing since those components travel with different speeds and hence get out of step. In the nonlinear case, the linear dispersion is cancelled out by the nonlinear effects.) These are now termed solitary waves, or solitons. There is a wonderful story of John Scott Russell who, in 1834, followed such a wave on horseback as it traveled along the Union Canal near Edinburgh. The paper by Porter, Zabusky, Hu, and Campbell has a wonderful photogra
ph of surfers on the famous solitary waves traveling along the Severn River in England.
In the FPU calculations, the emphasis was on modal energies rather than actual total displacement shapes for the system (only one of the nine graphs deals with the spatial effects). Thus the surprising recurrence phenomena detected by Fermi, Pasta, and Ulam was not linked by them to waves like solitons, which provide a natural explanation for such phenomena.
The area of soliton physics has grown enormously with many wonderful results and explanations of physical effects driven by nonlinearities. In fact, some of this work might well have made it on to a list of important calculations.
12.4.3 A Mysterious Lady
When the use of the MANIAC computer is mentioned in the FPU report, there is a footnote which reads:
We thank Miss Mary Tsingou for efficient coding of the problems and for running the computations on the Los Alamos MANIAC machine.7
You might think, well that is nice to give a thank-you to an assistant; but was that really enough? Remember, this work was done in the 1950s when electronic computers were in their infancy. To convert a mathematical problem into an algorithm suitable for a machine to follow, then code that algorithm to implement it, and finally to run it on a computer was a major task. Problems that would have been considered immensely difficult to program in 1950 are today tackled by writing a few lines of program in MAPLE or some other computing package, so it may be hard to appreciate just what the early pioneers of computer experiments achieved. (Figure 2 in Thierry Dauxois's article about Mary Tsingou shows the intricate and extensive flow diagram used to program the computer to tackle the FPU problem. Also, according to Dauxois, “Mary Tsingou and John Pasta were the first ones to create graphics on the computer, when they considered the problem with an explosion and visualized it on an oscilloscope.”8) We saw above that a graphical presentation was the most suitable way to present the FPU results. Clearly, Mary Tsingou was more than just some everyday programmer; it seems likely that she was an essential and integral part of the team leading to the FPU report. Was this another case of a woman in science missing out on some well deserved credit?
Mary Tsingou obtained a bachelor of science degree in 1951 and a master of science in mathematics in 1955. In that era of the Korean War there was a shortage of young men for positions at Los Alamos, and so Mary was part of a group of young women hired to do hand calculations. Obviously she gained entry to the new world of electronic computers and demonstrated a talent for programming. In 1958, she married and became Mary Tingsou Menzel, and she later published scientific papers under that name. More details of this story are in Dauxois's article, and like him, I wonder if really to be fair we should not talk about the FPUT problem.
12.5 DISCOVERING CHAOS
In the preface I wrote that the calculations discussed in this book “are mostly relatively simple (at least in principle, if not always in practice) but still of major significance.” In chapter 1 I gave Parson Malthus's calculation of population growth as an example. The calculation itself is almost trivial, but the message contained in the output is profound. It is now time to come full circle and see how time has added to Malthus's work. The calculation I want to tell you about is almost as simple as the original one, and like its predecessor, it has had profound consequences, leading to what might well be called a revolution in the way we think about science and its mathematical description. Because the calculation itself is so important, I will deal with it in detail, and luckily, I believe, in this case everyone can follow those details.
12.5.1 Revisiting Malthus
In chapter 1 I introduced the idea of time periods and population values at the end of each period. Thus pn is the population at the end of time period n; it could also be the population of generation n when the growth of nonoverlapping populations is being considered. Malthus was interested in the human population, but the growth of populations of all sorts of animals and plants is of major interest today when we seem to be putting so much pressure on our environment. The Malthusian growth problem can be expressed in the simple formula that says the population in one time period is just the growth rate constant g times the preceding or breeding population:
This equation has a simple solution telling us how the population changes from a give initial value p0:
When g < 1, there is no growth and the population dies out. When g > 1, the population increases without limit. (The special case g = 1 has pn = p0 for all time.) But no population can grow forever, and as Malthus pointed out, the resources for supporting it will come under more and more pressure. Eventually lack of resources—which could mean food, hunting territory, breeding sites, and so on—will cause a decline in the growth of the population. The calculation discussed in this section concerns the mathematical modeling of that process and its implications.
The simplest model of growth limitation may be introduced by defining P to be the maximum population that can be supported by the system under consideration. As the population grows closer to P, the pressures (on food and other resources) will be such that the growth rate g will be reduced. As a simple expression of that we make the change:
Now λ is the new growth rate parameter. When the population is small the growth rate is close to λ, but as the population increases and becomes closer to the possible maximum, the effective growth rate decreases. For example, if pn = ¼P, the effective growth rate is ¾λ; but if the population increases to be ¾P, the growth rate falls to be ¼λ.
With this change to the description of the growth rate, equation (12.12) for the evolution of the population becomes
The second form of the equation (obtained from the first by dividing both sides by P) tells us about the population as a fraction of the possible maximum population P. If we call that fraction x, the new population evolution equation may be written in the standard form found in the literature:
This equation is so important that it has a name: the logistic equation. It is important because it proves to be mathematically interesting and because it also has applications in many areas of biology, social sciences, physics, and engineering. It may be hard to accept, but that simple equation (12.14) reveals a whole new world of population behaviors and challenges some of the most cherished assumptions at the heart of science itself.
12.5.2 Exploring the Logistic Equation
There is no closed form solution for the logistic equation, equation (12.14), as there was for simple Malthusian growth in equation (12.13). The strategy must be to explore the results generated by the logistic equation numerically. That is the “evangelical plea” made by Robert M. May in his seminal 1976 paper “Simple Mathematical Models with Very Complicated Dynamics.” He urged that “people be introduced to…[the logistic equation]…early in their mathematical education. This equation can be studied phenomenologically by iterating it on a calculator or even by hand.”9
When I was a young graduate student in Sydney in the 1960s, Bob May was a brilliant and inspiring associate professor working on plasma physics. He subsequently changed to biology and has become a world-renowned figure in the study of ecology. He has held professorships at Princeton and in England, became the president of the Royal Society, was chief scientific adviser to the UK government, and is now Baron May of Oxford, where he is a professor. Lord Robert May's 1976 paper is a wonderful exposition of the topic of this section written by one of the founders of the theory.
Suppose we begin by taking the growth constant λ in equation (12.14) to be 2, so that for small populations, the growth rate is almost Malthus's population doubling. For example, for the case of an initial population one-tenth of the possible maximum, we set x0 = 0.1 and equation (12.14) gives x1 = 0.18; the population has almost doubled. If we continue to calculate x2, x3, and so on, we discover that the population increases as 0.2952, 0.4161, 0.4859, 0.4996…, but it does not reach the maximum possible value of 1. In fact, the population eventually settles down to the value x = ½. Try substituting xn = ½ in equation (12.14
) with λ = 2 and you will find it returns xn+1 = ½; the population stays the same.
If we start with a different value for x0, the population still finishes up being ½ and then stays at that value. So it is easy to see from the equation that a population of 50 percent of the maximum is what is called a steady-state solution. We can interpret this solution as a balance of population and resources. Try a larger population, xn = 0.6 say, and you will find that the population is pulled down by resources pressures to xn+1 = 0.48. Everything balances for xn = 0.5. I will call this the equilibrium population xeq.
Decreasing the growth constant λ from 2 to 1.5 reduces the equilibrium population to xeq = ⅓ while increasing λ to 2.5 increases the equilibrium population to xeq = 0.6. In fact, putting xn+1 = xn = xeq into equation (12.14) gives a simple formula for the equilibrium population: xeq = (1 – 1/λ).
This is all very satisfying. The logistic equation does indeed model the effect of resource pressures on population growth and tells us that the population and resources are in balance for an equilibrium population which has a simple dependence on the growth rate constant λ. This seems like a good example of mathematical modeling producing numbers that make sense to us for a given physical situation.
12.5.3 A Surprise
Many animal and plant populations “explode,” so we might expect some cases to be modeled by a growth constant λ well beyond 2. We know that there is a general equilibrium value xeq = (1 – 1/λ), and we have seen that once the population gets to that value it remains there. However, there is a surprise in store for when we move to λ > 3: if the population deviates even just a little from xeq it now oscillates around that value while moving further and further away from it. The steady state is now an unstable state. What happens next is quite strange.