The+Fermi+-+Pasta+-+Ulam+-+Tsingou+chain

= =
 * DM2 The Fermi - Pasta - Ulam - Tsingou chain **

= Introduction = Consider a system which evolves deterministically (there is no noise) and observe it for a long time. You can //e.g.// think of hard spheres on a billiard, or of particles interacting on a lattice, or more generally of a system described by a Hamiltonian dynamics math \dot p_n = - \frac{\partial\mathcal H}{\partial q_n}\quad math math \dot q_n = \frac{\partial\mathcal H}{\partial p_n} \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad(1) math
 * Don't hesitate to ask questions and make remarks on this wiki page.**

=
(Here, // q n // (//t//) represents a position and //p// //n // (//t//) a momentum). The initial condition is // q n // (0), //p// //n // (0) //.// The total energy i s conserved: it is constant in time: ====== math E=\mathcal H(q_n(0),p_n(0)) \,. math One may hope that in the long time and large system size limits, one recovers the results of statistical physics: in other words, we expect that the long-time average of an observable coincides with the ensemble average of statistical physics. In this case, energy is conserved, the ensemble is microcanonical. For this ensemble we expect equipartition: the probability density for the system to have visited a point ( // q n //, //p// //n // ) in the phase space is uniform on the energy shell //E :// math E=\mathcal H(q_n,p_n) math One prerequisite for this property to hold is that the system has to explore the whole energy surface. There are simple counter-examples to this expectation: one is given by a chain of **harmonically** coupled particles on a lattice (//i.e.// particles connected to their neighbors with elastic springs). Roughly speaking, the number of constants of motion in this system is high enough to forbid the dynamics to explore the whole energy surface. = The harmonic chain = = Fourier representation = We start by studying a discrete elastic chain: L+1 particles are interacting on a one dimensional lattice (see figure). The particles interact with their neighbors at left and at right with a **harmonic** potential. The distance between the rest position of the particle //n// and its actual position is denoted by //u//// n // (//n//=0...//L//). We consider the case with fixed boundary conditions and with the particles initially at rest: math u_0=u_L=0, \quad \dot u_0=\dot u_L=0, \quad v_n(0)=0 math The Hamiltonian is math \mathcal H(\{u_n,v_n\}) = \frac 12 \sum_{n=0}^{L-1} \Big[ v_n^2 + (u_{n+1}-u_n)^2\Big] math where the velocity is the moment conjugated to the position (taking unit mass). Using equation (1), show that the equations of motion are: math \ddot u_n = u_{n+1}+u_{n-1}-2u_n \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad(2) math The problem is easier to study in the Fourier space: one defines the Fourier modes //k// = 0,1...//L// − 1 math \hat u_k = \sqrt{\tfrac 2L}\sum_{n=0}^{L-1} u_n \sin\frac{nk\pi}{L}\quad ; \quad\quad\quad \hat v_k = \sqrt{\tfrac 2L}\sum_{n=0}^{L-1} v_n \sin\frac{nk\pi}{L} math which are compatible with the boundary conditions. The inverse Fourier transform is given by math u_n = \sqrt{\tfrac 2L}\sum_{k=0}^{L-1} \hat u_k \sin\frac{nk\pi}{L}\quad ; \quad\quad\quad v_n = \sqrt{\tfrac 2L}\sum_{k=0}^{L-1} \hat v_k \sin\frac{nk\pi}{L} math The corresponding expression of the Hamiltonian is math \mathcal H(\{\hat u_k,\hat v_k\}) = \sum_{k=0}^{L-1} E_k = \sum_{k=0}^{L-1} \frac 12 \Big[ \hat v_k^2 + \omega_k^2 \hat u_k^2\Big] \quad\text{where}\quad \omega_k=2\sin\frac{k\pi}{2L} math
 * Count the modes of strictly positive energy for the system described above.
 * At // t // =0, we prepare the system at energy // E // . If the equipartition of energy applies, what is the average energy associated to each mode k?

Evolution in time
At t=0, we prepare the system in the Fourier mode k=1: math u_n(t=0) ~=~ L\sin\frac{n\pi}{L} \quad\quad v_n(t=0) ~=~ 0 math

and the energy is concentrated in the first mode, of energy math E(t=0)=E_1=L^3 \sin^2 \frac{\pi}{2L} math We may wonder whether the energy stays in this mode or decays on the others modes, as predicted by statistical physics. We now study this question numerically, by discretizing the equation of evolution (2) with a time step dt. code format="python" import pylab,math,numpy Ntime = 1000 L = 32 dt = 0.001 def HARMONIC(u,n): return u[n+1]+u[n-1]-2.*u[n] def FT(e,k): return sum([e[n]*math.sin(n*k*math.pi/L) for n in range(1,L)])/math.sqrt(L/2.) u = [0]+[L*math.sin(n*math.pi/L) for n in range(1,L)]+[0] v = [0. for n in range(L+1)] mode1 = [] for itime in range(Ntime): oldu = u[:] oldv = v[:] v = [0]+[oldv[n]+dt*HARMONIC(oldu,n) for n in range(1,L)]+[0] u = [0]+[oldu[n]+dt*v[n] for n in range(1,L)]+[0] k = 1; en1 = FT(v,k)**2/2.+2.*(math.sin(k*math.pi/(2*L)))**2*FT(u,k)**2 #k = 2; en2 = FT(v,k)**2/2.+2.*(math.sin(k*math.pi/(2*L)))**2*FT(u,k)**2 mode1.append(en1) pylab.title('The harmonic chain') pylab.xlabel('time') pylab.ylabel('E_1, energy of mode 1') time = [dt*itime for itime in range(Ntime)] exact = [L*(L*math.sin(math.pi/(2*L)))**2 for itime in range(Ntime)] pylab.plot(time,exact,'b-') pylab.plot(time,mode1,'r-') pylab.axis([0,Ntime*dt,75.,82.]) pylab.show code
 * 1) initial condition for the position u and velocity v
 * 1) dynamics with fixed boundary conditions: u[0] = u[L-1] = 0

Note that one could use the [|Fast Fourier Transform] algorithm to compute The Fourier Transform. Moreover, we could also use more efficient numerical integration such as the [|Leapfrog method], instead of the modified [|Euler] one.
 * Check the precision of the algorithm when decreasing the time step dt.
 * Try the initial condition u=[0]+[L/2.-abs(L/2.-n) for n in range(1,L)]+[0] . Plot some snapshots of the displacement //u//(//t//) . Plot the time evolution of the energy //E// // k // for modes //k//=1 to 4. Comment.
 * Questions on Python language:
 * 1) What is the difference between a list (e.g. [1,2,3]), a tuple (e.g. (1,2,4))?
 * 2) Define a list a (e.g. a=[1,2,2]). Define b=a. Modify a. What happens to b?
 * 3) If a and b are two lists, what is the difference between the commands a.append(b), a.extend(b) and a+b?

More generically, //integrable systems// also break ergodicity. It was thought that non-integrable systems would on the contrary be ergodic, until... = A counter-example by Fermi, Pasta, Ulam and Tsingou =

System and ergodicity
In 1955, Enrico Fermi, John Pasta, Stanislaw Ulam and Mary Tsingou proposed to add a slightly **non-harmonic** interaction, that would allow the mixing of the energy between the modes so as to recover the microcanonical prediction of statistical physics. The equations of motions that they considered is math \ddot u_n = u_{n+1}+u_{n-1}-2u_n + \alpha \Big[(u_{n+1}-u_n)^2-(u_{n}-u_{n-1})^2\Big] \quad \text{with} \; \alpha \ll 1 math The numerical study of this equation was performed with one of the very first computers, the [|Maniac I]. To repeat the same experiment, you can use the following force code format="python" def FPU(u,n): return u[n+1]+u[n-1]-2.*u[n]+alpha*((u[n+1]-u[n])**2-(u[n]-u[n-1])**2) code instead of HARMONIC(u,n). We start with the following parameters: code format="python" Ntime = 10000 L = 32 alpha = 0.3 dt = 0.2 code and this initial state code format="python" u = [0]+[math.sin(n*math.pi/L) for n in range(1,L)]+[0] code
 * 1) What do you observe? Explain why Fermi, Pasta, Ulam and Tsingou were happy.
 * 2) Take a larger number of time steps Ntime=44000. What do you really conclude about ergodicity?
 * 3) Compare different values of α.
 * 4) Plot the shape of the profile //u// //n// at different times on a same graph, slightly shifted so as to visualize the evolution.

Study of solitons (bonus; difficult)
This complicated evolution corresponds to the propagation of many undeforming profiles, called "solitons". To determine the possible shape of a single soliton, one first takes a macroscopic limit by defining continuous variables math x=\varepsilon n \quad \text{and} \quad\tau=\varepsilon, \quad \text{with} \varepsilon \ll 1 math

math u_n(t)=\varepsilon \phi(\varepsilon n, \epsilon t) math Show that at second order the field φ(//x//,τ) verifies the equation math \partial_\tau^2\phi = \partial_x^2\phi+\bigg[2\alpha\partial_x\phi\partial_x^2\phi +\frac 1{12} \partial_x^4\phi\bigg]\varepsilon^2 math math \partial_\tau^2 \phi math and show that math \partial_\tau\Phi +\bigg[\alpha\Phi\partial_X\Phi +\frac 1{24} \partial_X^3\phi\bigg]\varepsilon^2=0 \quad \text{where} \Phi = \partial_x \phi math math \Phi(X,t)= C_0 \Phi_0\big(C_1(X-v\tau)\big) math where the invariant shape is math \Phi_0(\bar X)=\frac 1{\cosh^2\bar X} math
 * One sets a field φ(//x//,τ) from the definition
 * To study a travelling wave, one sets // X = x − τ //. Find the equation for φ(//X//,τ).
 * Assuming that variations in time are slow, neglect the term
 * This equation is called the Korteweg-de Vries equation. Find a traveling wave solution, //i.e.// of the following form:
 * Come back to the initial variable φ(//x//,τ) (Hint: when integrating from Φ to φ, ensure that the profile φ goes to zero at infinity).
 * Test numerically the stability of the obtained soliton. Don't forget to compute the initial velocity correctly.
 * How do two colliding solitons of opposite velocities interact? Illustrate this fact numerically.

= References =
 * The original article by Fermi, Pasta and Ulam.
 * Scholarpedia article on the Fermi-Pasta-Ulam-Tsingou problem.
 * Mathematical aspects of the Fermi-Pasta-Ulam-Tsingou problem.
 * An historical perspective.

url}?f=print|print this page