Dissipative time-stepping scheme for first order in time system

Question

When solving semi-discrete equations (originating from finite element models, for example), which are second-order in time of the form \begin{equation} M\ddot d + C\dot d + Kd = F, \end{equation} where $d$ is the solution vector, $M$, $C$ and $K$ are matrices, and $F$ is a vector, one can make use of methods that damp out spurious high-frequency oscillations such as damped versions of the Newmark method, HHT-$\alpha$, etc.

If one wishes to solve instead a system of the form \begin{equation} M\dot d + Kd = F, \end{equation} the obvious choice would be to use a generalized trapezoidal method. However, I am looking for a method that exhibits damping out of spurious high-frequency oscillations, as in the second-order case. I don't need a highly accurate method, but it would preferably be an explicit one.

Answer 1

There is an excellent family of time integration schemes that fit your description called Generalized Single Step Single Solve (GS4). The original work on the implicit methods for first order systems can be found in [1].

Here is the implicit algorithm:

\begin{equation} C \widetilde{\dot{d}} + K \widetilde{d} = \widetilde{F} \end{equation}

where the variables have been approximated between the $n$th and $(n+1)$th timestep as

\begin{align} & \widetilde{\dot{d}} = \dot{d}_n + \Lambda_6 W_1 ( \dot{d} _{n+1}- \dot{d} _n) \\ & \widetilde{d} = d_n + \Lambda_4 W_1\Delta t \dot{d}_n + \Lambda_5 W_2\Delta t( \dot{d} _{n+1}- \dot{d} _n) \\ & \widetilde{F} =(1-W_1)q_n+W_1 F_{n+1} \end{align}

Now one can solve for $ \Delta\dot{d} = \dot{d} _{n+1}- \dot{d} _n$ using

\begin{align} ( \Lambda_6 W_1 C &+ \Lambda_5 W_2 \Delta t K) \Delta \dot{d } \\ = & -C \dot{d} _n - K ( d _{n} +\Lambda_4 W_1 \Delta t \dot{d} _n) \\ &+ (1-W_1)F_n+W_1 F_{n+1} \end{align}

with the updates

\begin{equation} \dot{d} _{n+1} = \dot{d} _n + \Delta \dot{d } \end{equation} \begin{equation} d _{n+1} = d _{n} + \lambda_4 \Delta t \dot{d} _{n} + \lambda_5 \Delta t \Delta \dot{d } \end{equation}

where

\begin{align} &\Lambda_4 W_1 = \frac{1}{1+\rho_{\infty}} \nonumber \\ &\Lambda_5 W_2 = \frac{1}{(1+\rho_{\infty})(1+\rho_{\infty}^{s})}\nonumber \\ &\Lambda_6 W_1 = \frac{3+\rho_{\infty}+\rho_{\infty}^{s} - \rho_{\infty}\rho_{\infty}^{s}}{2(1+\rho_{\infty})(1+\rho_{\infty}^{s})} \\ &W_1 = \frac{1}{1+\rho_{\infty}} \nonumber \\ &\lambda_4 = 1, \quad \nonumber \lambda_5 = \frac{1}{1+\rho_{\infty}^{s}} \nonumber \end{align}

Admittedly, this may seem like a whole load of parameters and nonsense, but all parameters only depend upon two chosen values $(\rho_{\infty}, \rho_{\infty}^{s})$. The beauty is once you have it programmed you can choose from a whole family of algorithms just by choosing you values for $(\rho_{\infty}, \rho_{\infty}^{s})$. These parameters correspond to the eigenvalues of amplification matrix of the single DOF problem. Thus, you can choose the amount of damping (numerical dissipation) simply by choosing them. Some noteable choices: $(\rho_{\infty}, \rho_{\infty}^{s}) = (1,1)$ gives the Crank-Nicolson method (no damping, not for you) and $(\rho_{\infty}, \rho_{\infty}^{s}) = (0,0)$ gives an algorithm equivalent to the highly dissipative Gear's method (aka two-step backwards difference formula). Any choice will give you a second-order, unconditionally stable algorithm.

Note the restriction on your choices: $0 \leq \rho_{\infty}^{s} \leq \rho_{\infty} \leq 1$.

Now if you want an explicit algorithm, some algorithms have been developed using the same approach that led to the mess above. I don't think they have been published any place highly visible yet but early work can be found in a master's thesis here [2].

The easiest thing to do to get an explicit scheme, with the nice dissipation properties of the above algorithm, is to turn it into a predictor-corrector method. You, of course, lose the unconditional stability, but you will still have a second-order time integrator. To do so you can replace the $\widetilde{d}$ above with:

\begin{align} \widetilde{d} = d_n + \Lambda_4 W_1\Delta t \dot{d}_n \end{align}

and lump the $C$ matrix, then march to your hearts content. Everything else stays the same but the restrictions to the $\rho$ parameters above are lifted (they can be anything). The stability of the algorithm and the amount of dissipation still depends upon this choice. Hover around 0 and you should be okay.

[1] http://onlinelibrary.wiley.com/doi/10.1002/nme.3228/full

[2] http://conservancy.umn.edu/handle/11299/162393

Answer 2

These differential equations appear to be linear, so you should be able to solve them analytically. In order to do this you first have to write it as one matrix differential equation of the form,

$$ \dot{\vec{x}} = A\vec{x} + B \vec{u}. $$

This can be done by rewriting the differential equations to an expression for $\ddot{d}$,

$$ \ddot{d} = M^{-1} F - M^{-1} C \dot{d} - M^{-1} K d. $$

The vector $\vec{x}$ can be defined as $\begin{bmatrix}d & \dot{d}\end{bmatrix}^T$, such that $\dot{\vec{x}}$ is equal to $\begin{bmatrix}\dot{d} & \ddot{d}\end{bmatrix}^T$. By using this definition then the matrix $A$ can be written as,

$$ A = \begin{bmatrix} 0 & I \\ -M^{-1}K & -M^{-1}C \end{bmatrix}, $$

where $0$ and $I$ are the zero/null and identity matrix receptively, both of the same size as $M$, $C$ and $K$. And similar the matrix $B$ can be written as,

$$ B = \begin{bmatrix} 0 \\ M^{-1} \end{bmatrix}, $$

where $0$ is again the zero/null matrix of the same size as $M$, $C$ and $K$. And the vector $\vec{u}$ is equal to $F$.

The homogeneous or transient solution ($F=0$) can be found with,

$$ \vec{x}_t(t) = \sum{c_i\vec{v}_ie^{\lambda_it}}, $$

where $\lambda_i$ are eigenvalues of $A$, with $\vec{v}_i$ the corresponding eigenvectors and $c_i$ are constants which can be found with the help of the initial conditions.

The nonhomogeneous or steady state solution can be found by taking the Laplace transform,

$$ \mathcal{L}\{\dot{\vec{x}}_{ss}(t)\}(s) = A \mathcal{L}\{\vec{x}_{ss}(t)\}(s) + B \mathcal{L}\{\vec{u}(t)\}(s), $$

$$ s X_{ss}(s) = A X_{ss}(s) + B U(s), $$

$$ X_{ss}(s) = (s I - A)^{-1} B U(s), $$

$$ \vec{x}_{ss}(t) = \mathcal{L}^{-1}\{(s I - A)^{-1} B \mathcal{L}\{\vec{u}(t)\}(s)\}(t). $$

If the $\vec{u}$ is constant in time, then this simplifies down to,

$$ \vec{x}_{ss}(t) = -A^{-1} B \vec{u}. $$

The steady state solution can also be found by taking the Laplace transform of the original differential equations. The only difference between this and using $A$ is that the later also gives the derivative of $d_{ss}$ because you solve for $\vec{x}_{ss}(t)$. The total response can be found by adding the two solutions,

$$ \vec{x}(t) = \vec{x}_{ss}(t) + \vec{x}_{t}(t), $$

where the constants $c_i$ should only be solved for by using the initial or boundary conditions when $\vec{x}_{ss}(t)$ has also been added.