Skip to content
Snippets Groups Projects
Commit 8a70e1b4 authored by jung_42's avatar jung_42
Browse files

Update sheet3

parent be5cadc9
No related branches found
No related tags found
No related merge requests found
import numpy as np import numpy as np
from scipy.constants import Boltzmann as k_B from scipy.constants import Boltzmann as k_B
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from tqdm import tqdm
def verlet(force, x0, p0, m, dt, N): def verlet(force, x0, p0, m, dt, N):
assert(x0.shape == p0.shape) assert(x0.shape == p0.shape)
...@@ -11,7 +13,7 @@ def verlet(force, x0, p0, m, dt, N): ...@@ -11,7 +13,7 @@ def verlet(force, x0, p0, m, dt, N):
x[0] = x0 x[0] = x0
p[0] = p0 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 p[i] = p[i-1] + 1/2 * force(x[i-1]) * dt
x[i] = x[i-1] + p[i] / m * dt x[i] = x[i-1] + p[i] / m * dt
p[i] = p[i] + 1/2 * force(x[i]) * dt p[i] = p[i] + 1/2 * force(x[i]) * dt
...@@ -43,81 +45,87 @@ p0[-1] = -p0[0] ...@@ -43,81 +45,87 @@ p0[-1] = -p0[0]
dt = 0.1 dt = 0.1
t = np.arange(0, 300, dt) 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]): # for i in range(x.shape[1]):
plt.plot(t, x[:,i], label=f"Bead {i+1}") # plt.plot(t, x[:,i], label=f"Bead {i+1}")
plt.xlabel(r"$t / \tau$") # plt.xlabel(r"$t / \tau$")
plt.ylabel(r"$x / l$") # plt.ylabel(r"$x / l$")
plt.legend() # plt.legend()
plt.show() # plt.show()
# E = np.sum(p**2 / (2 * m), axis=1) + 1/2 * k * np.sum((x[:,1:] - x[:,:-1])**2, 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) # # K = np.sum(p**2 / (2 * m), axis=1)
# plt.plot(t, K) # # plt.plot(t, K)
# plt.show() # # plt.show()
# b) # # b)
for N in (2, 3, 11, 100): # for N in (2, 3, 11, 100):
m = np.full(N, 1) # m = np.full(N, 1)
I = np.arange(1, N+1) # I = np.arange(1, N+1)
x0 = l * I # x0 = l * I
p0 = np.zeros(x0.shape) # p0 = np.zeros(x0.shape)
p0[0] = 1/2 * np.sqrt(m[0] * eps) # p0[0] = 1/2 * np.sqrt(m[0] * eps)
p0[-1] = -p0[0] # p0[-1] = -p0[0]
x, p = verlet(force, x0, p0, m, dt, len(t)) # x, p = verlet(force, x0, p0, m, dt, len(t))
K = np.sum(p**2 / (2 * m), axis=1) # K = np.sum(p**2 / (2 * m), axis=1)
K_l = np.average(K) # K_l = np.average(K)
T = 2 * K_l / N # T = 2 * K_l / N
print(f"<K> = {K_l}") # print(f"<K> = {K_l}")
print(f"k_B T = {T}") # print(f"k_B T = {T}")
# c) # # c)
for N in (2, 3, 11, 100): # for N in (2, 3, 11, 100):
for i in range(3): # for i in range(3):
m = np.full(N, 1) # m = np.full(N, 1)
I = np.arange(1, N+1) # I = np.arange(1, N+1)
x0 = l * I # 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)
x, p = verlet(force, x0, p0, m, dt, len(t)) # x, p = verlet(force, x0, p0, m, dt, len(t))
R_l = np.abs(x[:,0] - x[:,-1]) # 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.plot(t, R_l)
plt.xlabel(r"$t / \tau$") # plt.xlabel(r"$t / \tau$")
plt.ylabel(r"$\overline{R_l} / l$") # plt.ylabel(r"$\overline{R_l} / l$")
plt.show() # plt.show()
# d) # d)
def d(x1, p1, x2, p2): def d(x1, p1, x2, p2):
return np.sqrt(1 / (6 * x1.shape[0]) * np.sum(k * (x1 - x2)**2 + (p1 - p2)**2 / m)) 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) N_vals = [10,20]
for N in N_vals:
for N in (10, 20):
m = np.full(N, 1) m = np.full(N, 1)
I = np.arange(1, N+1) I = np.arange(1, N+1)
x0 = l * I 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)) 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) 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)]) 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)}") print(f"N = {N}: d_min = {np.min(ds)}")
plt.figure(figsize=(48, 6))
plt.xlabel(r"$t / \tau$") plt.xlabel(r"$t / \tau$")
plt.ylabel(r"$d / \sqrt{k l^2}$") plt.ylabel(r"$d / \sqrt{k l^2}$")
plt.plot(ts, ds) plt.plot(ts, ds)
plt.show() # plt.tight_layout()
plt.savefig(f'4c_N={N}.png')
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment