In 1953-54 Enrico Fermi, John Pasta and Stanislaw Ulam used the MANIAC I computer of the Los Alamos laboratory for one of the first numerical experiments in history: the simulation of a chain of oscillators with a weak non-linear coupling. The programs were written and run by Mary Tsingou, now recognised as a co-author: for this reason one also speaks of the FPUT problem^[1].

The original paper discussing the problem examined here, which first appeared in report LA-1940 (May 1955) of the Los Alamos National Laboratory, was published after Fermi's death in the collection of his works.

Normal modes in a non-linear system

• As illustrated in more detail on the dedicated page, the normal modes of a linear system are always independent: when one mode is excited it does not exchange energy with the other normal modes, and the motion of the system remains stationary in that mode.
• However, when the system is non-linear, the superposition principle no longer applies, and the normal modes can interact, exchanging energy. Then, even if we could initially excite a single mode, the energy given to the system could transfer to other modes over time. The motion of the system is no longer stationary.

The Fermi, Pasta and Ulam problem

The question, then, was: in the presence of a weak non-linear coupling, which modes will be excited after a long enough time?

The answer Fermi, Pasta and Ulam expected was based on the equipartition theorem: energy is shared equally among all the degrees of freedom of the system.

• The simplest physical system on which to test this conjecture is the chain of oscillators, into which, however, a weak non-linear coupling is introduced. This problem does not admit a simple analytical solution like its linear counterpart, and only a numerical simulation can reveal its dynamics.

• You can try the simulation yourself with the program at the bottom of this page. A program similar to this one ran in 1953-54 on the MANIAC I computer of the Los Alamos laboratory, simulating chains of 16, 32 and 64 oscillators, and revealed to Fermi, Pasta and Ulam that their conjecture was wrong.

• Letting the virtual system evolve in time, one first observes a redistribution of energy: initially exciting mode 1, modes 2, 3 and 4 become active in turn. However, if a long enough time is allowed to pass (157 periods of mode 1), the energy returns almost completely to mode 1, and there is a periodicity in the way the different modes are excited.

• Fermi, Pasta and Ulam had thus shown that a small non-linear coupling is not enough to obtain equipartition of energy.

Check your understanding

Simulation program

The following Matlab program is published in Thierry Dauxois, Michel Peyrard, Stefano Ruffo, The Fermi-Pasta-Ulam "numerical experiment": history and pedagogical perspectives, European Journal of Physics 26, S3-S11 (2005), and reproduced here courtesy of Stefano Ruffo.

clear
N=32;		% Number of particles must be a power of 2
TMAX=10000;	% tmax   10000
DT=20;  	% Delta t
alpha=0.5;	% Nonlinear parameter
a=1;		% Amplitude of the initial condition
tspan=[0:DT:TMAX];
options=odeset('Reltol',1e-4,'OutputFcn','odeplot','OutputSel',[1,2,N]);
% Test different tolerances, changing Reltol

for I=1:N,
b(I)=a*sin(pi*I/(N+1));    b(I+N)=0;    % FPU initial condition
%b(I)=a*sin(pi*N*I/(N+1)); b(I+N)=0;    % Zabusky-Deem initial condition
%k=0.4;
%b(I)=-0.5/alpha*log((1+exp(2*k*(I-1-N/4)))/(1+exp(2*k*(I-N/4))))+0.5/alpha*log((1+exp(2*k*(I-1-3*N/4)))/(1+exp(2*k*(I-3*N/4)))); % Kink anti-kink initial condition
%b(I+N)= 0.5*sinh(k)/alpha*(1/(exp(-2*k*(I-1-  N/4))+1)-1/(exp(-2*k*(I-  N/4))+1) )-0.5*sinh(k)/alpha*(1/(exp(-2*k*(I-1-3*N/4))+1)-1/(exp(-2*k*(I-3*N/4))+1) );
omegak2(I)=4*(sin(pi*I/2/N))^2;		% Frequencies of the linear modes
end
[T,Y]=ode45('fpu1',tspan,b',options,N);	% Time integration
for IT=1:(TMAX/DT),
  TIME(IT)=IT*DT*sqrt(omegak2(1))/2/pi;	% Time iteration loop
  YX(IT,1:N+1)=[0 Y(IT,1:N   )];
  YV(IT,1:N+1)=[0 Y(IT,N+1:2*N   )];
  sinXFFT(IT,:)= sqrt(2/(N+1))*0.5*imag(fft([YX(IT,1:N+1) 0  -YX(IT,N+1:-1:2)]));
  sinVFFT(IT,:)= sqrt(2/(N+1))*0.5*imag(fft([YV(IT,1:N+1) 0  -YV(IT,N+1:-1:2)]));
  Energy(IT,1:N)=(omegak2(1:N).*(sinXFFT(IT,2:N+1).^2)+sinVFFT(IT,2:N+1).^2)/2;
  for J=2:N-1,  		% Space loop  %YY(IT,N-1)=Y(IT,N);
    DifY(IT,J)=Y(IT,J+1)-Y(IT,J);
    EE(IT,J)=(Y(IT,N+J)^2+(Y(IT,J+1)-Y(IT,J))^2/2+(Y(IT,J)-Y(IT,J-1))^2/2+alpha/3*((Y(IT,J+1)-Y(IT,J))^3+(Y(IT,J)-Y(IT,J-1))^3))/2;
end
  EE(IT,1)=Y(IT,N+1)^2/2+(Y(IT,2)-Y(IT,1))^2/4+(Y(IT,1))^2/4+alpha/6*((Y(IT,2)-Y(IT,1))^3+(Y(IT,1))^3);
  EE(IT,N)=Y(IT,N+N)^2/2+(Y(IT,N))^2/4 +(Y(IT,N)-Y(IT,N-1))^2/4+alpha/6*(-Y(IT,N)^3+(Y(IT,N)-Y(IT,N-1))^3);
end

surf(EE) ;

plot(TIME,Energy(:,1),TIME, Energy(:,2),TIME, Energy(:,3),TIME, Energy(:,4))
xlabel('\omega_1t/2\pi');ylabel('E_k');

surf(DifY); 	% representation of the space derivative field to show the soliton dynamics

The program calls the function fpu1, which computes the equations of motion and must be saved in the file fpu1.m in the same folder (from the same paper, with the value of alpha of the main program and N passed as an argument):

function dy=fpu1(t,y,N)
alpha=0.5;	% Nonlinear parameter: same value as in the main program
D(N+1)=y(2)-2*y(1)+alpha*((y(2)-y(1))^2-y(1)^2);	D(1)=y(N+1);
D(2*N)=y(N-1)-2*y(N)+alpha*(y(N)^2-(y(N)-y(N-1))^2);	D(N)=y(2*N);
for I=2:N-1,
D(N+I)=y(I+1)+y(I-1)-2*y(I)+alpha*((y(I+1)-y(I))^2-(y(I)-y(I-1))^2);
D(I)=y(N+I);
end
dy=D';

Notes