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.

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]: