Python Lax-Friedrichs and central scheme programming

Question

Good day :)

I found this exercise:

Consider $q_t + aq_x = 0$ with a constant velocity $a$.

Write a computer program that implements this equation using both the central scheme and also the Lax-Friedrichs scheme.

Use periodic boundary conditions. The initial data is a box followed by a Gaussian. Check the instabiliy of the central scheme for varaious CFL numbers. Check the dissipation of the Lax-Friedrichs scheme for large times.

Tryed to do it with Python and this is my code:

import numpy as np
import matplotlib.pyplot as plt

N = 100
t_min = 0
t_max = 1
k = (t_max-t_min)/N
x_min = 0
x_max = 1
a = 1
h = (x_max-x_min)/N
alpha = a*k/(2*h)
x = np.arange(x_min, x_max+2*h, h)
u_0 = np.exp(-225*(x-0.3)**2)
u_0[np.where((x>=0.6) & (x <= 0.8))] = 1.0
u_lf = u_0.copy()
u_cs = u_0.copy()
u_num1 = u_0.copy()
u_num2 = u_0.copy()
t = 10000

plt.clf()

for j in range(N):
    u_num1[j] = (u_lf[j+1]+u_lf[j-1])/2 - alpha*(u_lf[j+1]-u_lf[j-1])
    u_num2[j] = u_cs[j] - alpha*(u_cs[j+1]-u_cs[j-1])
u_lf = u_num1.copy()
u_cs = u_num2.copy()

u_ex = np.exp(-225*(x-0.3-(a*t)%t)**2)
u_ex[np.where((x>=0.6) & (x <= 0.8))] = 1.0

u_ex=u_ex%len(u_ex)
u_lf=u_lf%len(u_lf)
u_cs=u_cs%len(u_cs)

plt.plot(x, u_ex, 'r-', label="Exact solution", fillstyle='none')
plt.plot(x, u_lf, 'o', label="Lax-Friedrichs", fillstyle='none')
plt.plot(x, u_cs, '^', label="Central scheme", fillstyle='none')
plt.axis((0, 1, -.5, 1.5))
plt.legend(loc='lower right')
plt.suptitle("t = %1d" % (t))

plt.show()

I think there are some errors which I don't see, since the central scheme is perfectly stable and the LF scheme is not dissipating for any time.

Answer 1

With the above choice of parameters, we have $a k/h = 1$, which is a very particular case for which some numerical schemes provide the exact solution. This property is sometimes called the unit CFL condition, see Exercise 4.2 p. 85 of the monograph by LeVeque (2002). Besides this feature, OP uses the same number of points $N$ both in space and time (there is no reason to do so in practice). Moreover, OP forgot to implement the periodic boundary conditions $$ \begin{aligned} Q_{0}^{n} &= Q_{N-1}^{n} \\ Q_{N}^{n} &= Q_{1}^{n} \end{aligned} $$ see this post for instance, and other related posts on this site.