Update UB4

UB4/UB4.py

 import numpy as np
 import matplotlib.pyplot  as plt
 from tqdm import trange
+from typing import Callable
+
 
 # Ex1: Numerical integrators 
 
@@ -17,12 +19,12 @@ def implicit_euler(f, y0: np.ndarray, t: float, dt: float, eps = 1e-5, **kwargs)
     y[0] = y0
 
     for i in trange(1, no_of_steps):
-        y_0 = y[i-1] + dt * f(y[i-1]) # predicted state obtained through EE
-        y_1 = y[i-1] + dt * f(y_0)
+        y_0 = y[i-1] + dt * f(y[i-1], **kwargs) # predicted state obtained through EE
+        y_1 = y[i-1] + dt * f(y_0, **kwargs)
 
-        while abs(y_1 - y_0) > eps:
+        while np.linalg.norm(y_1 - y_0) > eps:
             y_0 = y_1
-            y_1 = y[i-1] + dt * f(y_0)
+            y_1 = y[i-1] + dt * f(y_0, **kwargs)
         y[i] = y_1
 
     return y
@@ -58,7 +60,7 @@ def implicit_midpoint(f, y0: np.ndarray, t: float, dt: float, eps = 1e-5, **kwar
         y_0 = y[i-1] + dt * f(y[i-1] + dt/2 * f(y[i-1], **kwargs), **kwargs) # predicted state obtained through IM
         y_1 = y[i-1] + dt * f( 1/2 * (y[i-1] + y_0) ,**kwargs)
         
-        while abs(y_1 - y_0) > eps:
+        while np.linalg.norm(y_1 - y_0) > eps:
             y_0 = y_1
             y_1 = y[i-1] + dt * f( 1/2 * (y[i-1] + y_0) ,**kwargs)
         y[i] = y_1
@@ -102,17 +104,17 @@ exact_solution = lambda t : np.exp(lamda * t) * y0
 ### t = 10, dt = 1e-2 <-- Strange behaviour !
 # They all behave similarly. They do not decay as rapidly as ES.
 
-for i in range(len(dt)):
-    plt.figure()
-    y_implicit_euler = implicit_euler(f, y0, t[1], dt[i])
-    plt.plot(implicit_euler(f, y0, t[1], dt[i]), '--', label = "implicit euler (IE)")
-    plt.plot(explicit_euler(f, y0, t[1], dt[i]), 'x', label = "explicit euler (EE)")
-    plt.plot(implicit_midpoint(f, 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[1]}, dt={dt[i]}')
-    plt.show()
+# for i in range(len(dt)):
+#     plt.figure()
+#     y_implicit_euler = implicit_euler(f, y0, t[1], dt[i])
+#     plt.plot(implicit_euler(f, y0, t[1], dt[i]), '--', label = "implicit euler (IE)")
+#     plt.plot(explicit_euler(f, y0, t[1], dt[i]), 'x', label = "explicit euler (EE)")
+#     plt.plot(implicit_midpoint(f, 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[1]}, dt={dt[i]}')
+#     plt.show()
 
 #### Observation
 ### t = 100, dt = 0.9  <-- Strange behaviour ! 
@@ -124,50 +126,212 @@ for i in range(len(dt)):
 
 
 # Ex3: 
-# def f(y: np.ndarray, k: float) -> np.ndarray:
-#     assert y.shape[0] == 3, 'y has the wrong dimension. It should be 3'
+def f(y: np.ndarray, k: float) -> np.ndarray:
+    assert y.shape[0] == 3, 'y has the wrong dimension. It should be 3'
 
-#     A = np.array([[-1, 1, 0],
-#                     [1, -1-k, k],
-#                     [0, k, -k]])
+    A = np.array([[-1, 1, 0], # S1 -> S2
+                  [1, -1-k, k], # S2 -> S1, S2 -> S3
+                  [0, k, -k]]) # S3 -> S2
     
-#     return np.dot(A, y)
+    return np.dot(A, y)
+
+dt = [0.5, 1e-1, 1e-2]
+t = 10
+y0 = np.array([1,0,0])
+
+### short time
+# fig = plt.figure()
+# for dt_i in dt:
+#     y = implicit_euler(f, y0, t, dt_i, k=1)
+#     plt.plot(y[:,0], '--',label = 'State 1')
+#     plt.plot(y[:,1], 'x', label = 'State 2')
+#     plt.plot(y[:,2], '>', label = 'State 3')
+#     plt.legend()
+#     plt.title(f"IE, dt = {dt_i}")
+#     plt.show()
+
+
+# fig = plt.figure()
+# for dt_i in dt:
+#     y = explicit_euler(f, y0, t, dt_i, k=1)
+#     plt.plot(y[:,0], '--', label = 'State 1')
+#     plt.plot(y[:,1], 'x', label = 'State 2')
+#     plt.plot(y[:,2], '>', label = 'State 3')
+#     plt.legend()
+#     plt.title(f"EE, dt = {dt_i}")
+#     plt.show()
+
+
+# fig = plt.figure()
+# for dt_i in dt:
+#     y = implicit_midpoint(f, y0, t, dt_i, k=1)
+#     plt.plot(y[:,0], '--', label = 'State 1')
+#     plt.plot(y[:,1], 'x', label = 'State 2')
+#     plt.plot(y[:,2], '>', label = 'State 3')
+#     plt.legend()
+#     plt.title(f"IM, dt = {dt_i}")
+#     plt.show()
+
 
-# dt = [1, 1e-1, 1e-2, 1e-3]
-# t = [10,100]
-# y0 = np.array([1,0,0])
+#### Observation
+# Diverging behaviour for dt in [0.9, 0.7] for EE
+# Funny behaviour fot dt in [0.9, 0.45] for IE (only 1 point)
+# Funny behaviour for dt in [0.9, 0.7] for IM (only 1 point)
 
-# # short time
+### k = 1e-3
 # fig = plt.figure()
 # for dt_i in dt:
-#     y = implicit_euler(f, y0, t[0], dt_i, k=1)
+#     y = implicit_euler(f, y0, t, dt_i, k=1e-3)
 #     plt.plot(y[:,0], '--',label = 'State 1')
 #     plt.plot(y[:,1], 'x', label = 'State 2')
 #     plt.plot(y[:,2], '>', label = 'State 3')
 #     plt.legend()
-#     plt.title(f'{integrators[0]}, dt = {dt_i}')
+#     plt.title(f"IE, dt = {dt_i}")
 #     plt.show()
 
 
 # fig = plt.figure()
 # for dt_i in dt:
-#     y = explicit_euler(f, y0, t[0], dt_i, k=1)
+#     y = explicit_euler(f, y0, t, dt_i, k=1e-3)
 #     plt.plot(y[:,0], '--', label = 'State 1')
 #     plt.plot(y[:,1], 'x', label = 'State 2')
 #     plt.plot(y[:,2], '>', label = 'State 3')
 #     plt.legend()
-#     plt.title(f'{integrators[1]}, dt = {dt_i}')
+#     plt.title(f"EE, dt = {dt_i}")
 #     plt.show()
 
 
 # fig = plt.figure()
 # for dt_i in dt:
-#     y = implicit_midpoint(f, y0, t[0], dt_i, k=1)
+#     y = implicit_midpoint(f, y0, t, dt_i, k=1e-3)
 #     plt.plot(y[:,0], '--', label = 'State 1')
 #     plt.plot(y[:,1], 'x', label = 'State 2')
 #     plt.plot(y[:,2], '>', label = 'State 3')
 #     plt.legend()
-#     plt.title(f'{integrators[2]}, dt = {dt_i}')
+#     plt.title(f"IM, dt = {dt_i}")
 #     plt.show()
 
 
+### different init cond
+y0 = np.array([0.5,0.2,0.3])
+
+### k = 1e-3
+# fig = plt.figure()
+# for dt_i in dt:
+#     y = implicit_euler(f, y0, t, dt_i, k=1e-3)
+#     plt.plot(y[:,0], '--',label = 'State 1')
+#     plt.plot(y[:,1], 'x', label = 'State 2')
+#     plt.plot(y[:,2], '>', label = 'State 3')
+#     plt.legend()
+#     plt.title(f"IE, dt = {dt_i}")
+#     plt.show()
+
+
+# fig = plt.figure()
+# for dt_i in dt:
+#     y = explicit_euler(f, y0, t, dt_i, k=1e-3)
+#     plt.plot(y[:,0], '--', label = 'State 1')
+#     plt.plot(y[:,1], 'x', label = 'State 2')
+#     plt.plot(y[:,2], '>', label = 'State 3')
+#     plt.legend()
+#     plt.title(f"EE, dt = {dt_i}")
+#     plt.show()
+
+
+# fig = plt.figure()
+# for dt_i in dt:
+#     y = implicit_midpoint(f, y0, t, dt_i, k=1e-3)
+#     plt.plot(y[:,0], '--', label = 'State 1')
+#     plt.plot(y[:,1], 'x', label = 'State 2')
+#     plt.plot(y[:,2], '>', label = 'State 3')
+#     plt.legend()
+#     plt.title(f"IM, dt = {dt_i}")
+#     plt.show()
+
+
+# Ex4: HO
+## A. Integrator
+# 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
+
+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)
+
+    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_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
+
+
+