A critical assessment of two-stage composite time integration schemes with a unified set of time approximations

Abstract

A simple unified set of displacement and velocity approximations is presented to express various two-stage composite time integration schemes. Based on the unified approximations, two novel sets of optimized parameters are newly proposed for the implicit composite schemes to enhance the capability of conserving total energy. Two special cases of the unified approximations are also considered to overcome some shortcomings of the existing schemes. Besides, the newly proposed unified set of approximations can include many of the existing composite time integration schemes. To be specific, both implicit and explicit composite schemes can be expressed by using the unified set of approximations. Thus, both implicit and explicit types of composite schemes can be selected from the unified set by simply changing algorithmic parameters. To demonstrate advantageous features of the proposed unified set of approximations, various numerical examples are solved, and results are analyzed.

Keywords:
Direct time integration; structural dynamics; composite time scheme; Kim and Reddy method; generalized composite time integration algorithm; Bathe method

Graphical Abstract

1 Introduction

For decades, direct time integration schemes have been widely used for the effective analysis of linear and nonlinear structural dynamics [¹, ², ³, ⁴, ⁵, ⁶, ⁷, ⁸, ⁹, ¹⁰, ¹¹]. Recently, the concept of subdividing a complete time step (or time interval) into two sub-steps is frequently employed for the development of direct time integration schemes with enhanced numerical performances. Time integration schemes developed based on this strategy are often called the composite time integration scheme. Occasionally, they are also called the two-stage time integration scheme. Implicit and explicit schemes of this type are known to possess improved numerical performances when compared with the conventional single-stage time integration schemes [⁸, ¹², ¹³, ¹⁴, ¹⁵, ¹⁶, ¹⁷, ¹⁸, ¹⁹]. In composite schemes, two different kinds of schemes are employed in the first and second sub-steps, respectively. Then, a complete scheme is obtained by properly combining the two different schemes. Even though various composite schemes are designed to solve the equation of structural dynamics numerically, they have never been expressed in a unified form.

In a sense, direct time integration schemes are essential tools for the efficient and effective analysis of structural dynamics. Since the computational aspect of structural dynamics is heavily dependent on the characteristic of a given problem, it is very difficult to choose a perfect scheme that works nicely for all cases. For example, the long-term analysis of highly nonlinear dynamic problems may require a non-dissipative higher-order accurate time integration scheme to conserve the total energy of the system while minimizing the period error [¹⁴, ²⁰, ²¹, ²², ²³]. In this case, higher-order accurate explicit time integration schemes may become more efficient than implicit schemes, because matrix factorizations and iterative nonlinear solution finding procedures can be avoided if the mass matrix is diagonal in explicit schemes. On the other hand, the analysis of wave propagation and impact problems can be done more effectively by using dissipative implicit time integration schemes and large time steps. In general, dissipative implicit methods can eliminate a multiple number of spurious high-frequency modes more effectively than explicit methods, while dissipative explicit methods require less computational efforts.

Although categorizations of time integration schemes and their general numerical characteristics are relatively well studied in the literature [⁴, ⁷, ¹⁴, ²⁰, ²², ²⁴, ²⁵, ²⁶, ²⁷, ²⁸, ²⁹, ³⁰], it is still very difficult to assert that a certain type of time integration scheme is the most suitable for all problems. For example, many studies state that implicit schemes are more accurate than explicit schemes due to their unconditional stability. Strictly speaking, this generalization is incorrect. In fact, implicit methods are unconditionally stable, and they allow the use of time steps that are greater than the critical time step of explicit schemes. If time steps of implicit schemes are greater than the critical time step of explicit schemes, it can be said that implicit schemes are more accurate than explicit schemes. However, if a time step that is smaller than the critical time step of explicit schemes is used in both implicit and explicit schemes, than explicit schemes may give more accurate numerical solutions. This can easily be illustrated by using the trapezoidal rule and the central difference method. Considering these discussions, an optimal time integration scheme should be selected carefully based on the purpose of the analysis and the characteristic of a given problem.

Composite time integration schemes for structural dynamics have been developed based on various numerical methods and techniques. For this reason, they may have quite different numerical characteristics. In addition to various numerical characteristics, understanding computational structures and algorithmic parameters of various composite schemes developed based on different mathematical frameworks can be a time-consuming process for analysts. Considering this difficulty, a concise unified form that can include various implicit and explicit composite schemes under its umbrella could be of benefit to researchers whose main interest is focused on the analysis of various structural problems by using the existing numerical techniques. If a unified set of time approximation with different sets of optimized algorithmic parameters is properly implemented on computers, then various schemes would become more easily accessible to users with different research purposes. With a unified set of time approximations, various time integration schemes can be expressed and an optimal scheme can be selected by simply changing algorithmic parameters. The Newmark method [², ³¹, ³², ³³] is a good example of using unified expressions with adjustable algorithmic parameters. In the Newmark method, for example, the trapezoidal rule, the linear acceleration method, and the central different method are obtained by simply changing values of parameters.

The purpose of this article is to provide a novel unified set of time approximations which can include various implicit and explicit composite time integration schemes. Since the aspect of transient analyses heavily dependent on characteristics of dynamics problems, it is very difficult to find a scheme that works perfectly for every situation. With a unified mathematical framework, however, more effective analyses can be conducted by selecting an optimal scheme. To this end, a novel unified set of time approximations is considered for the displacement and velocity vectors. In the unified set of approximations, collocation parameters are considered to control locations of time points where the displacement and velocity vectors are approximated, and weighting parameters are used to truncate second- and higher-order time derivatives of the approximated variables in a systematic manner. These algorithmic parameters are precisely optimized to achieve improved numerical performances.

In this article, two novel sets of optimized parameters are proposed for the implicit composite schemes to enhance the total energy conserving capability. The implicit composite schemes with the newly optimized parameters can preserve the total energy of dynamic systems almost perfectly when compared to other cases. Besides, several interesting cases of implicit schemes are also considered to overcome some drawbacks of the exiting implicit schemes. One of them is designed to have fixed values of the splitting parameters (i.e., the collocation parameters), but different algorithmic parameters are used to overcome some disadvantages causes by varying the splitting ratio beyond the range of the current time interval. Another special family of the implicit schemes does not require the factorization of the mass matrix, because the initial acceleration vector is not required to start the procedure. In this family, the computational time of large linear system can be greatly reduced due to the absence of the factorization of the mass matrix.

In addition to various implicit composite schemes, many of the recently developed explicit composite schemes are also expressed in the form of the proposed unified set of time approximations. To illustrate that some schemes may become more effective than other schemes depending on the purpose of the analysis and the characteristic of a given problem, implicit and explicit composite schemes are precisely compared with each other by using various problems of structural dynamics, which was not done before for the composite type time schemes. Through the analysis of numerical experiments, specific guide lines of choosing a proper time integration scheme are provided. Based on these discussions, an optimal scheme can be selected from the proposed unified set of time approximations by simply changing algorithmic parameters.

2 A unified set of time approximations

To express various implicit composite schemes, a unified set of discrete equations is used for the time approximations of the displacement and velocity vectors over the time interval $t_s\le t\le t_s+\Delta t$, where $t_s$ and $\Delta t$ are the beginning of the time interval and the size of the time interval, respectively. $\Delta t$ is also called the size of the time step or just the time step.

To make the composite schemes more intuitive and easier to use, the velocity and displacement vectors are approximated by truncating their second- and higher-order time derivatives in a systematic manner. At the time point $t=t_s+{\tau }_i\Delta t$, the velocity and displacement vectors are approximated by using the weighted sums of the velocity and acceleration vectors, respectively as

${}^i{\dot{\mathbf{u}}}_{t_s+{\tau }_i\Delta t}={\dot{\mathbf{u}}}_{t_s}+\left({\tau }_i\Delta t\right)\sum _{j=0}^i\left({\beta }_{ij}{}^j{\ddot{\mathbf{u}}}_{t_s+{\tau }_j\Delta t}\right)\mathrm{f}\mathrm{o}\mathrm{r}i=1,2,3$ (1)

${}^i{\mathbf{u}}_{t_s+{\tau }_i\Delta t}={\mathbf{u}}_{t_s}+\left({\tau }_i\Delta t\right)\sum _{j=0}^i\left({\alpha }_{ij}{}^j{\dot{\mathbf{u}}}_{t_s+{\tau }_j\Delta t}\right)\mathrm{f}\mathrm{o}\mathrm{r}i=1,2,3$ (2)

where $i$ of ${}^i{\dot{\mathbf{u}}}_{t_s+{\tau }_i\Delta t}$ and ${}^i{\mathbf{u}}_{t_s+{\tau }_i\Delta t}$ denotes that the approximations belong to the $i$th stage, ${\tau }_i$ is the collocation parameter that determines the $i$th time point within the time interval, and ${\alpha }_{ij}$ and ${\beta }_{ij}$ are algorithmic weighting parameters for the velocity and acceleration vectors, respectively. It is noted that ${\tau }_0$ and ${\tau }_3$ should always be chosen as

${\tau }_0=0,{\tau }_3=1$ (3)

In the current approach, ${\tau }_1$ and ${\tau }_2$ can be chosen as

${\tau }_1\in \left[0,1\right],{\tau }_2\in \left[\mathrm{0,1}\right]$ (4)

${\tau }_1$ of the proposed approximations plays the role of the splitting ratio of the existing composite schemes. In the existing composite schemes [¹²,¹⁴,¹⁵,³⁴], the splitting ratio cannot become zero or one, but the proposed time approximations do not have this restriction. In Eqs. (1)-(2), the weighting parameters ${\alpha }_{ij}$ and ${\beta }_{ij}$ should also satisfy

$\sum _{j=0}^i{\alpha }_{ij}=\sum _{j=0}^i{\beta }_{ij}=1\mathrm{f}\mathrm{o}\mathrm{r}i=\mathrm{1,2},3$ (5)

and

${\beta }_{33}=0$ (6)

With the restriction on ${\beta }_{33}$ given in Eq. (6), the unified set of approximations does not include the unknown acceleration vector in the third stage (i.e., for $i=3$), and the variables of the third stage can simply be computed without a major matrix computation. By using the properties given in Eqs. (5)-(6), one of the weighting parameters of the $i$th stage can be stated in terms of other weighting parameters of the $i$th stage. Since ${\beta }_{33}$ is always zero, ${}^3{\ddot{\mathbf{u}}}_{t_s+{\tau }_3\Delta t}$ is not included in the approximations given in Eqs. (1)-(2). ${}^1{\ddot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}$ and ${}^2{\ddot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}$ are the two unknown acceleration vectors to be determined by solving two sets of algebraic equations. If ${}^1{\ddot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}$ and ${}^2{\ddot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}$ are determined, then the displacement and velocity vectors can also be computed according to Eqs. (1)-(2). To determine these acceleration vectors, two sets of the dynamic equilibrium equations are required. The dynamic equilibrium equation at $t=t_s+{\tau }_i\Delta t$ is expressed as

$\mathbf{M}{}^i{\ddot{\mathbf{u}}}_{t_s+{\tau }_i\Delta t}=\mathbf{f}\left({}^i{\mathbf{u}}_{t_s+{\tau }_i\Delta t},{}^i{\dot{\mathbf{u}}}_{t_s+{\tau }_i\Delta t},t_s+{\tau }_i\Delta t\right)\mathrm{f}\mathrm{o}\mathrm{r}i=1,2$ (7)

where $\mathbf{M}$ is the mass matrix, and $\mathbf{f}$ is the total force vector. In linear structural dynamics, $\mathbf{f}$ is often expressed as

$\mathbf{f}\left({}^i{\mathbf{u}}_{t_s+{\tau }_i\Delta t},{}^i{\dot{\mathbf{u}}}_{t_s+{\tau }_i\Delta t},t_s+{\tau }_i\Delta t\right)=\mathbf{q}\left(t_s+{\tau }_i\Delta t\right)-\mathbf{C}{}^i{\dot{\mathbf{u}}}_{t_s+{\tau }_i\Delta t}-\mathbf{K}{}^i{\mathbf{u}}_{t_s+{\tau }_i\Delta t}$ (8)

where $\mathbf{q}$ is the vector of externally applied force, $\mathbf{C}$ is the damping matrix, $\mathbf{K}$ is the stiffness matrix. To determine the two unknown acceleration vector, Eq. (7) should be solved by using Eqs. (1)-(2). Again, it should be emphasized that the approximations given in Eqs. (1)-(2) do not contain ${}^3{\ddot{\mathbf{u}}}_{t_s+{\tau }_3\Delta t}$. Thus, only two equation solving procedures are required to find ${}^1{\ddot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}$ and ${}^2{\ddot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}$ in this type of composite scheme.

This composite scheme can be summarized as follows. In the first sub-step, ${}^1{\ddot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}$ is computed by using

${}^1{\dot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}={\dot{\mathbf{u}}}_{t_s}+\left({\tau }_1\Delta t\right)\left({\beta }_{10}{\ddot{\mathbf{u}}}_{t_s}+{\beta }_{11}{}^1{\ddot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}\right)$ (9)

${}^1{\mathbf{u}}_{t_s+{\tau }_1\Delta t}={\mathbf{u}}_{t_s}+\left({\tau }_1\Delta t\right)\left({\alpha }_{10}{\dot{\mathbf{u}}}_{t_s}+{\alpha }_{11}{}^1{\dot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}\right)$ (10)

$\mathbf{M}{}^1{\ddot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}=\mathbf{f}\left({}^1{\mathbf{u}}_{t_s+{\tau }_1\Delta t},{}^1{\dot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t},t_s+{\tau }_1\Delta t\right)$ (11)

In the second sub-step, ${}^2{\ddot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}$ is computed by using

${}^2{\dot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}={\dot{\mathbf{u}}}_{t_s}+\left({\tau }_2\Delta t\right)\left({\beta }_{20}{\ddot{\mathbf{u}}}_{t_s}+{\beta }_{21}{}^1{\ddot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}+{\beta }_{22}{}^2{\ddot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}\right)$ (12)

${}^2{\mathbf{u}}_{t_s+{\tau }_2\Delta t}={\mathbf{u}}_{t_s}+\left({\tau }_2\Delta t\right)\left({\alpha }_{20}{\dot{\mathbf{u}}}_{t_s}+{\alpha }_{21}{}^1{\dot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}+{\alpha }_{22}{}^2{\dot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}\right)$ (13)

$\mathbf{M}{}^2{\ddot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}=\mathbf{f}\left({}^2{\mathbf{u}}_{t_s+{\tau }_2\Delta t},{}^2{\dot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t},t_s+{\tau }_2\Delta t\right)$ (14)

By using the ${}^1{\ddot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}$ and ${}^2{\ddot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}$, ${}^3{\dot{\mathbf{u}}}_{t_s+{\tau }_3\Delta t}$ and ${}^3{\mathbf{u}}_{t_s+{\tau }_3\Delta t}$ are updated as

${}^3{\dot{\mathbf{u}}}_{t_s+{\tau }_3\Delta t}={\dot{\mathbf{u}}}_{t_s}+\left({\tau }_3\Delta t\right)\left({\beta }_{30}{\ddot{\mathbf{u}}}_{t_s}+{\beta }_{31}{}^1{\ddot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}+{\beta }_{32}{}^2{\ddot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}+{\beta }_{33}{}^3{\ddot{\mathbf{u}}}_{t_s+{\tau }_3\Delta t}\right)$ (15)

${}^3{\mathbf{u}}_{t_s+{\tau }_3\Delta t}={\mathbf{u}}_{t_s}+\left({\tau }_3\Delta t\right)\left({\alpha }_{30}{\dot{\mathbf{u}}}_{t_s}+{\alpha }_{31}{}^1{\dot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}+{\alpha }_{32}{}^2{\dot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}+{\alpha }_{33}{}^3{\dot{\mathbf{u}}}_{t_s+{\tau }_3\Delta t}\right)$ (16)

where ${\tau }_3=1$, ${\beta }_{33}=0$, and the computation of ${}^3{\ddot{\mathbf{u}}}_{t_s+{\tau }_3\Delta t}$ is not required. To advance another step, ${\ddot{\mathbf{u}}}_{t_s}$, ${\dot{\mathbf{u}}}_{t_s}$, and ${\mathbf{u}}_{t_s}$ of the next time step are updated by using ${}^2{\ddot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}$, ${}^3{\dot{\mathbf{u}}}_{t_s+{\tau }_3\Delta t}$, and ${}^3{\mathbf{u}}_{t_s+{\tau }_3\Delta t}$, respectively. Then, procedures given in Eqs. (9)-(16) can be repeated again. It is noted that various different kinds of time integration schemes can be obtained by determining the algorithmic parameters (i.e., ${\tau }_i$, ${\alpha }_{ij}$ and ${\beta }_{ij}$) included in Eqs. (9)-(16).

3 Implicit schemes

To develop improved implicit composite time schemes with preferable numerical characteristics by using Eqs. (9)-(16), the algorithmic parameters ${\tau }_i$, ${\alpha }_{ij}$ and ${\beta }_{ij}$ should be determined carefully. To optimize these algorithmic parameters, the single-degree-of-freedom problem is used [⁴,⁵,²⁵]. The single-degree-of-freedom problem is expressed as

$\ddot{u}\left(t\right)+2\xi \omega \dot{u}\left(t\right)+{\omega }^{2}u\left(t\right)=q\left(t\right)$ (17)

where $u$ is the displacement, $q$ is the externally applied force, $\omega$ is the natural frequency, and $\xi$ is the damping ratio. The initial conditions are given by

$u\left(0\right)=u_0,\dot{u}\left(0\right)={\dot{u}}_0$ (18)

By setting $\mathbf{M}=1$, $\mathbf{C}=2\xi \omega$, $\mathbf{K}={\omega }^2$, $\mathbf{u}=u$, $\mathbf{q}=q=0$, and $t_s=0$, the newly proposed approximations given in Eqs. (9)-(16) can be applied to Eqs. (17)-(18). Application of Eqs. (9)-(16) to Eqs. (17)-(18) gives

${\mathbf{x}}_{\Delta t}=\mathbf{A}{\mathbf{x}}_0$ (19)

where ${\mathbf{x}}_{\Delta t}=\{u_{{\tau }_3\Delta t},{\dot{u}}_{{\tau }_3\Delta t},{\ddot{u}}_{{\tau }_2\Delta t}{\}}^T$, ${\mathbf{x}}_0=\{u_0,{\dot{u}}_0,{\ddot{u}}_0{\}}^T$, and $\mathbf{A}$ is called the amplification matrix. In this study, important numerical characteristics, such as the accuracy, the stability, and the algorithmic dissipation, are also investigated by using the single-degree-of-freedom problem given in Eqs. (17)-(18) and the amplification matrix given in Eq. (19). For example, the local truncation error [⁵] of a scheme is defined by

${\tau }_e=\frac{\sum _{i=0}^3{\left(-1\right)}^i{A}_{i}u\left(t_s+(1-i)\Delta t\right)}{{\Delta t}^2}$ (20)

where $A_0=1$, $A_1$ is the trace of $\mathbf{A}$, $A_2$ is the sum of the principal minors of $\mathbf{A}$, and $A_3$ is the determinant of $\mathbf{A}$. An algorithm is $k$th-order accurate if

${\tau }_e=O\left(\Delta t^k\right)$ (21)

is satisfied. In time integration schemes, the spectral radius of $\mathbf{A}$ is defined by

$\rho \left(\mathbf{A}\right)=\mathrm{m}\mathrm{a}\mathrm{x}\left(\left|{\lambda }_1\right|,\left|{\lambda }_2\right|,\left|{\lambda }_3\right|\right)$ (22)

where ${\lambda }_i$ is the $i$th eigenvalue of $\mathbf{A}$. The spectral radius defined in Eq. (22) is a direct measurement of the stability and the numerical dissipation. It is noted that $\rho \left(\mathbf{A}\right)$ is a function of $\Omega =\omega \Delta t$. A scheme is stable if

$\rho \left(\Omega \right)\le 1$ (23)

is satisfied. If Eq. (23) is satisfied for all $\Delta t\in [0,\mathrm{\infty })$, a scheme is unconditionally stable for linear problems. In implicit schemes, the level of numerical dissipation is often measured by the ultimate spectral radius. The ultimate spectral radius is defined by

${\rho }_{\mathrm{\infty }}=\underset{\Delta t\to \mathrm{\infty }}{\text{lim}}\rho \left(\mathbf{A}\right)$ (24)

By stating one of the algorithmic parameters of a scheme in terms of ${\rho }_{\mathrm{\infty }}$ according to Eq. (24), desired level of numerical dissipation in the high-frequency limit can be prescribed through ${\rho }_{\mathrm{\infty }}$. The most dissipative case is obtained when ${\rho }_{\mathrm{\infty }}$ is set as zero, and less dissipative cases are obtained as ${\rho }_{\mathrm{\infty }}$ approaches to one. When ${\rho }_{\mathrm{\infty }}=1$ is used, the non-dissipative case is obtained. The non-dissipative case is usually suitable for most of general analyses. Thus, the case of ${\rho }_{\mathrm{\infty }}=1$ is considered standard. The most dissipative case obtained with ${\rho }_{\mathrm{\infty }}=0$ is also called the asymptotic annihilating case. The case of ${\rho }_{\mathrm{\infty }}=0$ is useful for the elimination of the spurious high-frequency modes.

3.1 The first family of implicit schemes

In the first family of implicit schemes, ${}^3{\dot{\mathbf{u}}}_{t_s+{\tau }_3\Delta t}$ and ${}^3{\mathbf{u}}_{t_s+{\tau }_3\Delta t}$ given in Eqs. (15)-(16) are treated as ${}^2{\dot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}$ and ${}^2{\mathbf{u}}_{t_s+{\tau }_2\Delta t}$ given in Eqs. (12)-(13) to eliminate the spurious root of the characteristic polynomial equation of the amplification matrix $\mathbf{A}$ given in Eq. (19). This condition can be imposed as

$\begin{array}{r}{\beta }_{3j}={\beta }_{2j},{\alpha }_{3j}={\alpha }_{2j}\mathrm{f}\mathrm{o}\mathrm{r}j=0,1,2\\ {\alpha }_{33}=0\end{array}$ (25)

With the relations given in Eq. (25), the amplification matrix have two eigenvalues. To make two eigenvalues of the amplification matrix complex conjugate, ${\beta }_{ij}$ and ${\alpha }_{ij}$ can be related as

$\begin{array}{r}{\beta }_{ij}={\alpha }_{ij}\mathrm{f}\mathrm{o}\mathrm{r}i=\mathrm{1,2},3\mathrm{a}\mathrm{n}\mathrm{d}j=\mathrm{0,1},\mathrm{2,3}\end{array}$ (26)

To obtain ${\ddot{\mathbf{u}}}_{t_s+\Delta t}$, ${\dot{\mathbf{u}}}_{t_s+\Delta t}$, and ${\mathbf{u}}_{t_s+\Delta t}$ from ${}^2{\ddot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}$, ${}^3{\dot{\mathbf{u}}}_{t_s+{\tau }_3\Delta t}$, and ${}^3{\mathbf{u}}_{t_s+{\tau }_3\Delta t}$, ${\tau }_2$ and ${\tau }_3$ should be chosen as

$\begin{array}{r}{\tau }_2={\tau }_3=1\end{array}$ (27)

From Eq. (5), ${\alpha }_{10}$ and ${\alpha }_{20}$ can be determined as

$\begin{array}{r}{\alpha }_{10}=1-{\alpha }_{11}\\ {\alpha }_{20}=1-{\alpha }_{21}-{\alpha }_{22}\end{array}$ (28)

When the parameters are determined according to Eqs. (25)-(28), the new approximations given in Eqs. (11)-(16) can be simplified as

${}^1{\dot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}={\dot{\mathbf{u}}}_{t_s}+\left({\tau }_1\Delta t\right)\left(\left(1-{\alpha }_{11}\right){\ddot{\mathbf{u}}}_{t_s}+{\alpha }_{11}{}^1{\ddot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}\right)$ (29a)

${}^1{\mathbf{u}}_{t_s+{\tau }_1\Delta t}={\mathbf{u}}_{t_s}+\left({\tau }_1\Delta t\right)\left(\left(1-{\alpha }_{11}\right){\dot{\mathbf{u}}}_{t_s}+{\alpha }_{11}{}^1{\dot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}\right)$ (29b)

$\mathbf{M}{}^1{\ddot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}=\mathbf{f}\left({}^1{\mathbf{u}}_{t_s+{\tau }_1\Delta t},{}^1{\dot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t},t_s+{\tau }_1\Delta t\right)$ (29c)

${}^2{\dot{\mathbf{u}}}_{t_s+\Delta t}={\dot{\mathbf{u}}}_{t_s}+\Delta t\left(\left(1-{\alpha }_{21}-{\alpha }_{22}\right){\ddot{\mathbf{u}}}_{t_s}+{\alpha }_{21}{}^1{\ddot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}+{\alpha }_{22}{}^2{\ddot{\mathbf{u}}}_{t_s+\Delta t}\right)$ (29d)

${}^2{\mathbf{u}}_{t_s+\Delta t}={\mathbf{u}}_{t_s}+\Delta t\left(\left(1-{\alpha }_{21}-{\alpha }_{22}\right){\dot{\mathbf{u}}}_{t_s}+{\alpha }_{21}{}^1{\dot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}+{\alpha }_{22}{}^2{\dot{\mathbf{u}}}_{t_s+\Delta t}\right)$ (29e)

$\mathbf{M}{}^2{\ddot{\mathbf{u}}}_{t_s+\Delta t}=\mathbf{f}\left({}^2{\mathbf{u}}_{t_s+\Delta t},{}^2{\dot{\mathbf{u}}}_{t_s+\Delta t},t_s+\Delta t\right)$ (29f)

${}^3{\dot{\mathbf{u}}}_{t_s+\Delta t}={}^2{\dot{\mathbf{u}}}_{t_s+\Delta t}$ (29g)

$\begin{array}{r}{}^3{\mathbf{u}}_{t_s+\Delta t}={}^2{\mathbf{u}}_{t_s+\Delta t}\end{array}$ (29h)

To guarantee second-order accuracy, ${\alpha }_{21}$ should be chosen to satisfy the condition given in Eq. (21) for $k=2$. From this condition, ${\alpha }_{21}$ is computed as

$\begin{array}{r}{\alpha }_{21}=-\frac{2{\alpha }_{22}-1}{2{\tau }_1}\end{array}$ (30)

With the relation given in Eq. (30), the ultimate spectral radius ${\rho }_{\mathrm{\infty }}$ defined in Eq. (24) is computed as

$\begin{array}{r}{\rho }_{\mathrm{\infty }}=\frac{2{\alpha }_{11}{\alpha }_{22}{\tau }_1-2{\alpha }_{11}{\tau }_1-2{\alpha }_{22}+1}{2{\alpha }_{11}{\alpha }_{22}{\tau }_1}\end{array}$ (31)

For the first family of implicit schemes, the level of numerical dissipation in the high-frequency limit can be determined by stating one of the remaining algorithmic parameters in terms of ${\rho }_{\mathrm{\infty }}$ according to Eq. (31). When ${\alpha }_{22}\ne {\alpha }_{11}$, ${\alpha }_{22}$ can be determined from Eq. (31) as

$\begin{array}{r}{\alpha }_{22}=-\frac{2{\alpha }_{11}{\tau }_1-1}{2\left({\alpha }_{11}{\tau }_1{\rho }_{\mathrm{\infty }}-{\alpha }_{11}{\tau }_1+1\right)}\end{array}$ (32)

For the dissipation control, ${\rho }_{\mathrm{\infty }}$ can be chosen from zero to one. If ${\rho }_{\mathrm{\infty }}$ is chosen in the range of $0\le {\rho }_{\mathrm{\infty }}\le 1$, the implicit scheme is unconditionally stable for linear problems. The choice of ${\rho }_{\mathrm{\infty }}=0$ gives the most dissipative case, and the choice of ${\rho }_{\mathrm{\infty }}=1$ gives the non-dissipative case. It should be noted that the relation given in Eq. (32) and the complete set of algorithmic parameter presented as the first general form in Table 1 can include most of the existing implicit composite schemes as its special cases.

To obtained an identical effective system matrix for the first and second sub-steps in linear problems, ${\alpha }_{11}$ can be chosen as

$\begin{array}{r}{\alpha }_{11}={\alpha }_{22}\end{array}$ (33)

When ${\alpha }_{11}$ is chosen according to Eq. (33), ${\alpha }_{22}$ can be determined from Eq. (31) as

$\begin{array}{r}{\alpha }_{22}=-\frac{{\tau }_1+1-\sqrt{{{\tau }_1}^2+2{\tau }_1{\rho }_{\mathrm{\infty }}+1}}{2{\tau }_1({\rho }_{\mathrm{\infty }}-1)}\end{array}$ (34)

The algorithmic parameters of various implicit composite time integration are presented in Table 1. Recently, the cases 1-1 and 1-2 presented in Table 1 are proposed by Kim and Reddy [²³]. The cases 1-1 and 1-2 can handle the externally applied forces as simply as the Newmark method does, because they use one identical time point (i.e., $t=t_s+\Delta t$) for the first and second dynamic equilibriums. The case 1-3 is equivalent to the scheme proposed by Kim and Reddy [³⁵]. The case 1-4 is equivalent to the scheme proposed by Kim and Choi [¹⁵]. It should be noted that the case 1-3 can include the standard Bathe method [⁸,³⁶] and the ${\rho }_{\mathrm{\infty }}$-Bathe method [³⁴]. Thus, the case 1-3 can be considered as a generalized expression for the existing implicit composite schemes. Other than the cases presented in Table 1, many other schemes can also be developed by using the first general form given in Table 1.

Table 1
Algorithmic parameters of the first family of implicit schemes.

	the first general form	case 1-1	case 1-2	case 1-3	case 1-4
${\tau }_0$	0	0	0	0	0
${\tau }_1$	adjustable	1	1	adjustable	$\frac{1}{2}$
${\tau }_2$	1	1	1	1	1
${\tau }_3$	1	1	1	1	1
${\alpha }_{10}$	$1-{\alpha }_{11}$	$1-{\alpha }_{11}$	$\frac{2{\rho }_{\mathrm{\infty }}-\sqrt{2{\rho }_{\mathrm{\infty }}+2}}{2({\rho }_{\mathrm{\infty }}-1)}$	$\frac{1}{2}$	$\frac{{\rho }_{\mathrm{\infty }}+1-\sqrt{2{\rho }_{\mathrm{\infty }}+2}}{{\rho }_{\mathrm{\infty }}-1}$
${\alpha }_{11}$	adjustable	adjustable	$\frac{-2+\sqrt{2{\rho }_{\mathrm{\infty }}+2}}{2({\rho }_{\mathrm{\infty }}-1)}$	$\frac{1}{2}$	$\frac{-2+\sqrt{2{\rho }_{\mathrm{\infty }}+2}}{{\rho }_{\mathrm{\infty }}-1}$
${\alpha }_{20}$	$\frac{2{\alpha }_{11}{\tau }_1{\rho }_{\mathrm{\infty }}-{\alpha }_{11}{\rho }_{\mathrm{\infty }}-{\alpha }_{11}+1}{2\left({\alpha }_{11}{\tau }_1{\rho }_{\mathrm{\infty }}-{\alpha }_{11}{\tau }_1+1\right)}$	$\frac{1}{2}$	$\frac{1}{2}$	$\frac{2{\tau }_1{\rho }_{\mathrm{\infty }}-{\rho }_{\mathrm{\infty }}+1}{2({\tau }_1{\rho }_{\mathrm{\infty }}+2-{\tau }_1)}$	$\frac{-2+\sqrt{2{\rho }_{\mathrm{\infty }}+2}}{2({\rho }_{\mathrm{\infty }}-1)}$
${\alpha }_{21}$	$\frac{({\rho }_{\mathrm{\infty }}+1){\alpha }_{11}}{2\left({\alpha }_{11}{\tau }_1{\rho }_{\mathrm{\infty }}-{\alpha }_{11}{\tau }_1+1\right)}$	$\frac{({\rho }_{\mathrm{\infty }}+1){\alpha }_{11}}{2\left({\alpha }_{11}{\rho }_{\mathrm{\infty }}-{\alpha }_{11}+1\right)}$	$\frac{{\rho }_{\mathrm{\infty }}+1-\sqrt{2{\rho }_{\mathrm{\infty }}+2}}{2({\rho }_{\mathrm{\infty }}-1)}$	$\frac{{\rho }_{\mathrm{\infty }}+1}{2({\tau }_1{\rho }_{\mathrm{\infty }}+2-{\tau }_1)}$	$\frac{{\rho }_{\mathrm{\infty }}+1-\sqrt{2{\rho }_{\mathrm{\infty }}+2}}{{\rho }_{\mathrm{\infty }}-1}$
${\alpha }_{22}$	$-\frac{2{\alpha }_{11}{\tau }_1-1}{2\left({\alpha }_{11}{\tau }_1{\rho }_{\mathrm{\infty }}-{\alpha }_{11}{\tau }_1+1\right)}$	$\frac{-2{\alpha }_{11}+1}{2\left({\alpha }_{11}{\rho }_{\mathrm{\infty }}-{\alpha }_{11}+1\right)}$	$\frac{-2+\sqrt{2{\rho }_{\mathrm{\infty }}+2}}{2({\rho }_{\mathrm{\infty }}-1)}$	$\frac{1-{\tau }_1}{{\tau }_1{\rho }_{\mathrm{\infty }}+2-{\tau }_1}$	$\frac{-2+\sqrt{2{\rho }_{\mathrm{\infty }}+2}}{2({\rho }_{\mathrm{\infty }}-1)}$
${\alpha }_{30}$	${\alpha }_{20}$	${\alpha }_{20}$	${\alpha }_{20}$	${\alpha }_{20}$	${\alpha }_{20}$
${\alpha }_{31}$	${\alpha }_{21}$	${\alpha }_{21}$	${\alpha }_{21}$	${\alpha }_{21}$	${\alpha }_{21}$
${\alpha }_{32}$	${\alpha }_{22}$	${\alpha }_{22}$	${\alpha }_{22}$	${\alpha }_{22}$	${\alpha }_{22}$
${\alpha }_{33}$	0	0	0	0	0
${\beta }_{ij}$	${\alpha }_{ij}$

In the existing composite time integration schemes, the size of the first sub-step is determined through the algorithmic parameter called the splitting ratio, while the size of the second sub-step is the same as a complete time interval. In the unified form, these existing schemes are expressed as the case 1-3 given in Table 1. The splitting ratio is defined as the ratio of the size of the first sub-step ($\tau \Delta t$) to the complete time step ($\Delta t$), and ${\tau }_1$ plays exactly the same role in the case 1-3. In the existing time integration methods, the splitting ratio is adjusted to control some spectral properties such as the level of numerical dissipation in the low-frequency range. In fact, the same spectral properties can also be obtained through alternative algorithmic parameters in the proposed unified set. In the case 1-1 given in Table 1, for example, ${\alpha }_{11}$ may be adjusted to have the effect which is the same as adjusting the splitting ratio in the existing schemes.

Figure 1
(a) Spectral radius of the case 1-1 for varying values of $\Delta t/T$. (b) Period error of the case 1-1 for varying values of $\Delta t/T$.

As shown in Fig. 1 (a) with the case 1-1, the level of numerical dissipation in the high-frequency limit can be adjusted through ${\rho }_{\mathrm{\infty }}$, and the level of numerical dissipation in the low-frequency range can also be adjusted through ${\alpha }_{11}$. For the case 1-1, ${\alpha }_{11}$ can be chosen in the range of $\frac{1}{4}\le {\alpha }_{11}<\frac{1}{2}$ or $\frac{1}{2}<{\alpha }_{11}<1$. When ${\alpha }_{11}=\frac{1}{4}$, the spectral properties of the case 1-1 become similar to those of the case 1-3 with ${\tau }_1=\frac{1}{2}$. Like the choice of ${\tau }_1=\frac{1}{2}$ is the standard setting in the existing implicit composite schemes, the choice of ${\alpha }_{11}=\frac{1}{4}$ is considered standard in the case 1-1. As ${\alpha }_{11}$ approaches to $\frac{1}{2}$, weaker numerical damping is introduced into the low-frequency range. As ${\alpha }_{11}$ approaches to $\frac{1}{2}$, the period and damping errors of the low-frequency range become almost identical to those of the trapezoidal rule with $\Delta t$ while a specified level of numerical dissipation can be imposed in the high-frequency limit through ${\rho }_{\mathrm{\infty }}$. On the other hand, stronger numerical damping is introduced into the low-frequency range when ${\alpha }_{11}$ is chosen in the range of $\frac{1}{2}<{\alpha }_{11}<1$. As ${\alpha }_{11}$ approaches to one, strong numerical damping is introduced into the low-frequency range. However, both period and damping errors increase as ${\alpha }_{11}$ approaches to one. Useful discussions regarding the role of ${\tau }_1$ are presented in Refs. [¹²,³⁴,³⁵].

3.2 The second family of implicit schemes

In the second family of implicit schemes, ${}^3{\dot{\mathbf{u}}}_{t_s+{\tau }_3\Delta t}$ and ${}^3{\mathbf{u}}_{t_s+{\tau }_3\Delta t}$ given in Eqs. (15)-(16) are not treated as ${}^2{\dot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}$ and ${}^2{\mathbf{u}}_{t_s+{\tau }_2\Delta t}$ given in Eqs. (12)-(13), which is different from the case of the first implicit family. To eliminate the spurious root, on the other hand, the acceleration vector of the previous step (i.e., ${\ddot{\mathbf{u}}}_{t_s}$) is excluded in the time approximations. For this reason, ${\alpha }_{i0}$ and ${\beta }_{i0}$ are chosen as

$\begin{array}{r}{\alpha }_{i0}={\beta }_{i0}=0\mathrm{f}\mathrm{o}\mathrm{r}i=1,2,3\end{array}$ (35)

To make two eigenvalues of the amplification matrix $\mathbf{A}$ become complex conjugate, ${\alpha }_{ij}$ and ${\beta }_{ij}$ can be related according to Eq. (26), and ${\alpha }_{33}$ is chosen as

$\begin{array}{r}{\alpha }_{33}=0\end{array}$ (36)

From Eqs. (5) and (35), ${\alpha }_{11}$, ${\alpha }_{22}$, and ${\alpha }_{32}$ of the second family can be determined as

$\begin{array}{r}{\alpha }_{11}=1\\ {\alpha }_{22}=1-{\alpha }_{21}\\ {\alpha }_{32}=1-{\alpha }_{31}\end{array}$ (37)

Since ${\ddot{\mathbf{u}}}_{t_s}$ is not participating in the approximations, the amplification matrix of the second family of implicit schemes does not have the spurious root. In addition, the computation of the initial acceleration vector ${\ddot{\mathbf{u}}}_0={\mathbf{M}}^{-1}\mathbf{q}\left({\mathbf{u}}_0,{\dot{\mathbf{u}}}_0,0\right)$ is not required in the second family of implicit schemes. This kind of time integration schemes is called a self-starting scheme. However, most of the existing time integration schemes require the computation of ${\ddot{\mathbf{u}}}_0={\mathbf{M}}^{-1}\mathbf{q}\left({\mathbf{u}}_0,{\dot{\mathbf{u}}}_0,0\right)$ to start the procedure. If the effective coefficient matrices of the first and second sub-steps are constructed identically in linear analyses, the second family of implicit schemes requires only one matrix factorization throughout the entire computation, which is not achievable in the existing schemes including the Newmark method.

This property may be useful for the analysis of big linear systems. The conditions given in Eqs. (26)-(27) are also used in the second family. When the parameters are determined according to Eqs. (26)-(27), (35) and (36)-(37), the new method given in Eqs. (11)-(16) can be simplified as

${}^1{\dot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}={\dot{\mathbf{u}}}_{t_s}+\left({\tau }_1\Delta t\right)\left({}^1{\ddot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}\right)$ (38a)

${}^1{\mathbf{u}}_{t_s+{\tau }_1\Delta t}={\mathbf{u}}_{t_s}+\left({\tau }_1\Delta t\right)\left({}^1{\dot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}\right)$ (38b)

$\mathbf{M}{}^1{\ddot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}=\mathbf{f}\left({}^1{\mathbf{u}}_{t_s+{\tau }_1\Delta t},{}^1{\dot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t},t_s+{\tau }_1\Delta t\right)$ (38c)

${}^2{\dot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}={\dot{\mathbf{u}}}_{t_s}+\left({\tau }_2\Delta t\right)\left({\alpha }_{21}{}^1{\ddot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}+(1-{\alpha }_{21}){}^2{\ddot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}\right)$ (38d)

${}^2{\mathbf{u}}_{t_s+{\tau }_2\Delta t}={\mathbf{u}}_{t_s}+\left({\tau }_2\Delta t\right)\left({\alpha }_{21}{}^1{\dot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}+(1-{\alpha }_{21}){}^2{\dot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}\right)$ (38e)

$\mathbf{M}{}^2{\ddot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}=\mathbf{f}\left({}^2{\mathbf{u}}_{t_s+{\tau }_2\Delta t},{}^2{\dot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t},t_s+{\tau }_2\Delta t\right)$ (38f)

${}^3{\dot{\mathbf{u}}}_{t_s+\Delta t}={\dot{\mathbf{u}}}_{t_s}+\Delta t\left({\alpha }_{31}{}^1{\ddot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}+(1-{\alpha }_{31}){}^2{\ddot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}\right)$ (38g)

$\begin{array}{r}{}^3{\mathbf{u}}_{t_s+\Delta t}={\mathbf{u}}_{t_s}+\Delta t\left({\alpha }_{31}{}^1{\dot{\mathbf{u}}}_{t_s+{\tau }_1\Delta t}+(1-{\alpha }_{31}){}^2{\dot{\mathbf{u}}}_{t_s+{\tau }_2\Delta t}\right)\end{array}$ (38h)

To guarantee second-order accuracy, ${\alpha }_{31}$ can be determined to satisfy the condition given in Eq. (21) with $k=2$. Form this condition ${\alpha }_{31}$ is computed as

$\begin{array}{r}{\alpha }_{31}=\frac{2{\tau }_2-1}{2({\tau }_2-{\tau }_1)}\end{array}$ (39)

For the second family of implicit schemes, the ultimate spectral radius ${\rho }_{\mathrm{\infty }}$ defined in Eq. (24) is computed as

$\begin{array}{r}{\rho }_{\mathrm{\infty }}=\frac{2{\alpha }_{21}{\tau }_1{\tau }_2+2{\tau }_2-2{\tau }_1{\tau }_2-2{\alpha }_{21}{\tau }_2-1+2{\tau }_1}{2{\tau }_1{\tau }_2({\alpha }_{21}-1)}\end{array}$ (40)

Like the case of the first family of implicit schemes, the level of numerical dissipation in the high-frequency limit can be determined by stating one of the remaining algorithmic parameters in terms of ${\rho }_{\mathrm{\infty }}$ according to Eq. (40). From Eq. (40), ${\alpha }_{21}$ can be determined as

$\begin{array}{r}{\alpha }_{21}=\frac{2{\tau }_1{\tau }_2{\rho }_{\mathrm{\infty }}-2{\tau }_1{\tau }_2+2{\tau }_2+2{\tau }_1-1}{2({\tau }_1{\rho }_{\mathrm{\infty }}-{\tau }_1+1){\tau }_2}\end{array}$ (41)

For the dissipation control, ${\rho }_{\mathrm{\infty }}$ can be chosen from zero to one. To obtained an identical effective system matrix for the first and second sub-steps in linear problems, ${\tau }_1$ can be chosen as

$\begin{array}{r}{\tau }_1=-{\tau }_2({\alpha }_{21}-1)\end{array}$ (42)

In this case, ${\alpha }_{21}$ can be determined from Eq. (40) as

$\begin{array}{r}{\alpha }_{21}=\frac{2{\tau }_2{\rho }_{\mathrm{\infty }}-2{\tau }_2+2-\sqrt{2{\rho }_{\mathrm{\infty }}+2}}{2({\rho }_{\mathrm{\infty }}-1){\tau }_2}\end{array}$ (43)

When Eq. (43) is used, only one matrix factorization is required in linear analyses. This is a huge advantage when compared with the existing single and composite implicit schemes.

Like case of the first family of implicit schemes, the second family is summarized in Table 2. The second general form given in Table 2 is newly developed. The cases 2-1 and 2-2 are spectrally identical to the scheme developed by the author [¹²]. The second general form given in Table 2 can also be used to develop various new schemes in this family.

Table 2
Algorithmic parameters of the second family of implicit schemes.

	the second general form	case 2-1	case 2-2
${\tau }_0$	0	0	0
${\tau }_1$	adjustable	adjustable	$\frac{-2+\sqrt{2{\rho }_{\mathrm{\infty }}+2}}{2({\rho }_{\mathrm{\infty }}-1)}$
${\tau }_2$	adjustable	1	1
${\tau }_3$	1	1	1
${\alpha }_{10}$	0	0	0
${\alpha }_{11}$	1	1	1
${\alpha }_{20}$	0	0	0
${\alpha }_{21}$	$1-\frac{-2{\tau }_1+1}{2\left({\tau }_1{\rho }_{\mathrm{\infty }}-{\tau }_1+1\right){\tau }_2}$	$1-\frac{-2{\tau }_1+1}{2\left({\tau }_1{\rho }_{\mathrm{\infty }}-{\tau }_1+1\right)}$	$1-\frac{-2+\sqrt{2{\rho }_{\mathrm{\infty }}+2}}{2({\rho }_{\mathrm{\infty }}-1)}$
${\alpha }_{22}$	$\frac{-2{\tau }_1+1}{2\left({\tau }_1{\rho }_{\mathrm{\infty }}-{\tau }_1+1\right){\tau }_2}$	$\frac{-2{\tau }_1+1}{2\left({\tau }_1{\rho }_{\mathrm{\infty }}-{\tau }_1+1\right)}$	$\frac{-2+\sqrt{2{\rho }_{\mathrm{\infty }}+2}}{2({\rho }_{\mathrm{\infty }}-1)}$
${\alpha }_{30}$	0	0	0
${\alpha }_{31}$	$1-\frac{2{\tau }_1-1}{2({\tau }_1-{\tau }_2)}$	$1-\frac{2{\tau }_1-1}{2({\tau }_1-1)}$	$1-\frac{{\rho }_{\mathrm{\infty }}+1-\sqrt{2{\rho }_{\mathrm{\infty }}+2}}{2{\rho }_{\mathrm{\infty }}-\sqrt{2{\rho }_{\mathrm{\infty }}+2}}$
${\alpha }_{32}$	$\frac{2{\tau }_1-1}{2({\tau }_1-{\tau }_2)}$	$\frac{2{\tau }_1-1}{2({\tau }_1-1)}$	$\frac{{\rho }_{\mathrm{\infty }}+1-\sqrt{2{\rho }_{\mathrm{\infty }}+2}}{2{\rho }_{\mathrm{\infty }}-\sqrt{2{\rho }_{\mathrm{\infty }}+2}}$
${\alpha }_{33}$	0	0	0
${\beta }_{ij}$	${\alpha }_{ij}$

3.3 Novel sets of parameters for the conservation of total energy

In this section, two novel sets of optimized parameters are presented for the implicit schemes to enhance the total energy conserving capability. In fact, two general forms of the first and second families of implicit schemes given in Tables 1 and 2 can also be used to enhance the total energy conserving capability. To find optimal values of the undetermined parameters, nonlinear conservative dynamics systems are used. It is also noted, the newly proposed sets of parameters do not increase any computational costs at all.

The simple nonlinear pendulum problem is frequently used for the test of time schemes [¹⁷], because its exact solution is computable at discrete time points based on the conservation of total energy [¹¹,³⁷]. In current study, on the other hand, we use the problem for the optimization of the extra algorithmic parameters. The governing equation of the simple nonlinear pendulum problem is given by

$\ddot{\theta }\left(t\right)+{\omega }^2\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \left(t\right)\right)=0$ (44)

and the initial conditions are

$\theta \left(0\right)={\theta }_0,\dot{\theta }\left(0\right)={\dot{\theta }}_0$ (45)

where, $\omega =\sqrt{g/L}$, $\theta \left(t\right)$ is the angle between the rod and the vertical line, $g$ is the gravitational constant, and $L$ is the length of the rod. Dimensionless values $g=1.0$ and $L=1.0$ are used in this study. Two sets of the initial conditions are used to synthesize highly nonlinear cases. With arbitrary initial displacement and velocity, the total energy of the pendulum can be stated as

$E_t=\frac{{\left({\dot{\theta }}_0\right)}^2}{2}-\mathrm{c}\mathrm{o}\mathrm{s}\left({\theta }_0\right)$ (46)

To enhance the total energy conserving capability in the implicit schemes, we impose a conditions of

$\left|E_t-E_n\right|=O\left(\Delta t^{k+1}\right)$ (47)

where $k$ is the targeted order of the convergence rate of the total energy error, and $E_n$ is the total energy which is obtained by using two general forms of the first and second families of implicit schemes given in Tables 1 and 2 at time $t=\Delta t$. To satisfy the condition given in Eq. (47), ${\alpha }_{11}$ and ${\tau }_1$ of the general form given in Table 1 can be chosen as

${\alpha }_{11}=-\frac{1}{{\tau }_1\left(3{\tau }_1{\rho }_{\mathrm{\infty }}+3{\tau }_1-2{\rho }_{\mathrm{\infty }}-4\right)}$ (48)

${\tau }_1=\frac{1}{2}$ (49)

With the algorithmic parameters given in Eqs. (48)-(49), the total energy error of the first general form given in Table 1 is computed as

$\left|E_t-E_n\right|=-\frac{\left({\rho }_{\mathrm{\infty }}-1\right)}{36{\left({\rho }_{\mathrm{\infty }}+5\right)}^2}\Lambda {\Delta t}^4+O\left({\Delta t}^5\right)$ (50)

where $\Lambda$ is given by

$\begin{array}{rr}\Lambda =& -2({\rho }_{\mathrm{\infty }}+1){\text{cos}}^3\left({\theta }_0\right)+({\rho }_{\mathrm{\infty }}-3){{\dot{\theta }}_0}^2{\text{cos}}^2\left({\theta }_0\right)\\ & +2({\rho }_{\mathrm{\infty }}+2)\text{cos}\left({\theta }_0\right)+({\rho }_{\mathrm{\infty }}+5){{\dot{\theta }}_0}^2\end{array}$ (51)

From Eq. (50), it can be shown that the choice of ${\rho }_{\mathrm{\infty }}=1$ (the non-dissipative case) gives $\left|E_t-E_n\right|=O\left(\Delta t^5\right)$. Since the condition given in Eq. (50) is satisfied between the exact initial conditions and the solutions at $t=\Delta t$, the convergence rate of the total energy error at arbitrary time $t=t_s+\Delta t$ becomes one order lower than that of Eq. (50). Thus, the convergence rate of the total energy error with the choice of ${\rho }_{\mathrm{\infty }}=1$ becomes fourth-order. Similarly, ${\tau }_1$ and ${\tau }_2$ of the general form of the second family given in Table 2 can be chosen as

${\tau }_1=\frac{3{\tau }_2-2}{3\left(2{\tau }_2-1\right)}$ (52)

${\tau }_2=\frac{\sqrt{3}+3}{6}$ (53)

With the algorithmic parameters given in Eqs. (52)-(53), the total energy error of the second general form given in Table 2 for the case ${\rho }_{\mathrm{\infty }}=1$ is computed as

$\left|E_t-E_n\right|=-\frac{{\dot{\theta }}_0}{4320}\Lambda {\Delta t}^5+O\left({\Delta t}^6\right)$ (54)

where $\Lambda$ is given by

$\begin{array}{rr}\Lambda =& -30(\sqrt{3}-2){\text{cos}}^2\left({\theta }_0\right)+10(\sqrt{3}-1){{\dot{\theta }}_0}^2\text{cos}\left({\theta }_0\right)+{{\dot{\theta }}_0}^4+30(\sqrt{3}-2)\end{array}$ (55)

If the parameters given in Eqs. (48) - (49) and Eqs. (52) - (53) are used for the first and second general forms, respectively, the total energy conserving capability can be enhanced noticeably. For $0\le {\rho }_{\mathrm{\infty }}\le 1$, the non-dissipative and various dissipative cases are obtained. It is noted that the convergence rate of the existing implicit composite schemes is only second-order.

${\theta }_0=0$ is used, and two different values are considered for ${\dot{\theta }}_0$. ${\dot{\theta }}_0=1.999999238456499$ and ${\dot{\theta }}_0=2.000000761543501$ are used to synthesize two highly nonlinear cases [¹⁷]. The minimum total energy required for the pendulum to make a complete rotation about the pivot point is 1.0. With ${\theta }_0=0$ and ${\dot{\theta }}_0=1.999999238456499$, the total energy of the first case becomes about $E=0.999998476913288$, and the pendulum should not make a complete rotation. For this case, the period is about $T=33.7210$.

With ${\theta }_0=0$ and ${\dot{\theta }}_0=2.000000761543501$, the total energy of the second case becomes about $E=1.000001523087292$, and the pendulum passes the peak points slowly. For this case, the period becomes $T=16.8605$. In this particular problem, the ability of keeping the total energy is very important. It is noted that only the non-dissipative case (${\rho }_{\mathrm{\infty }}=1$) is considered to satisfy the condition of $\left|E_t-E_n\right|=O\left(\Delta t^5\right)$.

The numerical results presented in Figs. 2 and 3 (the first and second sets of the initial conditions respectively) are in a good agreement with the enhanced energy conserving capability of the new scheme shown in Fig. 4. As shown in Fig. 2 (a), the solution of the new second-order implicit scheme (red line) is superposing the exact reference solution (black dashed line) while the solution of the KR and NB methods (green line, the case 1-3) shows deviations from the reference solution. The Kim method (yellow line, the case 2-2) also gives very accurate predictions when the smaller time step ($\Delta t=T/2500$) is used, but the KR and NB methods (the case 1-3) cannot. As shown in Fig. 3, similar tendencies are also observed for the second set of the initial conditions.

Figure 2
Angle of oscillating pendulum. (a) Solution of various schemes with $\Delta t=T/500$. (b) Solution of various schemes with $\Delta t=T/5000$. T=33.7210. Non-dissipative cases (${\rho }_{\infty }=1$) are used.

Figure 3
Angle of rotating pendulum. (a) Solution of various schemes with $\Delta t=T/250$. (b) Solution of various schemes with $\Delta t=T/2500$. T=16.8605. Non-dissipative cases (${\rho }_{\infty }=1$) are used.

Figure 4
Convergence rate of total energy measured at T/4. $\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}=\left|E_n-E_t\right|$. (a) Oscillating case. (b) Rotating case.

Two additional nonlinear SDOF problems are also used to show that the enhanced energy conserving capability of the implicit composite time schemes is not limited to the single pendulum problem. The softening and hardening elastic spring problems presented in [³⁸,³⁹] are solved by using various methods. These problems have also been used to test the total energy conserving capability of various existing time schemes. Since no energy dissipation mechanism is considered in the problems, total energy should also be conserved [²²]. The softening problem is described by

$\begin{array}{r}\ddot{u}\left(t\right)+s\text{tanh}\left(u\left(t\right)\right)=0\end{array}$ (56a)

$\begin{array}{r}u\left(0\right)=u_0,\dot{u}\left(0\right)=v_0\end{array}$ (56b)

where $s>0.0$. In this case, $s=100.0$, $u_0=4.0$, and $v_0=0.0$ are also used [³⁸]. The period of the problem is computed as $T=1.141876$[¹¹,²²]. The motion of the hardening elastic spring is stated as

$\begin{array}{r}\ddot{u}\left(t\right)+s_{1}u\left(t\right)\left(1+s_2{u\left(t\right)}^2\right)=0\end{array}$ (57a)

$\begin{array}{r}u\left(0\right)=u_0,\dot{u}\left(0\right)=v_0\end{array}$ (57b)

where $s_1=100.0$, $s_2=10.0$, $u_0=1.5$, and $v_0=0.0$ [³⁸]. The period of the problem is computed as $T=0.151532$. In the additional examples, as shown in Fig. 5, the total energy convergence rate of the proposed scheme is also observed as fourth-order.

Figure 5
Convergence rate of total energy measured at T/4. $\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}=\left|E_n-E_t\right|$. (a) The softening spring $\ddot{u}+s\text{tanh}\left(u\right)=0$, $s=100.0$, $u_0=4.0$, and $v_0=0.0$. (b) The hardening spring $\ddot{u}+s_{1}u\left(1+s_2{u}^2\right)=0$, $s_1=100.0$, $s_2=10.0$, $u_0=1.5$, and $v_0=0.0.$

4 Explicit schemes

Interestingly, several improved explicit composite schemes can also be developed based on the approximations given in Eqs. (9)-(16). Unlike the case of the implicit schemes, explicit schemes do not use the conditions given in Eqs. (25) and (26). For the development of explicit methods, the displacement and velocity vectors are explicitly approximated through

$\begin{array}{r}{\beta }_{11}={\beta }_{22}={\beta }_{33}=0\end{array}$ (58)

With the condition given in Eq. (58), new explicit methods can also be developed by using the same procedures used in the development of the implicit schemes. For example, ${\tau }_0$, ${\tau }_1$, ${\tau }_2$, and ${\tau }_3$ can be chosen as

$\begin{array}{r}{\tau }_0=0,{\tau }_1={\tau }_2={\tau }_3=1\end{array}$ (59)

which is the same as the cases 1-1 and 1-2 as shown in Table 1. The values of ${\tau }_i$ given in Eq. (59) may make handling of the external force vector simpler than other cases. To make one of the three eigenvalues of the amplification matrix $\mathbf{A}$ become zero for the undamped case (i.e., $\xi =0$), ${\alpha }_{30}$, ${\alpha }_{31}$ and ${\alpha }_{33}$ can be determined according to

$\begin{array}{r}{\alpha }_{30}={\alpha }_{20},{\alpha }_{31}={\alpha }_{21},{\alpha }_{33}=0\end{array}$ (60)

Then, to satisfy ${\tau }_e=O\left(\Delta t^3\right)$, ${\alpha }_{20}$, ${\alpha }_{21}$ and ${\beta }_{31}$ should be chosen as

$\begin{array}{r}{\alpha }_{20}={\beta }_{30}=\frac{1}{2},{\alpha }_{21}=\frac{6{\beta }_{20}-5}{12({\beta }_{20}-1)},{\beta }_{31}=\frac{12{\beta }_{20}-7}{12({\beta }_{20}-1)}\end{array}$ (61)

and ${\beta }_{20}$ is used to control the level of the numerical dissipation. In this family of explicit schemes, two of the eigenvalues become complex conjugate, and the point that two eigenvalues change from real to complex is called the bifurcation point. The spectral radius at the bifurcation point is frequently used to specify the level of the numerical dissipation in explicit schemes, because the ultimate spectral radius is not available in explicit methods due to the conditional stability. The spectral radius at the bifurcation point is defined by

$\begin{array}{r}{\rho }_b=\rho \left(\mathbf{A}\left({\Delta t}_b\right)\right)\end{array}$ (62)

where ${\Delta t}_b$ is the time step where the bifurcation of two roots occurs. In explicit schemes, the level of numerical dissipation is often measured by the spectral radius at the bifurcation. By stating one of the algorithmic parameters of a scheme in terms of ${\rho }_b$ according to Eq. (62), the desired level of numerical dissipation at the bifurcation point is prescribed through ${\rho }_b$. The most dissipative case is obtained when ${\rho }_b$ is zero, and less dissipative cases are obtained as ${\rho }_b$ approaches to one. When ${\rho }_b$ is one, the non-dissipative case is obtained. In this particular family of explicit schemes, ${\beta }_{20}$ is determined by

$\begin{array}{r}{\beta }_{20}=\frac{5{{\rho }_b}^2+71{\rho }_b+38-5\sqrt{-3{{\rho }_b}^4+15{{\rho }_b}^2+18{\rho }_b+6}}{48(2{\rho }_b+1)}\end{array}$ (63)

where $0\le {\rho }_b\le 1$. When the case of ${\rho }_b=1$ is used to determine ${\beta }_{20}$ according to Eq. (63), the non-dissipative case is obtained. This case is recommended as a standard setting for general analyses. On the other hand, more dissipative cases are obtained as ${\rho }_b$ changes from one to zero, which is similar to the role of the ultimate spectral radius ${\rho }_{\mathrm{\infty }}$ in the implicit schemes.

In summary, the algorithmic parameters of the improved explicit composite schemes are presented in Table 3. With the algorithmic parameters given in Table 3, the level of numerical dissipation is controlled by specifying ${\rho }_b$. Like cases of the implicit schemes, the choice of ${\rho }_b=1$ gives the non-dissipative cases, and more dissipative cases are obtained as ${\rho }_b$ approaches to zero. Here, ${\rho }_b$ is the value of the spectral radius at the bifurcation point. $\Delta t_c$ denotes the critical time step. The critical time step is defined as the largest time step that satisfies the condition $\rho \left(\mathbf{A}\right)\le 1.0$.

In Table 3, three different families of explicit schemes with dissipation control capabilities are presented. The case 3-1 is spectrally identical to the recent scheme developed by the author and Reddy [²³], the case 3-2 is spectrally identical to the scheme presented in Ref. [¹³]. The case 3-3 is spectrally identical to the scheme developed by the author [¹⁸]. The explicit scheme proposed by Soares [⁴⁰] also uses two stages like the explicit composite schemes, and its spectral properties are almost identical to the case 3-2 for the undamped case (i.e., $\xi =0$). However, it becomes only a conditionally stable implicit scheme for damped cases (i.e., $\xi \ne 0$) and velocity dependent problems. Besides, the Soares method is not extended to general nonlinear problems.

Table 3
Algorithmic parameters of the explicit schemes with adjustable numerical dissipation.

	case 3-1	case 3-2	case 3-3
${\tau }_0$	0	0	0
${\tau }_1$	1	$\frac{-2+\sqrt{2+2{\rho }_b}}{{\rho }_b-1}$	$\frac{1}{2}$
${\tau }_2$	1	1	$\frac{1}{2}$
${\tau }_3$	1	1	1
${\alpha }_{10}$	$\frac{1}{2}$	$\frac{1}{2}$	0
${\alpha }_{11}$	$\frac{1}{2}$	$\frac{1}{2}$	1
${\alpha }_{20}$	$\frac{1}{2}$	$-\frac{1}{2({\tau }_1-2)}$	0
${\alpha }_{21}$	$\frac{6{\beta }_{20}-5}{12({\beta }_{20}-1)}$	$-\frac{1}{2({\tau }_1-2)}$	$\frac{2}{3}$
${\alpha }_{22}$	$-\frac{1}{12({\beta }_{20}-1)}$	$\frac{{\tau }_1-1}{{\tau }_1-2}$	$\frac{1}{3}$
${\alpha }_{30}$	$\frac{1}{2}$	$-\frac{1}{2({\tau }_1-2)}$	0
${\alpha }_{31}$	$\frac{6{\beta }_{20}-5}{12({\beta }_{20}-1)}$	$-\frac{1}{2({\tau }_1-2)}$	$\frac{1-{\alpha }_{32}}{2}$
${\alpha }_{32}$	$-\frac{1}{12({\beta }_{20}-1)}$	$\frac{{\tau }_1-1}{{\tau }_1-2}$	$\frac{2({\rho }_b+1)({{\rho }_b}^2-2{\rho }_b-5+2\sqrt{-3{{\rho }_b}^2+6{\rho }_b+6})}{({\rho }_b-1{)}^3}$
${\alpha }_{33}$	0	0	$\frac{1-{\alpha }_{32}}{2}$
${\beta }_{10}$	1	1	0
${\beta }_{11}$	0	0	0
${\beta }_{20}$	$\frac{5{{\rho }_b}^2+71{\rho }_b+38-5\sqrt{-3{{\rho }_b}^4+15{{\rho }_b}^2+18{\rho }_b+6}}{48(2{\rho }_b+1)}$	$\frac{{\tau }_1}{2}$	0
${\beta }_{21}$	$1-{\beta }_{20}$	$\frac{2-{\tau }_1}{2}$	1
${\beta }_{22}$	0	0	0
${\beta }_{30}$	$\frac{1}{2}$	$-\frac{{{\tau }_1}^2-3{\tau }_1+1}{2{\tau }_1}$	0
${\beta }_{31}$	$\frac{12{\beta }_{20}-7}{12({\beta }_{20}-1)}$	$-\frac{{\tau }_1-1}{2{\tau }_1}$	0
${\beta }_{32}$	$-\frac{6{\beta }_{20}-1}{12({\beta }_{20}-1)}$	$\frac{{\tau }_1}{2}$	1
${\beta }_{33}$	0	0	0

In these schemes, the level of numerical dissipation at the bifurcation point is determined by specifying the value of ${\rho }_b$. Unlike the cases 3-1 and 3-2, the case 3-3 does not require the computation of the initial acceleration vector to start the scheme. This property is similar to the cases of 2-1, 2-2, and 2-3. However, the advantage of the self-starting capability of the case 3-3 is not that huge when compared with the cases of the implicit schemes, because the explicit schemes inevitably require the result of the factorization of the mass matrix in all steps. Thus, omitting the computation of ${\ddot{\mathbf{u}}}_0={\mathbf{M}}^{-1}\mathbf{q}\left({\mathbf{u}}_0,{\dot{\mathbf{u}}}_0,0\right)$ does not give a noticeable computational advantage. However, the case 3-3 shows improved accuracy when applied to velocity-independent linear and nonlinear problems. Among the three explicit families, the case 3-1 is the most accurate. The key values of ${\rho }_b$, $\Delta t_b/T$ and $\Delta t_c/T$ of the cases 3-1, 3-2, and 3-3 are presented in Tables 4 and 5, where $T$ is the shortest natural period of a given dynamic system.

Table 4
${\beta }_{20}$(only for the case 3-1), ${\alpha }_{32}$(only for the case 3-3), $\Delta t_b$ and $\Delta t_c$ of the cases 3-1 and 3-3 for varying values of ${\rho }_b$, where $T=2\pi$ is the period.

${\mathit{\rho }}_{\mathit{b}}$	${\mathit{\beta }}_{20}$ (for the case 3-1 only)	${\mathit{\alpha }}_{32}$ (for the case 3-3 only)	$\mathit{\Delta }{\mathit{t}}_{\mathit{b}}/\mathit{T}$	$\mathit{\Delta }{\mathit{t}}_{\mathit{c}}/\mathit{T}$
0.0	0.536511	0.202041	0.525428	0.551329
0.1	0.539104	0.191927	0.530881	0.551329
0.2	0.542188	0.179747	0.535519	0.551329
0.3	0.545748	0.165482	0.539443	0.551329
0.4	0.549771	0.149088	0.542728	0.551329
0.5	0.554249	0.130500	0.545429	0.551329
0.6	0.559174	0.109610	0.547590	0.551329
0.7	0.564544	0.086299	0.549242	0.551329
0.8	0.570357	0.060403	0.550406	0.551329
0.9	0.576619	0.031720	0.551099	0.551329
1.0	0.583333	0.000000	0.551329	0.551329

Table 5
${\tau }_1$, $\Delta t_b$ and $\Delta t_c$ of the explicit scheme of the case 3-2 for varying values of ${\rho }_b$, where $T=2\pi$ is the period.

${\mathit{\rho }}_{\mathit{b}}$	${\mathit{\tau }}_1$ (for the case 3-2 only)	$\mathit{\Delta }{\mathit{t}}_{\mathit{b}}/\mathit{T}$	$\mathit{\Delta }{\mathit{t}}_{\mathit{c}}/\mathit{T}$
0.0	0.585786	0.543390	0.568311
0.1	0.574178	0.554375	0.573862
0.2	0.563508	0.564873	0.579795
0.3	0.553641	0.574938	0.586050
0.4	0.544467	0.584627	0.592585
0.5	0.535898	0.593976	0.599381
0.6	0.527864	0.603016	0.606406
0.7	0.520304	0.611776	0.613651
0.8	0.513167	0.620284	0.621106
0.9	0.506411	0.628560	0.628764
1.0	0.500000	0.636620	0.636620

When the dissipation control capability is unnecessary, all of the algorithmic parameters can be used to enhance the accuracy. The schemes presented in Table 6 are developed to maximize the accuracy. As shown in Table 6, these schemes do not have adjustable free parameters. Since they are moderately dissipative, they are not suitable for long-term analyses where the total energy of a given system should be conserved. To use these explicit schemes for long-term analyses, the time step should be small enough not decease the total energy noticeably. When $\Delta t\le T/100$ is guaranteed, for example, the numerical dissipation of these explicit schemes is almost negligible. For long-term analyses, the cases 3-1 and 3-3 with ${\rho }_b=1.0$ are recommended. However, the case 4-2 and 4-3 are fourth-order accurate for velocity independent problems and third-order accurate for velocity dependent problems, and more accurate solutions can be obtained with large time steps. The case 4-1 is spectrally equivalent to the scheme developed in Ref. [¹⁷], and the cases 4-2 and 4-3 are equivalent to the schemes developed in Ref. [²³].

The cases 3-1, 3-2, and 3-3 can change the level of numerical dissipation by specifying the value of the spectral radius at the bifurcation point. As shown in Fig. 6 (a), the cases 3-1 and 3-2 can specify the level of numerical dissipation through ${\rho }_b$. Interestingly, the case 3-1 has very small period error when compared with the case 3-2 as shown in Fig. 6 (b), but the case 3-2 has slightly higher stability limit as shown in Fig. 6 (a). The spectral radius and the period error of the higher-order accurate schemes are presented in Fig. 7. As shown in Fig. 7, the period error of the case 4-2 and 4-3 is very small. Unlike the cases 3-1, 3-2, and 3-3, the case 4-2 and 4-3 is third-order accurate for both linear and nonlinear problems.

Table 6
Algorithmic parameters of the explicit schemes without enhanced accuracy.

	case 4-1	case 4-2	case 4-3
${\tau }_0$	0	0	0
${\tau }_1$	$\frac{1}{3}$	$\frac{1}{3}$	$\frac{1}{2}$
${\tau }_2$	1	1	1
${\tau }_3$	1	1	1
${\alpha }_{10}$	$\frac{1}{2}$	$\frac{1}{2}$	$\frac{1}{2}$
${\alpha }_{11}$	$\frac{1}{2}$	$\frac{1}{2}$	$\frac{1}{2}$
${\alpha }_{20}$	$\frac{1}{6}$	0	$\frac{2}{9}$
${\alpha }_{21}$	$\frac{1}{2}$	$\frac{3}{4}$	$\frac{5}{9}$
${\alpha }_{22}$	$\frac{1}{3}$	$\frac{1}{4}$	$\frac{2}{9}$
${\alpha }_{30}$	$\frac{1}{6}$	$\frac{1}{12}$	$\frac{1}{3}$
${\alpha }_{31}$	$\frac{1}{2}$	$\frac{5}{8}$	$\frac{1}{3}$
${\alpha }_{32}$	$\frac{1}{6}$	$\frac{1}{8}$	0
${\alpha }_{33}$	$\frac{1}{6}$	$\frac{1}{6}$	$\frac{1}{3}$
${\beta }_{10}$	1	1	1
${\beta }_{11}$	0	0	0
${\beta }_{20}$	$-\frac{1}{2}$	-1	$-\frac{1}{2}$
${\beta }_{21}$	$\frac{3}{2}$	2	$\frac{3}{2}$
${\beta }_{22}$	0	0	0
${\beta }_{30}$	0	0	$\frac{1}{6}$
${\beta }_{31}$	$\frac{3}{4}$	$\frac{3}{4}$	$\frac{2}{3}$
${\beta }_{32}$	$\frac{1}{4}$	$\frac{1}{4}$	$\frac{1}{6}$
${\beta }_{33}$	0	0	0

Figure 6
(a) Spectral radius of the cases 3-1 and 3-2 for varying values of $\Delta t/T$. (b) Period error of the cases 3-1 and 3-2 for varying values of $\Delta t/T$.

Figure 7
(a) Spectral radius of the cases 4-2 and 4-2 for varying values of $\Delta t/T$. (b) Period error of the cases 4-2 and 4-3 for varying values of $\Delta t/T$.

5 Discussion

In this section, various benchmark problems are solved by using time schemes and their results are compared. Through numerical tests, general advantages of using composite schemes, the high-frequency filtering capability, the use the splitting ratio, and the computational cost are discussed in detail.

5.1 General advantages of composite schemes

To demonstrate general advantages of using the composite schemes, $u_0={\dot{u}}_0=0.0$, $\xi =0.05$, $\omega =2\pi$, $q\left(t\right)=10.0\mathrm{s}\mathrm{i}\mathrm{n}\left(t\right)+5.0\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{t}{3}\right)$ are used for the single-degree-of-freedom problem given in Eqs. (17)- (18), and the problem is numerically solved by using the composite schemes and the implicit and explicit generalized-$\alpha$ methods [5,26]. Then, numerical results are compared with each other. To guarantee enough accuracy in each scheme, $\Delta t=0.01$ is used for all composite schemes, and $\Delta t=0.005$ (i.e., a half time step) is used for the implicit and explicit generalized-$\alpha$ methods.

In general, the computation cost per time step required in implicit composite time schemes is twice higher than that of the implicit generalized-$\alpha$ method. As observed in the results presented in Fig. 8 (a), the maximum accuracy level of the non-dissipative cases of implicit composite schemes (${\rho }_{\mathrm{\infty }}=1.0$) is almost the same as the cases of using the trapezoidal rule and the generalized-$\alpha$ method (${\rho }_{\mathrm{\infty }}=1.0$) twice with a half time step. On the other hand, the asymptotic annihilating cases (${\rho }_{\mathrm{\infty }}=0.0$) of the implicit composite schemes are presenting much smaller errors than the asymptotic annihilating cases (${\rho }_{\mathrm{\infty }}=0.0$) of the generalized-$\alpha$ method as shown in Fig. 8 (b).

The explicit composite schemes are also tested by using the same problem. As shown in Fig. 9, the non-dissipative and the most dissipative cases of the explicit composite schemes are presenting much smaller errors. The most dissipative case of the explicit generalized-$\alpha$ method gives less accurate solutions when compared with its non-dissipative case. However, the most dissipative cases of the explicit composite schemes give more accurate solutions when compared with their non-dissipative cases. As shown in Fig. 6 (b), the period error of the explicit composite methods decreases as ${\rho }_{\mathrm{\infty }}$ approaches to zero. For this particular problem, the case 4-2 is presenting almost negligible errors as shown in Fig. 10.

As shown in this test, the composite schemes can give more accurate solutions than using the improved single-stage schemes twice with a half step. Thus, more efficient computation is possible by using the composite schemes with large time steps. Interestingly, all of the explicit composite schemes are presenting much smaller errors than the implicit schemes. Thus, using the explicit composite schemes for the analysis of general dynamic problems may be not only more efficient but also more accurate.

Figure 8
Errors of implicit schemes. (a) Non-dissipative cases (${\rho }_{\mathrm{\infty }}=1.0$). (b) The most dissipative case (${\rho }_{\mathrm{\infty }}=0.0$). The implicit generalized-$\alpha$ method used a half time step to make its computational effort equivalent to the implicit composite schemes.

Figure 9
Errors of explicit schemes. (a) Non-dissipative cases (${\rho }_{\mathrm{\infty }}=1.0$). (b) The most dissipative case (${\rho }_{\mathrm{\infty }}=0.0$). The explicit generalized-$\alpha$ method used a half time step to make its computational effort equivalent to the implicit composite schemes.

Figure 10
Errors of explicit schemes without the explicit generalized-$\alpha$ method. (a) Non-dissipative cases (${\rho }_{\mathrm{\infty }}=1.0$). (b) The most dissipative case (${\rho }_{\mathrm{\infty }}=0.0$).

5.2 High-frequency filtering capabilities

A poor discretization of spatial domains of structural problems often introduces the spurious high-frequency modes. Due to the presence of the spurious high-frequency modes, numerical solutions of transient analyses may be seriously contaminated. However, refining the spatial discretization may dramatically increase the computational cost in some cases. To fix this problems more practically, the combination of lower-order elements and the algorithmic dissipation of time integration schemes is frequently used.

Figure 11
Description of the spring-mass problem [¹⁴, ⁴¹].

In the literature [¹², ¹⁴, ⁴⁰, ⁴¹], the high-frequency filtering capability of dissipative time integration schemes was often tested by using the spring-mass problem. In the spring-mass problem presented in Fig. 11, one of the two springs is very stiff when compared with the other one. In Ref. [⁴¹], the spring constant of the stiff spring is chosen as $k_2={10.0}^p$, where $p=7$. $m_1=0.0$, $m_2=1.0$, $m_3=1.0$, ${\omega }_p=6/5$ and zero initial conditions are used. The spring constant of the soft spring is $k_3=1.0$. With $p=7$, two natural frequencies of the system are given by ${\omega }_1=3162.27781828$ and ${\omega }_2=0.999999949999$. In this particular problems, the highest natural frequency of the system is more than three thousand times greater than the lowest natural frequency. The asymptotic annihilating case (${\rho }_{\mathrm{\infty }}=0$) of the implicit composite schemes can effectively filter out the high frequency of the system as presented in Refs. [¹², ⁴¹].

In this study, on the other hand, a different case is considered. To make the difference of the highest and lowest frequencies small, the value of $p$ is chosen as $2$. In the explicit case 3-1 with ${\rho }_b=0.0$, $\Delta t=0.328481363056$ is used to introduce the maximum numerical dissipation into the high frequency mode. In the implicit case 1-3 and the implicit generalized-$\alpha$ method with ${\rho }_{\mathrm{\infty }}=0.0$, three different values are used for $\Delta t$ to verify that the implicit composite scheme is not suitable for this particular case [⁴²].

With $p=2$, two natural frequencies of the system are computed as ${\omega }_1=10.0503730777$ and ${\omega }_2=0.994987939525$. With ${\omega }_1$ and ${\omega }_2$, two natural periods are computed as $T_1=2\pi /{\omega }_1=0.625169360244$ and $T_2=2\pi /{\omega }_2=6.31483564532$. The bifurcation point of the explicit scheme of the case 3-1 with ${\rho }_b=0.0$ is computed as $\frac{\Delta t_b}{T}=0.52542780235992$, where $T$ is the natural period. When $\Delta t$ is chosen as $\frac{\Delta t_b}{T}\times T_1$, the maximum amount of numerical damping is introduced into the mode associated with ${\omega }_1$.

With $\Delta t=0.52542780235992\times 0.625169360244=0.328481363056$, the explicit scheme can filter out the high-frequency mode of this problem with a fast rate. The spectral radius is a direct measurement of the amount of numerical dissipation for a given value of $\Delta t$. As shown in Fig. 12, the case 3-1 (explicit) can introduce the maximum numerical dissipation into the mode associated with ${\omega }_1$, while providing enough accuracy for the mode associated with ${\omega }_2$. On the other hand, the case 1-3 (implicit) can introduce only a moderate amount of numerical dissipation in the high-frequency mode. This explains the reason that the explicit case 3-1 can give better high-frequency filtering as shown in Fig. 13 when $\Delta t=0.328481363056$ is used in both cases. Like the explicit case 3-1, the explicit generalized-$\alpha$ method also gives good high-frequency filtering capability when $\Delta t=0.140712$ is used.

Figure 12
Illustration of levels of spectral radii for $T_1$ and $T_2$ when $\Delta t$ is chosen as $\frac{\Delta t_b}{T}\times T_1$.

Figure 13
Acceleration of the second mass. (a) The cases 3-1 (explicit) with ${\rho }_b=0.0$ and $\Delta t=0.328481$. (b) The cases 1-4(implicit) with ${\rho }_{\mathrm{\infty }}=0.0$ and $\Delta t=0.328481$.

To improve the results of the implicit case 1-3 and the implicit generalized-$\alpha$ method for this case, larger time steps are employed. When $\Delta t$ is chosen as $0.656962726112$, the result is slightly improved as shown in Fig. 14 (a), but the quality of the numerical solution is not that accurate when compared with the result of the case 3-1 presented in Fig. 13. When $\Delta t$ is chosen as $1.313925452224$ to introduce even strong numerical dissipation, the overall accuracy of the numerical solution decreases significantly as shown in Fig. 14 (b). Increasing the size of the time step to introduce strong numerical damping into high-frequency mode also accompanies decease of the accuracy for the low-frequency mode.

Unlike the test conducted with the spring-mass problem, the implicit schemes may become more effective than the explicit schemes for the problems that have multiple spurious high-frequency modes [¹²]. The excited elastic bar problem given in Fig. 15 is used to simulate a situation that multiple spurious high-frequency modes are influencing the numerical solution. The governing equation of the excited elastic bar problem given in Fig. 15 is expressed as

$\rho A\frac{{\partial }^{2}u}{\partial t^2}-\frac{\partial }{\partial x}\left(EA\frac{\partial u}{\partial x}\right)=f\left(x,t\right),0\le x\le L,t\ge 0$ (64)

and the initial and boundary conditions are given by

$u(0,t)=0,{\left.EA\frac{\partial u}{\partial x}\right|}_{x=L}=Q_r=\mathrm{10,000},u(x,0)=0,\dot{u}(x,0)=0$ (65)

For the spatial discretization of the problem given in Fig. 15, one thousand uniform linear elements are used. After the spatial discretization, $\mathbf{M}$, $\mathbf{C}$, $\mathbf{K}$, and $\mathbf{q}$ of this problems are given by

$\begin{array}{r}\mathbf{K}=\frac{EA}{h}{\left[\begin{array}{ccccc}2& -1& & & \\ -1& 2& -1& & \\ & \ddots & \ddots & \ddots & \\ & & -1& 2& -1\\ & & & -1& 1\end{array}\right]}_{N\times N},\\ \mathbf{M}=\frac{\rho Ah}{6}{\left[\begin{array}{ccccc}4& 1& & & \\ 1& 4& 1& & \\ & \ddots & \ddots & \ddots & \\ & & 1& 4& 1\\ & & & 1& 2\end{array}\right]}_{N\times N},\mathbf{q}\left(t\right)={\left\{\begin{array}{c}\\ \\ \\ \\ 10,000\end{array}\right\}}_N\end{array}$(66)

where $E=30.0\times {10.0}^6$, $A=1.0$, $\rho =0.00073$, $L=200.0$, $h=L/N$, and $N=\mathrm{1,000}$.

Figure 14
Acceleration of the second mass. (a) The cases 1-4 (implicit) with ${\rho }_{\mathrm{\infty }}=0.0$ and $\Delta t=0.328481\times 2$. (b) The cases 1-4(implicit) with ${\rho }_{\mathrm{\infty }}=0.0$ and $\Delta t=0.328481\times 4$.

Figure 15
Description of the excited elastic bar problem [²³,³¹].

As shown in Fig. 16, the numerical solution of the central difference method shows severe spurious oscillations when compared with the reference solution. The central difference method cannot eliminate spurious oscillations because it is non-dissipative. On the other hand, two of the asymptotic annihilating cases are presenting the numerical solutions with less spurious oscillations. Among the two asymptotic annihilating cases, the numerical solution of the implicit case 1-1 is more accurate when compared with the explicit case 3-1. This can be explained by the fact that the asymptotic annihilating case of the implicit scheme has the combination of the unconditional stability and the numerical dissipation distributed over a wide frequency range. Simply, the implicit case can eliminated spurious modes distributed over a wide frequency range more efficiently. In this kind of situation, dissipative implicit schemes are more advantageous than dissipative explicit schemes in contrast to the previous example.

Figure 16
(a) Numerical solutions of the excited elastic bar problem obtained from implicit and explicit methods. (b) Enlargement of A-A.

5.3 Splitting ratio and alternative parameters

As shown in Table 1, ${\tau }_1$ of the case 1-1 and 1-2 is always one. The setting of ${\tau }_1=1$ is not allowed in the existing composite time integration schemes, but the proposed unified set of time approximations admits this special setting. The choice of ${\tau }_1=1$ gives a computational advantage when compared with the existing implicit composite methods. With the setting of the algorithmic parameters given in the case 1-1 and 1-2, handling the external force vector becomes as simple as the Newmark method or the central difference method. This can be advantageous for problems where the external force vector is recorded discontinuously with uniform time intervals.

The proposed unified set of time approximations can also include most of the existing composite schemes [³⁴,³⁵] with the setting of algorithmic parameters given in the case 1-3. In the existing implicit composite methods, on the other hand, ${\tau }_1$ is adjusted to synthesize preferable numerical characteristics, such as the level of numerical dissipation in the low-frequency range [³⁴,³⁵]. Then, the time point of the first dynamic equilibrium given in Eq. (11) also changes accordingly.

In Ref. [³⁴], it has been stated that the choice of ${\tau }_1=1.99$ and ${\rho }_{\mathrm{\infty }}=0.0$ can improve the quality of the numerical solutions of the excited elastic bar problem by introducing additional numerical dissipation into the low-frequency range when compared with the standard case obtained by ${\tau }_1=0.5$ and ${\rho }_{\mathrm{\infty }}=0.0$. However, the choice of ${\tau }_1=1.99$ can cause a serious problems in some situations. Basically, the choice of ${\tau }_1=1.99$ worsens the accuracy. To remedy this, very small time step (due to the use of CFL=0.1) should be employed [³⁴]. When smaller time steps are used, the computational time also increases. To support this discussion, two of the implicit cases are presented together by using different CFL numbers as shown in Figs. 17 and 18. Among the three cases presented in Fig. 18, the case 1-1 with the CFL=1.0 gives the most accurate results with much faster computation time.

To verify that the case 1-1 with ${\alpha }_{11}=0.995$ and ${\rho }_{\mathrm{\infty }}=0.0$ can give the same high-frequency filtering capability of the case 1-2 with ${\tau }_1=1.99$ and ${\rho }_{\mathrm{\infty }}=0.0$, the excited elastic bar problem given in Fig. 15 is used again. As shown in Fig 17, the numerical solutions of the cases 1-1 and 1-3 are identical. It can be said that using the case 1-1 with a proper value of ${\alpha }_{11}$ is safer than varying ${\tau }_1$ beyond the range of the current time interval in the case 1-3.

Figure 17
(a) Numerical solutions of the excited elastic bar problem obtained from implicit schemes with CFL=0.1. (b) Enlargement of A-A. Algorithmic parameters of the cases 1-1 and 1-3 are adjusted to introduce additions numerical damping, and smaller time steps are used accordingly.

Figure 18
(a) Numerical solutions of the excited elastic bar problem obtained from various schemes. (b) Enlargement of A-A. Algorithmic parameters are adjusted to introduce additions numerical damping, and smaller time steps are used accordingly.

To show a possible misuse of the splitting ratio in the existing implicit composite schemes more clearly, $u_0={\dot{u}}_0=0.0$, $\xi =0.05$, and $\omega =2\pi$ are used for the single-degree-of-freedom problem given in Eqs. (17)-(18). Then, $q\left(t\right)=10.0$ is considered for $0.0\le t\le 1.0$, and $q\left(t\right)=0.0$ is considered for $t>1.0$ as shown in Fig. 19. For the case 1-3, $t_1=1.01$ and ${\rho }_{\infty }=0.0$ are used. In the case 1-1, however, ${\alpha }_{11}=0.505$ and ${\rho }_{\infty }=0.0$ are used to introduce equivalent amounts of numerical damping in the low-frequency range when compares with the case 1-3 with $t_1=1.01$ and ${\rho }_{\infty }=0.0$. To guarantee enough accuracy,$\Delta t=0.05$ is used in both cases.

In this particular numerical experiment, the numerical solutions of the existing implicit scheme (the case 1-3) shows noticeable errors at $t=1.0$ as shown in Fig. 20. To compute the variables at $t=1.0$, two dynamic equilibrium equations at $t=1.0005$ and $t=1.0$ are used in the first and second sub-steps, respectively. The external forces of the first and second sub-steps are computed as $q\left(1.0005\right)=0.0$ and $q\left(1.0\right)=1.0$, respectively. As shown in Fig. 19, the first sub-step reflects the effect of the external force beyond the range of the current time step (i.e, $q\left(t\right)=0.0$ for $t>1.0$), which is not desirable in this particular case. On the other hand, the numerical solutions of the case 1-1 do not present abruptly increased errors at $t=1.0$, because $q\left(1.0\right)=1.0$ is used for the first and second sub-steps.

Figure 19
Description of external loading.

Figure 20
Comparison of numerical solutions obtained of the single-degree-of-freedom problem with the discontinuous external force. $q\left(t\right)=10.0$ for $0.0\le t\le 1.0$ and $q\left(t\right)=0.0$ for $t>1.0$. (a) The velocity. (b) The acceleration.

5.4 Computational cost

Among the implicit schemes, the second general form, the case 2-1, and the case 2-2 do not require the calculation ${\ddot{\mathbf{u}}}_0={\mathbf{M}}^{-1}\mathbf{f}\left({\mathbf{u}}_0,{\dot{\mathbf{u}}}_0,0\right)$. For this reason, the computational time of linear analyses can be reduced in these schemes when compare with the other implicit schemes. Especially, the case 2-2 requires only one matrix factorization because the effective system matrices of the first and second sub-steps are identical in linear analyses. Even the Newmark family requires the factorization of two system matrices (the mass matrix and effective system matrix) in general linear problems. For nonlinear analyses, this particular advantage of the second general form, the case 2-1, and the case 2-2 becomes negligible, because the factorization of the mass matrix becomes a very small portion of the entire computation in large and complex nonlinear problems.

Unlike the implicit composite schemes, the explicit composite schemes are not only more efficient but also more accurate when the same size of the time step is used within the stability limit. This point is already illustrated through the results of the single-degree-of-freedom problem given in Figs. 8 and 10. For this reason, some nonlinear analyses can be conducted more efficiently with explicit schemes. In addition to the superior accuracy, explicit schemes do not require iterative nonlinear solution finding procedures. The unknown acceleration vector is directly computable as

${}^i{\ddot{\mathbf{u}}}_{t_s+{\tau }_i\Delta t}={\mathbf{M}}^{-1}\mathbf{f}\left({}^i{\mathbf{u}}_{t_s+{\tau }_i\Delta t},{}^i{\dot{\mathbf{u}}}_{t_s+{\tau }_i\Delta t},t_s+{\tau }_i\Delta t\right)$ (67)

because ${}^i{\mathbf{u}}_{t_s+{\tau }_i\Delta t}$ and ${}^i{\dot{\mathbf{u}}}_{t_s+{\tau }_i\Delta t}$ are completely explicit. Probably, this feature is the most significant advantage of explicit schemes when compared with implicit schemes. Unlike explicit schemes, the unknown vectors of nonlinear problems are not directly computable in implicit schemes. In each time step, implicit schemes require several times of iterative solutions finding procedures such as the Picard and Newton-Raphson iterative methods [³²,⁴³], and each iteration requires constructions and factorizations of new effective system matrices, which may increase overall computational time dramatically for nonlinear problems of big sizes. If analyses should be conducted over a long duration of time, using implicit schemes for nonlinear problems is not a good strategy in general.

On the other hand, using implicit schemes may also accomplish a required computational goal faster with less computational effort in some situations. For example, many transient analyses of structural members require only the dominant low-frequency solutions for a short duration of time. In this case, dissipative implicit schemes with a large time step can finish the analyses by numerically eliminating the high-frequency modes. Nonlinear transient analyses of thin structural members are good examples [³¹,³⁵]. On the other hand, explicit schemes require a very small time step that satisfies the stability condition normally dictated by the unimportant highest frequency mode. Thus, it is very difficult to conclude that a certain type of time scheme is always more efficient than others. More details regarding the computational time and cost of the composite time integration schemes can be found in Refs. [¹²,¹⁸,²³].

6 Concluding remarks

The unified set of time approximations presented in this article can include many of the existing composite schemes, and some novel implicit composite schemes can be newly developed base on it. The two general forms of the implicit families proposed in this article may be used for the development of different composite schemes. For example, two novel sets of optimized parameters are presented for the implicit schemes to enhance the total energy conserving capability. To demonstrate general advantages of the composite time integration schemes, numerical performances of the composite schemes are compared with those of the implicit and explicit generalized-$\alpha$ methods.

Optimized algorithmic parameters for various implicit and explicit schemes are summarized and presented. By using these parameters, an optimal scheme can be selected by simply changing values of the algorithmic parameters. In numerical examples, it was shown that some schemes might become more effective than others depending on characteristics of problems and goals of analyses. For this reason, being able to choose different composite schemes from the unified set of time approximations can be more advantageous than working with a specialized time scheme. It is expected that the proposed unified set of time approximations and optimized algorithmic parameters would give more precise options when dealing with different types of challenging problems in structural dynamics.

Acknowledgments

The author deeply appreciates the support from the Republic of Korea Army. The second author acknowledges partial support from the NSF Grant, award number (FAIN): 1952873.

References

• ^[1]
  J. C. Houbolt. A recurrence matrix solution for the dynamic response of aircraft in gusts. NACA Report No. 1010, 1951.
• ^[2]
  N. M. Newmark. A method of computation for structural dynamics. Journal of the Engineering Mechanics Division, 85(3):67-94, 1959.
• ^[3]
  K. C. Park. An improved stiffly stable method for direct integration of nonlinear structural dynamic equations. Journal of Applied Mechanics, 42(2):464-470, 1975.
• ^[4]
  H. M. Hilber, T. J. R. Hughes, and R. L. Taylor. Improved numerical dissipation for time integration algorithms in structural dynamics. Earthquake Engineering & Structural Dynamics, 5(3):283-292, 1977.
• ^[5]
  J. Chung and G. M. Hulbert. A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-alpha method. Journal of Applied Mechanics, 60:271-275, 1993.
• ^[6]
  T. C. Fung. Unconditionally stable higher-order accurate Hermitian time finite elements. International Journal for Numerical Methods in Engineering, 39(20):3475-3495, 1996.
• ^[7]
  T. C. Fung and S. K. Chow. Solving non-linear problems by complex time step methods. Communications in numerical methods in engineering, 18(4):287-303, 2002.
• ^[8]
  K. J. Bathe and M. M. I. Baig. On a composite implicit time integration procedure for nonlinear dynamics. Computers & Structures, 83(31):2513-2524, 2005.
• ^[9]
  W. Kim, S. Park, and J. N. Reddy. A cross weighted-residual time integration scheme for structural dynamics. International Journal of Structural Stability and Dynamics, 14(06):1450023, 2014.
• ^[10]
  Delfim Soares. An implicit family of time marching procedures with adaptive dissipation control. Applied Mathematical Modelling, 40(4):3325-3341, 2016.
• ^[11]
  W. Kim and J. N. Reddy. A new family of higher-order time integration algorithms for the analysis of structural dynamics. Journal of Applied Mechanics, 84(7):071008, 2017.
• ^[12]
  W. Kim. An improved implicit method with dissipation control capability: The simple generalized composite time integration algorithm. Applied Mathematical Modelling, 2020.
• ^[13]
  G. Noh and K. J. Bathe. An explicit time integration scheme for the analysis of wave propagations. Computers & structures, 129:178-193, 2013.
• ^[14]
  W. Kim and J. N. Reddy. An improved time integration algorithm: A collocation time finite element approach. International Journal of Structural Stability and Dynamics, 17(02):1750024, 2017.
• ^[15]
  W. Kim and S. Y. Choi. An improved implicit time integration algorithm: The generalized composite time integration algorithm. Computers & Structures, 196:341-354, 2018.
• ^[16]
  W. Kim and J. H. Lee. An improved explicit time integration method for linear and nonlinear structural dynamics. Computers & Structures, 206:42-53, 2018.
• ^[17]
  W. Kim. An accurate two stage explicit time integration scheme for structural dynamics and various dynamic problems. International Journal for Numerical Methods in Engineering, 120:1-28, 2019.
• ^[18]
  W. Kim. A new family of two-stage explicit time integration methods with dissipation control capability for structural dynamics. Engineering Structures, 195:358-372, 2019.
• ^[19]
  J. Zhang. A-stable two-step time integration methods with controllable numerical dissipation for structural dynamics. International Journal for Numerical Methods in Engineering, 121(1):54-92, 2020.
• ^[20]
  Wooram Kim and J. N. Reddy. Effective higher-order time integration algorithms for the analysis of linear structural dynamics. Journal of Applied Mechanics, 84(7):071009, 2017.
• ^[21]
  W. Kim. Higher-order explicit time integration methods for numerical analyses of structural dynamics. Latin American Journal of Solids and Structures, 16:2-29, 2019.
• ^[22]
  W. Kim and Jin Ho Lee. A comparative study of two families of higher-order accurate time integration algorithms. International Journal of Computational Methods, 2019.
• ^[23]
  W. Kim and J.N. Reddy. A novel family of two-stage implicit time integration schemes for structural dynamics. International Journal of Computational Methods, 2020.
• ^[24]
  H. M. Hilber. Analysis and design of numerical integration methods in structural dynamics. PhD thesis, University of California Berkeley, 1976.
• ^[25]
  T. J. R. Hughes. Analysis of transient algorithms with particular reference to stability behavior. Computational methods for transient analysis(A 84-29160 12-64). Amsterdam, North-Holland, 1983, pages 67-155, 1983.
• ^[26]
  G. M. Hulbert and J. Chung. Explicit time integration algorithms for structural dynamics with optimal numerical dissipation. Computer Methods in Applied Mechanics and Engineering, 137(2):175-188, 1996.
• ^[27]
  C. Semler, W. C. Gentleman, and M. P. Païdoussis. Numerical solutions of second order implicit non-linear ordinary differential equations. Journal of Sound and Vibration, 195(4):553-574, 1996.
• ^[28]
  K. J. Bathe. Conserving energy and momentum in nonlinear dynamics: a simple implicit time integration scheme. Computers & structures, 85(7-8):437-445, 2007.
• ^[29]
  W. Kim. A simple explicit single step time integration algorithm for structural dynamics. International Journal for Numerical Methods in Engineering, 119:383-403, 2019.
• ^[30]
  M. Drolia, M. S. Mohamed, O. Laghrouche, M. Seaid, and A. E. Kacimi. Explicit time integration with lumped mass matrix for enriched finite elements solution of time domain wave problems. Applied Mathematical Modelling, 77:1273-1293, 2020.
• ^[31]
  J. N. Reddy. An Introduction to the Finite Element Method. McGraw-Hill, New York, NY, 4th edition, 2019.
• ^[32]
  K. J. Bathe. Finite Element Procedures. Klaus-Jurgen Bathe, 2006.
• ^[33]
  T. J. R. Hughes. The finite element method: linear static and dynamic finite element analysis. Courier Corporation, 2012.
• ^[34]
  G. Noh and K. J. Bathe. The bathe time integration method with controllable spectral radius: The ρ∞-bathe method. Computers & Structures, 212:299-310, 2019.
• ^[35]
  W. Kim and J. N. Reddy. An improved time integration algorithm: A collocation time finite element approach. International Journal of Structural Stability and Dynamics, 17(02):1750024, 2017.
• ^[36]
  M. M. I. Baig and K. J. Bathe. On direct time integration in large deformation dynamic analysis. In 3rd MIT Conference on Computational Fluid and Solid Mechanics, pages 1044-1047, 2005.
• ^[37]
  T. C. Fung. Solving initial value problems by differential quadrature method part 2: second-and higher-order equations. International Journal for Numerical Methods in Engineering, 50(6):1429-1454, 2001.
• ^[38]
  Y. M. Xie. An assessment of time integration schemes for non-linear dynamic equations. Journal of Sound and Vibration, 192(1):321-331, 1996.
• ^[39]
  J. Liu and X. Wang. An assessment of the differential quadrature time integration scheme for nonlinear dynamic equations. Journal of Sound and Vibration, 314(1):246-253, 2008.
• ^[40]
  D. Soares. A novel family of explicit time marching techniques for structural dynamics and wave propagation models. Computer Methods in Applied Mechanics and Engineering, 311:838-855, 2016.
• ^[41]
  K. J. Bathe and G. Noh. Insight into an implicit time integration scheme for structural dynamics. Computers & Structures, 98:1-6, 2012.
• ^[42]
  D. Soares. An explicit family of time marching procedures with adaptive dissipation control. Inter-national Journal for Numerical Methods in Engineering, 100(3):165-182, 2014.
• ^[43]
  J. N. Reddy. An Introduction to Nonlinear Finite Element Analysis: with applications to heat transfer, fluid mechanics, and solid mechanics. Oxford University Press, 2nd. edition, 2015.

Edited by

Editor:

Marcílio Alves.

Publication Dates

• Publication in this collection
  03 May 2021
• Date of issue
  2021

History

• Received
  26 Feb 2021
• Reviewed
  10 Mar 2021
• Accepted
  15 Mar 2021
• Published
  15 Mar 2021