# Why is my simulation of a first-order wave equation not stable?

According to the equation

$$\frac{\partial y}{\partial t} = -a\frac{\partial y}{\partial x}$$

I simulated this in python. I used center differentiation, and I determined step size based on Von-Neumann stability criterion for a linear equation $\Delta t = \frac{\Delta x^2}{2a}$ .

This is what I get:https://gyazo.com/86024db42f71a6b5cb34eca7eb0d115f

import os
import numpy as np

import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.animation as animation

alpha = 1
seconds = 1000
dx = 1
dt = .1/alpha*dx**2
iter = int(np.round(seconds/dt))
print(iter, dt)

x = np.linspace(0,1000,1000)
# z = np.linspace(0,990, 10, dtype= int)
y = np.zeros([iter,1000])
y[0,:] = 50*np.sin(x[:]*np.pi/1000)
dy_t = np.zeros(1000)

for t in range(0,iter-1):

for x in range(1, 998):

dy_t[x] = -alpha*(y[t, x+1]-y[t,x-1])/(2*dx)

for x in range(1,998):
y[t+1,x] = y[t,x]+dy_t[x]*dt

# animation
fig, ax = plt.subplots()
line, = ax.plot(y[0,:])
def init():
line.set_ydata([np.nan] * 1000)
return line,

def animate(i):
line.set_ydata(y[i,:])
return line,

ani = animation.FuncAnimation(
fig, animate, range(0, iter), init_func=init, interval=5, blit=True, save_count=50)
plt.show()

• For one thing, your problem does not look well-posed to me, similar to this question: scicomp.stackexchange.com/questions/5425. Jun 23, 2018 at 14:41
• @BillGreene I was about to make that comment as an answer, but since you have already pointed it out you should instead? An alternative time-stepper would be RK4, I guess. Jun 23, 2018 at 21:00
• OK, I replaced my previous comment with an answer and included a link to a Strang video, which I like, on this topic. Jun 23, 2018 at 21:36

Numerical solution of the advection equation with centered differences in space and forward Euler in time is unconditionally unstable. So the behavior you are seeing is expected.

Here is a nice lecture by Gil Strang ( MIT 18.086) where he discusses the instability of this method. He also shows how the simple centered difference method can be modified to yield the Lax-Friedrichs and Lax-Wendroff methods which are conditionally stable.

I am going to answer this assuming you really want to use central difference scheme in space and forward difference in time.

Consider the problem $$u_t + a u_x = 0, \quad x \in (x_L,x_R), \quad t \in [0,T]$$ and a numerical scheme $$u_j^{n+1} = H_j(u^n; a, \Delta t, h)$$ The scheme is said to be weakly stable if there exists $C = C(T)$ such that $$\| u^n \| \le C \| u^0 \|, \qquad n \Delta t \le T$$ An equivalent definition is $$\| u^{n+1} \| \le ( 1 + C \Delta t) \| u^n \|$$ If $C=1$ then we say it is stable (strongly stable). The FTCS scheme $$\frac{u_j^{n+1} - u_j^n}{\Delta t} + \frac{u_{j+1}^n - u_{j-1}^n}{2h} = 0$$ is never strongly stable. It is weakly stable under the condition $$\Delta t \le \left( \frac{h}{a} \right)^2$$ This is not a useful scheme: (1) The time step is too small. (2) It is not useful for long time integration.

You seem to have some bugs in your code which you have to fix yourself. I give my own implementation below where time step is chosen as $$\Delta t = CFL \left( \frac{h}{a} \right)^2$$ Run this with Python2; for Python3 you have to make some changes.

import numpy as np
import matplotlib.pyplot as plt

def smooth(x):
return np.sin(2*np.pi*x)

def update_ftcs(nu, u):
unew = np.empty_like(u)
unew = u + 0.5*nu*(u[-2] - u)
unew[1:-1] = u[1:-1] + 0.5*nu*(u[0:-2] - u[2:])
unew[-1] = unew
return unew

def solve():
scheme     = 'FTCS'
N          = 100
cfl        = 0.95
Tf         = 1.0
uinit      = smooth
xmin, xmax = 0.0, 1.0
a          = 1.0

h = (xmax - xmin)/N
dt= cfl * (h / np.abs(a))**2
nu= a * dt / h

x = np.linspace(xmin, xmax, N+1)
u = uinit(x)

fig = plt.figure()
line1, = ax.plot(x, u, 'r')
line2, = ax.plot(x, u, 'b')
ax.set_xlabel('x'); ax.set_ylabel('u')
plt.legend(('Numerical','Exact'))
plt.title('N='+str(N)+', CFL='+str(cfl)+', Scheme='+scheme)
plt.draw(); plt.pause(0.1)
wait = raw_input("Press enter to continue ")

t, it = 0.0, 0
while t < Tf:
u = update_ftcs(nu, u)
t += dt; it += 1
if it%100 == 0:
line1.set_ydata(u)
line2.set_ydata(uinit(x-a*t))
plt.draw(); plt.pause(0.1)
plt.show()

solve()