From e7695548c941e95657d6fc75b9f0755e91560fd8 Mon Sep 17 00:00:00 2001
From: nguyed99 <nguyed99@zedat.fu-berlin.de>
Date: Tue, 12 Dec 2023 11:51:07 +0100
Subject: [PATCH] Add UB5
---
UB5/UB5.py | 44 ++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 44 insertions(+)
create mode 100644 UB5/UB5.py
diff --git a/UB5/UB5.py b/UB5/UB5.py
new file mode 100644
index 0000000..7ce0279
--- /dev/null
+++ b/UB5/UB5.py
@@ -0,0 +1,44 @@
+"""
+ODE
+\partial c / \partial t - D \partial**2 / \partial x**2 c = 0
+
+Central finite differences in space
+Explicit Euler discretization in time
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+
+T = 50.0
+
+def time_space_integrator(c0: np.ndarray, dt: float, dx: float, no_time_steps: int, D: float = 1) -> np.ndarray:
+
+ no_of_sites = c0.shape[0]
+ t = np.zeros(no_time_steps)
+ c_t = np.zeros((no_time_steps, no_of_sites))
+ c_t[0] = c0
+
+ Id = np.identity(no_of_sites)
+ B = 2 * np.identity(no_of_sites) + (-1) * np.diag([1] * (no_of_sites-1), k = -1) + (-1) * np.diag([1] * (no_of_sites-1), k = 1)
+ B[-1][0] = -1
+ B[0][-1] = -1
+ M = Id - dt * (D / dx**2) * B
+ for i in range(1, no_time_steps):
+ t[i] += i*dt
+ c_t[i] = np.dot(M, c_t[i-1])
+
+ return t, c_t
+
+c0 = np.zeros(51)
+dx = 5e-2
+dt = 1e-1
+C0 = 3.14
+c0[51//2] = C0 * (1/dx)
+
+t, c_t = time_space_integrator(c0, dt, dx, int(T/dt), D=0.05)
+plt.figure()
+for i in range(0, int(T/dt), 100):
+ print(i)
+ plt.plot(c_t[i], label = f"t = {t[i]}")
+plt.legend()
+plt.show()
+
--
GitLab