Skip to content
Snippets Groups Projects
Commit f979ba20 authored by jasLogic's avatar jasLogic
Browse files

Add sheet4.py

parent 6d86de64
No related branches found
No related tags found
No related merge requests found
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 4.1
# a)
N = 50
k = 1 # ?
m = np.full(N, 1) # ?
M = np.sum(m)
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, N+1)
r0 = l * I
p0 = np.random.normal(0, np.sqrt(m[0] * eps), N)
dt = 0.1
dn = int(1 / dt)
t_max = 5000 # tau
t = np.arange(0, t_max, dt)
r = np.zeros((len(t), N))
p = np.zeros((len(t), N))
r[:dn], p[:dn] = verlet(force, r0, p0, m, dt, dn)
for i in range(1, t_max):
p[i * dn - 1] = np.random.normal(0, np.sqrt(m[0] * eps), N)
r[i * dn:(i+1) * dn], p[i * dn:(i+1) * dn] = verlet(force, r[i * dn - 1], p[i * dn - 1], m, dt, dn)
r_cm = 1/M * np.sum(m * r, axis=1)
R_g2 = 1/M * np.sum(m * (r.T - r_cm).T**2, axis=1)
R_l = np.abs(r[:,-1] - r[:,0])
E_pot = 1/2 * k * np.sum((r[:,1:] - r[:,:-1])**2, axis=1)
plt.plot(t, np.sqrt(R_g2))
plt.title(r"$R_g(t)$")
plt.xlabel(r"$t / \tau$")
plt.ylabel(r"$R_g / l$")
plt.show()
plt.plot(t, R_l)
plt.title(r"$R_l(t)$")
plt.xlabel(r"$t / \tau$")
plt.ylabel(r"$R_l / l$")
plt.show()
plt.plot(t, E_pot / N)
plt.title(r"$E_{pot} / N$")
plt.xlabel(r"$t / \tau$")
plt.ylabel(r"$E_{pot} / \epsilon$")
plt.show()
# b)
t_eq = 1500
print(f"<R_g²> = {np.mean(R_g2[t_eq * dn:])} ± {np.std(R_g2[t_eq * dn:])} l")
print(f"<R_l²> = {np.mean(R_l[t_eq * dn:]**2)} ± {np.std(R_l[t_eq * dn:]**2)} l")
print(f"<E_pot / N> = {np.mean(E_pot[t_eq * dn:] / N)} ± {np.std(E_pot[t_eq * dn:] / N)} ϵ")
for _ in range(4):
r[:dn], p[:dn] = verlet(force, r0, p0, m, dt, dn)
for i in range(1, t_max):
p[i * dn - 1] = np.random.normal(0, np.sqrt(m[0] * eps), N)
r[i * dn:(i+1) * dn], p[i * dn:(i+1) * dn] = verlet(force, r[i * dn - 1], p[i * dn - 1], m, dt, dn)
r_cm = 1/M * np.sum(m * r, axis=1)
R_g2 = 1/M * np.sum(m * (r.T - r_cm).T**2, axis=1)
R_l = np.abs(r[:,-1] - r[:,0])
E_pot = 1/2 * k * np.sum((r[:,1:] - r[:,:-1])**2, axis=1)
print(f"<R_g²> = {np.mean(R_g2[t_eq * dn:])} ± {np.std(R_g2[t_eq * dn:])} l")
print(f"<R_l²> = {np.mean(R_l[t_eq * dn:]**2)} ± {np.std(R_l[t_eq * dn:]**2)} l")
print(f"<E_pot / N> = {np.mean(E_pot[t_eq * dn:] / N)} ± {np.std(E_pot[t_eq * dn:] / N)} ϵ")
# c)
Ns = np.array([10, 20, 50, 100, 200])
R_g2_N = np.zeros(Ns.shape)
for j, (N, t_max, t_eq) in enumerate(zip(Ns, (5000, 5000, 5000, 8000, 12000), (1500, 1500, 1500, 4000, 8000))):
m = np.full(N, 1)
M = np.sum(m)
t = np.arange(0, t_max, dt)
I = np.arange(1, N+1)
r0 = l * I
p0 = np.random.normal(0, np.sqrt(m[0] * eps), N)
r = np.zeros((len(t), N))
p = np.zeros((len(t), N))
r[:dn], p[:dn] = verlet(force, r0, p0, m, dt, dn)
for i in range(1, t_max):
p[i * dn - 1] = np.random.normal(0, np.sqrt(m[0] * eps), N)
r[i * dn:(i+1) * dn], p[i * dn:(i+1) * dn] = verlet(force, r[i * dn - 1], p[i * dn - 1], m, dt, dn)
r_cm = 1/M * np.sum(m * r, axis=1)
R_g2 = 1/M * np.sum(m * (r.T - r_cm).T**2, axis=1)
R_g2_N[j] = np.mean(R_g2[t_eq * dn:])
# print(f"<R_g²> = {np.mean(R_g2[t_eq * dn:])} ± {np.std(R_g2[t_eq * dn:])}")
# plt.plot(t, R_g2)
# plt.show()
# plt.loglog(Ns, R_g2_N)
plt.plot(np.log(Ns), np.log(R_g2_N))
plt.show()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment