%pylab inline
import numpy as np; np.random.seed(3)
import matplotlib.pyplot as plt
import seaborn as sns; sns.set(style="white", palette=sns.color_palette("Greys", n_colors=3))

Populating the interactive namespace from numpy and matplotlib

import io # these imports are needed to display .mp4 video
import base64
from IPython.display import HTML

def animate(traj, frame, scale=100): # store .png in output/ folder
    plt.figure()
    plt.axis("off")
    plt.xlim([-3, 3])
    plt.ylim([-15, 15])
    colors = cm.RdBu(np.append(np.linspace(0., 0.25, N/2), np.linspace(0.75, 1., N/2)))
    for k in range(N):
        plt.scatter(traj[2*k, frame], traj[2*k+1, frame], color=colors[k], s=2)
    plt.savefig("output/" + str(int(frame / scale)) + ".png")

a, b

The Hamiltonian $H = \frac{p^2}{2m} + V(q)$ governs the time evolution of the system. The changes in $p$ and $q$ are calculated from:

• $\dot{q} = \partial_p H$.
• $\dot{p} = -\partial_q H$.

In the given example we set the mass to one, thus $H = \frac12 p^2 + \frac12 q^2 \big( \frac12 q^2 - 1 \big)$. Calculate the derivates explicitly to construct $\phi$:

• $\dot{q} = p$.
• $\dot{p} = - q^3 + q$.

c, d

Introduction of $p=\dot{x}$, $q=x$ and $\lambda = t$ yields: $\begin{pmatrix} \dot{x} \\ \dot{p} \end{pmatrix} = \begin{pmatrix} p \\ \mu (1 - q^2) p - q - 3 \cos \lambda \end{pmatrix}$.

The harmonic oscillator $m\ddot{x} = -kx$ is rewritten by introducing $p$ and $g=\frac{k}{m}$: $\begin{pmatrix} \dot{x} \\ \dot{p} \end{pmatrix} = \begin{pmatrix} p \\ g x\end{pmatrix}$.

In order to determine the phase space trajectory of a gas with $N_A$ particles in $d$ dimensions, one would need some $\gamma_0 \in \mathbb{R}^{2 N_A d}$.

e, g

A first step is to discretize the equations of motion mentioned before. For a harmonic oscillator potential $V(q) = \frac12 kq^2$ one would get:

• $q_{t+1} = q_t + p_t\Delta t$.
• $p_{t+1} = p_t - kx_t\Delta t$.

The potential that was used to create the animation is $V(q) = \frac12 q^2 - 4 q^4 + \frac18 q^8$. Its derivative is: $V'(q) = q - 16q^3 + q^7$. Both are depicted in a below plot. To make things more interesting we also add a local distortion to the force, ie. $F(q) = -V'(q) + \varepsilon$ where $\varepsilon \sim \mathcal{N}(0,\sigma^2)$. Set $\sigma^2 = \frac{N}{150}$.

def V(x):
    return .5 * (x ** 2) - 4 * (x ** 4) + .125 * (x ** 8)
def F(x, sigma=1e-6):
    return -(x - 16 * (q ** 3) + q ** 7) + np.random.normal(0, sigma ** 2, x.shape[0])

q = np.linspace(-2.5, 2.5, 200)
plt.fill_between(q, -100, V(q), lw=0, label="potential", color="darkblue", alpha=.2)
plt.plot(q, F(q), label="force", color="navy", alpha=.8)
plt.xlim(q[0], q[-1])
plt.ylim([-50, 100])
plt.xlabel("q")
plt.legend(loc="best");

Set the initial conditions to $q_0 = \vec{0}$ and $p_0 = (-\frac{N}{80},\dots,\frac{N}{80})$. Calculate the trajectories of $N=300$ particles.

N = 300
T = 10000
traj = np.zeros(shape=(2*N, T))
q = np.zeros(shape=N)
p = (np.arange(0, N) - (N-1)/2) * .025

For each point, generate $(q_{t+1} , p_{t+1})$ from $f(q_t,p_t)$.

def f(q, p, dt=2e-3, sigma=10.):
    return q + p * dt, p - ((q - 16 * (q ** 3) + (q ** 7)) + np.random.normal(0., sigma)) * dt

for t in range(T):
    for k in range(N):
        q[k], p[k] = f(q[k], p[k])
    traj[::2, t] = q
    traj[1::2, t] = p

%%capture
choose_scale = 10
for f in range(choose_scale, 3000, choose_scale):
    animate(traj, f, scale=choose_scale)

Display them via HTML5. Some observations:

• Over time, the points smear out in phase space. In the end, it seems to make no difference whether they started out with a negative (red) or positive (blue) momentum.
• The energy of each particle is finite, so the points never leave the window.
• There is no difference when animating backwards in time. All that happens is that blue is exchanged with red.
• For Van der Waals interactions, we need to calculate forces between particles. This is going to be a lot more demanding in terms of the required computational resources, as it involves the calculation of $\frac{N}{2}(N-1)$ terms.

video = io.open('video/output.mp4', 'r+b').read()
encoded = base64.b64encode(video)
HTML(data='''<video width="600" height="450" controls>
                <source src="data:video/mp4;base64,{0}" type="video/mp4" />
             </video>'''.format(encoded.decode('ascii')))