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}] $$
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
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}$$
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 $$
The amplitude of the modes is independent of the change($k$) and thus they do not decay at all.
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$.
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$.
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.
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)))
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$
ani1
$\Delta t = 0.5$
ani2
$\Delta t = 0.1$
ani3