import numpy as np from scipy.constants import Boltzmann as k_B import matplotlib.pyplot as plt from tqdm import tqdm 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 tqdm(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)) dt = 0.1 t = np.arange(0, 1*10e3, dt) N_vals = [10,20] for N in N_vals: 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) 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)) ts = t[int(10 / dt):] ds = np.zeros(ts.shape) for i in tqdm(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.figure(figsize=(48, 6)) plt.xlabel(r"$t / \tau$") plt.ylabel(r"$d / \sqrt{k l^2}$") plt.plot(ts, ds) # plt.tight_layout() plt.savefig(f'4c_N={N}.png')