Fermi-Pasta-Ulam problem

(1)

Fermi-Pasta-Ulam problem

Tristan Koskamp

Supervisors: J.S. Caux, D. Fioretto, M. Brockmann

ITFA

July 18, 2013

Abstract

The aim of this paper is to get more insight in the behaviour of the Fermi-Past-Ulam problem. We studied the behaviour of the FPU-chain and its thermalization. To obtain the results we wrote a program in C++ to simulate the chain for the α-model. Different initial conditions were used. The results for the α= 0, 25 model are in accordance with the findings of Fermi, Pasta and Ulam. For the α = 1 model we saw that there was chaotic behaviour as predicted with the KAM-theorem. For the model with the high initial energy, we saw a strong tendency towards equipartion of energy.

(2)

1 Introduction

In this paper we try to reproduce the numerical results from Fermi, Pasta and Ulam and investigate the (quasi-)periodicity and equipartition of energy of the Fermi-Pasta-Ulam problem1_{. In order to reproduce their results a code in C}++

was written. In the following sections the origin of the FPU will be discussed. Later the mechanics of a string of oscillators with a little perturbation will be explained and also the principals of equitation of energy were Fermi, Pasta and Ulam based their predictions on. The numerical techniques that are implemented in the C++ code will also be discussed. Examples of these techniques are the Runge-Kutta algorithm and the fast Fourier transform2_{. The goal of this research}

is the get more insight in the FPU-problem and its sometimes weird behaviour.

1_{From now on referred as FPU problem} 2_{From now on referred as FFT}

(4)

2 The Fermi-Pasta-Ulam problem

After the war Enrico Fermi becomes interested in the potential of the electronic computing machines. He thought that computer could contribute to the under-standing of physics problems, especially when problem are not analytically solv-able. In this case a computer model could give insight in the properties and per-haps solutions of the problem. Fermi also thought that future fundamental theories may involve nonlinear equations and operators.

He wanted to study a simple model with non-linear terms and look at its long time behaviour. The model that Fermi started with was a 1-D lattice, you could imagine it as a long chain of oscillators coupled by springs. In this model the springs had also a little perturbation in the form of a quadratic term, so the springs arent linear. Enrico Fermi, John Pasta and Stanislaw Ulam ran numerical simula-tions of this model and looked at the thermalization of this lattice. They expected that the energy of the initial state would slowly drift to other states and eventually the system would reach equipartion of energy.

On one night they left the experiment running by accident. The next morning they found something quite astonishing happens when the simulations runs for a long period of time. As opposed to their expectations equipartion of energy was not reached. They observed a quasi-periodic behavior, first the energy disperse into different models but comes back later in the initial mode for almost about 97%. When the simulation runs longer the energy in the initial comes back in the initial mode with an even higher accuracy, this is known as super recurrence. The results of the simulations were astonishing because there are strong arguments that equipartion also will be reached when there is a nonlinear interaction between particles. This paradoxical observations is known as the FPU problem[4].

(5)

2.1 Mechanics of N oscilators

In the FPU problem we study coupled oscillators. In the normal case when the springs only have linear terms we can use Hook’s law. If we set the constants k and m equal too 1 Hook’s law is

F(x)= ¨x = −x (1)

Fermi, Pasta and Ulam added a little perturbation to Hook’s law to in the form of a nonlinear term, −αx2, and obtained

F(x)= ¨x = −x − αx2 (2)

Figure 1: Chain of N oscillators

The displacement from the equilibrium position, x, of the n’th oscillator will be denoted by un. The equilibrium position of the oscillator is declared with the

dashed line in figure (1). If we now apply formula (2) to oscillator n as shown in figure (1) and only take into account its nearest neighbour left and right, we get the following equation for the unperturbated case after some rewriting.

¨un = un−1− 2un+ un+1 (3)

For the pertubated case we get

¨un= un−1− 2un+ un+1+ α[(un− un−1)2+ (un+1− un)2] (4)

This is the acceleration we will use for the numerical experiment. Fermi, Pasta and Ulam also used a cubic perturbation, βx3_{, but in this paper we won’t discuss}

(6)

2.2 Equipartition of Energy

The equipartition theorem is an important theorem in statistical mechanics and was used by Fermi, Pasta and Ulam to predict the outcome of their experiment. It states that at temperature T , the average energy of any quadratic degree of free-dom is 1₂kT. This means that every enery mode should be equally occupied after a long period of time. The theorem also states that all forms of energy for which the formula is a quadratic function of a coordinate or velocity component, the equipartition theorem holds. Degrees of freedom are for example the translational motion in the x, y and z directions. Other examples of degrees of freedom are vibration motion and elastic potential energy such as energy stored in a spring.

The equipartition theorem applies only to systems whose energy is in the form of quadratic degrees of freedom. This means it needs to be in the following form

E(q)= cq2 (5)

In this equation c is a constant and q is for example any coordinate or momentum variable like x or px. Examples of energy in the form of a quadratic degrees of

freedom are 1₂mv2x, the speed in the x-direction or 1₂kx2, the energy stored in a

spring. When we consider a system with on single degree of freedom and assume that it’s in equilibrium with a reservoir at temperature T . Further let the values for qbe discretely spaced, separated by small intervals∆q and we assume that ∆q is very small. Then the partition function for this system is

Z= X

q

e−βE(q) =X

q

e−βcq2 (6)

If we now multiply by ∆q inside the summation and divide by ∆q outside the summation we get

Z = 1

∆q X

e−βcq2 (7)

Because∆q is very small we can change the summation for an integral

Z = 1

∆q

Z ∞

−∞

e−βcq2dq (8)

We now make the substitution x = √βcq and dq = dx/√βc

Z = 1 ∆q 1 √ βc Z ∞ −∞ e−x2dx = 1 ∆q r π βc (9)

(7)

We now can calculate the average energy E with the following function when we put (9) into it E = −1 Z ∂Z ∂β = 1 2kT (10)

So every quadratic degree of freedom will eventually have an averaged energy of1₂kT. This is known as the equipartition theorem.[7]

In our model we start off with the simple configuration of an sine wave of wave-length N/2 and we want to see how this simple configuration changes over time. Due to the nonlinear forces the configuration will slowly change and the chain of oscillators will get into other states that are also described by sine waves. To make this visual the shape of the chain is analyzed with Fourier series. In our case the end points are fixed, so u0 = uN = 0. So only the complex part of the Fourier

series apply. The idea is to write the energy in terms of Fourier coefficients so we can see in which modes the energy is stored. The Fourier coefficients are

ak = N−1 X n=0 un(t) sin( πkn N ) (11)

The sum of the kinetic and the potential enery for the α-model is

Ekin+ Epot = N−1 X n=0 1 2˙u 2 n+ 1 2(un+1− un) 2₊ α 3(un+1− un) 3 ₍₁₂₎

When we neglect the contribution of the perturbation to the potential, this is only a small amount of the total energy, we can express the energy of the energy stored in every mode with the following equation.[4]

Ek = 1 2˙ak+ 2ak 2_sin2 (πk N) (13)

Fermi, Pasta and Ulam believed that also their system, with a small α, would reach equipartition. So all the energy modes Ek would have the same time averaged

en-ergy. When α is small we can treat the nonlinearity as a small perturbation and due the nonlinearity the different energy modes exchange energy. Also some analyti-cal arguments predicted equipartion of energy for a nonlinear system. But when they actually run the experiment they didn’t encounter equipartion of energy, as told before[2].

(8)

Almost at the same time as Fermi, Pasta and Ulam did their experiment, A.N. Kolmogorov made the first formulation of the KAM theorem. This could explain the results of the FPU experiment. The theorem states that all regular features of the dynamics are kept by integrable Hamiltonian systems subject to a small enough perturbation. This might be a solution to the FPU paradox. Probably the perturbation that Fermi, Pasta and Ulam used was to small in order to get the sys-tem to behave like they predicted and didn’t get equipartion of energy[5].

Another important discovery is the discovery of a stochasticity threshold. When the initial energy is large enough and the energy exceeds the stochasticity thresh-old the FPU-paradox will disappear and there will be equipartiton of energy. This discovery is done by Izrailev and Chirikov. For the explanation of this threshold Izrailev and Chirikov used the KAM theorem[1].

(9)

3 Numerical Techniques

There are a lot of numerical techniques to solve the differential equation stated in equation (4). In this section we will discuss the two techniques that are used to simulate the FPU problem. First the Verlet integration is derived and come with an algorithm that is implemented in the C++ code. Also the more accurate Runge Kutta algorithm will be discussed. This algorithm is an order more accurate than the Verlet integration. At the end the Fast Fourier Transform will be discussed. This is an algorithm to reduce the calculation time of the Fourier transform when you have to deal with large arrays.

3.1 Verlet Integration

One way to simulate the FPU problem is with use of Verlet integeration. This technique is pretty straight forward and has a fourth order error. In this section we will derive the technique step by step and at the end we have an algorithm to simulate the N coupled oscillators. First we start with a simple Taylor expansion of a function un(t) around the point t. Were u denotes the displacement of the

spring form its equilibrium position, t the time and n the n’th oscillator of the chain. If t0is the variable the expression becomes

un(t0)= un(t)+ ˙un(t) 1! (t0− t)+ ¨un(t0) 2! (t0− t) 2₊ ... un(t0) 3! (t0− t) 3_{+ ...} ₍₁₄₎

If we now plug in t+ dt for t0into equation (14) we get the following result

un(t+ dt) = un(t)+ ˙un(t) 1! dt+ ¨un(t) 2! dt 2₊ ... un(t) 3! dt 3_{+ O(dt}4₎ ₍₁₅₎

We also plug in t − dt for t0into expression (14) and then we get

un(t − dt)= un(t) − ˙un(t) 1! dt+ ¨un(t) 2! dt 2₋ ... un(t) 3! dt 3_{+ O(dt}4₎ ₍₁₆₎

The next step is to add the last two equations and we get the following expression after rearranging

un(t+ dt) = 2un(t) − un(t − dt)+ ¨un(t)dt2+ O(dt4) (17)

We now have an expression to calculate the displacement of the oscillators at time t + dt and we only have to know the deflection at time t and time t − dt. Also the acceleration at time t needs to be known, but this only depends on the displacement at time t, see equation (4). To implement this algorithm into the FPU simulation we assume that at t = 0 the acceleration is also zero (¨un(t) = 0).

This means that un(t) = un(t+ dt) at the start. With this initial condition we can

(10)

3.2 Runge-Kutta

Besides the method of Verlet integration you have a lot more methods to evalu-ate the problem numerical. One of these methods is Runge-Kutta method. This method is al little more complex to implement but is one order more accurate than the Verlet integration method. The standard Runge-Kutta method is simple to im-plement for first order differential equations but it also useful for second order differential equations although a little bit more work to implement. In the next section it is discussed how to the Runge-Kutta method is implemented for the FPU problem.

The standard Runge-Kutta method for first order differential equations is as fol-low. If we have the first order differential equations and it’s initial conditions defined

∂x

∂t = f(t, x) x(t0)= x0 (18)

We can calculate the next value for x, x(t+ dt), with the Runge-Kutta algorithm. In order to do this we need to calculate four variable - k1, k2, k3 and k4. They are

defined as follow k1 = f (t, x) (19) k2 = f (t + dt 2, x + k1 2) (20) k3 = f (t + dt 2, x + k2 2) (21) k4 = f (t + dt, x + k3) (22)

When all these k’s are calculated you can calculate x(t + dt) with the following function:

x(t+ dt) = x(t) + dt

2(k1+ 2k2+ 2k3+ k4) (23)

In our case we have to modify the algorithm a bit, because we are dealing with second order differential equations. The formula for the acceleration is

∂2_u n

∂2_t = un−1− 2un+ un+1+ α[(un− un−1) 2_{+ (u}

(11)

Now we are going to make a substitution ∂un ∂t =vn (25) ∂vn ∂t = F(t, u(t)) (26) ∂vn ∂t =un−1− 2un+ un+1+ α[(un− un−1)2+ (un+1− un)2] (27) (28)

With this we can use the Runge-Kutta method and apply it to second order dif-ferential equations. In our case an array of length N for j1, j2, j3 and j4 needs

to be created. This is due the fact that the acceleration depends not only on the displacement of the oscillator at place n but also on its neighbours. Also an array of length N for the variables k1, k2, k3and k4is needed. In our case j1, j2, j3and j4

are needed for the calculations of the velocity and the k1, k2, k3and k4are needed

for the calculations of the displacement. So the next thing to do is to build the following array for the Runge-Kutta algorithm

k10 = ˙u0 k11 = ˙u1 · · · k1N−2 = ˙uN−2 k1N−1 = ˙uN−1 j10 = ü0 j11 = ü1 · · · j1N−2 = üN−2 j1N−1 = üN−1 (29)

Now we need to recalculate all the displacement un of the oscillators so we can

calculate j2. This is because the acceleration of oscillator n also depends on the

(12)

the following way u20 = u1+ dt 2k10 u21 = u1+ dt 2k11 · · · u2N−2 = u1+ dt 2k1N−2 u2N−1 = u1+ dt 2k1N−1 (30)

The next step is to calculate al the values for k2 and j2. Because neither j1 nor k1

depends explicited on t the next step will be

k20 = ˙u0+ dt 2 j10 k21 = ˙u1+ dt 2 j11 · · · k2N−2 = ˙uN−2+ dt 2 j1N−2 k2N−1 = ˙uN−1+ dt 2 j1N−1 (31) ü2n = u2n−1 − 2u2n + u2n+1 + α[(u2n − u2n−1) 2_{+ (u} 2n+1 − u2n) 2_] ₍₃₂₎ j20 = ü20 j21 = ü21 · · · j2N−2 = ü2N−2 j2N−1 = ü2N−1 (33)

After we calculated all the values of k2 and j2 we can recalculated the

displace-ments u3n. With these displacements it’s possible to calculate all the values of j3

and with the values of k2 it’s possible to calculate all the values of k3. This

pro-cess is repeated till k4 and j4 are calculated. Once these values are obtained it’s

possible to calculate un(t+ dt) and vn(t+ dt) with the following formulas

un(t+ dt) = un(t)+ dt 6(k1n + 2k2n + 2k3n + k4n) (34) vn(t+ dt) = vn(t)+ dt 6( j1n + 2 j2n + 2 j3n + j4n) (35)

(13)

3.3 Fast Fourier Transform

The FFT is a great algorithm to reduce arithmetic operations with respect to the discrete Fourier transform3_{. It reduces the number of calculations from an order}

of N2 to an order of N log₂(N). For large N this saves tons of time. The only downside of this algorithm is that we have to use an array that has a length N that obeys the equations N = 2L _{were L} = 2, 3, 4, .... So N has to be a power of two.

The basic idea of the algorithm is to split the N-point DFT up into two (N/2)-point DFTs. This can be done a total of L times. Now we are going to discuss the exact way how the FFT works.

The DFT coefficient X(k) is defined by X(k)= 1 N N−1 X n=0 x(n)(e−i2π/N)kn = 1 N N−1 X n=0 x(n)Wkn (36)

Where k = 0, 1, 2, ...., N−1 and x(n) are the data points. If we have a data sequence that has in total N = 2L_{points. Then its possible to apply some manipulations to}

equation (36). Since N is a power of 2, we know that N/2 is an integer and we can combine samples that are separated by N/2 in the data sequence. So we get

X(k)= 1 N N/2−1 X n=0 [x(n)Wkn+ x(n + N/2)Wkn+N/2] (37) X(k)= 1 N N/2−1 X n=0 [x(n)+ x(n + N/2)WN/2]Wkn (38)

This result can be simplified because WkN/2can only take two values. Namely −1 and 1. This is because of the following

WkN/2 = exp(−i2π N kN 2 )= e −iπk _{= (−1)}k (39)

When k is even we see WkN/2 = 1. Let define the following functions

k = 2l l= 0, 1, 2, ..., N/2 − 1 (40)

g(n)= x(n) + x(n + N/2) (41)

(14)

If we use this to rewrite (37) we get X(2l)= 1 N N/2−1 X n=0 g(n)W2ln = 1 2 1 N/2 N/2−1 X n=0 g(n)(W2)ln (42)

When we look more closely to equation (42), we see that this is one-half times an (N/2)-point DFT because W2 = exp[−i2π/(N/2)] and the input is {g(0), g(1), g(2),...., g(N/2 − 1)}. We can conclude that for even values of k we can reduce the DFT in an (N/2)-point DFT.

Now we are going to look at the case when k is odd, so WkN/2 = −1. In this case we define

k= 2l + 1 l= 0, 1, 2, ..., N/2 − 1 (43)

y(n)= x(n) − x(n + N/2) (44)

If we use this to rewrite (37) we get

X(2l+ 1) = 1 N N/2−1 X n=0 y(n)W(2l+1)n= 1 2 1 N/2 N/2−1 X n=0 h(n)(W2)ln (45)

In this case h(n) = y(n)Wn_{. We can see now just as in the case of k is even, that}

equation (45) is one-half times an (N/2)-point DFT with the input {h(0), h(1), h(2),...., h(N/2 − 1)}. The parameter Wn _{in y(n)W}n _{is often referred as the twillde}

factor.

As an example we apply this method discussed above for an 8-point DFT. So in this case N = 8, we have eight inputs namely {x(0), x(1), x(2),...., x(7)}. To reduces the input by a factor two we need to add x(0) and x(4), x(1) and x(5), x(2) and x(6) and x(3) and x(7). For the odd values of k we need to subtract x(4) from x(0), x(5) from x(1), x(6) from x(2), x(7) from x(1) and multiply by the right twid-dlefactor W0, W1, W2 and W3. The inputs and outputs of the first stage described above are shown in table (1) and in figure (2) we can see how the outputs are build up. We also split the 4-point DFT’s up into four 2-point DFT. How this is done is shown in table (2) and figure (3). The last step is to calculate the 2-point DFT’s outputs, the inputs and outputs are shown in table (3) and figure (4). In the tables the normalization factor of 1/2 for every stage is left away in order to get a more organized table. At the end all the outputs in table (3) need to be multiplied with 1/8.

(15)

Figure 2: Reduction of an 8-point DFT to two 4-point DFT’s Input Output x(0) x(0)+ x(4) x(1) x(1)+ x(5) x(2) x(2)+ x(6) x(3) x(3)+ x(7) x(4) x(0) − x(4)W0 x(5) x(1) − x(5)W1 x(6) x(2) − x(6)W2 x(7) x(3) − x(7)W3

(16)

Figure 3: Reduction of an 8-point DFt to two 4-point DFT’s Input Output x(0)+ x(4) x(0)+ x(4) + x(2) + x(6) x(1)+ x(5) x(1)+ x(5) + x(3) + x(7) x(2)+ x(6) x(0)+ x(4) − [x(2) + x(6)] x(3)+ x(7) x(1)+ x(5) − [x(3) + x(7)] x(0) − x(4)W0 x(0) − x(4)W0+ [x(2) − x(6)W2] x(1) − x(5)W1 x(1) − x(5)W1+ [x(3) − x(7)W3] x(2) − x(6)W2 _{x(0) − x(4)W}0_{− [x(2) − x(6)W}2_] x(3) − x(7)W3 x(1) − x(5)W1− [x(3) − x(7)W3]W2 Table 2: Two 4-point DFT’s to four 2-point DFT’s

(17)

Input Output x(0)+ x(4) + x(2) + x(6) = g(0) g(0)+ g(1) = X(0) x(1)+ x(5) + x(3) + x(7) = g(1) g(0) − g(1)= X(4) x(0)+ x(4) − [x(2) + x(6)] = g(2) g(2)+ g(3) = X(2) x(1)+ x(5) − [x(3) + x(7)] = g(3) g(2) − g(3)= X(6) x(0) − x(4)W0+ [x(2) − x(6)W2_]= g(4) g(4) + g(5) = X(1) x(1) − x(5)W1+ [x(3) − x(7)W3]= g(5) g(4) − g(5) = X(5) x(0) − x(4)W0− [x(2) − x(6)W2]= g(6) g(6) + g(7) = X(3) x(1) − x(5)W1− [x(3) − x(7)W3_]_{= g(7) g(6) − g(7) = X(7)}

Table 3: Four 2-point DFT’s to final output

The FFT for an 8-point DFT is now ready. It’s pretty straight forward to extend this algorithm to larger N, if N is a power of two. The only problem we have now is that the output is in bit-reversed order (and we didn’t normalize the output by a factor 1/8). In order to get the output in the right order we need to bit-reverse the output. For the 8-point DFT this is done in table (4). The way to do it is to write your output in bit notations, mirror it and write it back to decimal notation.

Bit-reversed index Bit-reversed Bits Linear index

0 000 000 0 4 100 001 1 2 010 010 2 6 110 011 3 1 001 100 4 5 101 101 5 3 011 110 6 7 111 111 7

Table 4: Bit reversing

Now we are done with the FFT algorithm[3]. The C++ code for the FFT is taken from Numerical Recipes: The Art of Scientific Computing[6]. There is only one problem left, in our case the end points are fixed (u0 = uN = 0), so there are no

cosines. In this case we want to use the fast sine transform. It is a little more complex to implement then the FFT. We are now going to discuss the theoretical background of this algorithm.

(18)

We want to calculate the sine transform and it is given by F(k)= N−1 X n=0 x(n) sin(πnk N ) (46)

Where x(0), x(1),..., x(N − 1) are the data points and x(0)= 0. At the first glance it looks like the complex part of the FFT, but the term inside the sine differs by a factor of two in comparison with the FFT.

We can use the FFT algorithm to calculate the sine transform, but we need to modify the input first. The idea is to expand the data set into an odd function with

x(N)= 0 x(2N − n)= −x(n) n= 0, 1, ..., N − 1 (47)

The FFT of this doubled data set is

F(k)=

2N−1

X

n=0

x(n)e2πink/(2N) (48)

Now we make an substitution to rewrite half of this sum. This substitution is n0 = 2N − n. We can rewrite the summation of equation (48)

2N−1 X n=N x(n)e2πink/(2N) = N X n0₌₁ x(N − n0)e2πi(2N−n0)k/(2N) (49) = − N−1 X n0₌₀ x(n0)e−2πn0k/(2N) (50) Now rewrite equation (48) with equation (50) and we obtain

F(k) =

N−1

X

n=0

x(n)[e2πink/(2N)− e−2πink/(2N)] (51)

= 2i N−1 X n=0 x(n) sin(πnk N ) (52)

Except for the factor 2i we get the sine transform from the FFT for the extended function. The only problem is that we need to double the array length which reduces the efficiency of the program. With another trick we can correct this. We need to create an extra array yjthat has the following entries

y(0)= 0 (53)

y(n)= sin(nπ/N)[x(n) + x(N − j)] +1

(19)

This array has also length N, n is ranging from 1 till N − 1. If we use the FFT on this array we will get the following transform for the real part of the transform (R(k)) and the imaginary part of the transform (I(k))

R(k)= N−1 X n=0 y(x) cos(2πnk/N) (55) = N−1 X n=1 [x(n)+ x(N − n)] sin(nπ/N) cos(2πnk/N) (56) = N−1 X n=0 2x(n) sin(nπ/N) cos(2πnk/N) (57) = N−1 X n=0 x(n)[sin(2k+ 1)nπ N − sin (2k − 1)nπ N ] (58) = F(2k + 1) − F(2k − 1) (59) I(k)= N−1 X n=0 x(x) sin(2πnk/N) (60) = N−1 X n=1 [x(n) − x(N − n)]1 2sin(2πnk/N) (61) = F(2k) (62)

From these equations we can determine F(k)

F(2k)= I(k) F(2k+ 1) = F(2k − 1) + R(k) k= 0, ..., (N/2 − 1) (63)

We see that the even terms are determined directly, but for the odd terms we also need another term. For k= 0 we use the fact that F(1) = F(−1) and then F(1) can be calculated, in this case F(1) = 1₂R(0). Now we can caluclate all the values of F(k).

In the book Numerical Recipes they use a clever trick to store a real array into a complex array so the FFT only need to calculate a complex array of length N/2 instead of an complex array of length N with all zeroes for the complex parts. We won’t further discuss the detail of this clever way of calculating the FFT, but the resulting code is used in the program to simulate the FPU-problem.[6]

(20)

4 Results

In the next section the result of different intial conditions are analyzed. There are three initial conditions that are used to generate the plots. The first one is the famous sine wave of wavelength N/2, amplitude of 1 and α = 0, 25 were Fermi, Pasta and Ulam discovered that the energy would go back into the original mode after some time. The second initial condition is an sine wave of wavelength N/2, amplitude of 1 but α = 1 . The final inital condition is an sine wave of wavelength N/2 with an amplitude of 8 and and α = 0, 25. The main goal is to look at the recurrence of the of the initial state and look if there is any form of equipartion of energy. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 1 2 3 4 5 6 7 8 9 10 Ener gy (arb . unit) Time (×103_sec) Energymode 1 Energymode 2 Energymode 3 Energymode 4 Energymode 5 Energymode 6

Figure 5: Wavelength N/2, amplitude of 1, N = 32, α = 0.25 and dt = 0.01 In figure (5) the first six energy modes are plotted for the initial condition of a sine wave of wavelength N/2. We see that all the energy is in the first mode at t = 0 because this is the initial conditions. We see that after some time approximately t = 2 × 103 all the first six modes picked up energy and it’s looking like all the enrgy modes are coming closer and closer and the equipartion theorem holds. But then at t = 4.3×103_{almost all the energy is in the second mode and at t}= 9.4×103

(21)

In figure (6) we see how the chain of oscillators looks like at different times. The red dots show how the chain looks like at the initial condition. At t = 2000 we see in figure (5) that the chain is a superposition of several Fourier modes. For t = 4364 the only mode occupied is the second mode and at t = 8650 all the energy almost returned in the first mode. All the findings in figure (6) are in accordance with figure (5).

-1.5 -1 -0.5 0 0.5 1 1.5 5 10 15 20 25 30 Displacement (arb . unit) Oscillator n t = 0 t = 2000 t = 4364 t = 8650

(22)

The next thing to look at is if the recurrence in figure (5) is coming back multiple times. In figure (7) and figure (8) we plotted the energymodes for the same initial conditions only with a longer timescale. If we look at figure (7) we see that the energy comes back into the first mode, but each time less energy goes back into the first mode. It seems that equipartition of energy will be reached eventually. But when we take a glance at figure (8) we see that there is another recurrence. Eventually all the energy is back in the first mode and this recurrence repeats many times. This recurrence is known as super recurrence.

0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 1 2 3 4 5 6 7 8 9 10 Ener gy (arb . unit) Time (×104_sec) Energymode 1 Energymode 2 Energymode 3 Energymode 4 Energymode 5 Energymode 6

Figure 7: Wavelength N/2, amplitude of 1, N = 32, α = 0.25 and dt = 0.01

(23)

Let’s take a look at the time averaged energy for every mode. For the first twenty modes the time averaged energy is plotted in figure (9). The time averaged energy of the energy modes first diverge to each other, but at some point it stops and equipartion of energy is not reached. In figure (10) all the energy modes are plotted at the x-axis and the time averaged energy at the y-axis at t = 106. It is clear that equipartition is not reached for the model with 32 oscillators, α = 0, 25 and initial condition of a sine wave with wavelength N/2.

10−30 10−25 10−20 10−15 10−10 10−5 100 105 101 ₁₀2 ₁₀3 ₁₀4 ₁₀5 ₁₀6 T ime av eraged ener gy (arb . unit) Time (sec)

10−14 10−12 10−10 10−8 10−6 10−4 10−2 100 0 5 10 15 20 25 30 T ime av eraged ener gy (arb . unit) Ek

(24)

When α is set equal to 1 we see that the energy also comes back into the first mode as shown in figure (11). Also a second recursion is observed in figure (11) although this looks already more chaotic then for α equal to 0,25. When we plot the energy modes for a longer period of time, there isn’t such a nice pattern as for the α = 0, 25 model. Although there is a little bit of recurrence the pattern is also for longer timescales more chaotic and this might indicate that the time averaged energy for each mode lies closer together.

Figure 11: Wavelength N/2, amplitude of 1, N= 32, α = 1 and dt = 0.01

(25)

When we look at figure (13) we see indeed that the time averaged energies lay closer together if we compare it with the model were α is equal to 0,25. In figure (14) the time averaged energy for each mode is plotted at t = 105_{seconds. In this}

figure it is clear that the time averaged energies for all the energies modes lie closer together. The perturbation is probably too big to obey the KAM-theorem.

10−30 10−25 10−20 10−15 10−10 10−5 100 105 101 ₁₀2 ₁₀3 ₁₀4 ₁₀5 T ime av eraged ener gy (arb . unit) Time (sec)

Figure 13: Wavelength N/2, amplitude of 1, N= 32, α = 1 and dt = 0.01

10−4 10−3 10−2 10−1 100 0 5 10 15 20 25 30 T ime av eraged ener gy (arb . unit) Ek

(26)

For the last couple of plots the initial condition was again a sine wave with wave-length N/2, but this time the amplitude was 8 instead of 1 and α was set back to 0,25. If we look at figure (15) it is clear that the pattern is more chaotic then in the previous plots. When we look at (16) we see that there isn’t really recurrence. The energy is not going back in the first mode for more than 65%. This might indicate that there is equipartition after a long time.

0 10 20 30 40 50 60 70 80 0 1 2 3 4 5 6 7 8 9 10 Ener gy (arb . unit) Time (×103_sec) Energymode 1 Energymode 2 Energymode 3 Energymode 4 Energymode 5 Energymode 6

0 10 20 30 40 50 60 70 80 0 1 2 3 4 5 6 7 8 9 10 Ener gy (arb . unit) Time (×104_sec) Energymode 1 Energymode 2 Energymode 3 Energymode 4 Energymode 5 Energymode 6

(27)

In figure (17) and (18) we see that the time averaged energy of the energy modes lies closer together than the for the α= 0.25 and α = 1 models. This is something you might expect after seeing figure (15) and (16). The stochasticity threshold is exceeded because the initial energy was big enough[1] and therefore a tendency towards equipartion is reached.

10−30 10−25 10−20 10−15 10−10 10−5 100 105 101 ₁₀2 ₁₀3 ₁₀4 ₁₀5 T ime av eraged Ener gy (arb . unit) Time (sec)

100 101 0 5 10 15 20 25 30 T ime av eraged ener gy (arb . unit) Ek

(28)

10−14 10−12 10−10 10−8 10−6 10−4 10−2 100 102 0 5 10 15 20 25 30 T ime av eraged ener gy (arb . unit) Ek A= 1, α = 0.25 A= 1, α = 1 A= 8, α = 0.25

Figure 19: The time averaged energy for each energymode compared with the different initail conditions

In figure (19) all the time averaged energy of the different modes from the three models are plotted. We can clearly see that the model with the big initial amplitude is the closest to equipartiton of energy. The model with amplitude 1 and α equal to 0.25 is far from equipartiton. Also the model with amplitude 1 and α equal to 1 doesn’t really shows equipartiton but is defenitly more closer to it. The model with the most tendecy towerads equipartion of energy is the model with the big amplitude. The stochasticity threshold is exceeded in this model and equipartition of energy will be reached eventually.

(29)

5 Conclusion

The results for the α = 0.25 model are in accordance with the results that are obtained by Fermi, Pasta and Ulam. The energy is of the first model slowly flows away to other modes and comes back later in the first mode. For a long timescale (106_{sec) the pattern looks like a periodic pattern and doesnt look like chaos at all.}

If we look at the time averaged energy of the different modes we see that there is no equipartiton of energy and that the energy is mainly stored in the first few modes. This is also in accordance with the observations of Fermi, Pasta and ulam.

For the α = 1 model we see already for a rather small timescale less recurrence and a less orderd pattern in comparison with the α = 0.25 model. It seems that it is much more chaotic. When we look at the time averaged energy per mode we see that they lie closer together. But there isnt really equipartion of energy. In this model α is probably big enough to violate the KAM-theorem and thats why we see more chaotic behaviour.

The last model, with high initial energy, shows the most chaotic behaviour. For a relative short timescale of 104 seconds the behaviour seems to be completely chaotic. When plotted for a timescale of 105 _{we dont see any pattern. If we}

look at the time averaged energy we see that the energy modes dont differ more than a factor 10 and this model shows the most tendency to equipartition. The stochasticity threshold is exceeded and it looks like equipartition of energy will eventually be reached. This is in accordance with the observation of Izrailev and Chirikov.

(30)

References

[1] G. Benettin, A. Carati, L. Galgani, A. Giorgilli, The Fermi-Pasta-Ulam Prob-lem and the Metastability Perspective, Lecture Notes in Physics, 728, pp. 151-189, (2008)

[2] G. P. Berman, F. M. Izrailev, The Fermi-Pasta-Ulam problem: 50 years of progress, Choas 15, pp. 01504, (2005)

[3] D. F. Elliot, K. R. Raoi, FAST TRRANSFORMS: Algorithms, Analyses, Appli-cations(Academic press, 1982), pp. 58-98

[4] E. Fermi, J. Pasta, S. ” Ulam, Studies of Nonlinear Problems”, (Tech. Rep. LA-1940, Los Alamos National Laboratory, Los Almos, New Mexico, 1955) [5] R. Livi, A. Vulpiani, The Kolmogorov Legacy in Physics (Springer, 2003), pp

17-18

[6] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, Third Edition in C++ (Cambridge university press, 2007), pp. 600-623

[7] D. V. Schroeder, An Introduction to Thermal Physics (Addison Wesley Long-man, 1991), pp. 14-16, 238-240

Fermi-Pasta-Ulam problem