Hyperbolic Equation PDE (Python)

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 - x

k = -1.0
def odefunc(u, t):
ux = np.zeros(x.shape)
u = 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.xlabel('t')
plt.ylabel('u(x,t)')

plt.savefig('pde.png')

from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()

SX, ST = np.meshgrid(x, tspan)
ax.plot_surface(SX, ST, sol, cmap='jet')
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('u(x,t)')
ax.view_init(elev=15, azim=-100) # adjust view so it is easy to see
plt.savefig('pde-3d.png')

However the resulting graph is not a step function as it should be. What might be the problem here? • 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). Nov 12 '15 at 15:12
• @ChristianClason Thanks! I've just edited it with the complete code and the (painfully wrong) plot of the results. Thanks for you input. :) Nov 12 '15 at 15:27

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 - t
def odefunc(u, t):
ux = np.zeros(x.shape)
u = 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.

• 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! Nov 12 '15 at 20:46
• Looks like Lax-Wentroff is given for $\varepsilon = \frac{\Delta t}2$ Nov 13 '15 at 14:52
• Right you are, thanks for pointing out! (Should have written $\mathcal{O}(\Delta t)$.) I've corrected it now. Nov 13 '15 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)