Let's pretend we have a spatially discretized PDE of the following form:
\begin{align} \frac{\partial^2 \boldsymbol{u}^k}{\partial t^2} = D\boldsymbol{u}^k \end{align}
where $D$ can be any form for now and $k$ refers to this being the discretization at time $t_k$. Then let's suppose, for example, we use a central difference approximation for the time derivative. We would then arrive at the following:
\begin{align} \frac{\partial^2 \boldsymbol{u}^k}{\partial t^2} &= D\boldsymbol{u}^k \\ \frac{\boldsymbol{u}^{k+1} - 2\boldsymbol{u}^k + \boldsymbol{u}^{k-1}}{\Delta t^2} &= D\boldsymbol{u}^k \\ \boldsymbol{u}^{k+1} &= \left( 2I + \Delta t^2 D \right)\boldsymbol{u}^k - \boldsymbol{u}^{k-1}\\ \end{align}
Given the formulation above, how could one approach stability? First thoughts go to using Von Neumann stability analysis, but I have never seen it used in a vector-wise fashion such that one could take into account an arbitrary $D$ matrix. Any thoughts or references would be very useful.
Edit
The references provided in the comments were useful, but I found that only a few simple connections were needed to approach stability for this scenario. The key link was to take the expression above and cast it into a different difference equation with a more suitable form.
The more suitable case is by defining $\boldsymbol{w}^k = \left[(\boldsymbol{u}^{k})^{T}, (\boldsymbol{u}^{k-1})^{T}\right]^T$. We can than recast our difference equation into the form:
\begin{align} \boldsymbol{w}^{k+1} = G \boldsymbol{w}^k \end{align}
where
\begin{align} G = \begin{bmatrix} (2I + \Delta t^2 D) & -I \\ I & 0 \end{bmatrix} \end{align}
Then we know this system is stable if, given the set of eigenvalues $\lbrace \lambda_i \rbrace$ for $G$, the following holds:
\begin{align} \max_{i} \left| \lambda_i\right| \lt 1 \end{align}
With these results, we can check whether some differential operator $D$ is going to be stable with a given time discretization, particularly if it uses older $\boldsymbol{u}$ states.
