Alternatives to von neumann stability analysis for finite difference methods

Question

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?

Answer 1

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.

Answer 2

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.