Why is my simulation of a first order wave eqaution not stable?

Question

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()

Answer 1

One direct problem that you have is that as a first order ODE there should be only one boundary condition in $x$ direction. You impose $2$, $x(t,0)=0$ which is sensible, and $x(t,L)=0$ which leads to large difference quotients as the incoming wave is cut off to zero. Replace the last condition with a backward difference quotient of order 2,

for t in range(0,iter-1):
    dy_t[1:-1] = -alpha*(y[t, 2:]-y[t,:-2])/(2*dx)
    dy_t[-1]   = -alpha*(3*y[t,-1]-4*y[t,-2]+y[t,-3])/(2*dx)

    y[t+1,1:] = y[t,1:]+dy_t[1:]*dt

and the start of the animation looks much more reasonable.

The inevitable fact is that the Euler method in time direction is of order $1$. The balancing condition that you cite only ensures that both time and space discretization give about the same contribution to the discretization error. The error still accumulates as $O(t\cdot Δt)$, for larger times it is in first order proportional to $\sim \frac1L(e^{Lt}-1)⋅Δt$ where $L$ is the Lipschitz constant of $\dot y=f(y)$, so in the main part $L=\frac1{|a|⋅Δx}$ which is $L=1$ with your data. This exponentially exploding error is very well visible in the animation, even in the corrected version.