Alternatives to von neumann stability analysis for finite difference methods

I'm working on solving the coupled one-dimensional poroelasticity equations (biot's model), given as:

$$-(\lambda+ 2\mu) \frac{\partial^2 u}{\partial x^2} + \frac{\partial p}{\partial x} = 0$$ $$\frac{\partial}{\partial t} \left[ \gamma p + \frac{\partial u}{\partial x}\right] -\frac{\kappa}{\eta}\left[\frac{\partial^2 p}{\partial x^2}\right] =q(x,t)$$ on the domain $\Omega=(0,1)$ and with the boundary conditions:

$p=0, (\lambda + 2\mu)\frac{\partial u}{\partial x}=-u_0$ at $x=0$ and $u=0, \frac{\partial p}{\partial x} = 0$ at $x=1$.

I discretized these equations using a centered finite difference scheme:

$$(\lambda + 2\mu)\frac{u^{t+1}_{i+1}-2u^{t+1}_i+u^{t+1}_{i-1}}{\Delta x^2} + \frac{p^{t+1}_{i+1}-p^{t+1}_{i-1}}{2\Delta x} = 0$$ $$\gamma \frac{p^{t+1}_i-p^t_i}{\Delta t} + \frac{u^{t+1}_{i+1}-u^{t+1}_{i-1}}{2\Delta x \Delta t} -\left[ \frac{u^t_{i+1}-u^t_{i-1}}{2\Delta x \Delta t}\right] - \frac{\kappa}{\eta}\left[ \frac{p^{t+1}_{i+1} -2p^{t+1}_i + p^{t+1}_{i-1}}{\Delta x^2}\right]=q^{t+1}_i$$

I'm currently working out the details of the scheme's convergence by analyzing its consistency and stability. The consistency part seems fairly straightforward to me, but I'm already foreseeing some difficulties with the stability analysis. First of all, there are two variables and two equations. Secondly, there is also a mixed spatiotemporal derivative term in the second equation. I'm familiar with von neumann stability analysis and can see that it will be very tough to establish stability with this method. Are there any alternatives to von neumann analysis that I can use?

• If you don't feel comfortable doing the analysis with a system of equations, just differentiate the first equation with respect to $t$ and the second one with respect to $x$. Then use equality of mixed partial derivatives to eliminate $u$. – David Ketcheson May 24 '13 at 16:02
• @DavidKetcheson: Interesting. In essence, you're suggesting that I could reduce the system to a single variable and conduct the standard von neumann analysis on $p$ without any loss of generality to $u$? – Paul May 24 '13 at 16:20
• It's the same problem, whether you write it as a system or a scalar PDE. – David Ketcheson May 24 '13 at 17:16

If you substitute, at least for your analysis, $\frac{\partial u}{\partial x}$ by $u_x$, you can write your system as $$\begin{bmatrix} 0 & 0 \\ I & I \end{bmatrix} \frac{d}{dt} \begin{bmatrix} p_h(t) \\ u_{x,h}(t) \end{bmatrix}+\begin{bmatrix} -\partial_h & \partial_h \\ -\Delta_h & 0 \end{bmatrix}\begin{bmatrix} p_h(t) \\ u_{x,h}(t) \end{bmatrix}= \begin{bmatrix} q_h(t) \\ 0 \end{bmatrix} \quad (*)$$ where all constants are set to $1$ and where the subscript ${}_h$ refers to the space discretization both of the variables and the differential operators. Your scheme is then obtained by approximating $\frac{d}{dt}$ via implicit Euler.

Now the differential-algebraic (DAE) structure is evident. For the variables there are both differential (in time) and algebraic equations.

If you can show that $\begin{bmatrix} -\partial_h & \partial_h \\ I & I \end{bmatrix}$ is invertible, cf. this preprint [p. 3] and the edit below, than the DAE is of index 1 or strangeness-free and implicit Euler is known to be convergent, see Theorem 5.12 in this book. (Disclaimer: This book is not freely available and written by my PhD supervisor)

With this approach, you maybe get around the stability analysis.

For a direct proof of $L^2$ stability, I would try to use Equation $(*)$ to apply von Neumann stability analysis using the eigenfunctions of $\Delta_h$ and investigating the effect of $\partial_h$ onto the eigenfunctions.

However, if stability cannot be established for $(*)$, it doesn't mean your scheme is not convergent -- because of the substitution of $u\leftarrow u_x$. Generally speaking, one can expect stability for schemes approximating the actual variables, rather than for schemes approximating their derivatives.

APPENDIX: A DAE is said to be index 1, if it can be transformed into an ODE without differentiating the equations.

Say, the DAE is of the form $$\begin{bmatrix} E_1 \\ 0 \end{bmatrix}\dot y + \begin{bmatrix} A_1 \\ A_2 \end{bmatrix} y = f.$$ Then invertibility of $\begin{bmatrix} E_1 \\ A_2 \end{bmatrix}$ implies, that there is a variable transform $\tilde y \to y$ that eventually swaps the columns of the coefficients so that $\begin{bmatrix} E_1 \\ A_2 \end{bmatrix} \to \begin{bmatrix} \tilde E_{11} & \tilde E_{12} \\ \tilde A_{21} & \tilde A_{22} \end{bmatrix}$ with $\tilde A_{22}$ invertible (full rank property of $A_2$) and $\tilde A_{11} - \tilde E_{12} \tilde A_{22}^{-1} \tilde A_{21}$ invertible (the Schur complement).

For the system $(*)$ this means that the algebraic part defined with $A_2:=[-\partial_h ~ \partial_h]$ can be used to solve for a part $\tilde y_2$ of $(p_h,u_{x,h})$. Then, one can eliminate $\frac{d}{dt}\tilde y_2$ from the differential part (the second block line in $(*)$), to obtain an ODE for the remaining variables.

• This is a very interesting technique. I looked at the paper you referenced, and I'm curious how you concluded that $$\left[\begin{array}{cc}-\partial_h &\partial_h\\ I & I \end{array}\right]$$ must be invertible. Which theorem did you apply? – Paul May 25 '13 at 1:26
• @Paul I didn't find a theorem for reference, so I will insert the arguments into my answer... – Jan May 25 '13 at 12:24

I am not familiar with the equations given here, but I remember learning another method for checking the stability of a numerical scheme in my coursework. It is known as Modified Equation analysis.

Here is a good reference for that,

http://193.146.160.29/gtb/sod/usu/\$UBUG/repositorio/10291890_Warming.pdf


In the above reference, the connection between stability theory based on Modified Equation analysis and Von Neumann stability analysis is established.

After a bit of online search, I came across following references,

This paper discusses Finite difference modeling of Biot's poroelastic equations at seismic frequencies. It has a section on stability of numerical scheme as well.

This paper presents a solution strategy of decoupling the coupled system, and checking the stability of numerical scheme.

• I have not performed the modified equation analysis on above equations, but as the question asked for alternatives to Von Neumann analysis, I wrote the above answer. It is quite possible that it does not answer the question. But somebody might find the listed references useful in his/her work. – Subodh May 23 '13 at 6:37
• Thank you for the reference! I can see that the form needed in your Modified Equation Analysis paper doesn't quite fit the equations I'm using, but it's quite intriguing to learn new analysis techniques! – Paul May 25 '13 at 1:29