One version of the advection equation is
$$\partial_t q(x, t) = c \, \partial_x q(x, t)$$Let's simulate this numerically.
Discretizing this with finite differences (first order in time and second order centered in space) leads to:
$$ q_i^ {n+1} = q_i^ {n} + \frac{c \Delta t}{2 \Delta x} (q_{i+1}^n - q_{i-1}^n) $$On the edges, I used a first order discretization (I know, that's kind of silly...).
Let's implement the initial, conditions and the main loop:
import numpy as np
x = np.linspace(0, 1, num=151)
dx = (x[1] - x[0]) / x.size
dt = 0.00001
c = -1.
Courant = abs(c) * dt / dx
print(f"Courant number: {Courant}")
def f(x, x0=0.5):
return np.exp(-30 * (x - x0)**2)
q0 = f(x)
def step(q_curr, c, dt, dx):
q_next = q_curr.copy()
q_next[1:-1] += c * dt / (2 * dx) * (q_curr[2:] - q_curr[:-2])
q_next[0] += c * dt / dx * (q_curr[1] - q_curr[0])
q_next[-1] += c * dt / dx * (q_curr[-1] - q_curr[-2])
return q_next
def step_loop(q_curr, c, dt, dx):
q_next = np.zeros_like(q_curr)
for i in range(1, x.size - 1):
q_next[i] = q_curr[i] + c * dt / (2 * dx) * (q_curr[i+1] - q_curr[i-1])
# boundaries
q_next[0] = q_curr[0] + c * dt / dx * (q_curr[1] - q_curr[0])
q_next[-1] = q_curr[-1] + c * dt / dx * (q_curr[-1] - q_curr[-2])
return q_next
Now, let's run the simulation.
q_curr = q0.copy()
q_next = np.zeros_like(q_curr)
snapshots = {}
for i in range(700):
t = i * dt
if i % 10 == 0:
snapshots[t] = q_curr.copy()
q_next = step_loop(q_curr, c, dt, dx)
q_curr, q_next = q_next, q_curr
len(snapshots)
Finally, let's do an animated plot.
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
fig, ax = plt.subplots()
line, = ax.plot([], [], label='numerical')
line2, = ax.plot(x, q0, label='initial')
def setup():
ax.set_xlim((x.min(), x.max()))
ax.set_ylim((-0.2, 1.2))
ax.legend(loc='upper right')
def update(frame):
key = list(snapshots.keys())[frame]
ydata = snapshots[key]
line.set_data(x, ydata)
ax.set_title(f'time: {key:.2e}')
return line, line2
anim = FuncAnimation(fig, update, frames=range(len(snapshots)), init_func=setup)
plt.close(fig)
HTML(anim.to_jshtml())
It appears that the scheme is not stable. It seems to be stable at first but the solutions grow and explode!
To highlight this feature a little better, let's plot the maximum and minimum value over time.
fig, ax = plt.subplots()
ax.plot(list(snapshots.keys()), [vals.max() for vals in snapshots.values()], label='max')
ax.plot(list(snapshots.keys()), [vals.min() for vals in snapshots.values()], label='min')
ax.legend()
What appears in this graph is that an exponential explosion happens at longer propagation times.
Why is that?
A good explanation by Langtangen is provided here: http://hplgit.github.io/fdm-book/doc/pub/book/sphinx/._book012.html#simplest-scheme-forward-in-time-centered-in-space
How does the explanation go? We consider a plane wave solution of the equation.