from IPython.display import HTML
s="""<h1>Leapfrog, Lax-Friedrichs and Lax-Wendroff Schemes (Sinusoidal Input)</h1></br><div id="toc"></div>
<script src="https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js" defer></script>""";
h=HTML(s); h
We want 2 periods in 1m, i.e. wavelength = $ 1 \over 2 $ so the wavenumber $k = {{2 \pi} \over \lambda} = {{2 \pi} \over 0.5} = 4 \pi$
t = 0s
and x = 0m
to x = 1m
$\quad \Rightarrow \quad$ $u = 1.0*sin(4 \pi (x - ct))$x = 0m
$\Rightarrow$ u = 0m/s
x = 5m
$\Rightarrow$ u = 0m/s
The Wave Equation does have a temporal component, so we use nt
nt
= ? (number of temporal points)nx
= 51 (number of x spatial points)tmax
= 1.0xmax
= 5.0c
= 1.0Formulate nt
based on the CFL Condition sigma = c (dt / dx) = 1
$\forall (x, t) \quad x = 0 \rightarrow u(0,t) = 0 $
$\forall (x, t) \quad x = 4 \rightarrow u(5,t) = 0 $
This is simply the wave travelling in the x-direction
$$ u_i^{n+1} = u_{i-1}^n $$We must limit the number of timesteps used for the analytical solution, proportional to the value of $\sigma$ in order to compare with the numerical solution. Simply limit $t_{max}$ to $\sigma \times t_{max}$
#Constants
xmax = 5.0
nx = 51
dx = xmax/(nx-1)
tmax = 1.0
sigma = 1.0
dt = (dx / sigma)
nt = int(tmax/dt + 1)
#Boundary Conditions
for n between 0 and nt
u(0,n)=1
u(nx-1,n)=0
#Initial Conditions
for i between 1 and nx-2
if(x[i] < 1.0)
u[i,0] = 1.0*sin(4 * pi * (x[i]))
else
u[i,0] = 0
#Iteration
for n between 0 and nt-1
for i between 1 and nx-2
u[i,n+1] = u[i,n]-sigma*(u[i,n]-u[i-1,n])
def initial_and_boundary_conditions(xmax, nx, nt):
# Initialise data structures
import numpy as np
u = np.zeros((nx,nt))
# Boundary conditions
u[0,:] = 0.0
u[nx-1,:] = 0.0
# X Loop
x = np.zeros(nx)
x = np.linspace(0.0,xmax,nx)
quarter = int((nx-1)/4)
i=0
# Initial conditions
u[i:quarter,0] = 1.0*sin(4 * pi * (x[i:quarter]))
return u, x
def analytical(sigma, nx, tmax, xmax, c):
# Increments
# dt = tmax/(nt-1)
# dx = (c * dt) / sigma
dx = xmax/(nx-1)
dt = ((dx*sigma) / c)
nt = int(sigma*tmax/dt + 1)
# Initial and Boundary Conditions
u, x = initial_and_boundary_conditions(xmax, nx, nt)
i = 1
# Loop - must loop as NumPy cannot be used for dependency reasons
for n in range(0,nt-1):
u[i:nx-1, n+1] = u[i-1:nx-2, n]
return u, x, nt
def upwind_convection(sigma, nx, tmax, xmax, c):
"""
Returns the velocity field and distance for 1D linear convection
"""
# Increments
dx = xmax/(nx-1)
dt = ((dx * sigma) / c)
nt = int(tmax/dt + 1)
# Initial and Boundary Conditions
u, x = initial_and_boundary_conditions(xmax, nx, nt)
i = 1
# Loop - must loop as NumPy cannot be used for dependency reasons
for n in range(0,nt-1):
u[i:nx-1, n+1] = u[i:nx-1, n]-sigma*(u[i:nx-1, n]-u[i-1:nx-2, n])
return u, x, nt
def plot(u,x,NT,u_analytical, x_analytical, NT2):
"""
Plots the 1D velocity field
"""
import matplotlib.pyplot as plt
import matplotlib.cm as cm
plt.figure()
ax=plt.subplot(111)
colour=iter(cm.rainbow(np.linspace(0,100,NT)))
ax.plot(x,u[:,0],linestyle='-',c='b',label='n='+str(0))
ax.plot(x,u[:,NT-1],linestyle='-',c='r',label='n='+str(NT-1))
ax.plot(x_analytical,u_analytical[:,NT2-1],linestyle='--',c='k',label='n='+str(NT2-1)+' analytical')
box=ax.get_position()
ax.set_position([box.x0, box.y0, box.width*2,box.height*2])
ax.legend( bbox_to_anchor=(1.02,1), loc=2)
plt.ylim([-1.5,1.5])
plt.xlim([0.0,3.0])
plt.xlabel('x (m)')
plt.ylabel('u (m/s)')
plt.show()
u00, x00, nt00 = analytical(1.0, 501, 1.0, 5.0, 1.0)
u0, x0, nt0 = upwind_convection(1.0, 501, 1.0, 5.0, 1.0)
plot(u0,x0,nt0, u00, x00,nt00)
u11, x11, nt11 = analytical(0.5, 501, 1.0, 5.0, 1.0)
u1, x1, nt1 = upwind_convection(0.5, 501, 1.0, 5.0, 1.0)
plot(u1,x1,nt1,u11,x11,nt11)
u22, x22, nt22 = analytical(0.25, 501, 1.0, 5.0, 1.0)
u2, x2, nt2 = upwind_convection(0.25, 501, 1.0, 5.0, 1.0)
plot(u2,x2,nt2,u22, x22, nt22)
For the one first derivative of velocity in time: 2nd order CD in time
For the one first derivative of velocity in space: 2nd order CD in space
Use upwind to initialise the leapfrog scheme
#Constants
xmax = 4.0
nx = 51
dx = xmax/(nx-1)
tmax = 1.0
sigma = 1.0
dt = (dx / sigma)
nt = int(tmax/dt + 1)
#Boundary Conditions
for n between 0 and nt
u(0,n)=1
u(nx-1,n)=0
#Initial Conditions
for i between 1 and nx-2
if(i < (nx-1)/2)
u(i,0) = 1
else
u(i,0) = 0
#Initialise using Upwind
n = 0
for i between 1 and nx-2
u(i,n+1) = u(i,n)-c*(dt/dx)*(u(i,n)-u(i-1,n)
#Proceed Using Leapfrog
for n between 1 and nt-1
for i between 1 and nx-2
u(i,n+1) = u(i,n-1)-c*(dt/dx)*(u(i+1,n)-u(i-1,n)
def leapfrog_convection(sigma, nx, tmax, xmax, c):
"""
Returns the velocity field and distance for 1D linear convection
"""
# Increments
dx = xmax/(nx-1)
dt = ((dx * sigma) / c)
nt = int(tmax/dt + 1)
# Initial and Boundary Conditions
u, x = initial_and_boundary_conditions(xmax, nx, nt)
# Initialise using Upwind
n = 0
i = 1
u[i:nx-1, n+1] = u[i:nx-1, n]-sigma*(u[i:nx-1, n]-u[i-1:nx-2, n])
# Proceed using Leapfrog
for n in range(1, nt-1):
u[i:nx-1, n+1] = u[i:nx-1, n-1]-sigma*(u[i+1:nx, n]-u[i-1:nx-2, n])
return u, x, nt
u44, x44, nt44 = analytical(0.5, 501, 1.0, 5.0, 1.0)
u4, x4, nt4 = leapfrog_convection(0.5, 501, 1.0, 5.0, 1.0)
plot(u4,x4,nt4,u44, x44, nt44)
Oscillatory, behind the wave, more accurate than Lax-Friedrichs
#Constants
xmax = 4.0
nx = 51
dx = xmax/(nx-1)
tmax = 1.0
sigma = 1.0
dt = (dx / sigma)
nt = int(tmax/dt + 1)
#Boundary Conditions
for n between 0 and nt
u(0,n)=1
u(nx-1,n)=0
#Initial Conditions
for i between 1 and nx-2
if(i < (nx-1)/2)
u(i,0) = 1
else
u(i,0) = 0
#Iteration
for n between 0 and nt-1
for i between 1 and nx-2
u(i,n+1) = 0.5*(u(i-1,n)+u(i+1,n))-0.5*sigma*(u(i+1,n)-u(i-1,n)
def lax_friedrichs_convection(sigma, nx, tmax, xmax, c):
"""
Returns the velocity field and distance for 1D linear convection
"""
# Increments
dx = xmax/(nx-1)
dt = ((dx * sigma) / c)
nt = int(tmax/dt + 1)
# Initial and Boundary Conditions
u, x = initial_and_boundary_conditions(xmax, nx, nt)
i = 1
# Loop - must loop as NumPy cannot be used for dependency reasons
for n in range(0,nt-1):
u[i:nx-1, n+1] = 0.5*(u[i-1:nx-2, n]+u[i+1:nx, n])-0.5*sigma*(u[i+1:nx, n]-u[i-1:nx-2, n])
return u, x, nt
u55, x55, nt55 = analytical(0.5, 501, 1.0, 5.0, 1.0)
u5, x5, nt5 = lax_friedrichs_convection(0.5, 501, 1.0, 5.0, 1.0)
plot(u5,x5,nt5,u55, x55, nt55)
def lax_wendroff_convection(sigma, nx, tmax, xmax, c):
"""
Returns the velocity field and distance for 1D linear convection
"""
# Increments
dx = xmax/(nx-1)
dt = ((dx * sigma) / c)
nt = int(tmax/dt + 1)
# Initial and Boundary Conditions
u, x = initial_and_boundary_conditions(xmax, nx, nt)
i = 1
# Loop - must loop as NumPy cannot be used for dependency reasons
for n in range(0,nt-1):
u[i:nx-1, n+1] = (u[i:nx-1, n] - 0.5*sigma*(u[i+1:nx, n] - u[i-1:nx-2, n]) +
0.5*(sigma**2)*(u[i+1:nx, n] - 2.0*u[i:nx-1, n] + u[i-1:nx-2, n]))
return u, x, nt
u66, x66, nt66 = analytical(0.5, 501, 1.0, 5.0, 1.0)
u6, x6, nt6 = lax_wendroff_convection(0.5, 501, 1.0, 5.0, 1.0)
plot(u6,x6,nt6,u66, x66, nt66)
This more accurately represents the sinusoid
However, there is an oscillatory response
Initial condition had a non-smooth slope, which means Lax-Wendroff cannot follow this initial change
def initial_and_boundary_conditions_2(xmax, nx, nt):
# Initialise data structures
import numpy as np
u = np.zeros((nx,nt))
# Boundary conditions
u[0,:] = 0.0
u[nx-1,:] = 0.0
# X Loop
x = np.zeros(nx)
x = np.linspace(0.0,xmax,nx)
quarter = int((nx-1)/4)
i=0
# Initial conditions
u[i:quarter,0] = 1.0*sin(12.0 * pi * (x[i:quarter]))
return u, x
def analytical_2(sigma, nx, tmax, xmax, c):
# Increments
# dt = tmax/(nt-1)
# dx = (c * dt) / sigma
dx = xmax/(nx-1)
dt = ((dx*sigma) / c)
nt = int(sigma*tmax/dt + 1)
# Initial and Boundary Conditions
u, x = initial_and_boundary_conditions_2(xmax, nx, nt)
i = 1
# Loop - must loop as NumPy cannot be used for dependency reasons
for n in range(0,nt-1):
u[i:nx-1, n+1] = u[i-1:nx-2, n]
return u, x, nt
def upwind_convection_2(sigma, nx, tmax, xmax, c):
"""
Returns the velocity field and distance for 1D linear convection
"""
# Increments
dx = xmax/(nx-1)
dt = ((dx * sigma) / c)
nt = int(tmax/dt + 1)
# Initial and Boundary Conditions
u, x = initial_and_boundary_conditions_2(xmax, nx, nt)
i = 1
# Loop - must loop as NumPy cannot be used for dependency reasons
for n in range(0,nt-1):
u[i:nx-1, n+1] = u[i:nx-1, n]-sigma*(u[i:nx-1, n]-u[i-1:nx-2, n])
return u, x, nt
def plot_2(u,x,NT,u_analytical, x_analytical, NT2):
"""
Plots the 1D velocity field
"""
import matplotlib.pyplot as plt
import matplotlib.cm as cm
plt.figure()
ax=plt.subplot(111)
colour=iter(cm.rainbow(np.linspace(0,100,NT)))
ax.plot(x,u[:,0],linestyle='-',c='b',label='n='+str(0))
ax.plot(x,u[:,NT-1],linestyle='-',c='r',label='n='+str(NT-1))
ax.plot(x_analytical,u_analytical[:,NT2-1],linestyle='--',c='k',label='n='+str(NT2-1)+' analytical')
box=ax.get_position()
ax.set_position([box.x0, box.y0, box.width*2,box.height*2])
ax.legend( bbox_to_anchor=(1.02,1), loc=2)
plt.ylim([-1.5,1.5])
plt.xlim([0.0,3.0])
plt.xlabel('x (m)')
plt.ylabel('u (m/s)')
plt.show()
u0000, x0000, nt0000 = analytical_2(1.0, 501, 1.0, 5.0, 1.0)
u000, x000, nt000 = upwind_convection_2(0.5, 501, 1.0, 5.0, 1.0)
plot_2(u000, x000, nt000,u0000, x0000, nt0000)
def lax_friedrichs_convection_2(sigma, nx, tmax, xmax, c):
"""
Returns the velocity field and distance for 1D linear convection
"""
# Increments
dx = xmax/(nx-1)
dt = ((dx * sigma) / c)
nt = int(tmax/dt + 1)
# Initial and Boundary Conditions
u, x = initial_and_boundary_conditions_2(xmax, nx, nt)
i = 1
# Loop - must loop as NumPy cannot be used for dependency reasons
for n in range(0,nt-1):
u[i:nx-1, n+1] = 0.5*(u[i-1:nx-2, n]+u[i+1:nx, n])-0.5*sigma*(u[i+1:nx, n]-u[i-1:nx-2, n])
return u, x, nt
u1111, x1111, nt1111 = analytical_2(1.0, 501, 1.0, 5.0, 1.0)
u111, x111, nt111 = lax_friedrichs_convection_2(0.5, 501, 1.0, 5.0, 1.0)
plot_2(u111, x111, nt111 ,u1111, x1111, nt1111)
def lax_wendroff_convection_2(sigma, nx, tmax, xmax, c):
"""
Returns the velocity field and distance for 1D linear convection
"""
# Increments
dx = xmax/(nx-1)
dt = ((dx * sigma) / c)
nt = int(tmax/dt + 1)
# Initial and Boundary Conditions
u, x = initial_and_boundary_conditions_2(xmax, nx, nt)
i = 1
# Loop - must loop as NumPy cannot be used for dependency reasons
for n in range(0,nt-1):
u[i:nx-1, n+1] = (u[i:nx-1, n] - 0.5*sigma*(u[i+1:nx, n] - u[i-1:nx-2, n]) +
0.5*(sigma**2)*(u[i+1:nx, n] - 2.0*u[i:nx-1, n] + u[i-1:nx-2, n]))
return u, x, nt
u2222, x2222, nt2222 = analytical_2(1.0, 501, 1.0, 5.0, 1.0)
u222, x222, nt222 = lax_wendroff_convection_2(0.5, 501, 1.0, 5.0, 1.0)
plot_2(u222, x222, nt222, u2222, x2222, nt2222)
def leapfrog_convection_2(sigma, nx, tmax, xmax, c):
"""
Returns the velocity field and distance for 1D linear convection
"""
# Increments
dx = xmax/(nx-1)
dt = ((dx * sigma) / c)
nt = int(tmax/dt + 1)
# Initial and Boundary Conditions
u, x = initial_and_boundary_conditions_2(xmax, nx, nt)
# Initialise using Upwind
n = 0
i = 1
u[i:nx-1, n+1] = u[i:nx-1, n]-sigma*(u[i:nx-1, n]-u[i-1:nx-2, n])
# Proceed using Leapfrog
for n in range(1, nt-1):
u[i:nx-1, n+1] = u[i:nx-1, n-1]-sigma*(u[i+1:nx, n]-u[i-1:nx-2, n])
return u, x, nt
u3333, x3333, nt3333 = analytical_2(1.0, 501, 1.0, 5.0, 1.0)
u333, x333, nt333 = leapfrog_convection_2(0.5, 501, 1.0, 5.0, 1.0)
plot_2(u333, x333, nt333, u3333, x3333, nt3333)