sheet3.py

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')