Copyright (C) 2010-2020 Luke Olson
Copyright (C) 2020 Andreas Kloeckner
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
Here you will set up the problem for $$ u_t + c u_x = 0$$ with periodic BC on the interval [0,1]
c = 1.0
T = 1.0 / c # end time
dx
will be the grid spacing in the $x$-directionx
will be the grid coordinatesxx
will be really fine grid coordinatesnx = 82
x = np.linspace(0, 1, nx, endpoint=False)
dx = x[1] - x[0]
xx = np.linspace(0, 1, 1000, endpoint=False)
Now define an initial condition:
def f(x):
u = np.zeros(x.shape)
u[np.intersect1d(np.where(x>0.4), np.where(x<0.6))] = 1.0
return u
plt.plot(x, f(x), lw=3, clip_on=False)
[<matplotlib.lines.Line2D at 0x7fafe3a977c0>]
Have spatial grid. Now we need a time step. So define a ratio parameter $\lambda$. Let $$ \Delta t = \Delta x \frac{\lambda}{c}$$
lmbda = 0.93 * np.sign(c)
dt = dx * lmbda / abs(c)
nt = int(T/dt)
print('T = %g' % T)
print('tsteps = %d' % nt)
print(' dx = %g' % dx)
print(' dt = %g' % dt)
print('lambda = %g' % lmbda)
T = 1 tsteps = 88 dx = 0.0121951 dt = 0.0113415 lambda = 0.93
Now make an index list, called $J$, so that we can access $J+1$ and $J-1$ easily
J = np.arange(0, nx - 1) # all vertices
Jm1 = np.roll(J, 1)
Jp1 = np.roll(J, -1)
import time
method = 'ETBS'
plotit = True
u = f(x)
fig = plt.figure(figsize=(10,10))
plt.title('u vs x')
line1, = plt.plot(x, u, lw=3, clip_on=False)
line2, = plt.plot(x, u, lw=3, clip_on=False)
def timestepper(n):
if method == 'ETBS':
u[J] = u[J] - lmbda * (u[J] - u[Jm1])
if method == "ETFS":
u[J] = u[J] - lmbda * (u[Jp1] - u[J])
if method == "ETCS":
u[J] = u[J] - lmbda * (1.0 / 2.0) * (u[Jp1] - u[Jm1])
uex = f((xx - c * (n+1) * dt) % 1.0)
line1.set_data(xx, uex)
line2.set_data(x, u)
return line1, line2
from matplotlib import animation
from IPython.display import HTML
ani = animation.FuncAnimation(fig, timestepper, frames=nt, interval=30)
html = HTML(ani.to_jshtml())
plt.clf()
html