4
$\begingroup$

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()
$\endgroup$
  • 2
    $\begingroup$ For one thing, your problem does not look well-posed to me, similar to this question: scicomp.stackexchange.com/questions/5425. $\endgroup$ – Kirill Jun 23 '18 at 14:41
  • $\begingroup$ @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. $\endgroup$ – Reid.Atcheson Jun 23 '18 at 21:00
  • $\begingroup$ OK, I replaced my previous comment with an answer and included a link to a Strang video, which I like, on this topic. $\endgroup$ – Bill Greene Jun 23 '18 at 21:36
8
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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()
$\endgroup$

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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