# Explicit Scheme Stability

Where, $\left[\text{A}\right]$ is amplification matrix. A spectral analysis of this matrix highlights the stability of the integration scheme.

The numerical algorithm is stable if and only if the radius spectral of $\left[\text{A}\right]$ is less than unity. In other words, when the module of all eigen values of $\left[\text{A}\right]$ are smaller than unity, the numerical stability is ensured.

Then, the equations are developed for the systems with or without damping. ^{1}

- ${A}_{1}=\frac{1}{2}tr\left[\text{A}\right]=\frac{1}{2}\left({A}_{11}+{A}_{22}\right)$
- ${A}_{2}=\mathrm{det}\left[\text{A}\right]={A}_{11}{A}_{22}-{A}_{12}{A}_{21}$

- Roots are real and one of them is equal to 1: You then have:
(6) $$1-2{A}_{1}+{A}_{2}=0$$This yields:(7) $$\begin{array}{l}{A}_{2}=2{A}_{1}-1\\ {\lambda}_{1}=1\\ {\lambda}_{2}={A}_{2}\end{array}$$The corresponding part of the boundary of the stability domain is the segment analytically defined by $1-2{A}_{1}+{A}_{2}=0$ and $-1\le {A}_{2}\le 1$ .

- Roots are real and one of them is equal to -1: You then have:
(8) $$1+2{A}_{1}+{A}_{2}=0$$This yields:(9) $$\begin{array}{l}{A}_{2}=-2{A}_{1}-1\\ {\lambda}_{1}=-1\\ {\lambda}_{2}=-{A}_{2}\end{array}$$In this case, the corresponding part of the boundary is the segment given by $1+2{A}_{1}+{A}_{2}=0$ and $-1\le {A}_{2}\le 1$ .

- Roots are complex conjugate: Their modulus is equal to 1. You then have, using ${\lambda}_{1,2}={e}^{\pm i\alpha}$ :
(10) $$\begin{array}{l}0={e}^{2i\alpha}-2{A}_{1}{e}^{i\alpha}+{A}_{2}\\ =\left(\mathrm{cos}2\alpha -2{A}_{1}\mathrm{cos}\alpha +{A}_{2}\right)+i\left(\mathrm{sin}2\alpha -2{A}_{1}\mathrm{sin}\alpha \right)\\ =\left(2\mathrm{cos}\alpha \left(\mathrm{cos}\alpha -{A}_{1}\right)+{A}_{2}-1\right)+i\left(2\mathrm{sin}\alpha \left(\mathrm{cos}\alpha -{A}_{1}\right)\right)\end{array}$$This yields:(11) $$\begin{array}{l}2\mathrm{cos}\alpha \left(\mathrm{cos}\alpha -{A}_{1}\right)+{A}_{2}-1=0\\ 2\mathrm{sin}\alpha \left(\mathrm{cos}\alpha -{A}_{1}\right)=0\end{array}$$Since $\mathrm{sin}\alpha \ne 0$ , you obtain:(12) $$\begin{array}{l}{A}_{1}=\mathrm{cos}\alpha \\ {A}_{2}=1\end{array}$$The corresponding part of the boundary is thus the segment given by ${A}_{2}=1$ and $-1\le {A}_{2}\le 1$ .

The 3 segments introduced above define a closed contour. Point ${A}_{1}={A}_{2}=0$ is located inside this contour and in this case, $\rho \left(\text{\hspace{0.05em}}\left[\text{A}\right]\text{\hspace{0.05em}}\right)=0$ . Since $\rho \left(\text{\hspace{0.05em}}\left[\text{A}\right]\text{\hspace{0.05em}}\right)$ varies continuously with respect to ${A}_{1}$ and ${A}_{2}$ , you can conclude that the stability domain corresponds to the interior of the contour. To precisely define the stability domain, you must also have points leading to double eigen value of modulus 1, that is, the intersections between the parabola ${A}_{1}^{2}={A}_{2}$ and the boundary of the domain. This corresponds to Points $\left({A}_{1},{A}_{2}\right)=\left(\pm 1,1\right)$ .You can analytically summarize the description of the stability by means of the following two sets of conditions:(13) $$\begin{array}{l}\left(1\right)\text{\hspace{1em}}-\frac{\left({A}_{2}+1\right)}{2}\le {A}_{1}\le \frac{\left({A}_{2}+1\right)}{2}\text{\hspace{0.17em}},\text{\hspace{1em}}-1\le {A}_{2}<1\\ \left(2\right)\text{\hspace{1em}}-1<{A}_{1}<1\text{\hspace{0.17em}},\text{\hspace{1em}}{A}_{2}=1\end{array}$$

## Numerical Stability of Undamped Systems

Where, $\omega =\sqrt{\frac{k}{m}}$ is the angular frequency of the considered mode.

## Numerical Stability with Viscous Damping: Velocities at Time Steps

This yields ${A}_{1}=\frac{1-\frac{{\omega}^{2}\text{\Delta}{t}^{2}}{2}}{1+\frac{c\text{\Delta}t}{2m}}$ and ${A}_{2}=\frac{1-\frac{c\text{\Delta}t}{2m}}{1+\frac{c\text{\Delta}t}{2m}}$ .

$-1\le \frac{1-\frac{c\text{\Delta}t}{2m}}{1+\frac{c\text{\Delta}t}{2m}}<1$

You find the same condition as in the undamped case, which echoes a conclusion given in
^{1}. You may yet remark that damping has changed the strict
inequality into a large inequality, preventing from weak instability due to a double eigen
value of modulus unity.

It is important to note that the relation Equation 29 is obtained by using the expression Equation 25 to compute nodal velocities at time steps. However, in an explicit scheme generally the mid-step velocities ${\dot{u}}^{n+\frac{1}{2}}$ and ${\dot{u}}^{n-\frac{1}{2}}$ are used. This will be studied in the next section.

## Numerical Stability with Viscous Damping: Velocities at Mid Steps

You have in this case ${A}_{1}=1-\frac{{\omega}^{2}\text{\Delta}{t}^{2}}{2}-\frac{c\text{\Delta}t}{2m}$ and ${A}_{2}=1-\frac{c\text{\Delta}t}{m}$ .

Therefore, the critical time step depends not only to $\omega $ but also to the mass and the damping. However, the critical time step depends only to $\omega $ when using the exact velocities to compute the viscous forces as described in the previous section.

## Numerical Stability with Rayleigh Damping

- $\left[\text{C}\right]$
- Viscous damping matrix of the system
- $\left[\text{M}\right]$
- Mass matrix
- $\left[\text{K}\right]$
- Stiffness matrix

- $m$
- Modal mass
- $c$
- Associated modal damping
- $k$
- Nodal stiffness

In this case, ${A}_{1}=\frac{1-\frac{{\omega}^{2}\text{\Delta}{t}^{2}}{2}-\frac{\beta {\omega}^{2}\text{\Delta}t}{2}}{1+\frac{\alpha \text{\Delta}t}{2}}$ and ${A}_{2}=\frac{1-\frac{\alpha \text{\Delta}t}{2}-\beta {\omega}^{2}\text{\Delta}t}{1+\frac{\alpha \text{\Delta}t}{2}}$ .

## Example: Critical Time Step for a Mass-Spring System

$M\ddot{X}+KX=0$ | (a) |

The element matrix expressions are given as:

$\left[M\right]=m\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$ ; $\left[K\right]=k\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]$

$X=Cos(\omega t-\alpha )\left[I\right]\Psi $ | (b) |

$\left(K-{\omega}^{2}M\right)\Psi Cos(\omega t-\alpha )=0$ | (c) |

$\mathrm{det}\left(K-{\omega}^{2}M\right)=0$ => ${\omega}^{2}=\frac{2k}{m}$ | (d) |

Assuming the following numerical values $m=1$ and $k=10$ , you have $\omega =\sqrt{\frac{2k}{m}}=4.472136$ .

The critical time step of the system is given by Equation 20:

$\text{\Delta}t\le \frac{2}{\omega}\text{\hspace{1em}}\Rightarrow \text{\hspace{1em}}\text{\Delta}t\le \frac{2}{4.472136}\text{\hspace{1em}}\Rightarrow \text{\hspace{1em}}\text{\Delta}t\le 0.4472$

## Example: Critical Time Step for Dropping Body

A dropping body is studied in Body Drop Example with analytical and numerical approaches. As shown in Figure 5, the numerical results are closed to the analytical solution if you use a step-by-step time discretization method with a constant time step $\text{\Delta}t=0.1$ . This value is less than the critical time step obtained by:

which may be computed as:

$\omega =\sqrt{\frac{k}{m}}=\sqrt{20}\Rightarrow \text{\hspace{1em}}\text{\Delta}{t}_{cr}=\frac{2}{4.472136}\text{\hspace{1em}}\Rightarrow \text{\hspace{1em}}\text{\Delta}{t}_{cr}=0.4472$

^{1}for more details on the computation of nodal time step).

^{1}Hughes T.J.R., “Analysis of Transient Algorithms with Particular Reference to Stability Behavior”, Computational Methods for Transient Analysis, eds. T. Belytschko and T.J.R. Hugues, 67-155, North Holland, 1983.