How to correctly plot Order of accuracy for different finite difference schemes

Question

I have implemented Upwind, Lax, Lax-Wendroff, Leapfrog and macCormak method for the linear advection equation with Dirichlet boundary conditions. I am trying to create the order of accuracy plots for these methods but I receive a constant lines for them. Not sure what is the issue.

def u0(ig):
    x = np.arange(0, ig, 1)
    u = .5*(1+np.tanh(250*(x-20)))

    return u

def exact(x, t):
    u = .5*(1+np.tanh(250*((x-.5*t)-20)))
    u[-1]=0
    u[0]=0
    return u

def error(analytical, numerical):
    return np.max(np.abs(numerical[-1, 1:-1]-analytical[1:-1]))

def upwind(Nt, Nx, dx, dt):
    solution = np.zeros((Nt, Nx))
    solution[0, :] = u0(Nx)
    solution[:,0] = 0
    solution[:, -1] = 0
    for n in range(0,Nt-1):
        solution[n+1,1:Nx-1] = solution[n,1:Nx-1] - dt/dx*(  max(.5,0)*(solution[n,1:Nx-1] - solution[n,0:Nx-2] + min(.5,0)*(solution[n,2:Nx] 
        - solution[n,1:Nx-1] )))
    
    return solution

def lax(Nt, Nx, dx, dt):
    solution = np.zeros((Nt, Nx))
    solution[0, :] = u0(Nx)
    solution[:,0] = 0
    solution[:,-1] = 0
    for n in range(0, Nt-1):
        solution[n+1,1:Nx-1] = 0.5*(solution[n,2:Nx] + solution[n,0:Nx-2] ) - 0.5*.5*dt/dx*( solution[n,2:Nx] - solution[n,0:Nx-2] )

    return solution

def wendroff(Nt, Nx, dx, dt):
    c = .5
    solution = np.zeros((Nt, Nx))
    solution[0, :] = u0(Nx)
    solution[:,0] = 0
    solution[:,-1] = 0
    for n in range(0, Nt-1):
        solution[n+1,1:Nx-1] = solution[n, 1:Nx-1] - 0.5*c*dt/dx*( solution[n,2:Nx] - solution[n,0:Nx-2] ) 
        + 0.5 * (c*dt/dx)**2*(solution[n,2:Nx] - 2*solution[n,1:Nx-1] + solution[n,0:Nx-2] )
    
    return solution

def leapfrog(Nt, Nx, dx, dt):
    c=.5
    solution = np.zeros((Nt, Nx))
    solution[0,:] = u0(Nx)
    solution[:, 0] = 0
    solution[:, -1] = 0
    solution[1,1:Nx-1] = solution[0,1:Nx-1] - dt/dx*(max(c,0)*(solution[0,1:Nx-1] - solution[0,0:Nx-2] + min(c,0)*(solution[0,2:Nx] - solution[0,1:Nx-1])))
    for n in range(1, Nt-1):
        solution[n+1,1:Nx-1] = solution[n-1, 1:Nx-1] - c*dt/dx*( solution[n,2:Nx] - solution[n,0:Nx-2])
    
    return solution

def macCormack(Nt, Nx, dx, dt):
    a = (.5 *dt)/dx
    solution = np.zeros((Nt, Nx))
    solution[0,:] = u0(Nx)
    solution[:, 0] = 0
    solution[:, -1] = 0
    up = solution.copy()
    for n in range(0, Nt-1):
        up[n+1, :-1] = solution[n, :-1] - a*(solution[n ,1:]-solution[n, :-1])
        solution[n+1, 1:] = .5*(solution[n, 1:]+up[n, 1:] -  a*(up[n, 1:]-up[n, :-1]))
    
    return solution

Nxs = [41, 100, 1000]
dxs = [1, .5, 1/10]
dts = [1, .5, 1/10]

errors = {}
errors['upwind'] = []
errors['lax'] = []
errors['wendroff'] = []
errors['leapfrog'] = []
errors['maccormak'] = []

for i in range(len(Nxs)):
    x = np.arange(0, Nxs[i], 1)
    ue = exact(x, Nts[i])
    sol_upwind_study = upwind(Nts[i], Nxs[i], dxs[i], dts[i])
    errors['upwind'].append(error(ue, sol_upwind_study))
    sol_lax_study = lax(Nts[i], Nxs[i], dxs[i], dts[i])
    errors['lax'].append(error(ue, sol_lax_study))
    sol_wendroff_study = wendroff(Nts[i], Nxs[i], dxs[i], dts[i])
    errors['wendroff'].append(error(ue, sol_wendroff_study))
    sol_leapfrog_study = leapfrog(Nts[i], Nxs[i], dxs[i], dts[i])
    errors['leapfrog'].append(error(ue, sol_leapfrog_study))
    sol_maccormak_study = macCormack(Nts[i], Nxs[i], dxs[i], dts[i])
    errors['maccormak'].append(error(ue, sol_maccormak_study))

plt.figure(figsize=(16, 9))
plt.plot(Nxs, errors['upwind'], 'x-', label='Upwind scheme', ms=12)
plt.plot(Nxs, errors['wendroff'], 'o-', label='Lax-Wendroff scheme')
plt.plot(Nxs, errors['lax'], 'x-', label='Lax-Friedrich scheme', ms=12)
plt.plot(Nxs, errors['leapfrog'], 'x-', label='Leapfrog scheme', ms=12)
plt.plot(Nxs, errors['maccormak'], 'x-', label='Maccormak scheme', ms=12)
#plt.xscale('log')
plt.yscale('log')
plt.legend(fontsize=14)
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)
plt.xlabel('x', fontsize=18)
plt.ylabel('u(x, t=1)', fontsize=18)
plt.grid()
plt.title(fr'The norm of the error compared to the exact solution', fontsize=20)
plt.show()

This create the following plot:

Answer 1

Here are a few issues I noticed:

1. Your boundary conditions you've chosen in your numerical solutions are inconsistent with the analytical solution. You set the right boundary to 0, however it should really be exact(x_right, t), which until the transition point of the hyperbolic tangent moves to the right is about 1. Ideally this shouldn't matter since you are upwinding to the right, but you include the BC nodes in your error calculation.

2. Convergence rates are often specified as $O(\Delta x^N)$ for a reason: they depend on the local resolution. You keep $\Delta x$ fixed but keep increasing the domain size. Note that as you decrease $\Delta x$ you might run into stability limits for $\Delta t$. Try choosing a $\Delta t$ which will be stable for the smallest $\Delta x$ you plan to test.

edit:

looking at your updated code, you're still computing the exact solution using the domain [0, Nx] with dx = 1. Your initial conditions for the solvers also still initialize $u_0$ with this domain, and then solve it with a completely different $\Delta x$. You're also evaluating your exact solution at time t=Nts[i] (whatever Nts is, it's not in your code segment). However, you're evaluating all of your numerical solutions at Nts[i] * dts[i]. You are also varying $\Delta t$ simultaneously with mesh resolution, so you have no way to tell if the convergence rate is due to temporal discretization error vs. spatial discretization error.

Plot your solution over the exact solution as well as the difference betwen the two and you should be able to see any obvious discrepancies relatively easily.

The last suggestion I have is I don't think the max norm is a particularly good method for measuring convergence. If the numerical solution predicts the speed is slightly faster or slower than the exact solution that may show up as a rather large difference in the max norm, which other norms like an L-1 or L-2 norm wouldn't see. I'm unsure if this change will fix your convergence rates, though I have run into this type of problem in the past and almost exclusively use the L-2 norm when I can.