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
radius = 10
f = np.zeros(I)
J = int(np.floor(radius/dx))
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.
