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)))
    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):
    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.xlabel('x', fontsize=18)
plt.ylabel('u(x, t=1)', fontsize=18)
plt.title(fr'The norm of the error compared to the exact solution', fontsize=20)

This create the following plot: enter image description here

  • $\begingroup$ Is the problem with the visualization or with the data? What happens if you output the solution? What happens if you output the numbers you want to visualize? $\endgroup$ Sep 28, 2022 at 22:48
  • 1
    $\begingroup$ Good question. When plotting the solution, it looks good. Data I want to visualise is strange. As I increase the mesh, I should see smaller error but instead, I get the same errors. $\endgroup$
    – nyaki
    Sep 28, 2022 at 23:03
  • $\begingroup$ Well, so then there is where you should look: How the error is computed. $\endgroup$ Sep 29, 2022 at 16:33

1 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.


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.

  • $\begingroup$ I corrected my code w.r.t your first point. I also made some correction to the steps but now some of them are constant and for some the error increase, even though they should be good under the CFL condition. $\endgroup$
    – nyaki
    Sep 29, 2022 at 14:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.