Why is my simulation of a first-order wave equation 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()

For one thing, your problem does not look well-posed to me, similar to this question: https://scicomp.stackexchange.com/questions/5425. — Kirill, Jun 23 '18 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. — Reid.Atcheson, Jun 23 '18 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. — Bill Greene, Jun 23 '18 at 21:36

score 8 · Answer 1 · answered Jun 23 '18 at 21:34

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.

cfdlab · Answer 2 · 2018-06-24T04:53:15.620

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[0] = u[0] + 0.5*nu*(u[-2] - u[1])
    unew[1:-1] = u[1:-1] + 0.5*nu*(u[0:-2] - u[2:])
    unew[-1] = unew[0]
    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()
    ax = fig.add_subplot(111)
    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()

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

2 Answers2