Fermi-Pasta-Ulam Problem

Set-up

Considering particles connected with springs along one dimension, we investigate the motion of particles.
The springs have a non-linear property.
The mass of the particles is the same $m$, $y_n$ is displacement from an equilibrium state.
Let a force from growth ($\Delta$) of the springs $F=-k(\Delta+\alpha\Delta^2)$ where $k$ is a spring constant and $\alpha$ is a nonlinear constant.
The equation of motion is written as: $$\frac{{{d^2}{y_n}}}{{d{t^2}}} = {\left( {{y_{n + 1}} - 2{y_n} + {y_{n - 1}}} \right) + \alpha \left[ {{{\left( {{y_{n + 1}} - {y_n}} \right)}^2} - {{\left( {{y_n} - {y_{n - 1}}} \right)}^2}} \right]} .$$

Here, we normalized $t$ by $\sqrt{k/m}$.
We put a fixed boundary ($y_0=0,\,\, y_{N+1}=0$) where $N$ is the number of particles.

Linear Springs

When $\alpha=0$, the natural frequencies are obtained as $$\omega_l=2\sin\frac{l\pi}{2(N+1)}.$$

In this linear case, $y_n(t)$ can be written as \begin{gathered} {y_n}(t) = \sqrt {\frac{2}{{N + 1}}} \sum\limits_{l = 1}^N {{q_l}(t)\sin \frac{{nl\pi }}{{N + 1}}} \hfill \\ {q_l}(t) = \sqrt {\frac{2}{{N + 1}}} \sum\limits_{n = 1}^N {{y_n}(t)\sin \frac{{nl\pi }}{{N + 1}}}. \hfill \\ \end{gathered}

Here $q_l$ is a Fourier component.

In [1]:
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
from IPython import display
import scipy.fftpack as sf
import matplotlib.animation as animation
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5

Simulation setup

  • $N=32$
  • $\alpha=1/4$
  • $\Delta t=0.1$
  • Initial condition: $y_n=\sin\frac{n\pi}{N+1}, \frac{dy_n}{dt}=0, \,\, (n=1,\cdots, N)$
In [2]:
dt=0.1
nt=105000
n=32
n=n+2  #add the boundaries
alpha=0.25
jj=np.arange(n)
yp=np.sin(jj*np.pi/(n-1))#np.random.rand(n)#np.zeros(n)
yp[0]=0; yp[n-1]=0
vp=np.zeros(n)
yphst=np.zeros((nt,n))
for it in range(nt):
    vp[1:n-1]=vp[1:n-1]\
              +dt*(yp[2:n]-2*yp[1:n-1]+yp[0:n-2])\
              +dt*alpha*((yp[2:n]-yp[1:n-1])**2-(yp[1:n-1]-yp[0:n-2])**2)
    yp[1:n-1]=yp[1:n-1]+dt*vp[1:n-1]
    
    yphst[it,:]=yp
In [3]:
fig=plt.figure()
ims=[]
for it in range(0,nt,200):
    im=plt.plot(yphst[it,:],'k.-')
    plt.xlabel('# of particles')
    plt.ylabel('y')
    ims.append(im)
    
plt.rcParams['animation.html'] = 'html5'
ani = animation.ArtistAnimation(fig, ims)
plt.close()
In [4]:
ani
Out[4]:
In [5]:
T1=2*np.pi/2/np.sin(np.pi/2/(n-1))
plt.imshow(yphst,origin='lower',aspect='auto',cmap='bwr',extent=([0,n,0,nt*dt/T1]))
plt.xlabel('# of particles')
plt.ylabel(r'$t/T_p$')
clb=plt.colorbar()
clb.set_label('y')
plt.show()

Time-Evolution of Energy

The energy of the mode of the number $l$ is given $${E_l}(t) = \frac{1}{2}\left( {{{\left( {\frac{{d{q_l}}}{{dt}}} \right)}^2} + \omega _l^2q_l^2} \right).$$

In [6]:
jj=np.arange(n-2)+1
qphst=np.zeros((nt,n-2))
for ifr in range(n-2):
    csin=np.sin(jj*(ifr+1)*np.pi/(n-1))
    qphst[:,ifr]=np.sum(yphst[:,1:n-1]*csin,axis=1)

ephst=np.zeros((nt-1,n-2))
wp=2*np.sin(jj*np.pi/2/(n-1))
for it in range(nt-1):
    ephst[it,:]=0.5*(((qphst[it+1,:]-qphst[it,:])/dt)**2+wp**2*qphst[it,:]**2)
In [7]:
fig=plt.figure()
plt.plot(np.linspace(0,dt*nt/T1,nt-1),np.abs(ephst[:,0]),label='mode 1')
plt.plot(np.linspace(0,dt*nt/T1,nt-1),np.abs(ephst[:,1]),label='mode 2')
plt.plot(np.linspace(0,dt*nt/T1,nt-1),np.abs(ephst[:,2]),label='mode 3')
plt.plot(np.linspace(0,dt*nt/T1,nt-1),np.abs(ephst[:,3]),label='mode 4')
plt.legend(bbox_to_anchor=(0.9, 0.8, 0.5, .100))
plt.xlabel(r'$t/T_p$')
plt.ylabel(r'$E_l$')
plt.show()