Stiff "hell in the cell" championship: RK4 vs Samarskii's AOS scheme
The following matlab scripts illustrate the performance of the considered integrators for the 1D semi-discrete heat equation.
This example is of a 1D diffusion equation discretized in space using a central difference scheme. The time integration is done using the fourth-order Runge--Kutta method (explicit, conditional stability) and Samarskii's additive operator splitting method (implicit unconditionally stable).
To run the demo you have to execute the following at the MATLAB command prompt:
tdiffusion1d_rk4_vs_samarskii (DiffusionCoeff, TotalTime, NumberOfTimeSteps, NumberOfElements)
Examples:
tdiffusion1d_rk4_vs_samarskii (1e-3, 5, 200, 300) %% stable regime of RK4
tdiffusion1d_rk4_vs_samarskii (1e-3, 5, 200, 334) %% stability limit of RK4
tdiffusion1d_rk4_vs_samarskii (1e-3, 5, 200, 335) %% perturbations amplify for RK4
Time integration
The momentum and the stabilized mass balance equations discretized in space by the FEM leads to the following system of equations,
\begin{gather} \label{eq_mom_bal} \mathbf{M}\dot{\mathsf{V}} + \mathbf{K} \mathsf{V} + \mathbf{Q} \mathsf{P} = \mathsf{F} \\ \label{eq_mas_bal} (\mathbf{Q}^{\text{T}} - \mathbf{S}_\mathsf{qv}) \mathsf{V} - \mathbf{S}_\mathsf{qp} \mathsf{P} = \mathsf{F}_\mathsf{q} \end{gather}where, the matrices and vectors with sub-scripts are the ones that appear due to the stabilization of the mass balance equation. Let \(\mathsf{A}:= \dot{\mathsf{V}}\) denote the vector of unknown nodal accelerations. The momentum balance equation, i.e. Eq. (\ref{eq_mom_bal}), is a dynamic equation. On the contrary, the stabilized mass balance equation, i.e. Eq. (\ref{eq_mas_bal}), is a quasi-static equation.
The well-known Newmark algorithm for the numerical integration of the momentum balance equation is defined by,
\begin{gather} \label{eq_newmark_x} \mathsf{X}^{n+1} = \mathsf{X}^{n} + {\Delta t} \mathsf{V}^{n} + \frac{\Delta t^{2}}{2} \left[ (1-2\beta)\mathsf{A}^{n} + 2\beta \mathsf{A}^{n+1} \right] \\ \label{eq_newmark_v} \mathsf{V}^{n+1} = \mathsf{V}^{n} + {\Delta t} \left[ (1-\theta) \mathsf{A}^{n} + \theta \mathsf{A}^{n+1} \right] \\ \label{eq_newmark_a} \mathbf{M}\mathsf{A}^{n+1} + \mathbf{K}^{n+1} \mathsf{V}^{n+1} + \mathbf{Q}^{n+1} \mathsf{P}^{n+1} = \mathsf{F}^{n+1} \end{gather}Thus, the Newmark algorithm employs only two time levels and two parameters \(\beta, \theta\) that control the accuracy and damping properties of the scheme. Further, the dynamic equilibrium equation is satisfied at the end of every time step, cf. Eq. (\ref{eq_newmark_a}). Note that in the Lagrangian description of the kinematics, the mass matrix \(\mathbf{M}\) is an invariant entity. Likewise, if the force vector \(\mathsf{F}\) has contributions only from some body force, then even it is invariant.
The Bossak–Newmark algorithm which introduces an additional parameter \(\alpha\) for controlling the damping properties of the Newmark's algorithm, is defined by,
\begin{gather} \label{eq_bossak_x} \mathsf{X}^{n+1} = \mathsf{X}^{n} + {\Delta t} \mathsf{V}^{n} + \frac{\Delta t^{2}}{2} \left[ (1-2\beta)\mathsf{A}^{n} + 2\beta \mathsf{A}^{n+1} \right] \\ \label{eq_bossak_v} \mathsf{V}^{n+1} = \mathsf{V}^{n} + {\Delta t} \left[ (1-\theta) \mathsf{A}^{n} + \theta \mathsf{A}^{n+1} \right] \\ \label{eq_bossak_a} \mathbf{M}\left[ (1-\alpha)\mathsf{A}^{n+1} + \alpha \mathsf{A}^{n} \right] + \mathbf{K}^{n+1} \mathsf{V}^{n+1} + \mathbf{Q}^{n+1} \mathsf{P}^{n+1} = \mathsf{F}^{n+1} \end{gather}The only change in the above equations with respect to the original Newmark algorithm is in the satisfaction of an \(\alpha\)-modified dynamic equilibrium equation, i.e. Eq. (\ref{eq_bossak_a}). Using the Bossak–Newmark algorithm the balance equations given in Eqs. (\ref{eq_mom_bal}) and (\ref{eq_mas_bal}) are integrated in time (\(\{\mathsf{X}^{n}, \mathsf{V}^{n}, \mathsf{A}^{n}, \mathsf{P}^{n}\} \rightarrow \{\mathsf{X}^{n+1}, \mathsf{V}^{n+1}, \mathsf{A}^{n+1}, \mathsf{P}^{n+1}\}\)) as follows.
\begin{gather} \label{eq_bossak_step1a} \begin{bmatrix} \dfrac{(1-\alpha)}{\Delta t} \mathbf{M} + \theta \mathbf{K}^{n+1} & \theta \mathbf{Q}^{n+1} \\ \theta (\mathbf{Q}^{\text{T}} - \mathbf{S}_\mathsf{qv})^{n+1} & - \theta \mathbf{S}_\mathsf{qp}^{n+1} \end{bmatrix} \begin{Bmatrix} \mathsf{V}^{n+1} \\ \mathsf{P}^{n+1} \end{Bmatrix} = \begin{Bmatrix} \theta \mathsf{F}^{n+1} + \mathsf{R}^{n} \\ \theta \mathsf{F}_\mathsf{q}^{n+1} \end{Bmatrix} \\ \label{eq_bossak_step1b} \mathsf{R}^{n} := \frac{(1-\alpha)}{\Delta t} \mathbf{M} \mathsf{V}^{n} + (1-\alpha-\theta) \mathbf{M}\mathbf{A}^{n} \\ \label{eq_bossak_step2} \mathsf{A}^{n+1} = \frac{1}{\theta \Delta t} \left( \mathsf{V}^{n+1} - \mathsf{V}^{n} \right) - \frac{(1-\theta)}{\theta} \mathsf{A}^{n} \\ \label{eq_bossak_step3} \mathsf{X}^{n+1} = \mathsf{X}^{n} + {\Delta t} \mathsf{V}^{n} + \frac{\Delta t^{2}}{2} \left[ (1-2\beta)\mathsf{A}^{n} + 2\beta \mathsf{A}^{n+1} \right] \end{gather}Making the choice \(\alpha = 0\) in the above equations we get the system of equations integrated in time using the standard Newmark algorithm.