I'm trying to solve the following first order hyperbolic PDE problem using method of lines:

Hyperbolic Equation: $u_t = -u_x$

with initial condition: $u(0,x) = 0, 0 < x < 1$

boundary condition: $u(t,0) = 1, t \ge 1$

The solution should be a step function to the right with velocity $1$. I'm using centered finite difference to get an approximation of $u_x$.

Following the code in this tutorial, my code becomes

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from scipy.fftpack import diff as psdiff

N = 100 #no. of mesh points
L = 1.0
x = np.linspace(0, L, N) #mesh points xi, 0 < xi < 1
h = x[1] - x[0]

k = -1.0
def odefunc(u, t):
    ux = np.zeros(x.shape)
    u[0] = 1 # boundary condition
    for i in range(1,N-1):
        ux[i] = float(u[i+1] - u[i-1])/(2*h) 
    #   ux[i] = float(u[i] - u[i-1])/h 
    dudt = -ux
    return dudt

init = np.zeros(x.shape, np.float) #initial condition
tspan = np.linspace(0.0, 2.0, N)
sol = odeint(odefunc, init, tspan, mxstep=5000)

for i in range(0, len(tspan), 2):
    plt.plot(x, sol[N-1], label='t={0:1.2f}'.format(tspan[i]))

plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))

plt.subplots_adjust(top=0.89, right=0.77)

from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

SX, ST = np.meshgrid(x, tspan)
ax.plot_surface(SX, ST, sol, cmap='jet')
ax.view_init(elev=15, azim=-100) # adjust view so it is easy to see

However the resulting graph is not a step function as it should be. What might be the problem here?

The resulting graph is here

  • 1
    $\begingroup$ Welcome to Scicomp.SE! Please include a complete code, and ideally also a plot showing your results. The answer (which I'll leave to the experts on hyperbolic equations) is that your discretization is not stable (see Randy LeVeque's book on Finite Difference Methods, Chapter 10). $\endgroup$ Nov 12, 2015 at 15:12
  • $\begingroup$ @ChristianClason Thanks! I've just edited it with the complete code and the (painfully wrong) plot of the results. Thanks for you input. :) $\endgroup$
    – meraxes
    Nov 12, 2015 at 15:27

2 Answers 2


First-order hyperbolic equations model conservation laws; as the alternative name "transport equations" suggests, they transport information along so-called "characteristic curves" with a finite speed of propagation.

For the simple model equation $u_t + u_x = 0$, you can show (e.g., by separation of variables) that the solution should be of the form $$u(t,x) = u(0,x-t),$$ where the missing information (where $x-t<0$) is taken from the boundary conditions. As you correctly write, with the data in your question you'd expect a step profile moving to the right with constant speed $1$: $$u(t,x) = \begin{cases} 1 & x\leq t,\\0& x>t. \end{cases}$$

Here you can already see that something should be going wrong if you naively use a finite difference approximation: This solution is not differentiable (it's not even continuous), and yet you try to approximate its derivative using a finite difference quotient. So you need to be a bit more clever.

One way is to slightly modify the equation so that it does have a differentiable solution, for example by considering $u_t +u_x = \varepsilon u_{xx}$ for $\varepsilon>0$ very small, and apply a finite difference method to that. It turns out that taking $\varepsilon = \frac{\Delta t}{2}$ (half the discrete time step size) is a good idea; this is called the Lax-Wendroff method: $$U^{n+1}_{j} = U^n_j - \frac{\Delta t}{2\Delta x}(U^n_{j+1}-U^n_{j-1}) + \frac{\Delta t^2}{2\Delta x^2}(U^n_{j+1}-2U^n_j + U^n_{j-1}).$$ You can show that this is second-order accurate and stable as long as $\Delta t\leq \Delta x$. (This is discussed in Chapter 10 of Randy LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM 2007, together with other schemes).

So if you replace your odefunc routine by

t = np.linspace(0.0, 2.0, N)
k = t[1] - t[0]
def odefunc(u, t):
    ux = np.zeros(x.shape)
    u[0] = 1 # boundary condition
    for i in range(1,N-1):
        ux[i] = (u[i+1] - u[i-1])/(2*h) 
        ux[i] -= k*(u[i+1] - 2*u[i] + u[i-1])/(2*h**2)
    dudt = -ux
    return dudt

(note that I've renamed tspan to t for simplicity), you should get much better results.

  • $\begingroup$ Great! Looks much smoother now when I applied the Lax-Wendroff method. :) On a side note, I also seem to get nice results when I use upwind scheme. Thanks again for all your input! $\endgroup$
    – meraxes
    Nov 12, 2015 at 20:46
  • $\begingroup$ Looks like Lax-Wentroff is given for $\varepsilon = \frac{\Delta t}2$ $\endgroup$
    – Steve
    Nov 13, 2015 at 14:52
  • $\begingroup$ Right you are, thanks for pointing out! (Should have written $\mathcal{O}(\Delta t)$.) I've corrected it now. $\endgroup$ Nov 13, 2015 at 14:53

Don't have the rep to comment, but you're missing some brackets in odefunc

for i in range(1,N):
    ux[i] = float(u[i+1] - u[i-1])/(2*h)
#               parentheses added  ^   ^ 
#   ux[i] = float(u[i] - u[i-1])/h 
  • $\begingroup$ Good catch; it should also be range(1,N-1). Then it works (as badly as) expected. $\endgroup$ Nov 12, 2015 at 18:29
  • $\begingroup$ Compare my line with yours: I've added the parentheses $\endgroup$
    – Steve
    Nov 12, 2015 at 19:06
  • $\begingroup$ Should that read parentheses in my answer? In British english we call them all brackets. $\endgroup$
    – Steve
    Nov 12, 2015 at 19:09
  • $\begingroup$ @Steve Oh right! Thanks for pointing that out. :) Oh, sorry. From where I'm from, we call theme parentheses and refer to [] as brackets. Anyway, thanks again for pointing that bug! The graph looks ever weirder now though, haha. $\endgroup$
    – meraxes
    Nov 12, 2015 at 19:10
  • $\begingroup$ Ah, the ol' fixing a bug makes it look worse trick. You might want to try with some smoother initial conditions and see how it fares with them. Good luck $\endgroup$
    – Steve
    Nov 12, 2015 at 19:14

Your Answer

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

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