Hyperbolic Equation PDE (Python)

Question

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.xlabel('t')
plt.ylabel('u(x,t)')

plt.subplots_adjust(top=0.89, right=0.77)
plt.savefig('pde.png')

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

Answer 1

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.

Answer 2

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