Update sheet 10

PS10_Jung/PS10.py

 import numpy as np
+from matplotlib import pyplot as plt
+
+plt.rcParams.update({
+    "text.usetex": True,
+    "font.family": "sans-serif",
+    "font.sans-serif": "Helvetica",
+})
 
 # a - Implementation based on Leimkuhler & Matthews
 gamma = 1
@@ -28,8 +35,8 @@ 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_hat = np.exp(-gamma*dt) * p[i] + xi_2 * r/np.sqrt(m)
-        x[i] = x[i] + (dt/2) * p_hat/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
 
     return x, p
@@ -37,4 +44,54 @@ def BAOAB(force, x0, p0, m, dt, N):
 r0 = np.zeros((N, d))
 p0 = np.zeros((N, d))
 
-r, p = BAOAB(force=force, x0=r0, p0=p0, m=m, dt=dt, N=N)
\ No newline at end of file
+r, p = BAOAB(force=force, x0=r0, p0=p0, m=m, dt=dt, N=N)
+
+# b
+# # kin and pot
+
+E_kin = np.sum(p**2, axis=(1,2)) / (2*m*N)
+E_pot = (m/(2*N)) * omega_0**2 * np.linalg.norm(r, axis=(1,2))**2
+
+# plt.figure()
+# plt.plot(E_kin, label=r'$E_{kin}$')
+# plt.plot(E_pot, label=r'$E_{pot}$')
+# plt.plot(np.full(fill_value=np.average(E_kin),shape=len(t)), label=r'$E_{kin\_av}$', linestyle='--')
+# plt.plot(np.full(fill_value=np.average(E_pot),shape=len(t)), label=r'$E_{kin\_pot}$', linestyle='dashed')
+# plt.vlines(50, min(np.min(E_kin), np.min(E_pot)), max(np.max(E_kin), np.max(E_pot)), linestyles='dashed', colors='purple')
+# plt.xlabel(r'$t / \gamma^{-1}$')
+# plt.ylabel(r'$E / k_BT$')
+# plt.grid(True)
+# plt.legend(loc='lower right')
+# plt.tight_layout()
+# plt.savefig('10.1_b.png')
+
+# # 3D-gauss dist
+# r_x = r[51:,:,0].flatten()
+# r_y = r[51:,:,1].flatten()
+# r_z = r[51:,:,2].flatten()
+# p_x = p[51:,:,0].flatten()
+# p_y = p[51:,:,1].flatten()
+# p_z = p[51:,:,2].flatten()
+
+
+# fig = plt.figure()
+# ax = fig.add_subplot(projection='3d')
+# ax.scatter(r_x, r_y, r_z, s=10, facecolors='none', edgecolors='r')
+# plt.savefig('norm_dist_x.png')
+# coords = ['x','y','z']
+# for i in range(3):
+#     r_i =r[51:,:,i].flatten()
+#     hist_r, r_range = np.histogram(r_i, bins = len(r_i),density=True)
+#     plt.figure()
+#     plt.plot(hist_r, label=f'r_{coords[i]}')
+#     plt.legend()
+#     plt.tight_layout()
+#     plt.savefig(f'r_{coords[i]}_dist.png')
+
+#     p_i =r[51:,:,i].flatten()
+#     hist_p, p_range = np.histogram(p_i, bins = len(p_i),density=True)
+#     plt.figure()
+#     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