In [1]:
import numpy as np
import matplotlib.pyplot as plt

Advection equation:

$$ \partial_t u(x, t) + a \partial_x u(x, t) = 0 $$

Boundary conditions: 0

Initial conditions: $ u(x, 0) = \psi(x) $

Analytical solution: $$ u(x, t) = \psi(x - at) $$

In [5]:
L = 1
a = 1

xx = np.linspace(0, L, 1000)
tt = np.linspace(0, 10, 100)

def adv_eq_solution(x, t, a, ini=lambda x: 0, bnd=lambda t: 0):
    return np.where(a*t <= x, ini(x-a*t), bnd(t-x/a))

def my_ini(x):
    return np.where((x >= 0.1) & (x < 0.3), 1, 0)

uu = adv_eq_solution(xx, 0.5, a, ini=my_ini)

plt.plot(xx, uu)
plt.show()

Forward difference (time), central difference (space)¶

In [21]:
Nx = 50
dt = 0.001

dx = L/Nx
Nt = 1000
Tmax = dt*Nt

xnum = np.linspace(0, L, Nx)
unum = np.zeros((Nx, Nt))
unum[:, 0] = my_ini(xnum)
$$ \dfrac{u(x_i, t_{j+1}) - u(x_i, t_j)}{\Delta t} + a \dfrac{u(x_{i+1}, t_j) - u(x_{i-1}, t_j)}{2\Delta x} = 0 $$
$$ u(x_i, t_{j+1}) = u(x_i, t_j) - \dfrac{a \Delta t}{2\Delta x} \left[u(x_{i+1}, t_j) - u(x_{i-1}, t_j) \right] $$$$ C = \dfrac{a \Delta t}{2\Delta x} $$
In [22]:
C = a*dt/(2*dx)
print('C = ', C)

C =  0.025
In [23]:
for j in range(0, Nt - 1):
    for i in range(1, Nx-1):
        unum[i, j+1] = unum[i,j] - C*(unum[i+1,j] - unum[i-1,j])
In [28]:
tstep = 100
t1 = tstep*dt

uu = adv_eq_solution(xx, t1, a, ini=my_ini)

plt.plot(xx, uu, 'b-', xnum, unum[:, tstep], 'ro-')
plt.show()

Forward difference (time), left difference (space): upwind scheme¶

$$ \dfrac{u(x_i, t_{j+1}) - u(x_i, t_j)}{\Delta t} + a \dfrac{u(x_{i}, t_j) - u(x_{i-1}, t_j)}{\Delta x} = 0 $$
$$ u(x_i, t_{j+1}) = u(x_i, t_j) - \dfrac{a \Delta t}{\Delta x} \left[u(x_{i}, t_j) - u(x_{i-1}, t_j) \right] $$$$ C_1 = \dfrac{a \Delta t}{\Delta x} = 2C $$
In [49]:
Nx = 500
dt = 0.001

dx = L/Nx
Nt = 1000
Tmax = dt*Nt

xnum = np.linspace(0, L, Nx)
unum = np.zeros((Nx, Nt))
unum[:, 0] = my_ini(xnum)

C1 = a*dt/(dx)
print('C1 = ', C1)

for j in range(0, Nt - 1):
    for i in range(1, Nx-1):
        unum[i, j+1] = unum[i,j] - C1*(unum[i,j] - unum[i-1,j])

C1 =  0.5
In [50]:
tstep = 500
t1 = tstep*dt

uu = adv_eq_solution(xx, t1, a, ini=my_ini)

plt.plot(xx, uu, 'b-', xnum, unum[:, tstep], 'ro-')
plt.title('Time: {}'.format(t1))
plt.show()
In [ ]: