From b7dd52c76e3850eee2fe3654f0e658eca65c79ff Mon Sep 17 00:00:00 2001
From: nguyed99 <nguyed99@zedat.fu-berlin.de>
Date: Wed, 22 Nov 2023 14:30:46 +0100
Subject: [PATCH] Update UB4
---
UB4/UB4.py | 64 ++++++++++++++++++++++++++++--------------------------
1 file changed, 33 insertions(+), 31 deletions(-)
diff --git a/UB4/UB4.py b/UB4/UB4.py
index 4c6ed98..1d81dfc 100644
--- a/UB4/UB4.py
+++ b/UB4/UB4.py
@@ -1,9 +1,10 @@
import numpy as np
import matplotlib.pyplot as plt
+from tqdm import trange
# Ex1: Numerical integrators
-def implicit_euler(f, y0: np.ndarray, t: float, dt: float, **kwargs) -> np.ndarray:
+def implicit_euler(f, y0: np.ndarray, t: float, dt: float, eps = 1e-5, **kwargs) -> np.ndarray:
"""
A-stable
:param f: function to be integrated
@@ -14,10 +15,15 @@ def implicit_euler(f, y0: np.ndarray, t: float, dt: float, **kwargs) -> np.ndarr
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] + f(y[i-1], **kwargs) * dt
- y[i] = y[i-1] + f(y[i], **kwargs) * dt
+
+ 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)
+
+ while abs(y_1 - y_0) > eps:
+ y_0 = y_1
+ y_1 = y[i-1] + dt * f(y_0)
+ y[i] = y_1
return y
@@ -32,14 +38,13 @@ def explicit_euler(f, y0: np.ndarray, t: float, dt: float, **kwargs) -> np.ndarr
y = np.zeros((no_of_steps, *y0.shape))
y[0] = y0
- for i in range(1, no_of_steps):
+ for i in trange(1, no_of_steps):
y[i] = y[i-1] + f(y[i-1], **kwargs) * dt
return y
-def implicit_midpoint(f, y0: np.ndarray, t: float, dt: float, **kwargs) -> np.ndarray:
+def implicit_midpoint(f, y0: np.ndarray, t: float, dt: float, eps = 1e-5, **kwargs) -> np.ndarray:
"""
- Not L-stable, doesn't decay properly - oscillation
:param f: function to be integrated
:param y0: initial value
:param t: time interval
@@ -48,18 +53,24 @@ def implicit_midpoint(f, y0: np.ndarray, t: float, dt: float, **kwargs) -> np.nd
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 * 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)
-
+
+ for i in trange(1, no_of_steps):
+ 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:
+ y_0 = y_1
+ y_1 = y[i-1] + dt * f( 1/2 * (y[i-1] + y_0) ,**kwargs)
+ y[i] = y_1
+
return y
+
# Ex2: dy/dt = lambda * y, y(0) = 1, lambda = -1
-y0 = np.array([1])
+y0 = np.array(1)
lamda = -1
t = [10, 100]
-dt = [1.5, 1, 1e-1, 1e-2]
+dt = [0.9, 1e-1, 1e-2]
f = lambda y: lamda * y
exact_solution = lambda t : np.exp(lamda * t) * y0
@@ -84,14 +95,8 @@ exact_solution = lambda t : np.exp(lamda * t) * y0
#### 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 = 0.9 <-- Strange behaviour !
+# IM decas as rapidly as ES, IE not as fast, EE too fast
### 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 !
@@ -110,17 +115,14 @@ for i in range(len(dt)):
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 !
+### t = 100, dt = 0.9 <-- Strange behaviour !
+# Same as in the case t = 10, dt = 0.9
+### t = 100, dt = 0.1 <-- Strange behaviour !
+# Same as in the case t = 100, dt = 0.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'
--
GitLab