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.


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}


\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.

  • 1
    $\begingroup$ Do you not just need the eigenvalues of $(2I + \Delta t^2 D)$? $\endgroup$
    – user20857
    Jun 6, 2017 at 6:27
  • $\begingroup$ @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. $\endgroup$
    – spektr
    Jun 6, 2017 at 6:30
  • 2
    $\begingroup$ @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. $\endgroup$
    – GoHokies
    Jun 6, 2017 at 9:15
  • $\begingroup$ @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. $\endgroup$
    – spektr
    Jun 6, 2017 at 14:06
  • $\begingroup$ 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. $\endgroup$
    – VorKir
    Jun 6, 2017 at 16:50

1 Answer 1


What you're looking for is a discussion of absolute stability of multistep methods for ODEs. See for example chapter 7 of LeVeque's book on finite difference methods .

Once you know the region of absolute stability for your time discretization, it's just a matter of checking that the eigenvalues of D (scaled by the time step) are inside that region. This is typical method-of-lines stability analysis.

The approach you've mentioned in your edit is equivalent, but in a sense less flexible since it doesn't separate the analysis of the time and space discretizations so clearly. Also, the approach I'm recommending can be done entirely with pencil and paper.

  • $\begingroup$ You're definitely right. After I did my edit, I ended up stumbling on stability for multistep schemes and seeing how you could obtain the stability region just like in method of lines. This will certainly be the approach I take in the future after some further study. I also agree that it is nice to have separation of time and space when looking at stability. $\endgroup$
    – spektr
    Jun 7, 2017 at 3:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.