In [179]:
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib import animation
from matplotlib.animation import FuncAnimation
import scipy.constants as spc
from IPython.core.display import HTML
from IPython.display import display, Math, Latex
from scipy.ndimage import convolve, generate_binary_structure
%matplotlib inline                  
%config InlineBackend.figure_format = 'retina'
plt.rcParams.update(
    {
    'figure.constrained_layout.use': False,
    'mathtext.fontset':'cm',
    'font.family':'serif', 
    'font.sans-serif':'Times New Roman',
    'figure.titlesize':22,
    "figure.dpi":90, 
    'savefig.dpi':90,
    'axes.titlesize':20, 
    'axes.labelsize':18, 
    'axes.titley': 1.0, 
    'axes.titlepad': 5.0,
    'axes.edgecolor':'black', 
    'axes.grid': False, 
    'grid.alpha': .5,
    'xtick.labelsize':14,'ytick.labelsize':14,
    'xtick.major.size':6,'ytick.major.size':6,
    'xtick.major.width':1.25, 'ytick.major.width':1.25, 
    'xtick.direction':'inout','ytick.direction':'inout',
    'xtick.top':False, 'ytick.right':False,
    'legend.title_fontsize':14, 'legend.fontsize':14,
    'legend.borderaxespad': 1, 'legend.borderpad': 0.5,
    'legend.framealpha': 1,
    'legend.handleheight': 0.5, 'legend.handlelength': 2.0, 'legend.handletextpad':0.5,
    'legend.labelspacing': 0.25,
    'legend.fancybox':False,
    'legend.edgecolor': '0',
    'legend.frameon': True,
    'legend.markerscale': 1.25,
    'animation.embed_limit':2**128,
    'animation.html': 'jshtml'
    }
)
pi = spc.pi

1D Fermi Pasta Ulam Tsingou Lattice¶

History¶

One year before his death, Enrico Fermi had to make a very frustrating choice. Deep within Los Alamos National Laboratory, he found himself face-to-face with a $1000 \mathrm{lbs}$ hunk of untapped computational power, MANIAC I. He needed to sniff out a problem in Physics that was worthy of numerical simulation.

Fermi chose a problem rooted in alarming mechanical simplicity that posed an even more alarming threat to an idea fundamental to thermodynamics — the Ergodic Hypothesis, an assumption at the bedrock of statistical mechanics as it was. Ergodicity in a system posits an equivalence between averaging a process as it evolves over long periods of time, and averaging the process over a statistical ensemble at a single instant in time. Simply put, the Ergodic Hypothesis dramatically diminishes the problem of simulating a system forever into the incomparably more feasible task of averaging over a plethora of independent realization of the same system.

Fermi wondered whether Ergodicity — antecedent to the Equipartition Theorem — was a valid assumption to make at all. If it proved to be invalid upon simulation, the study of thermodynamics would have been a scientific transgression dragged through centuries of erroneous development. Maybe there was some humor in that, and maybe Fermi just needed a good laugh.

With the programming talent of Mary Tsingou and expertise of John Pasta and Stanislaw Ulam, Fermi's team studied a 1D lattice of point-masses coupled to nearest neighbors with an interaction potential containing the usual elastic potential energy and its familiar linear restoring force, but also the very next term in the Taylor-expanded potential — one that delivers a directionally-biased force quadratic in the separation.

Once the system was placed in a normal mode, Fermi hypothesized that as the system advances through a substantial number of cycles, all other modes would, in time, be explored with equal likelihood. the influence of initially excited normal modes would fade and eventually lead to thermalization.

The team let the system run overnight, and to their surprise, it snuck back into the initial configuration like a teenager stumbling home at 3 a.m., greeting Mom and Dad with a smile that belied a boozed-up night of partying. That's their kid all right, but who knows where they've been all night. What's going on here?

A recent breakthrough.

The Discrete $\alpha$-FPUT Lattice¶

Here we consider a chain of $N$ particles with:

  • Mass $m$ $\left[\mathrm{kg}\right]$
  • Equilibrium spacing $\ell$ $\left[\mathrm{m}\right]$
  • Linear spring constant $\kappa$ $\left[\dfrac{\mathrm{kg}}{\mathrm{s^2}}\right]$
  • Non-linear coupling parameter $\alpha$ $\left[\dfrac{1}{\mathrm{m}}\right]$

The displacement and momentum of the $n^\mathrm{th}$ particle along the chain are denoted $q_{n}$ and $p_{n}$. Accordingly, we can express the Hamiltonian of the entire system as: \begin{align} H = \sum_{j=0}^{N} \frac{1}{2m}p_{j}^{2} + \frac{\kappa}{2}(q_{j+1} - q_{j})^2 + \frac{\alpha\kappa}{3}(q_{j+1} - q_{j})^3 \end{align}

which leads to the following canonical equations of motion for the $n^\mathrm{th}$ particle: \begin{align} \dot{q}_{n} &= \frac{p_{n}}{m}\\ \dot{p}_{n} &= \kappa\left(q_{n+1} + q_{n-1} - 2 q_{n}\right) \Big(1 + \alpha \left(q_{n+1} - q_{n-1}\right)\Big) \end{align}

The linear frequency of the system (when $\alpha = 0$) is $\displaystyle \omega^2 = \frac{\kappa}{m}$.

Dimensionless variables are introduced for time, displacement and momenta coordinates: \begin{align} t &\longrightarrow {t}\cdot \frac{1}{\omega}\\ q_n &\longrightarrow {q}_n\cdot \ell \\ p_n &\longrightarrow {p}_n\cdot m\omega\ell \end{align}

The non-linear coefficient $\alpha$ has dimensions of reciprocal length. In the equations of motion, we see that the non-linearity extends over the region between the left and right nearest neighbors of the particle at $n$. Hence, we define the dimensionless parameter $\varepsilon$ as: $$\varepsilon = 2\alpha\ell$$ to characterize the strength of the non-linear interaction between nearest neighbors.

The dimensionless equations of motion are: \begin{align} \dot{q}_{n} &= p_n\\ \dot{p}_{n} &= \left( q_{n+1} + q_{n-1} - 2 q_{n} \right) \Big(1 + \frac{\varepsilon}{2}\left(q_{n+1} - q_{n-1}\right) \Big) \end{align}

We scale the Hamiltonian as $\displaystyle H \longrightarrow H \cdot \frac{1}{2} m\omega^2\ell^2$, so the dimensionless Hamiltonian in terms of the dimensionless variables is: \begin{align} {H}_{\alpha} = \sum_{j=0}^{N} {p}_{j}^{2}+({q}_{j+1}-{q}_{j})^2 +\frac{\varepsilon}{3}({q}_{j+1}-{q}_{j})^3 \end{align}

In [352]:
def fput_solve(q0, p0, N = 32, epsilon = 1, dt = 0.25, steps = 10_000):
    """ Solve for the dynamics of the discrete alpha-FPUT lattice.
        Args:
            - q0(np.array) : initial displacement profile.
            - p0(np.array) : initial momentum profile.
        Kwargs:
            - N (int) : Number of active particles in the lattice (excluding the fixed ends).
            - epsilon (float) : non-linear coupling strength of cubic potential.
            - dt (float) : time step.
            - steps (int) : Number of time steps.
        Returns:
            t, q, p (dimensionless arrays for simulation time, displacement profiles and momentum profiles).
    """
    def F(x):
        d1 = np.roll(x,-1)+np.roll(x,1)-2*x
        d2 = np.roll(x,-1)-np.roll(x,1)
        d1[0] = d1[-1] = d2[0] = d2[-1] = 0.0        
        force = d1*(1 + 0.5*epsilon*d2)
        force[0] = force[-1] = 0.0
        return force
    
    t = np.arange(steps) * dt
    q = np.zeros((steps+2, N+2))
    p = np.zeros((steps+2, N+2))
    
    q[0] = q[1] = q0
    p[0] = p[1] = p0

    for i in range(1,steps+1):
        # q[:,0] = q[:,-1] = 0
        # p[:,0] = p[:,-1] = 0
        p[i+1] = p[i] + dt * F(q[i])
        q[i+1] = q[i] + dt * p[i+1]
        q[:,0] = q[:,-1] = 0
        p[:,0] = p[:,-1] = 0
    t = t    
    q = q[1:-1]
    p = p[1:-1]
    return t, q, p

def dst(q):
    """Discrete Sine Transform"""
    N = q.shape[1]
    A = np.sqrt(2/(N-1))
    k = np.arange(N)
    n = k.reshape((N,1))
    sin_k = np.sin(pi*k*n/(N-1))
    Q = A*np.dot(q, sin_k)
    return Q

def get_energy(Q, P):
    Qmode, Pmode = Q.T, P.T
    N = Q.shape[1]
    E_k = np.zeros_like(Qmode)
    for k in range(N):
        freq_k = 2*np.sin(0.5*k*pi/(N-1))
        E_k[k] = (Pmode[k]**2 + freq_k**2 * Qmode[k]**2)/2
    return E_k.T

def init_mode(lattice, mode):
    N = len(lattice)
    A = np.sqrt(2/(N-1))
    return A * np.sin(pi*mode*lattice/(N-1))
In [385]:
N = 32
ε = 4
dt = 0.02
steps = 65_000
eq = np.arange(N+2)
q0 = init_mode(eq, 2)
p0 = 0*q0

t, q, p = fput_solve(q0, p0, N, ε, dt, steps)
Q, P = dst(q), dst(p)
E_k = get_energy(Q, P)
E_0 = sum(E_k[0])
In [386]:
%matplotlib inline                  
%config InlineBackend.figure_format = 'retina'
fig = plt.figure(figsize = (10,6), facecolor = 'white')
plt.plot(t/dt, E_k/E_0,'-',lw = 2)
plt.title('Energy Contributions from Normal Modes')
plt.xlabel('Simulation Steps ($n$)', labelpad = 20)
plt.ylabel(r'$\dfrac{\mathrm{E}_k}{E_{0}}$', rotation = 0, labelpad=20)
fig.tight_layout()
plt.show()
In [387]:
def anim1D(u, x, t, fps = 24, n_frames = 96):
    steps = int(len(t) - 1)
    frames = np.arange(0, steps, int(steps/n_frames))
    fig, ax = plt.subplots(1,1, figsize = (10,6))
    line, = ax.plot([], [], 'ko', mfc = 'w', mec = 'k')
    ax.set_xlim(x[0]-1, x[-1]+1)
    ax.set_ylim(1.25*u.min(),1.25*u.max())
    ax.axis('off')
    def animate(i):
        line.set_data(x, u[i])
        ax.set_title(r'${}\%$'.format(round(t[i]/t[-1]*100)))
        return line,
    plt.close()
    anim = FuncAnimation(fig, animate, interval = int(1e3/fps), frames = frames)
    return anim

anim1D(q, eq, t, fps = 30, n_frames = 600)
Out[387]: