3
$\begingroup$

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?

The resulting graph is here

$\endgroup$
  • 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$ – Christian Clason Nov 12 '15 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 '15 at 15:27
5
$\begingroup$

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.

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

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 
$\endgroup$
  • $\begingroup$ Good catch; it should also be range(1,N-1). Then it works (as badly as) expected. $\endgroup$ – Christian Clason Nov 12 '15 at 18:29
  • $\begingroup$ Compare my line with yours: I've added the parentheses $\endgroup$ – Steve Nov 12 '15 at 19:06
  • $\begingroup$ Should that read parentheses in my answer? In British english we call them all brackets. $\endgroup$ – Steve Nov 12 '15 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 '15 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 '15 at 19:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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