Skip to content
Snippets Groups Projects
Commit 596fe2c2 authored by ziskaj00's avatar ziskaj00
Browse files

Add problem 10.1 c

parent 51c2f6b2
No related branches found
No related tags found
No related merge requests found
......@@ -7,6 +7,8 @@ plt.rcParams.update({
"font.sans-serif": "Helvetica",
})
from tqdm import tqdm
# a - Implementation based on Leimkuhler & Matthews
gamma = 1
omega_0 = 1/3 * gamma
......@@ -15,8 +17,8 @@ d = 3
k_BT = 1
l = 1 # gamma**-1 * np.sqrt(k_BT / m)
N = 100
t_max = 1000
dt = 0.1
t_max = 1000
dt = 0.1
t = np.arange(0,t_max,dt)
def force(r: np.ndarray) -> np.ndarray:
......@@ -35,7 +37,7 @@ def BAOAB(force, x0, p0, m, dt, N):
r = np.random.normal(0,1,(N,d))
p[i] = p[i-1] + 1/2 * force(x[i-1]) * dt # same
x[i] = x[i-1] + p[i] / (2*m) * dt # diff
p[i] = np.exp(-gamma*dt) * p[i] + xi_2 * r/np.sqrt(m)
p[i] = np.exp(-gamma*dt) * p[i] + xi_2 * r * np.sqrt(m)
x[i] = x[i] + (dt/2) * p[i]/m
p[i] = p[i] + 1/2 * force(x[i]) * dt
......@@ -94,4 +96,23 @@ E_pot = (m/(2*N)) * omega_0**2 * np.linalg.norm(r, axis=(1,2))**2
# plt.plot(hist_p, label=f'p_{coords[i]}')
# plt.legend()
# plt.tight_layout()
# plt.savefig(f'p_{coords[i]}_dist.png')
\ No newline at end of file
# plt.savefig(f'p_{coords[i]}_dist.png')
# c)
r_f, p_f = BAOAB(lambda r: 0, r0, p0, m, dt, N)
dr2 = np.zeros(t.shape)
dr2_f = np.zeros(t.shape)
for i in tqdm(range(1, len(dr2))):
dr2[i] = np.average(np.linalg.norm(r[i:] - r[:-i], axis=2)**2)
dr2_f[i] = np.average(np.linalg.norm(r_f[i:] - r_f[:-i], axis=2)**2)
plt.loglog(t, dr2, label="harmonic")
plt.loglog(t, dr2_f, label="free")
plt.loglog(t, 3 * k_BT / m * t**2, ls="--", label=r"$3 k_B T t^2 / m$")
plt.loglog(t, 6 * k_BT * t / (m * gamma) - 6 * k_BT / (m * gamma**2) * (1 - np.exp(-gamma * t)), ls="-.", label="free, analytic")
plt.legend()
plt.xlabel(r"$\log t$")
plt.ylabel(r"$\log \delta r^2$")
plt.savefig("problem10.1c.png", dpi=300)
plt.show()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment