# Instability of discrete wave equation

I've been studying the wave equation, and simulating it in this discrete form:

$$u(\vec x, t+1) = u(\vec x, t) + v(\vec x, t)$$ $$v(\vec x, t+1) = v(\vec x, t) + c^2\nabla^2 u(\vec x, t+1)$$

where $u$ is the waving quantity, $v$ is its time derivative, $\vec x \in \mathbb Z^n$ is the position vector, $t$ is time, and $c$ is the wave speed. (Note that this is different from the Euler step method, which would have $t$ in place of $t+1$ on the right.)

$\nabla^2$ is the discrete Laplacian:

$$\nabla^2 u(\vec x) = \sum_{i=1}^n \big(u(\vec x+\vec 1_i) + u(\vec x-\vec 1_i) - 2u(\vec x)\big) = \sum_{\vec y\text{ is next to }\vec x} \big(u(\vec y) - u(\vec x)\big)$$

(The unit basis vector $\vec 1_i$ could instead be called $\hat e_i$ or something similar.)

It seems that the system is stable if $c^2 \leq \frac1n$, and unstable otherwise; any non-zero initial conditions lead to unbounded value of $u$ . Can this be proven?

I think $\nabla^2$ is sometimes defined with a factor of $1/n$ ; then the system would be stable if $c^2 \leq 1$, which seems more fundamental. But it doesn't help with the proof.

I could also add a viscosity term, $b\,\nabla^2 v(\vec x, t)$ . That doesn't completely stabilize it; I expect $b$ to appear with $c^2$ in an inequality.

This question looks related, but it's hard to tell what he's asking.

This question uses a numerical method equivalent to mine (in 1 dimension).

EDIT1

Though this question uses domain $\mathbb Z^n$, the wave equation is often used on a finite rectangular region of space. Such a situation can be modelled in the whole space (and give correct results) by using periodic initial conditions. For example, if a waving string has a loose end at $x=a>0$ and a fixed end at $x=0$, then the value of $u(x)$ can be extended outside of $[0,a]$ by the reflection $u(a+x)=u(a-x)$ and the inversion $u(-x)=-u(x)$ . See also fundamental domain.

With that in mind, I've found a 1D counter-example to the question's premise: $u(\frac12+x)=u(\frac12-x)$ , $u(\frac52+x)=u(\frac52-x)$ (and the same symmetry for $v(x)$). This is periodic (in space, not time), with period $4$ . The initial conditions $u(1)$, $u(2)$, $v(1)$, and $v(2)$ determine everything. The system is stable with $c^2<2$ instead of $c^2\leq1$ . (Actually, it's only stable if $v(1)+v(2)=0$; otherwise, $v(x,t)$ may be bounded, but $u(x,t)$ is unbounded.)

I assume that this type of counter-example would be removed by the following condition:

The initial conditions $u(\vec x,0)$ and $v(\vec x,0)$ must have compact support [must be zero everywhere except a finite number of points].

EDIT2

I rewrote the system as a second-order recursion by eliminating $v$, as in Mark's question that I linked to.

$$u(\vec x, t+1) = 2u(\vec x, t) - u(\vec x, t-1) + c^2\nabla^2 u(\vec x, t)$$ $$= 2(1-nc^2)u(\vec x, t) - u(\vec x, t-1) + c^2\sum_{i=1}^n \big(u(\vec x+\vec 1_i, t)+u(\vec x-\vec 1_i, t)\big)$$

The inequality $c^2 \leq 1/n$ is equivalent to $2(1-nc^2) \geq 0$ ... and the sum on the right is $2n$ times the average of the surrounding values of $u$. These facts look useful for an answer!

If you want to include viscosity (relative drag), here's what I used:

$$u(\vec x, t+1) = u(\vec x, t) + v(\vec x, t)$$ $$w(\vec x, t) = v(\vec x, t) + c^2\nabla^2 u(\vec x, t+1)$$ $$v(\vec x, t+1) = w(\vec x, t) + b\nabla^2 w(\vec x, t)$$

($w$ is just an intermediate velocity, before accounting for drag.) The corresponding differential equation is

$$\frac{\partial^2 u}{\partial t^2} = c^2\nabla^2 u + b\nabla^2\frac{\partial u}{\partial t}$$

EDIT3

For $n=2$, I tried using a hexagonal instead of square grid; each point has 3 neighbors instead of 4. This was stable with $c^2 \leq \frac23$ . (Again, $c^2$ is the coefficient of the sum of $u(\vec y)-u(\vec x)$ where $\vec y$ is next to $\vec x$ .) This can be rewritten as $3c^2 \leq 2$ .

This is consistent with the pattern for the square & cubic grids, that $2nc^2$ (the coefficient of the average of $u(\vec y)-u(\vec x)$ ) must be $\leq 2$ .

I also tried a triangular grid, where each point has 6 neighbors. This seemed to be stable with $c^2 \leq \frac49$ , or $6c^2 \leq \frac83 \neq 2$ . That doesn't fit the pattern...