Nel 1953-54 Enrico Fermi, John Pasta e Stanislaw Ulam usarono il calcolatore MANIAC I del laboratorio di Los Alamos per uno dei primi esperimenti numerici della storia: la simulazione di una catena di oscillatori con un debole accoppiamento non lineare. I programmi furono scritti ed eseguiti da Mary Tsingou, oggi riconosciuta come coautrice: per questo si parla anche di problema FPUT^[1].

L'articolo originale con la discussione del problema qui esaminato, originariamente apparso sul rapporto LA-1940 (maggio 1955) del Los Alamos National Laboratory, è stato pubblicato dopo la morte di Fermi nella raccolta delle sue opere.

Modi normali in un sistema non lineare

• Come illustrato meglio nella pagina dedicata, i modi normali di un sistema lineare sono sempre indipendenti, cioè quando un modo viene eccitato esso non scambia energia con gli altri modi normali, e il moto del sistema rimane stazionario in quel modo.
• Tuttavia, quando il sistema è non lineare, il principio di sovrapposizione non si applica più, e i modi normali possono interagire tra loro scambiandosi energia. Allora, anche supponendo di poter eccitare inizialmente un solo modo, l'energia fornita al sistema potrà trasferirsi, col passare del tempo, ad altri modi. Il moto del sistema non è più stazionario.

Il problema di Fermi, Pasta e Ulam

La domanda dunque era: in presenza di un debole accoppiamento non lineare quali modi saranno eccitati dopo un tempo abbastanza lungo?'

La risposta che Fermi, Pasta e Ulam si aspettavano era basata sul teorema di equipartizione dell'energia: l'energia si distribuisce in misura uguale tra tutti i gradi di libertà del sistema.

• Il sistema fisico più semplice su cui mettere alla prova questa congettura è la catena di oscillatori in cui però si introduca un debole accoppiamento non-lineare. Questo problema non ammette una soluzione analitica semplice come il suo analogo lineare, e solo una simulazione numerica può svelarne la dinamica.

• Potete provare voi stessi la simulazione con il programma che trovate in fondo a questa pagina. Un analogo di questo programma girava nel 1953-54 sul calcolatore MANIAC I del laboratorio di Los Alamos, simulando catene di 16, 32 e 64 oscillatori, e svelò a Fermi, Pasta e Ulam che la loro congettura era sbagliata.

• Lasciando evolvere nel tempo il sistema virtuale si osserva dapprima una ridistribuzione dell'energia: eccitando inizialmente il modo 1 si attivano successivamente i modi 2, 3, 4. Tuttavia, se si lascia passare un tempo sufficientemente lungo (157 periodi del modo 1), si osserva che l'energia ritorna quasi completamente nel modo 1, e che esiste una periodicità nel modo in cui i diversi modi vengono eccitati.

• Fermi, Pasta e Ulam avevano così dimostrato che, per ottenere l'equipartizione dell'energia, non è sufficiente un piccolo accoppiamento non lineare.

Verifica la comprensione

Programma per la simulazione

Il seguente programma Matlab è pubblicato 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), e qui riprodotto per cortesia di 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

Il programma chiama la funzione fpu1, che calcola le equazioni del moto e va salvata nel file fpu1.m nella stessa cartella (dallo stesso articolo, con il valore di alpha del programma principale e N passato come argomento):

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';

Note