# WENO Scheme for 1D linear advection equation

I have been stuck for a while now on a WENO scheme implementation and I have a few questions concerning the implementation, but also the theory behind it. For context: I'm currently working on numerical schemes for kinetic equations. Someone suggested to me to look at WENO schemes. However, I somehow am not able to implement it in a correct way. Therefore, as an exercise, I wanted to implement it for the 1D linear advection equation: \begin{align*} f_t+af_x=0 \end{align*} To do so, I have implemented 5th-order WENO with a simple Euler time-integration. Many WENO schemes use higher time-integrations, such as RK methods, but I think that this shouldn't play a role here. Indeed, in the specific case of this equation, RK-2 for example is merely an Euler method with time-step divided in half (if I'm not mistaken). Hence the schemes becomes: \begin{align*} f_i^{n+1}=f_i^n-a\frac{\Delta t}{\Delta x}\left( \hat{f}_{i+\frac{1}{2}}^n-\hat{f}_{i-\frac{1}{2}}^n\right) \end{align*} where: \begin{align*} &\hat{f}_{i+\frac{1}{2}}^n = w_1 \hat{f}_{i+\frac{1}{2}}^{n,(1)}+w_2 \hat{f}_{i+\frac{1}{2}}^{n,(2)}+w_3 \hat{f}_{i+\frac{1}{2}}^{n,(3)}\\ &\hat{f}_{i+\frac{1}{2}}^{n,(1)}=\frac{1}{6}\left(2f_{i-2}^n-7f_{i-1}^n+11f_{i}^n \right)\\ &\hat{f}_{i+\frac{1}{2}}^{n,(2)}=\frac{1}{6}\left(-f_{i-1}^n+5f_{i}^n+2f_{i+1}^n \right)\\ &\hat{f}_{i+\frac{1}{2}}^{n,(3)}=\frac{1}{6}\left(2f_{i}^n+5f_{i+1}^n-f_{i+2}^n \right) \end{align*} Now, my first question concerns the weights to compute the values at $$i+\frac{1}{2}$$ on each substencil. From what, I understand, these are the weights that one yields through polynomial interpolation. However, in some papers, I have read, that these weights approximate the value of the function at $$i+\frac{1}{2}$$ up to order 3. But when I do a Taylor expansion with these weights, I obtain something different (for example with the second formula): \begin{align*} &\frac{1}{6}\left( -f\left(x-\frac{3\Delta x}{2}\right)+5f\left(x-\frac{\Delta x}{2}\right) +2f\left(x+\frac{\Delta x}{2}\right)\right)\\ =&f(x)-\frac{1}{24}f''(x)\Delta x^2+\mathcal{O}(x^3) \end{align*} In fact, for the other formulas I obtain the same Taylor expansion. Am I missing something? Is this normal? Am I really bad at computing Taylor expansion? Or is this correct and what is then the interpretation of this second-order term?

Once I have implemented the code (which can be found below), even on very smooth initial data, after some time I obtain an instability (in the example below, I have used a Gaussian). In order, to track down the error, I have used the linear weights in the reconstruction, i.e. $$w_1,w_2,w_3=\frac{1}{10},\frac{6}{10},\frac{3}{10}$$ instead of the non-linear weights (which should anyhow be close to these values, since the data are smooth). Also, close to the boundaries, I have used simple UW values, in order to not be bothered with boundary conditions.

The parameters I have taken is $$a=1,T=10,L=100,\Delta x=0.1, \Delta t=0.05$$ (hence CFL $$0.5$$). And I obtain the following results ( Red, $$t=0$$; Green, $$t=5$$; Orange, $$t=6$$; Blue, $$t=7$$): I also include my (Python) code:

import numpy as np
import matplotlib.pyplot as plt

# Domain
L = 100
T = 10
# Parameters
a = 1
dx = 0.1
dt = dx/(2*a)
eps = 10**(-6)
# Discretization of the grid
nbsteps = int(np.floor(T/dt))
I = int(np.floor(L/dx))
# Setting initial data
f = np.zeros(I)
for i in range(2*J):
f[i] = (np.exp(-20*((i-J)/J)**2))
# Saving computed data
nbframes = 50
everynthframe = int(np.ceil(nbsteps/nbframes))
countframe = 0
f_data = np.zeros((nbframes,I))

for t in range(nbsteps):
# Saving data
if t%everynthframe==0:
time = countframe*everynthframe*dt
f_data[countframe,:] = f
countframe += 1
# Scheme for f
fhat = np.zeros(I)
f_x = np.zeros(I)
for i in range(2,I-2):
v1 = f[i-2]
v2 = f[i-1]
v3 = f[i]
v4 = f[i+1]
v5 = f[i+2]
S1 = (13/12*(v1-2*v2+v3)**2+1/4*(v1-4*v2+3*v3)**2)*a**2
S2 = (13/12*(v2-2*v3+v4)**2+1/4*(v2-v4)**2)*a**2
S3 = (13/12*(v3-2*v4+v5)**2+1/4*(3*v3-4*v4+v5)**2)*a**2
a1 = 1/(10*(eps+S1)**2)
a2 = 6/(10*(eps+S2)**2)
a3 = 3/(10*(eps+S3)**2)
w1 = a1/(a1+a2+a3)
w2 = a2/(a1+a2+a3)
w3 = a3/(a1+a2+a3)
#fhat[i]=fhat_{i+1/2}
#fhat[i] = w1*(v1/3-7*v2/6+11*v3/6)+w2*(-v2/6+5*v3/6+v4/3)+w3*(v3/3+5*v4/6-v5/6)
fhat[i] = 1/10*(v1/3-7*v2/6+11*v3/6)+6/10*(-v2/6+5*v3/6+v4/3)+3/10*(v3/3+5*v4/6-v5/6)
for i in range(3,I-2):
f_x[i] = (fhat[i]-fhat[i-1])/dx
# For cells close to the boundary, we use UW flux and NBC
i = 0
f_x[i] = 0
i = 1
f_x[i] = (-f[i-1]+f[i])/(dx)
i = 2
f_x[i] = (-f[i-1]+f[i])/(dx)
i = I-2
f_x[i] = (-f[i-1]+f[i])/(dx)
i = I-1
f_x[i] = (-f[i-1]+f[i])/(dx)
#Update f
f -= a*dt*f_x


My final question concerns another aspect of WENO schemes, which is not really related to the issue shown above. But especially in the application of WENO schemes to HJ Equations, I have noticed that in some cases, instead of computing $$\hat{f}_{i+\frac{1}{2}}^n$$, they computed $$\hat{f}_{x,i}^n$$ by taking as input values $$\frac{f_{i+k}-f_{i+k-1}}{\Delta x}$$ for $$k=-2,-1,0,1,2$$ (in the upwind case). So in that case, the scheme becomes: \begin{align*} f_i^{n+1}=f_i^n-a\Delta t \hat{f}_{x,i}^n \end{align*} In the case of the linear weights, this is obviously equivalent, because of linearity, but in the general case, these operations do not commute. Do you have an opinion, on which approach is better?

EDIT : I also include the standard WENO scheme (with non-linear weights). Here the scheme seems to be stable, but the profiles are everything, but satisfactory. • I wouldn't be so sure you have the right time stepper or time step size to ensure stability here. Mar 26 at 13:12
• Hmm, interesting! This is in fact a point, I'm not completely familiar with! I noticed that if I reduce my time step size, the instability appears later... But this isn't too surprising I guess. I have not tried any other time stepper. I still think my observation on RK-2 is valid. Furthermore, higher-order RK methods from what I understand have the property to be convex combinations of Euler methods (which is the key ingredient to show stability). Hence I'd find it very surprising that Euler methods lead to instability. But I'm not an expert, so perhaps you are right. Could you maybe develop? Mar 26 at 13:31
• @SpencerBryngelson You were in fact right! Thank you very much! I implemented the code with RK-3 and everything works wonderfully! And my comment on RK-2 was nonsense too! However, I still am a bit confused about the problem with Euler. If by any chance, you know why this method fails, I'd be very interested. All the best! Mar 26 at 15:10
• Can you leave the full code with graphs? May 8 at 12:33

As @SpencerBryngelson has pointed out, the scheme doesn't work, because of the Euler time integration method. One should rather use RK-3 (or maybe RK-2) (see for instance these excellent notes: https://www.brown.edu/research/projects/scientific-computing/sites/brown.edu.research.projects.scientific-computing/files/uploads/Essentially%20non-oscillatory%20and%20weighted.pdf). Furthermore, if anyone can explain why the Euler time integration method doesn't work, that'd be lovely. I also join the new code, for people that would have problems like me:

# Domain
L = 100
T = 70
# Parameters
a = 1
dx = 0.1
dt = dx/(2*a)
eps = 10**(-6)
# Discretization of the grid
nbsteps = int(np.floor(T/dt))
I = int(np.floor(L/dx))
# Setting initial data
f = np.zeros(I)
for i in range(2*J):
f[i] = (np.exp(-20*((i-J)/J)**2))
# Saving computed data
nbframes = 50
everynthframe = int(np.ceil(nbsteps/nbframes))
countframe = 0
f_data = np.zeros((nbframes,I))

def compute_f_x_pos(f,I,dx,eps,a):
fhat = np.zeros(I+1)
f_x = np.zeros(I)
for i in range(I+1):
#In the main case, v1=f[i-3],v2=f[i-2],v3=f[i-1],v4=f[i],v5=f[i+1]
if i==0:
v1 = f
elif i==1:
v1 = f
elif i==2:
v1 = f
else:
v1 = f[i-3]
if i==0:
v2 = f
elif i==1:
v2 = f
else:
v2 = f[i-2]
if i==0:
v3 = f
else:
v3 = f[i-1]
if i==I:
v4 = f[I-1]
else:
v4 = f[i]
if i==I:
v5 = f[I-2]
elif i==I-1:
v5 = f[I-1]
else:
v5 = f[i+1]
S1 = (13/12*(v1-2*v2+v3)**2+1/4*(v1-4*v2+3*v3)**2)*a**2
S2 = (13/12*(v2-2*v3+v4)**2+1/4*(v2-v4)**2)*a**2
S3 = (13/12*(v3-2*v4+v5)**2+1/4*(3*v3-4*v4+v5)**2)*a**2
a1 = 1/(10*(eps+S1)**2)
a2 = 6/(10*(eps+S2)**2)
a3 = 3/(10*(eps+S3)**2)
w1 = a1/(a1+a2+a3)
w2 = a2/(a1+a2+a3)
w3 = a3/(a1+a2+a3)
#fhat[i]=fhat_{i-1/2}
fhat[i] = w1*(v1/3-7*v2/6+11*v3/6)+w2*(-v2/6+5*v3/6+v4/3)+w3*(v3/3+5*v4/6-v5/6)
for i in range(I):
f_x[i] = (fhat[i+1]-fhat[i])/dx
return(f_x)

for t in range(nbsteps):
# Saving data
if t%everynthframe==0:
time = countframe*everynthframe*dt
f_data[countframe,:] = f
countframe += 1
# Scheme for f
f_x = compute_f_x_pos(f, I, dx, eps, a)
f1 = f-a*dt*f_x
f1_x = compute_f_x_pos(f1, I, dx, eps, a)
f2 = 3/4*f+1/4*f1-a*dt/4*f1_x
f2_x = compute_f_x_pos(f2, I, dx, eps, a)
f = 1/3*f+2/3*f2-a*dt*2/3*f2_x
$$$$
`
• Chi-Wang Shu is always right, this is the first law of high-order reconstruction schemes. Also, your comment about RK2 just being forward Euler with half the time step size is incorrect. Mar 26 at 17:24
• For why Euler is unstable, see here: dx.doi.org/10.1137/050637868 Mar 26 at 17:26

Turning my comments into an answer. Forward Euler is unstable for fifth-order-accurate WENO. Referenced here: dx.doi.org/10.1137/050637868. The key quote is in the abstract

... the combination of the widely used fifth-order WENO spatial discretization (WENO5) and the forward Euler time integration method is linearly unstable when numerically integrating hyperbolic conservation laws. Consequently it is not convergent. Furthermore we show that all two-stage, second-order explicit Runge–Kutta (ERK) methods are linearly unstable (and hence do not converge) when cou-pled with WENO5. We also show that all optimal first- and second-order strong-stability-preserving (SSP) ERK methods are linearly unstable when coupled with WENO5. Moreover the popular three-stage, third-order SSP(3,3) ERK method offers no linear stability advantage over non-SSP ERK methods, including ones with negative coefficients, when coupled with WENO5.

What makes a method stable, you ask?

For linear stability of an ERK time integration method coupled with WENO5, we show that it is sufficient that the classical linear stability region of the ERK method include a piece of the imaginary axis.