Update UB4

UB4/UB4.py

@@ -3,7 +3,7 @@ import matplotlib.pyplot  as plt
 
 # Ex1: Numerical integrators 
 
-def implicit_euler(f, y0: np.ndarray, t: float, dt: float) -> np.ndarray:
+def implicit_euler(f, y0: np.ndarray, t: float, dt: float, **kwargs) -> np.ndarray:
     """
     A-stable
     :param f: function to be integrated
@@ -16,12 +16,12 @@ def implicit_euler(f, y0: np.ndarray, t: float, dt: float) -> np.ndarray:
     y[0] = y0
     
     for i in range(1, no_of_steps):
-        y[i] = y[i-1] + f(y[i-1]) * dt
-        y[i] = y[i-1] + f(y[i]) * dt
+        y[i] = y[i-1] + f(y[i-1], **kwargs) * dt
+        y[i] = y[i-1] + f(y[i], **kwargs) * dt
 
     return y
 
-def explicit_euler(f, y0: np.ndarray, t: float, dt: float) -> np.ndarray:
+def explicit_euler(f, y0: np.ndarray, t: float, dt: float, **kwargs) -> np.ndarray:
     """
     :param f: function to be integrated
     :param y0: initial value
@@ -33,11 +33,11 @@ def explicit_euler(f, y0: np.ndarray, t: float, dt: float) -> np.ndarray:
     y[0] = y0
 
     for i in range(1, no_of_steps):
-        y[i] = y[i-1] + f(y[i-1]) * dt
+        y[i] = y[i-1] + f(y[i-1], **kwargs) * dt
     
     return y
 
-def implicit_midpoint(f, y0: np.ndarray, t: float, dt: float) -> np.ndarray:
+def implicit_midpoint(f, y0: np.ndarray, t: float, dt: float, **kwargs) -> np.ndarray:
     """
     Not L-stable, doesn't decay properly - oscillation
     :param f: function to be integrated
@@ -50,74 +50,122 @@ def implicit_midpoint(f, y0: np.ndarray, t: float, dt: float) -> np.ndarray:
     y[0] = y0
     
     for i in range(1, no_of_steps):
-        y[i] = y[i-1] + dt * f(y[i-1] + dt/2 * f(y[i-1]))
-        y[i] = y[i-1] + dt * f(1/2 *  (y[i-1] + y[i]))
+        y[i] = y[i-1] + dt * f(y[i-1] + dt/2 * f(y[i-1], **kwargs), **kwargs)
+        y[i] = y[i-1] + dt * f(1/2 *  (y[i-1] + y[i]), **kwargs)
     
     return y
 
 # Ex2: dy/dt = lambda * y, y(0) = 1, lambda = -1
 y0 = np.array([1])
 lamda = -1
-t = 100
-dt = 1e-2
+t = [10, 100]
+dt = [1.5, 1, 1e-1, 1e-2]
 f = lambda y: lamda * y
-y = [implicit_euler(f, y0, t, dt), explicit_euler(f, y0, t, dt), implicit_midpoint(f, y0, t ,dt)]
-integrators = ['implicit euler', 'explicit euler', 'implicit midpoint']
+exact_solution = lambda t : np.exp(lamda * t) * y0
 
-# for i in range(3):
+# A-stability: the method should produce a numerical solution that does not grow in magnitude for any dt
+# L-stability: stronger form of A-stability -> numerical solution decays (at least) as rapidly as the exact solution
+# Expectation: 
+## EE: A-stable (No), L-stable (No)
+## IE: A-stable (Yes), L-stable (Yes)
+## IM: A-stable (Yes), L-stable (No)
+
+# for i in range(len(dt)):
 #     plt.figure()
-#     plt.plot(y[i])
-#     plt.title(integrators[i])
+#     y_implicit_euler = implicit_euler(f, y0, t[0], dt[i])
+#     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()
 
-# Ex3: 
-def f(y: np.ndarray, k: float = 1) -> 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]])
-    
-    return np.dot(A, y)
-
-dt = [1, 1e-1, 1e-2, 1e-3]
-t = [10,100]
-y0 = np.array([1,0,0])
-
-# short time
-fig = plt.figure()
-for dt_i in dt:
-    y = implicit_euler(f, y0, t[0], dt_i)
-    plt.plot(y[:,0], '--',label = 'State 1')
-    plt.plot(y[:,1], 'x', label = 'State 2')
-    plt.plot(y[:,2], '>', label = 'State 3')
+#### For A-stability, dt < 1
+#### Observation
+### t = 10, dt = 1.5 <-- Strange behaviour !
+# IE grows and EE decays!
+# IM is robust
+
+### t = 10, dt = 1 <-- Strange behaviour !
+# IE gives a constant function (but it doesn't blow up as time progresses!).
+# However, this means it's not A-stable because it does not appropriately damp the solution. 
+# EE & IM decay faster than ES -> L-stable. EE gives a worse approximation in comparison to IM
+### t = 10, dt = 1e-1  <-- Strange behaviour ! 
+# They all behave similarly. They do not decay as rapidly as ES.
+### 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'{integrators[0]}, dt = {dt_i}')
+    plt.title(f't={t[1]}, dt={dt[i]}')
     plt.show()
 
+#### Observation
+### t = 100, dt = 1.5 <-- Strange behaviour !
+# IE blows up at the end.
+### t = 100, dt = 1 <-- Strange behaviour !
+# Same as in the case t = 10, dt = 1
+### t = 100, dt = 1e-1  <-- Strange behaviour ! 
+# Same as in the case t = 100, dt = 1e-1
+### t = 100, dt = 1e-2 <-- Strange behaviour !
+# Same as in the case t = 100, dt = 1e-2
 
 
 
+# Ex3: 
+# 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]])
+    
+#     return np.dot(A, y)
+
+# dt = [1, 1e-1, 1e-2, 1e-3]
+# t = [10,100]
+# y0 = np.array([1,0,0])
+
+# # short time
+# fig = plt.figure()
+# for dt_i in dt:
+#     y = implicit_euler(f, y0, t[0], 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'{integrators[0]}, dt = {dt_i}')
+#     plt.show()
 
-fig = plt.figure()
-for dt_i in dt:
-    y = explicit_euler(f, y0, t[0], dt_i)
-    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.show()
 
-fig = plt.figure()
-for dt_i in dt:
-    y = implicit_midpoint(f, y0, t[0], dt_i)
-    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.show()
+# fig = plt.figure()
+# for dt_i in dt:
+#     y = explicit_euler(f, y0, t[0], 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'{integrators[1]}, dt = {dt_i}')
+#     plt.show()
+
+
+# fig = plt.figure()
+# for dt_i in dt:
+#     y = implicit_midpoint(f, y0, t[0], 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'{integrators[2]}, dt = {dt_i}')
+#     plt.show()