# Why this error occurs in my code for Lax Wendroff?

I want to implement the Lax Wendroff method for a non linear advection equation which is

$$\frac{u_{i}^{n+1}-u_{i}^{n}}{t} + \frac{f(u_{i+1}^{n})-f(u_{i-1}^{n}) }{2h} -\frac{t}{2h} \left( F_{i+1/2}^{n} \frac{f(u_{i+1}^{n})-f(u_{i})^{n}}{h}) - F_{i-1/2}^{n} \frac{f(u_{i}^{n})-f(u_{i-1}^{n})}{h} \right)=0$$

for the problem:$$u_{t} +u u_{x} = 0$$ for $$t \in (0,1)$$,

$$x \in (0,1)$$

$$u(0,x) = \exp(-10(4x-1)^{2})$$

import numpy as np
def NLLaxWendroff(Nx,Nt,t_end):
Nx = Nx+1                         # Number of grid points
xmax = 1.                         # Domain limit to the right
xmin = 0.                         # Domain limit to the left
h = (xmax-xmin)/(Nx-1)           # Mesh size
x = np.arange(xmin,xmax,h)       # Discretized mesh

U =  np.exp(-10*(4*x-1)**2)       # Initial solution

dt = 1.0/Nt                        # Time step
t = 0.                            # Initial time

# upwind scheme (7.18) from Dr.Georgoulis Notes after rearrangement (solving for u_{i}^{n+1})
while (t <= t_end):
t = t+dt
v = t/(2*h)
a = (t**2)/(2*h)
Un = U
Um = np.roll(Un,1)
Up = np.roll(Un,-1)
F_plus  = 0.5 * (Up + Un)
F_minus = 0.5 * (Un + Um)
U =  Un - v * (Um-Up) + a * (F_plus  * Un-Um - F_minus *Um-Up)

plt.plot(x,U)
plt.title('t='+str(round(t_end,2)),fontsize=16)
plt.xlabel('x',fontsize=18)
plt.ylabel('u',fontsize=18)
plt.show()
NLLaxWendroff(100,100,0.6)


But the following error occurs:

<ipython-input-55-f6e78f745ea8>:24: RuntimeWarning: overflow encountered in multiply
U =  Un - (v * (Up-Um)) + a * (  (F_plus  * (Un-Um)) - (F_minus *(Up-Um))   )
<ipython-input-55-f6e78f745ea8>:24: RuntimeWarning: invalid value encountered in subtract
U =  Un - (v * (Up-Um)) + a * (  (F_plus  * (Un-Um)) - (F_minus *(Up-Um))   )
<ipython-input-55-f6e78f745ea8>:22: RuntimeWarning: invalid value encountered in subtract
F_plus  = 1/2 * (Un - Um)
<ipython-input-55-f6e78f745ea8>:23: RuntimeWarning: invalid value encountered in subtract
F_minus = 1/2 * (Un - Up)
<ipython-input-55-f6e78f745ea8>:24: RuntimeWarning: invalid value encountered in add
U =  Un - (v * (Up-Um)) + a * (  (F_plus  * (Un-Um)) - (F_minus *(Up-Um))   )

$$$$


I found a few issues with your implmentation:

• you should replace $$t$$ by $$dt$$ in your discrete equation (both in your code and in the question), otherwise it makes no sense !
• the second part of your equation (with $$a$$) seems incomplete or wrong, and there's a mix between $$F$$ and $$f$$. Moreover, you flux evaluation is wrong: Fplus should read 0.5*Up**2 for example.
• there are also missing brackets around Un-Um and Um-Up in the computation of U.

I've implemented the Lax-Wendroff scheme based on http://folk.ntnu.no/leifh/teaching/tkt4140/._main075.html and it seems to work fine.

The following code compares your implementation with the various fixes, and the new one. Both work and give identical results.

import numpy as np
import matplotlib.pyplot as plt

Nx=100;Nt=100;t_end=2.0
nChoice = 1

### LAX-WENDROFF SCHEME with periodic BCs
Nx = Nx+1                         # Number of grid points
xmax = 1.                         # Domain limit to the right
xmin = 0.                         # Domain limit to the left
h = (xmax-xmin)/(Nx-1)           # Mesh size
x = np.arange(xmin,xmax,h)       # Discretized mesh
U =  np.exp(-10*(4*x-1)**2)       # Initial solution
dt = 1.0/Nt                        # Time step
t = 0.                            # Initial time

aplus, aminus, Fplus, Fminus, Fn = (np.zeros(Nx-1) for i in range(5))
nsteps = 0
while (t <= t_end):
if np.any(np.isnan(U)):
raise Exception('NaNs detected --> exit')

v = dt/(2*h)
a = (dt**2)/(2*h**2)
Un = U
Um = np.roll(Un,1)
Up = np.roll(Un,-1)
# F_plus  = 0.5 * (Up + Un)
# F_minus = 0.5 * (Un + Um)
F_plus = 0.5*Up**2
F_minus = 0.5*Um**2
Fn = 0.5*Un**2
# U =  Un - v * (Um-Up) + a * (F_plus  * Un-Um - F_minus *Um-Up)
U =  Un - v*(F_plus-F_minus) + a*( 0.5*(Up+U)*(F_plus-Fn)- 0.5*(U+Um)*(Fn-F_minus))

else: # another implementation
aplus[:-1] = 0.5 * (U[:-1] + U[1:])
aplus[-1]  = 0.5 * (U[-1]  + U[0])
aminus[1:] = aplus[:-1]
aminus[0]   = aplus[-1]

Fn[:] = 0.5*U**2
Fplus[:-1] = Fn[1:] # or just use np.roll as you did
Fplus[-1]  = Fn[0]
Fminus[1:] = Fn[:-1]
Fminus[0]  = Fn[-1]

U[:] = U - dt/(2*h)*(Fplus-Fminus) + (dt**2)/(2*h**2)*(aplus*(Fplus-Fn) - aminus*(Fn-Fminus))

t = t+dt
nsteps+=1

# plot solution evolution every 10 time steps
if np.mod(nsteps,10)==1 or t>=t_end:
plt.figure()
plt.plot(x,U)
plt.title('t={:.2f}'.format(t),fontsize=16)
plt.xlabel('x',fontsize=18)
plt.ylabel('u',fontsize=18)
plt.grid()
plt.ylim(0,1)
plt.show()
$$$$

• thanks mate.You are awesome – user37062 Nov 28 '20 at 11:45