Stability of PDE Discretizations with Multistep Time Discretizations

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.

• Do you not just need the eigenvalues of $(2I + \Delta t^2 D)$? – Spencer Bryngelson Jun 6 '17 at 6:27
• @SpencerBryngelson I don't believe that is correct since that alone doesn't weigh in contributions from the $\boldsymbol{u}^{k-1}$ term. I have seen Von Neumann stability analysis used in the non vector case for situations like this where we are using past solution values, but I am just unsure how to handle this more general vector case. – spektr Jun 6 '17 at 6:30
• @C.Howard Have a look at section 4.6 of this book by Nick Trefethen. For vector finite-difference formulas, the Von Neumann condition - theorem 4.10 - (a necessary, but not sufficient condition for discrete stability) is formulated in terms of the spectral radius of the amplification matrix. There is a very interesting connection to pseudo-spectra (that I believe are one of the author's favourite research topics). See theorem 4.11 on page 181. – GoHokies Jun 6 '17 at 9:15
• @GoHokies wow that is an excellent reference. It looks like this reference is what I need to approach stability in this context. I will make an update once I have things sorted out. – spektr Jun 6 '17 at 14:06
• Maybe you can find something useful here, samarskii.ru/books/book2001_2.pdf, page 428-... But it assumes something about the matrix $D$. Actually, I don't understand the necessity of considering the general case - are you aiming to solve a particular kind of PDE's or developing general theory? The more general is the theory, the less strong are results often. – VorKir Jun 6 '17 at 16:50