diff --git a/UB4/UB4.py b/UB4/UB4.py
index 1d81dfccd2463fdb96a4d70c2fdc2629f144456c..a31f962baf1f002d6ca1d110c83097550a9d30c3 100644
--- a/UB4/UB4.py
+++ b/UB4/UB4.py
@@ -1,6 +1,8 @@
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
+
+
+