The wave equation 1D with constant density is defined as:
\begin{equation}
\frac{\partial^2 U}{\partial t ^2} = V^2 \frac{\partial^2 U}{\partial x ^2}
\label{eqa}
\end{equation}
And the implicit difference schema:
\begin{equation}
U_j^{n+1} = \left( \frac{\Delta t V_j^n}{\Delta s} \right) ^2 \left( U_{j+1}^{n+1} – 2 U_j^{n+1} + U_{j-1}^{n+1} \right) + 2 U_j^n – U_j^{n-1}
\label{eq:b}
\end{equation}
Where $\Delta x = \Delta s $.
Von Newman stability analysis.
Since $ u(x,t) = e^{i(wt+kx)} $ is a solution for the discrete schema above can be written as:
\begin{eqnarray}
u(x_j,t_n) &=& e^{i(wt_n+kx_j)} \nonumber \\
&=& e^{iwn\Delta t} e^{ikj\Delta s} \nonumber \\
&=& \epsilon^n e^{ikj\Delta s}
\label{eq:c}
\end{eqnarray}
Where $\epsilon = e^{iw \Delta t} $ and $ i = \sqrt{-1} $ .
To maintain stability we should make sure that $ \epsilon \leq 1 $ not growing exponentially with increasing time steps. So applying the previous equation in the implicit differences schema we can analyse the growth.
\begin{multline}
U_j^{n+1} = \left( \frac{\Delta t V_j^n}{\Delta s} \right) ^2 \left( U_{j+1}^{n+1} – 2 U_j^{n+1} + U_{j-1}^{n+1} \right) + 2 U_j^n – U_j^{n-1} \\
\epsilon^{n+1} e^{ikj\Delta s} = r^2 \left( \epsilon^{n+1} e^{ik(j+1)\Delta s} – 2 \epsilon^{n+1} e^{ikj\Delta s} + \epsilon^{n+1} e^{ik(j-1)\Delta s} \right) \\
+ 2 \epsilon^n e^{ikj\Delta s} – \epsilon^{n-1} e^{ikj\Delta s} \\
1 = r^2 (e^{ik\Delta s} -2 + e^{-ik\Delta s}) + 2 \epsilon^{-1} – \epsilon^{-2} \\
\epsilon^{-2} – 2 \epsilon^{-1} – r^2 (e^{ik\Delta s} -2 + e^{-ik\Delta s}) + 1 = 0 \\
\epsilon^{-2} – 2 \epsilon^{-1} – 4 r^2 \sin^2 \left( k \Delta s / 2 \right) + 1 = 0
\label{eq:d}
\end{multline}
Where $ r = \left( \frac{\Delta t V_j^n}{\Delta s} \right) $
At the previous equation we can substitute $ \phi = \epsilon^{-1} $ turning it to a second degree bellow:
\begin{eqnarray}
\phi^2 – 2 \phi -4 r^2 \sin^2 \left(k\Delta s/2 \right) +1 = 0 \nonumber \\
\phi^2 – 2 \phi + c = 0 \nonumber
\end{eqnarray}
With:
$$ c = 1 -4 r^2 \sin^2 \left(k\Delta s/2 \right) $$
Has roots $ \phi^{'} $ and $ \phi^{''} $ as bellow. Replacing $c$ also.
\begin{eqnarray}
&=& \frac{2 \pm \sqrt{4 – 4c}}{2} = 1 \pm \sqrt{1-c} \nonumber \\
&=& 1 \pm \sqrt{4 r^2 \sin^2 \left(k\Delta s/2 \right) } \nonumber \\
&=& 1 \pm 2 r \sin \left(k\Delta s/2 \right) \nonumber \\
\phi^{'} &=& 1 + 2 r \sin \left(k\Delta s/2 \right) \\
\phi^{''} &=& 1 – 2 r \sin \left(k\Delta s/2 \right)
\end{eqnarray}
Going back to $\epsilon$ we get :
\begin{eqnarray}
\epsilon^{'} &=& \frac{1}{1 + 2 r \sin \left(k\Delta s/2 \right)} \\
\epsilon^{''} &=& \frac{1}{1 – 2 r \sin \left(k\Delta s/2 \right) }
\end{eqnarray}
But with $r \sin \left(k\Delta s/2 \right) \to 0$ we get that $\epsilon^{''} \to \infty $
Should it not be unconditionally stable? What is wrong in here?
Update: After awesome answer by Jitse, fixing the signal we get:
\begin{eqnarray}
\phi^{'} &=& 1 + i 2 r \sin \left(k\Delta s/2 \right) \\
\phi^{''} &=& 1 – i 2 r \sin \left(k\Delta s/2 \right)
\end{eqnarray}
Analyzing the modulus, since $ \phi^{'} $ and $ \phi^{''} $ are conjugate pairs of the same complex number the modulus is the same for both.
We get:
$$ \| \phi^{'} \| = \| \phi^{''} \| = \sqrt{1+ 4 r^2 \sin ^2 \left(k\Delta s/2 \right)}$$
So since $ \| \phi \| \geq 1 \ \ \forall \ r, \ k,\ \Delta s $ than
$$ \| \epsilon \| = \frac{1}{ \sqrt{1 + 4 r ^2 \sin ^2 \left(k\Delta s/2 \right)}} $$
$$ \| \epsilon \| \leq 1 \ \ \forall \ r, \ k,\ \Delta s$$
Awesome, unconditionally stable!!
Best Answer
Your method for Neumann stability analysis is correct, and indeed I would expect the scheme to be unconditionally stable.
I do not understand the step where you introduce the sine. I think that $$ e^{ix} - 2 + e^{-ix} = -4\sin^2(\tfrac12x). $$ The minus sign seems to have been lost in your computation.