diff --git a/Jaslo/sheet3.py b/Jaslo/sheet3.py index 30d6bafba3139c5a044a4eabe588c7d814a69039..cc34dfbee7df69a037d8c18362c0cc65c63d91d7 100644 --- a/Jaslo/sheet3.py +++ b/Jaslo/sheet3.py @@ -1,6 +1,8 @@ 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) @@ -11,7 +13,7 @@ def verlet(force, x0, p0, m, dt, N): x[0] = x0 p[0] = p0 - for i in range(1, N): + 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 @@ -43,81 +45,87 @@ p0[-1] = -p0[0] dt = 0.1 t = np.arange(0, 300, dt) -x, p = verlet(force, x0, p0, m, dt, len(t)) +# x, p = verlet(force, x0, p0, m, dt, len(t)) -# fig, ax = plt.subplots() +# # fig, ax = plt.subplots() -for i in range(x.shape[1]): - plt.plot(t, x[:,i], label=f"Bead {i+1}") +# 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() +# 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) +# # 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() +# # 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] +# # 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 +# 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}") +# 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) +# # 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]) +# 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)}") +# 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() +# 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) -t = np.arange(0, 2e4, dt) - -for N in (10, 20): +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.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):] + ts = t[int(10 / dt):] ds = np.zeros(ts.shape) - for i in range(len(ds)): + 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.show() + # plt.tight_layout() + plt.savefig(f'4c_N={N}.png')