Add sheet3.py

Jaslo/sheet3.py

+import numpy as np
+from scipy.constants import Boltzmann as k_B
+import matplotlib.pyplot as plt
+
+def verlet(force, x0, p0, m, dt, N):
+    assert(x0.shape == p0.shape)
+
+    x = np.zeros((N, *x0.shape))
+    p = np.zeros((N, *p0.shape))
+
+    x[0] = x0
+    p[0] = p0
+
+    for i in range(1, N):
+        p[i] = p[i-1] + 1/2 * force(x[i-1]) * dt
+        x[i] = x[i-1] + p[i] / m * dt
+        p[i] = p[i] + 1/2 * force(x[i]) * dt
+
+    return x, p
+
+# Problem 3.1
+# a)
+k = 1 # ?
+m = np.full(11, 1) # ?
+l = 1 # ?
+eps = k * l**2
+tau = np.sqrt(m[0] / k)
+
+def force(x):
+    f = np.zeros(x.shape)
+
+    f[:-1] += -k * (x[:-1] - x[1:])
+    f[1:] += -k * (x[1:] - x[:-1])
+
+    return f
+
+I = np.arange(1, 12)
+x0 = l * I
+p0 = np.zeros(x0.shape)
+p0[0] = 1/2 * np.sqrt(m[0] * eps)
+p0[-1] = -p0[0]
+
+dt = 0.1
+t = np.arange(0, 300, dt)
+
+x, p = verlet(force, x0, p0, m, dt, len(t))
+
+# fig, ax = plt.subplots()
+
+for i in range(x.shape[1]):
+    plt.plot(t, x[:,i], label=f"Bead {i+1}")
+
+plt.xlabel(r"$t / \tau$")
+plt.ylabel(r"$x / l$")
+plt.legend()
+plt.show()
+
+# E = np.sum(p**2 / (2 * m), axis=1) + 1/2 * k * np.sum((x[:,1:] - x[:,:-1])**2, axis=1)
+# K = np.sum(p**2 / (2 * m), axis=1)
+
+# plt.plot(t, K)
+# plt.show()
+
+# b)
+for N in (2, 3, 11, 100):
+    m = np.full(N, 1)
+    I = np.arange(1, N+1)
+    x0 = l * I
+    p0 = np.zeros(x0.shape)
+    p0[0] = 1/2 * np.sqrt(m[0] * eps)
+    p0[-1] = -p0[0]
+
+    x, p = verlet(force, x0, p0, m, dt, len(t))
+    K = np.sum(p**2 / (2 * m), axis=1)
+    K_l = np.average(K)
+    T = 2 * K_l / N
+
+    print(f"<K> = {K_l}")
+    print(f"k_B T = {T}")
+
+# c)
+for N in (2, 3, 11, 100):
+    for i in range(3):
+        m = np.full(N, 1)
+        I = np.arange(1, N+1)
+        x0 = l * I
+        p0 = np.random.normal(0, np.sqrt(m[0] * eps), N)
+
+        x, p = verlet(force, x0, p0, m, dt, len(t))
+        R_l = np.abs(x[:,0] - x[:,-1])
+
+        print(f"N = {N}, i = {i}: R_l = {np.average(R_l)}")
+
+        plt.plot(t, R_l)
+        plt.xlabel(r"$t / \tau$")
+        plt.ylabel(r"$\overline{R_l} / l$")
+        plt.show()
+
+# d)
+def d(x1, p1, x2, p2):
+    return np.sqrt(1 / (6 * x1.shape[0]) * np.sum(k * (x1 - x2)**2 + (p1 - p2)**2 / m))
+
+t = np.arange(0, 2e4, dt)
+
+for N in (10, 20):
+    m = np.full(N, 1)
+    I = np.arange(1, N+1)
+    x0 = l * I
+    p0 = np.random.normal(0, np.sqrt(m[0] * eps), N)
+
+    x, p = verlet(force, x0, p0, m, dt, len(t))
+
+    ts = t[int(10 / dt):]
+    ds = np.zeros(ts.shape)
+    for i in range(len(ds)):
+        ds[i] = d(x[int(10 / dt) + i], p[int(10 / dt) + i], x[int(10 / dt)], p[int(10 / dt)])
+
+    print(f"N = {N}: d_min = {np.min(ds)}")
+
+    plt.xlabel(r"$t / \tau$")
+    plt.ylabel(r"$d / \sqrt{k l^2}$")
+    plt.plot(ts, ds)
+    plt.show()