Update UB4

UB4/UB4.py

@@ -67,6 +67,36 @@ def implicit_midpoint(f, y0: np.ndarray, t: float, dt: float, eps = 1e-5, **kwar
 
     return y
 
+def IE(y0, t, dt, lamda = -1):
+    no_of_steps = int(t // dt)
+    y = np.zeros((no_of_steps, *y0.shape))
+    y[0] = y0
+
+    for i in range(1, no_of_steps):
+        y[i] = 1 / (1 - dt*lamda) * y[i-1]
+    
+    return y
+
+def EE(y0, t, dt, lamda = -1):
+    no_of_steps = int(t // dt)
+    y = np.zeros((no_of_steps, *y0.shape))
+    y[0] = y0
+
+    for i in range(1, no_of_steps):
+        y[i] = y[i-1] + dt * lamda*y[i-1]
+    
+    return y
+
+def IM(y0, t, dt, lamda = -1):
+    no_of_steps = int(t // dt)
+    y = np.zeros((no_of_steps, *y0.shape))
+    y[0] = y0
+
+    for i in range(1, no_of_steps):
+        y[i] = (1 + lamda * dt/2) /(1 - lamda * dt/2) * y[i-1] 
+
+    return y
+
 
 # Ex2: dy/dt = lambda * y, y(0) = 1, lambda = -1
 y0 = np.array(1)
@@ -83,17 +113,19 @@ exact_solution = lambda t : np.exp(lamda * t) * y0
 ## IE: A-stable (Yes), L-stable (Yes)
 ## IM: A-stable (Yes), L-stable (No)
 
-# for i in range(len(dt)):
-#     plt.figure()
-#     y_implicit_euler = implicit_euler(f, y0, t[0], dt[i])
+for i in range(len(dt)):
+    plt.figure()
     # plt.plot(implicit_euler(f, y0, t[0], dt[i]), '--', label = "implicit euler (IE)")
     # plt.plot(explicit_euler(f, y0, t[0], dt[i]), 'x', label = "explicit euler (EE)")
     # plt.plot(implicit_midpoint(f, y0, t[0], dt[i]), '>', label = "implicit midpoint (IM)")
-#     time = np.arange(int(t[0] // dt[i]))
-#     plt.plot(exact_solution(time), '3', label="exact solution (ES)")
-#     plt.legend()
-#     plt.title(f't={t[0]}, dt={dt[i]}')
-#     plt.show()
+    plt.plot(IE(y0, t[1], dt[i]), '--', label = "implicit euler (IE)")
+    plt.plot(EE(y0, t[1], dt[i]), 'x', label = "explicit euler (EE)")
+    plt.plot(IM(y0, t[1], dt[i]), '>', label = "implicit midpoint (IM)")
+    time = np.arange(int(t[1] // dt[i]))
+    plt.plot(exact_solution(time), '3', label="exact solution (ES)")
+    plt.legend()
+    plt.title(f't={t[0]}, dt={dt[i]}')
+    plt.show()
 
 #### For A-stability, dt < 1
 #### Observation
@@ -254,84 +286,84 @@ y0 = np.array([0.5,0.2,0.3])
 # p_dot = - partial H / partial q = -q
 # q_dot = partial H / partial p =  p
 
-p_0 = np.array([0])
-q_0 = np.array([1])
-t = 100
-dt = 0.001
-
-def HO_IE(p_0: np.ndarray, q_0: np.ndarray, t: float, dt: float):
-    """
-    Integrator = implicit Euler
-    p_dot = - partial H / partial q = -q
-    q_dot = partial H / partial p = p
-    """
-    no_of_steps = int(t // dt)
-    p = np.zeros((no_of_steps, *p_0.shape))
-    q = np.zeros((no_of_steps, *q_0.shape))
-    p[0] = p_0
-    q[0] = q_0
-
-    for i in trange(1, no_of_steps):
-        p[i] = 1/(1 + dt**2) * (p[i-1] - dt * q[i-1])
-        q[i] = 1/(1 + dt**2) * (q[i-1] + dt * p[i-1])
-
-    return p, q
-
-def HO_EE(p_0: np.ndarray, q_0: np.ndarray, t: float, dt: float):
-    """
-    Integrator = explicit Euler
-    p_dot = - partial H / partial q = -q
-    q_dot = partial H / partial p = p
-    """
-    no_of_steps = int(t // dt)
-    p = np.zeros((no_of_steps, *p_0.shape))
-    q = np.zeros((no_of_steps, *q_0.shape))
-    p[0] = p_0
-    q[0] = q_0
-
-    for i in trange(1, no_of_steps):
-        p[i] = p[i-1] + dt * -q[i-1]
-        q[i] = q[i-1] + dt * p[i-1]
-    
-    return p, q
+# p_0 = np.array([0])
+# q_0 = np.array([1])
+# t = 100
+# dt = 0.001
 
-def HO_IM(p_0: np.ndarray, q_0: np.ndarray, t: float, dt: float):
-    """
-    Integrator = implicit midpoint
-    p_dot = - partial H / partial q = -q
-    q_dot = partial H / partial p = p
-    """
-    no_of_steps = int(t // dt)
-    p = np.zeros((no_of_steps, *p_0.shape))
-    q = np.zeros((no_of_steps, *q_0.shape))
-    p[0] = p_0
-    q[0] = q_0
+# def HO_IE(p_0: np.ndarray, q_0: np.ndarray, t: float, dt: float):
+#     """
+#     Integrator = implicit Euler
+#     p_dot = - partial H / partial q = -q
+#     q_dot = partial H / partial p = p
+#     """
+#     no_of_steps = int(t // dt)
+#     p = np.zeros((no_of_steps, *p_0.shape))
+#     q = np.zeros((no_of_steps, *q_0.shape))
+#     p[0] = p_0
+#     q[0] = q_0
+
+#     for i in trange(1, no_of_steps):
+#         p[i] = 1/(1 + dt**2) * (p[i-1] - dt * q[i-1])
+#         q[i] = 1/(1 + dt**2) * (q[i-1] + dt * p[i-1])
+
+#     return p, q
+
+# def HO_EE(p_0: np.ndarray, q_0: np.ndarray, t: float, dt: float):
+#     """
+#     Integrator = explicit Euler
+#     p_dot = - partial H / partial q = -q
+#     q_dot = partial H / partial p = p
+#     """
+#     no_of_steps = int(t // dt)
+#     p = np.zeros((no_of_steps, *p_0.shape))
+#     q = np.zeros((no_of_steps, *q_0.shape))
+#     p[0] = p_0
+#     q[0] = q_0
+
+#     for i in trange(1, no_of_steps):
+#         p[i] = p[i-1] + dt * -q[i-1]
+#         q[i] = q[i-1] + dt * p[i-1]
+    
+#     return p, q
+
+# def HO_IM(p_0: np.ndarray, q_0: np.ndarray, t: float, dt: float):
+#     """
+#     Integrator = implicit midpoint
+#     p_dot = - partial H / partial q = -q
+#     q_dot = partial H / partial p = p
+#     """
+#     no_of_steps = int(t // dt)
+#     p = np.zeros((no_of_steps, *p_0.shape))
+#     q = np.zeros((no_of_steps, *q_0.shape))
+#     p[0] = p_0
+#     q[0] = q_0
 
-    for i in trange(1, no_of_steps):
-        q[i] = (q[i-1] + dt*(p[i-1] - dt/4 * q[i-1])) / (1 + dt**2 / 4)
-        p[i] = (p[i-1] - dt*(q[i-1] - dt/4 * p[i-1])) / (1 + dt**2 / 4)
+#     for i in trange(1, no_of_steps):
+#         q[i] = (q[i-1] + dt*(p[i-1] - dt/4 * q[i-1])) / (1 + dt**2 / 4)
+#         p[i] = (p[i-1] - dt*(q[i-1] - dt/4 * p[i-1])) / (1 + dt**2 / 4)
 
-    return p, q
+#     return p, q
 
-p, q = HO_IE(p_0, q_0, t, dt)
-plt.figure()
-plt.plot(p, q)
-plt.title("Implicit Euler")
-plt.show()
-### Observation -> Could see that phase space volume increases
+# p, q = HO_IE(p_0, q_0, t, dt)
+# plt.figure()
+# plt.plot(p, q)
+# plt.title("Implicit Euler")
+# plt.show()
+# ### Observation -> Could see that phase space volume increases
 
-p, q = HO_EE(p_0, q_0, t, dt)
-plt.figure()
-plt.plot(p, q)
-plt.title("Explicit Euler")
-plt.show()
+# p, q = HO_EE(p_0, q_0, t, dt)
+# plt.figure()
+# plt.plot(p, q)
+# plt.title("Explicit Euler")
+# plt.show()
 
-p, q = HO_IM(p_0, q_0, t, dt)
-plt.figure()
-plt.plot(p, q)
-plt.title("Implicit midpoint")
-plt.show()
-### Observation -> Conserve phase space volume better than the other 2 methods
+# p, q = HO_IM(p_0, q_0, t, dt)
+# plt.figure()
+# plt.plot(p, q)
+# plt.title("Implicit midpoint")
+# plt.show()
+# ### Observation -> Conserve phase space volume better than the other 2 methods