How to derive the stability of time stepping schemes?

Question

This is more of a mathematical question but since we deal with this all the time in computational science, maybe it is relevant in this forum too.

I am an engineer and I am learning how to model the wave equation in Fenics. To obtain the fully discretised variational form in space and time, I have to implement a numerical method for evolution in time. I have been looking at Sympletic Euler, Crank Nicolson and Newmark Beta method.

In literature, it is reported that these methods are conditionally or unconditionally stable for certain values. Focussing on Newmark Beta method, it is reported that this method is unconditionally stable, if we choose the parameters 2β ≥ γ ≥ 1/2. Here beta and gamma are parameters in the Newmark scheme.

My question is: How do mathematicians derive this information about the stability of a method. Whether it is conditionally stable or unconditionally stable, and for which parameter values?

Can anyone kindly point me in the right direction/ document/ papers? I am not researching this topic so I don't want to dive too deep into the associated mathematical research. I am just curious if there are plain derivations for this.

Thanks a lot!

Answer 1

Another way you can find stability is to assume first you are solving the following model problem for complex $\lambda$:

$$\dot{x} = f(t,x) = \lambda x$$

with exact solution: $x(t) = x_0 e^{\lambda t}$

Then you define your scheme, let's start with Explicit Euler:

$$ x_{k+1} = x_{k} + h \; f(t,x_{k}) $$

where $h$ is the step size. Based on this scheme, you substitute the value for your simple equation:

$$ x_{k+1} = x_{k} + h \lambda x_{k} = (1 + h \lambda ) x_{k} $$

We then can say that an equation is stable if the following is true:

$$ \left| \frac{x_{k+1}}{x_{k}} \right| = \left| \sigma \right| \le 1 \;\;\;\;\;\forall k $$

Note that $\left| \cdot \right|$ is the complex modulus in this case. With this, we can find the ratio $\sigma$ is:

$$ \sigma = (1 + h\lambda_{Re}) + i h \lambda_{Im}$$

where $\lambda_{Re} = Re[\lambda]$ and $\lambda_{Im} = Im[\lambda]$. We can then work out and obtain the following:

$$ \left| \sigma \right|^2 \le 1 $$ $$ (1+h \lambda_{Re})^2 + h^2 \lambda_{Im}^2 \le 1$$ $$ 1 + 2 h \lambda_{Re} + h^2 (\lambda_{Re}^2 +\lambda_{Im}^2) \le 1$$ $$ 2 h Re[\lambda] + h^2 \left| \lambda \right|^2 \le 0 $$ $$ 2 Re[\lambda] + h \left| \lambda \right|^2 \le 0 $$ $$ h \le \frac{-2 Re[\lambda]}{\left| \lambda \right|^2} $$

Now first, given $h$ is positive, the term on the right hand side must be positive. This implies that $Re[\lambda] \lt 0 $ must be true for the problem to potentially be stable.

For purely real $\lambda$, you can simplify the expression to:

$$ h \le \frac{2 }{\left| \lambda \right|} $$

Now given you have discretized your problem to the following system equations:

$$ \dot{\textbf{u}} = A \textbf{u}$$

where $A$ is some matrix, specifically a differential operator when solving Partial Differential Equations, you can then theoretically transform this system using the eigensystem of $A$, $A=Q \Lambda Q^{-1}$, by doing the following:

$$ \dot{\textbf{u}} = A \textbf{u}$$ $$ \dot{\textbf{u}} = Q \Lambda Q^{-1} \textbf{u}$$ $$ Q^{-1} \dot{\textbf{u}} = \Lambda Q^{-1} \textbf{u}$$ Define $\textbf{z} = Q^{-1} \textbf{u}$ $$ \dot{\textbf{z}} = \Lambda \textbf{z}$$ or $$ \dot{z}_{j} = \lambda_{j} z_{j} \;\;\; \forall j$$

where $\lambda_{j}$ is the $j^{th}$ eigenvalue of $A$. Now that the system is diagonalized, you can find that the most limiting step size will be:

$$ h \le 2 \min_{j} \left(\frac{-Re[\lambda_{j}]}{\left| \lambda_{j} \right|^2}\right) $$

which will be based on the eigenvalue with the largest modulus in this case

Conclusion

Provided you have some scheme, you can derive/look up it's stability criteria relative to the model problem. After you have this stability criteria, you can apply it to problems by, if the scheme is conditionally stable, finding eigenvalues of your matrix operator $A$ and find the most limiting step size using the eigenvalues.

Answer 2

One way is via von Neumann stability analysis. This technique hinges on the fact that the Fourier transform can sometimes be regarded as an $L^2$ isometry. This and other techniques such as direct norm estimates are discussed in, for example, the first few chapters of:

Strikwerda, J. C. (2004). Finite difference schemes and partial differential equations (Second edition). Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM).

Answer 3

To show stability of a numerical scheme, as stated above, Von Neumann analysis is typically the way to go. It looks at a single Fourier mode of the scheme to find the amplification factor, and from that, one may deduce the sufficient conditions of the parameters such that the scheme is stable.

To wit, we show the Backward Euler scheme for the Cauchy problem for the diffusion equation is strongly $\ell_2$ stable, for any $\mu>0$, where $\mu = \frac{\Delta t}{\Delta x^2}$. Letting $w^{n+1}_l$ denote the approximation at time $\Delta t\cdot n$ and $\Delta x\cdot l$ in the domain, and define $h:=\Delta x$.

Let's assume that $w^n_l = e^{ilh\xi}$, that is $w^n_l$ is a single Fourier mode. Hence $w^{n+1}_l = g(\xi)w^n_l = g(\xi)e^{ilh\xi}$, and we seek $g(\xi)$, the amplification factor. So \begin{eqnarray*} w^n_l &=& -\mu(w^{n+1}_{l-1} + w^{n+1}_{l+1}) + (1+2\mu)w^{n+1}_{l}\\ &=& -\mu g(\xi)(e^{i(l-1)h\xi} + e^{i(l+1)h\xi}) + (1+2\mu)e^{ilh\xi}g(\xi)\\ &=& g(\xi)\left(1+2\mu(1-cos(h\xi))\right)e^{ilh\xi}\\ \end{eqnarray*}

Now we claim that the scheme is unconditionally stable. Hence $$ |g(\xi)| = |\frac{1}{1+2\mu(1-cos(h\xi))}|\leq 1 $$ or $$ 1 \leq |{1+2\mu(1-cos(h\xi))}| $$ We find $4\mu \geq 1+2\mu(1- cos(h\xi)) \geq 1$, conditional of $\mu$.

Similar analysis is carried out for other schemes to find the sufficient conditions of the parameters so that the scheme is stable. Specifically, the analysis for Crank Nicolson is more complicated since it has half a Forward Euler step and half a Backward Euler step, and is correspondingly stable for $\mu\leq 1$.

Answer 4

Basis of Stability

The stability of a method is derived on the basis of truncation error, an error that accumulates over each time step. In other words, the stability of a solution in space-time domain can be looked from a point of view of its dependence on a future time. Fully explicit schemes are known to be stable for a certain time step size because of the accumulation of truncation error over each time step. If the time step size is beyond the prescribed limit then a high truncation error would normally result in over estimating the final solution. However in implicit schemes the solution estimation depend on the current time step only, so they are independent of the time step size. This independence of time-step in an implicit scheme is gained by the introduction of numerical diffusion effects into the approximating equations, which causes

Explicit Scheme

We calculate the solution for a future time using the solution at the current time. Mathematically, if $c(t)$ is the solution at the current time and $c(t+1)$ is the solution at a future time, then for an explicit method we would have the following form of equation:

$$c(t+1) = f\left[c(t)\right]$$ where, $f\left[c(t)\right]$ is some known function of $c$ at current time $t$ and is known. Only unknown in this equation is $c(t+1)$ on left hand side.

Implicit Scheme

We calculate the solution by solving an equation which contains the solution at both the current time and future time. Using the same convention defined earlier, for an implicit scheme we would have the following equation: $$c(t+1) = c(t) + f[c(n+1)]$$

In this equation right hand side depends upon the quantity you are trying to calculate. Assuming $f(c)=Ac$, where $A$ is a matrix independent of $c$, then the implicit scheme requires solving the following system of linear equations:

$$(I-A).c(t+1) = c(t)$$ $$\Rightarrow c(t+1) = (I-A)^{-1}.c(t)$$ i.e. to solve for a solution at a future time you would need to perform matrix inversion.

Pros and Cons of Two Schemes

Explicit schemes are easy to program with only simple calculations performed at each timestep while implicit schemes are more difficult and can involve many iterations per timestep. Explicit schemes are more prone to be unstable and require very small timesteps compared to implicit schemes which are much stable and does not have any constraint on the size of a time step you can take.