Chapter 8 - Finite Differences: Partial Differential Equations

(8.1) Consider the $1D$ wave equation

$$\frac{\partial^2 u}{{\partial t}^2}=v^2\frac{\partial^2 u}{{\partial x}^2}\text{ }.$$

(a) Write down the straightforward finite-difference approximation.

Taking three consecutive points in time and space, we have $$\frac{(u^{n+1}_j-u^n_j)-(u^n_j - u^{n-1}_j)}{(\Delta t)(\Delta t)} = v^2\frac{(u^n_{j+1}-u^n_j)-(u^n_j-u^n_{j-1})}{(\Delta x)(\Delta x)} $$ Which gives the second order finite-difference approximation $$ u^{n+1}_j = 2u^n_j-u^{n-1}_j + v^2\frac{(\Delta t)^2}{(\Delta x)^2}[u^n_{j+1} -2u^n_j+u^n_{j-1}] $$

(b) What order of approximation is this in time and in space?

The equation is second order in space ($x$) and second order in time ($t$).
Useful link to find the order (spoiler: Taylor expansion of $u$) is the following
http://www.mathematik.uni-dortmund.de/~kuzmin/cfdintro/lecture4.pdf

(c) Use the con Neumann stability criterion to find the mode amplitudes.

Plugging in the approximation equation $u^n_j=A(k)^ne^{ijk\Delta x}$ we have $$ \begin{align} A^{n+1}e^{ijk\Delta x} &= 2A^ne^{ijk\Delta x} - A^{n-1}e^{ijk\Delta x} + v^2\frac{(\Delta t)^2}{(\Delta x)^2}[A^ne^{i(j+1)k\Delta x}-2A^ne^{ijk\Delta x}+A^ne^{i(j-2)k\Delta x}] \Rightarrow\\ A &= 2-\frac{1}{A} + v^2\frac{(\Delta t)^2}{(\Delta x)^2}[e^{ik\Delta x}-2+e^{-ik\Delta x}] \\ &= 2-\frac{1}{A} + v^2\frac{(\Delta t)^2}{(\Delta x)^2}[2\cos{k\Delta x}-2] \\ &= 2-\frac{1}{A} - 4v^2\frac{(\Delta t)^2}{(\Delta x)^2}\sin^2{\frac{k\Delta x}{2}} \end{align} $$ which results to the quadratic equation $$ A^2 + 2\beta A + 1=0 $$ where $\beta = 2v^2\frac{(\Delta t)^2}{(\Delta x)^2}\sin^2{\frac{k\Delta x}{2}} -1$.

Solving the equation, we get the two mode amplitudes: $$A = -\beta \pm \sqrt{\beta^2-1}$$

(d) Use this to find a condition on the velocity, time step, and space step for stability. (hint: consider the product of the two amplitude solutions)

Following the hint, the product of the solutions is $$A_1A_2 = \beta^2 - (\beta^2-1) = 1 $$. As the product is equal to $1$, then if one of the solutions is less than $1$ the other will have to be greater (and vice-versa) and thus the system will be unstable. Thus, both solutions should have a magnitude of $1$, i.e. they should be complex conjugates.

That is $A=-\beta \pm i\sqrt{-(\beta^2-1)}$$ $

Consequently, the only condition for stability is for the radicand to be positive, i.e. $b^2-1\leq0$, which yields $$\Big(2v^2\frac{(\Delta t)^2}{(\Delta x)^2}\sin^2{\frac{k\Delta x}{2}} -1\Big)^2 -1 \leq 0 $$

(e) Do different modes decay at different rates for the stable case?

The amplitude of the modes is independent of the change($k$) and thus they do not decay at all.

(g) If a damping term is added to the equation, assume that $u(x,t)=Ae^{i(kx-\omega t)}$ and find a relationship between $k$ and $\omega$, and simplify it for small $\gamma$. Comment on the relationship to the preceeding question.

The equation now is $$\frac{\partial^2 u}{\partial t^2}=v^2\frac{\partial^2 u}{\partial x^2}+\gamma\frac{\partial}{\partial t}\frac{\partial^2 u}{\partial x^2}$$

Inserting the solurion, we have $$ \begin{align} -\omega^2 &= -v^2k^2 + i\gamma k^2\omega \\ 0&=\omega^2 + i\gamma k^2 \omega -v^2k^2 \\ \omega &= \frac{-i\gamma k^2 \pm \sqrt{-\gamma^2k^4+4v^2k^2}}{2}\\ &\approx \frac{-i\gamma k^2 \pm \sqrt{4v^2k^2} }{2} \\ &\approx -\frac{i\gamma k^2}{2} \pm vk \end{align} $$

Plugging the modes back into the solution, we have $$y(x,t) = Ae^{ikx\pm vkt}e^{-\frac{i\gamma k^2}{2}t}$$ which corresponds to two waves travelling in opposite directions, damped exponentially in time, with a constant proportional to the square of the wavenumber $k$.

(f, h) Numerically solve the wave equation for the evolution from an initial condition with $u=0$ except for one nonzero node, and verify the stability criterion. For (h) include the damping term.

Taking inspiration (and some code, despite his Coffeescript implementation) from this rope project from Erik Ă–rjehag, I finally managed to write a simple and moldular three.js implementation of a physics engine. This will be super helpful for my final project (i.e. MD in the browser).

To have look at the code, click here: code

Essentially, the finite-difference approximation step happens in the this.apply(force) function call, where the wave constant $D$ is approximated by having nodes interconnected with springs having spring constant $k$.

Click the rope!

(8.2) Write a program to solve a 1D diffusion problem on a lattice of $500$ sites, with an initial condition of zero at all the sites, except the central site which starts at the value $1.0$. Take $D = \Delta x = 1$, and use fixed boundary conditions set equal to zero.

(a) Use the explicit finite difference scheme, and look at the behavior for $∆t =1, 0.5,$ and $0.1$. What step size is required by the Courant condition?

From (8.21) we can calculate the step size: $$\frac{2D\Delta t}{(\Delta x)^2} \leq 1 \Rightarrow \Delta t \le 0.5$$

The equation is solved below.

In [1]:
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib import rc
import numpy as np

def update_state_explicit(state, const):
    next_state = np.zeros_like(state)
    for j, u in enumerate(state):
        if j==0:
            next_state[j] = u + const*(state[j+1]-2*u + 0)
        elif j == len(state)-1:
            next_state[j] = u + const*(0-2*u + state[j-1])
        else:
            next_state[j] = u + const*(state[j+1]-2*u + state[j-1])
    return next_state

initial_state = np.zeros(501)
initial_state[251] = 1.
D = 1.
dx = 1.
dts = [1, .5, .1]
states = [[] for _ in dts]
for i, dt in enumerate(dts):
    states[i].append(initial_state)
    for s in range(300):
        states[i].append(update_state_explicit(states[i][-1], D*dt/(dx**2)))
In [2]:
fig, ax = plt.subplots()
ax.set_ylim(-1, 1)
ax.set_xlim(0, 501)
line, = ax.plot([], [], lw=2)
ax.grid()

def run0(i):
    # update the data
    line.set_data(np.arange(0,501), states[0][i])
    return line,


def run1(i):
    # update the data
    line.set_data(np.arange(0,501), states[1][i])
    return line,


def run2(i):
    # update the data
    line.set_data(np.arange(0,501), states[2][i])
    return line,

rc('animation', html='html5')
ani1 = animation.FuncAnimation(fig, run0, list(range(299)), blit=False, interval=50,repeat=True)
ani2 = animation.FuncAnimation(fig, run1, list(range(299)), blit=False, interval=50,repeat=True)
ani3 = animation.FuncAnimation(fig, run2, list(range(299)), blit=False, interval=50,repeat=True)

$\Delta t = 1.0$

In [3]:
ani1
Out[3]:

$\Delta t = 0.5$

In [4]:
ani2
Out[4]:

$\Delta t = 0.1$

In [5]:
ani3
Out[5]: