Skip to content
Snippets Groups Projects
Commit e7695548 authored by nguyed99's avatar nguyed99
Browse files

Add UB5

parent 6b62c70c
No related branches found
No related tags found
No related merge requests found
"""
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()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment