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?