Update for subprob 1.1.c

PS1_Jung/PS1.py

@@ -21,6 +21,7 @@ def heun(x0: float, p0: float, t: tuple, dt: float, f) -> np.array:
     :param f: force function
     """
     t = np.arange(t[0], t[1], dt)
+    dt=-dt
     x = np.zeros(len(t))
     x[0] = x0
 
@@ -39,6 +40,7 @@ def verlet(x0: float, p0: float, t: tuple, dt: float, f) -> np.array:
     :param t: (t0,tmax)
     """
     t = np.arange(t[0], t[1], dt)
+    dt = -dt
     x = np.zeros(len(t))
     x[0] = x0
 
@@ -68,49 +70,66 @@ if __name__ == '__main__':
 
     traj_heun_dt_1 = heun(X_0, P_0, t, dt_1, f)
     traj_heun_dt_2 = heun(X_0, P_0, t, dt_2, f)
-    traj_verlet_dt_1 = verlet(X_0, P_0, t, dt_1, f)
-    traj_verlet_dt_2 = verlet(X_0, P_0, t, dt_2, f)
+    # traj_verlet_dt_1 = verlet(X_0, P_0, t, dt_1, f)
+    # traj_verlet_dt_2 = verlet(X_0, P_0, t, dt_2, f)
 
 
-    # a) part 1 - x(t), p(t)
-    twoD_plot([traj_heun_dt_1[:2,:], traj_verlet_dt_1[:2,:]],['Heun', 'Verlet'],['t','x'], 'Heun_Verlet_dt=1e-3_x(t)', 'Heun_Verlet_dt=1e-3_x(t)', path )
-    twoD_plot([traj_heun_dt_2[:2,:], traj_verlet_dt_2[:2,:]],['Heun', 'Verlet'],['t','x'], 'Heun_Verlet_dt=1e-1_x(t)', 'Heun_Verlet_dt=1e-1_x(t)', path )
-    twoD_plot([traj_heun_dt_1[[0,2],:], traj_verlet_dt_1[[0,2],:]],['Heun', 'Verlet'],['t','p'], 'Heun_Verlet_dt=1e-3_p(t)', 'Heun_Verlet_dt=1e-3_p(t)', path )
-    twoD_plot([traj_heun_dt_2[[0,2],:], traj_verlet_dt_2[[0,2],:]],['Heun', 'Verlet'],['t','p'], 'Heun_Verlet_dt=1e-1_p(t)', 'Heun_Verlet_dt=1e-1_p(t)', path )
+    # # a) part 1 - x(t), p(t)
+    # twoD_plot([traj_heun_dt_1[:2,:], traj_verlet_dt_1[:2,:]],['Heun', 'Verlet'],['t','x'], 'Heun_Verlet_dt=1e-3_x(t)', 'Heun_Verlet_dt=1e-3_x(t)', path )
+    # twoD_plot([traj_heun_dt_2[:2,:], traj_verlet_dt_2[:2,:]],['Heun', 'Verlet'],['t','x'], 'Heun_Verlet_dt=1e-1_x(t)', 'Heun_Verlet_dt=1e-1_x(t)', path )
+    # twoD_plot([traj_heun_dt_1[[0,2],:], traj_verlet_dt_1[[0,2],:]],['Heun', 'Verlet'],['t','p'], 'Heun_Verlet_dt=1e-3_p(t)', 'Heun_Verlet_dt=1e-3_p(t)', path )
+    # twoD_plot([traj_heun_dt_2[[0,2],:], traj_verlet_dt_2[[0,2],:]],['Heun', 'Verlet'],['t','p'], 'Heun_Verlet_dt=1e-1_p(t)', 'Heun_Verlet_dt=1e-1_p(t)', path )
 
 
     # # a) part 2 - trajectory plot
     # phasespace_trajectory_plot(traj_heun_dt_1, 'Heun_dt=1e-3', 'Heun_dt=1e-3', path)
-    # phasespace_trajectory_plot(traj_heun_dt_2, 'Heun_dt=1e-1', 'Heun_dt=1e-1',path, fps=5)
+    # phasespace_trajectory_plot(traj_heun_dt_2, 'Heun_dt=1e-1', 'Heun_dt=1e-1',path)
     # phasespace_trajectory_plot(traj_verlet_dt_1, 'Verlet_dt=1e-3', 'Verlet_dt=1e-3', path)
     # phasespace_trajectory_plot(traj_verlet_dt_2, 'Verlet_dt=1e-1', 'Verlet_dt=1e-1', path)
 
-    # b) part 1 w/ E_tot = E_kin + E_pot = (p²)/2m + (1/2)*k*x²
-    E_heun_dt_1 = np.vstack((traj_heun_dt_1[0,:],(1/(2*M))*np.square(traj_heun_dt_1[2,:]) + (1/2)*K*np.square(traj_heun_dt_1[1,:])))
-    E_heun_dt_2 = np.vstack((traj_heun_dt_2[0,:],(1/(2*M))*np.square(traj_heun_dt_2[2,:]) + (1/2)*K*np.square(traj_heun_dt_2[1,:])))
-    E_verlet_dt_1 = np.vstack((traj_verlet_dt_1[0,:],(1/(2*M))*np.square(traj_verlet_dt_1[2,:]) + (1/2)*K*np.square(traj_verlet_dt_1[1,:])))
-    E_verlet_dt_2 = np.vstack((traj_verlet_dt_2[0,:],(1/(2*M))*np.square(traj_verlet_dt_2[2,:]) + (1/2)*K*np.square(traj_verlet_dt_2[1,:])))
-
-    # b) part 2
-    delta_E_heun_dt_1 = np.vstack((list(range(len(E_heun_dt_1[1,:])-1)),1/(dt_1)*np.diff(E_heun_dt_1[1,:])))
-    delta_E_heun_dt_2 = np.vstack((list(range(len(E_heun_dt_2[1,:])-1)),1/(dt_2)*np.diff(E_heun_dt_2[1,:])))
-    delta_E_verlet_dt_1 = np.vstack((list(range(len(E_verlet_dt_1[1,:])-1)),1/(dt_1)*np.diff(E_verlet_dt_1[1,:])))
-    delta_E_verlet_dt_2 = np.vstack((list(range(len(E_verlet_dt_2[1,:])-1)),1/(dt_2)*np.diff(E_verlet_dt_2[1,:])))
-
-    # twoD_plot(delta_E_heun_dt_1,('t', '$\Delta E$'), '$\Delta E$($\Delta t$)_heun_dt=1e-3','Delta_E(Delta_t)_heun_dt=1e-3',  path)
-    # twoD_plot(delta_E_heun_dt_2,('t', '$\Delta E$'), '$\Delta E$($\Delta t$)_heun_dt=1e-1','Delta_E(Delta_t)_heun_dt=1e-1',  path)
-    # twoD_plot(delta_E_verlet_dt_1,('t', '$\Delta E$'), '$\Delta E$($\Delta t$)_verlet_dt=1e-3','Delta_E(Delta_t)_verlet_dt=1e-3',  path)
-    # twoD_plot(delta_E_verlet_dt_2,('t', '$\Delta E$'), '$\Delta E$($\Delta t$)_verlet_dt=1e-1','Delta_E(Delta_t)_verlet_dt=1e-1',  path)
-
-    delta_E_heun_dt_1_sq = np.vstack((list(range(len(E_heun_dt_1[1,:])-1)),1/(dt_1**2)*np.diff(E_heun_dt_1[1,:])))
-    delta_E_heun_dt_2_sq = np.vstack((list(range(len(E_heun_dt_2[1,:])-1)),1/(dt_2**2)*np.diff(E_heun_dt_2[1,:])))
-    delta_E_verlet_dt_1_sq = np.vstack((list(range(len(E_verlet_dt_1[1,:])-1)),1/(dt_1**2)*np.diff(E_verlet_dt_1[1,:])))
-    delta_E_verlet_dt_2_sq = np.vstack((list(range(len(E_verlet_dt_2[1,:])-1)),1/(dt_2**2)*np.diff(E_verlet_dt_2[1,:])))
-
-    ## Problem 1.2
-    g = 4*np.pi**2
-    m_E = 3*10e-6
-    m_M = 3.69*10e-8
-    T_L = 27.3/365.3
-    R_EM = 2.57*10e-3
+    # # b) part 1 w/ E_tot = E_kin + E_pot = (p²)/2m + (1/2)*k*x²
+    # E_heun_dt_1 = np.vstack((traj_heun_dt_1[0,:],(1/(2*M))*np.square(traj_heun_dt_1[2,:]) + (1/2)*K*np.square(traj_heun_dt_1[1,:])))
+    # E_heun_dt_2 = np.vstack((traj_heun_dt_2[0,:],(1/(2*M))*np.square(traj_heun_dt_2[2,:]) + (1/2)*K*np.square(traj_heun_dt_2[1,:])))
+    # E_verlet_dt_1 = np.vstack((traj_verlet_dt_1[0,:],(1/(2*M))*np.square(traj_verlet_dt_1[2,:]) + (1/2)*K*np.square(traj_verlet_dt_1[1,:])))
+    # E_verlet_dt_2 = np.vstack((traj_verlet_dt_2[0,:],(1/(2*M))*np.square(traj_verlet_dt_2[2,:]) + (1/2)*K*np.square(traj_verlet_dt_2[1,:])))
 
+    # delta_E_heun_dt1_t = np.vstack((traj_heun_dt_1[0,:],np.insert(np.diff(E_heun_dt_1[1,:]),0,0)))
+    # delta_E_heun_dt2_t = np.vstack((traj_heun_dt_2[0,:],np.insert(np.diff(E_heun_dt_2[1,:]),0,0)))
+    # delta_E_verlet_dt1_t = np.vstack((traj_verlet_dt_1[0,:],np.insert(np.diff(E_verlet_dt_1[1,:]),0,0)))
+    # delta_E_verlet_dt2_t = np.vstack((traj_verlet_dt_2[0,:],np.insert(np.diff(E_verlet_dt_2[1,:]),0,0)))
+
+    # twoD_plot([delta_E_heun_dt1_t, delta_E_verlet_dt1_t], ['Heun', 'Verlet'], ['t','$\Delta E$'], 'Heun_Verlet_dt=1e-3_$\Delta E$(t)', 'Heun_Verlet_dt1_Delta_E(t)', path)
+    # twoD_plot([delta_E_heun_dt2_t, delta_E_verlet_dt2_t], ['Heun', 'Verlet'], ['t','$\Delta E$'], 'Heun_Verlet_dt=1e-1_$\Delta E$(t)', 'Heun_Verlet_dt2_Delta_E(t)', path)
+
+
+    # # b) part 2
+    # delta_E_heun_dt_1 = np.vstack((list(range(len(E_heun_dt_1[1,:])-1)),1/(dt_1)*np.diff(E_heun_dt_1[1,:])))
+    # delta_E_heun_dt_2 = np.vstack((list(range(len(E_heun_dt_2[1,:])-1)),1/(dt_2)*np.diff(E_heun_dt_2[1,:])))
+    # delta_E_verlet_dt_1 = np.vstack((list(range(len(E_verlet_dt_1[1,:])-1)),1/(dt_1)*np.diff(E_verlet_dt_1[1,:])))
+    # delta_E_verlet_dt_2 = np.vstack((list(range(len(E_verlet_dt_2[1,:])-1)),1/(dt_2)*np.diff(E_verlet_dt_2[1,:])))
+
+    # twoD_plot([delta_E_heun_dt_1, delta_E_verlet_dt_1], ['Heun', 'Verlet'], ['$\Delta t$','$\Delta E$'], '$\Delta E$($\Delta t$)dt=1e-3', 'Delta_E(Delta_t)_dt=1e-3', path)
+    # twoD_plot([delta_E_heun_dt_2, delta_E_verlet_dt_2], ['Heun', 'Verlet'], ['$\Delta t$','$\Delta E$'], '$\Delta E$($\Delta t$)dt=1e-1', 'Delta_E(Delta_t)_dt=1e-1', path)
+
+    # delta_E_heun_dt_1_sq = np.vstack((list(range(len(E_heun_dt_1[1,:])-1)),1/(dt_1**2)*np.diff(E_heun_dt_1[1,:])))
+    # delta_E_heun_dt_2_sq = np.vstack((list(range(len(E_heun_dt_2[1,:])-1)),1/(dt_2**2)*np.diff(E_heun_dt_2[1,:])))
+    # delta_E_verlet_dt_1_sq = np.vstack((list(range(len(E_verlet_dt_1[1,:])-1)),1/(dt_1**2)*np.diff(E_verlet_dt_1[1,:])))
+    # delta_E_verlet_dt_2_sq = np.vstack((list(range(len(E_verlet_dt_2[1,:])-1)),1/(dt_2**2)*np.diff(E_verlet_dt_2[1,:])))
+
+    # twoD_plot([delta_E_heun_dt_1_sq, delta_E_verlet_dt_1_sq], ['Heun', 'Verlet'], ['$\Delta t^{2}$','$\Delta E$'], '$\Delta E$($\Delta t^{2}$)dt=1e-3', 'Delta_E(Delta_t^2)_dt=1e-3', path)
+    # twoD_plot([delta_E_heun_dt_2_sq, delta_E_verlet_dt_2_sq], ['Heun', 'Verlet'], ['$\Delta t^{2}$','$\Delta E$'], '$\Delta E$($\Delta t^{2}$)dt=1e-1', 'Delta_E(Delta_t^2)_dt=1e-1', path)
+
+    # # c)
+    # (manually) time-inverted trajectories
+    x0 = 21372.097557449062 # from (correct) forward calculation
+    p0 = -589.1101949748515 # from (correct) forward calculation
+    traj_heun_dt_1_inv = heun(x0, p0, t, dt_1, f)
+    t1,x1,p1 = traj_heun_dt_1_inv
+    traj_heun_dt_1_inv= np.vstack((np.flip(t1),x1,p1))
+
+    traj_verlet_dt_1_inv = verlet(x0, p0, t, dt_1, f)
+    t2,x2,p2 = traj_verlet_dt_1_inv
+    traj_verlet_dt_1_inv= np.vstack((np.flip(t2),x2,p2))
+
+    phasespace_trajectory_plot(traj_heun_dt_1_inv, 'Heun_inv_dt=1e-3', 'Heun_inv_dt=1e-3', path)
+    phasespace_trajectory_plot(traj_verlet_dt_1_inv, 'Verlet_inv_dt=1e-3', 'Verlet_inv_dt=1e-3', path)

graphics.py

@@ -20,7 +20,8 @@ def twoD_plot(solutions: List[np.array], legend: List[str], axis: List[str], tit
         plt.plot(t,u, ls=linestyle[i], label=legend[i])
     plt.xlabel(x_axis); plt.ylabel(y_axis); plt.legend()
     plt.title(title)
-    plt.savefig(os.path.join(path, title + '.png'))
+    plt.tight_layout()
+    plt.savefig(os.path.join(path, fname + '.png'))
 
 
 def phasespace_trajectory_plot(solution: np.array, title: str, fname: str, path, fps=30, interval=500):