High-order accurate multi-sub-step implicit integration algorithms with dissipation control for second-order hyperbolic problems

Recently, some novel implicit integration methods, namely the so-called composite sub-step algorithms, have attracted the attention of some researchers. A composite two-sub-step algorithm was analyzed by Bathe [12] to successfully solve the strongly nonlinear problems that the non-dissipative trapezoidal rule fails to solve. The Bathe algorithm uses the composite two-sub-step technique. Firstly, the current integration interval $t\in[t_{n},~{}t_{n+1}]$ is split into two sub-intervals $[t_{n},~{}t_{n}+\gamma\varDelta\!t]$ and $[t_{n}+\gamma\varDelta\!t,~{}t_{n+1}]$, where $\gamma$ denotes the splitting ratio of sub-step size. Then, the non-dissipative trapezoidal rule and the three-point backward difference formula are used in the first and second sub-steps, respectively. Finally, the resulting scheme achieves L-stability, second-order accuracy, non-overshooting, and self-starting features. The unique free-parameter $\gamma$ adjusts numerical dissipation imposed in the low-frequency range. It has been shown in [13] that the Bathe algorithm [12] corresponds to the non-dissipative trapezoidal rule in either $\gamma=0$ or $1$, while the first-order backward Euler scheme is covered in $\gamma=2$. These facts also explain the reason why $\gamma$ can adjust numerical dissipation in the low-frequency range. When the parameter $\gamma$ gets close to either $0$ or $1$, the numerical behavior of the Bathe algorithm [12] naturally tends to that of the non-dissipative trapezoidal rule, so controlling dissipation in the low-frequency range. More importantly, the Bathe method [12, 14] corresponds essentially to the two-stage diagonally implicit Runge-Kutta (RK) method proposed by Bank et al. [15]. Hence, the composite sub-step methods, featuring consistent displacement and velocity updating formulas, essentially form part of the RK family. Over the past decades, some composite sub-step implicit methods developed in computational mechanics have predominantly belonged to the implicit RK schemes and are briefly viewed as follows.

The composite multi-sub-step implicit algorithms using the non-dissipative trapezoidal rule and more general multi-point backward difference scheme were constructed and analyzed in the literature [16, 17, 18, 19]. Besides, other integration algorithms without adopting the trapezoidal rule in the first sub-step were successfully developed in [20, 21, 22] and they can predict more accurate solutions for solving some dynamical problems. Two and three sub-steps have been widely used to construct composite sub-step algorithms since these sub-step schemes can be readily derived and analyzed. For example, another two Bathe algorithms, named as $\rho_{\infty}$-Bathe [17] and $\beta_{1}/\beta_{2}$-Bathe [23], are composite two-sub-step schemes, and they actually share the completely same algorithm structure [13]. The theoretical analysis [13] has shown that the second-order $\beta_{1}/\beta_{2}$-Bathe algorithm is algebraically identical to the $\rho_{\infty}$-Bathe algorithm, and the first-order $\beta_{1}/\beta_{2}$-Bathe algorithm is not competitive with other higher-order methods. The researches have clearly shown that the composite multi-sub-step implicit algorithms attain optimal spectral properties, namely the maximum numerical dissipation but the minimum relative period errors, when achieving identical effective stiffness matrices within each sub-step. Hence, the optimal composite two-sub-step implicit algorithm achieving dissipation control, second-order accuracy, and no overshoots has been proposed in [24]. Other composite two-sub-step schemes, such as $\rho_{\infty}$-Bathe [17], can cover it when achieving optimal spectral properties. Similarly, an optimal three-sub-step implicit method has been also considered in [18], where identical effective stiffness matrices and controllable algorithmic dissipation in the high-frequency range are achieved. All composite multi-sub-step implicit algorithms mentioned above are only second-order accurate.

High-order accurate integration algorithms often require much more computational cost than the common second-order ones. So far, there have been a great number of studies [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38] reporting high-order accurate integration algorithms. A preferred way to develop high-order schemes is to operate the first-order differential systems by converting the original second-order equations of motion. The studies [25, 29, 32, 31] used the time finite elements and the weighted residual approach to address the resulting first-order systems. These high-order schemes suffer mainly from two drawbacks. One is that the resulting integration algorithms cannot output acceleration responses, although they are directly self-starting. Moreover, for solving structures subjected to the external loads, the lack of external load analysis could cause the order reduction in displacement and velocity. For example, the third/fourth-order accurate schemes [25] from the extrapolated Galerkin time finite elements are only of second-order for solving forced vibrations. The other is how to compute external loads effectively and accurately. Because the load terms are often expressed as the integral form in the current time interval, how to compute these loads in the integral form is crucial to maintain the expected order of accuracy. It is a known fact [37, 27, 39] that the $(p+1)$th-order approximation for the second-order dynamical problems imposes $p$ unknown variables, and thus the order of accuracy can be optimized up to $(2p-1)$ and $2p$. In general, the developed high-order algorithms can be designed to have controllable dissipation capability. The high-order methods are generally $(2p-1)$th- and $2p$th-order accurate in the dissipative and non-dissipative cases, respectively. It should be emphasized that the high-order methods derived from the transformed first-order systems have to solve the equation system of dimension $pd\times pd$ where $p$ and $d$ denote the number of the involved variables and the degree-of-freedom of the spatially discretized model, respectively. In the work [37], the $p$th-order Lagrange polynomial and Gauss-Lobatto integration points are used to develop $(2p-1)$th- and $2p$th-order accurate algorithms, and a family [27] of high-order algorithms considers $p$th-order Hermite interpolation polynomials and higher-order time derivatives of the equilibrium equations. Note that these families [37, 27] of high-order algorithms belonging to the fully implicit RK family must solve the equation systems with higher dimensions per time step, so increasing computational costs. It should be pointed out that the effective stiffness matrices of these two methods [37, 27] are neither symmetric nor sparse. Some existing high-order implicit methods [40, 41, 31, 29, 42] suffer from these issues. Hence, some researchers [43, 31, 37] proposed various iterative solvers to solve the resulting higher dimensional equation systems. A recent overview of these high-order implicit algorithms can refer to the literature [39].

On the other hand, high-order integration algorithms [26, 33, 36] directly operating the second-order differential systems show some advantages over other high-order ones, but there are some shortcomings. One of the main advantages of these high-order implicit methods is avoiding solving higher dimensional equation systems per time step, so significantly reducing computational costs. For instance, the Tarnow and Simo technique [26] is often applied to the single-step non-dissipative schemes, such as the trapezoidal rule and the midpoint rule. When some single-step dissipative schemes are used in the Tarnow and Simo technique, the resulting methods suffer from the order reduction. Hence, the high-order schemes from the Tarnow and Simo technique are generally non-dissipative. If the Tarnow and Simo technique is applied to the midpoint rule, the resulting scheme also covers a directly self-starting, fourth-order accurate, and non-dissipative method [44]. A single-step high-order implicit method based on Pad$\acute{\text{e}}$ expansion was developed in [45]. The developed method [45] is always non-dissipative due to the use of diagonal Pad$\acute{\text{e}}$ approximation, and it involves the complicated calculations of external loads and complex operations. Soares developed a directly self-starting fourth-order integration scheme [46] for solving undamped models. It is reduced to be third- and second-order accurate for numerically and physically damped models, respectively. However, it is only conditionally stable, so its effective applications may be limited. Zhang et al. [33] adopted the composite sub-step technique to develop an implicit family of high-order algorithms with controllable high-frequency dissipation and acceptable computational cost, but these algorithms also present lower-order of accuracy for solving forced vibrations [36]. The second-order accurate $\rho_{\infty}$-Bathe algorithm [17] is reanalyzed to achieve third-order accuracy [47] via considering optimal load selections. The resulting third-order scheme [47, 48] cannot achieve a full range of numerical dissipation and the load information in the previous time step has to be used in the current step. The $\rho_{\infty}$-Bathe method has been revisited [48] to analyze the influence of the time splitting ratio $\gamma$ on accuracy. When a complex-valued $\gamma$ is adopted, the $\rho_{\infty}$-Bathe method can achieve a full range of dissipation control and the fourth-order accuracy is obtained in the non-dissipative case.

Adopting the composite sub-step technique, the authors have already developed an implicit family [36] of high-order algorithms which correspond essentially to the singly diagonally implicit Runge-Kutta methods (SDIRK). The developed algorithms are directly self-starting, avoiding both calling any starting procedures and computing the initial acceleration vector. In general, the directly self-starting $s$-sub-step implicit methods [36] can achieve $s$th-order accuracy, unconditional stability, no overshoots, and controllable high-frequency dissipation. As stressed by Li et al. [36], the directly self-starting high-order algorithms require using the second-order equations of motion once additionally to guarantee identical $s$th-order accuracy, which is computationally expensive. Otherwise, the methods [36] cannot output the same acceleration accuracy as those in displacement and velocity. In addition to the considerable computational effort required for the output of acceleration responses, the directly self-starting high-order methods also show poorer robustness for solving complicated dynamical problems due to losing the acceleration $\ddot{\mathbf{U}}_{n}$ in all sub-steps. Almost all directly self-starting algorithms are faced with these computational issues, and further studies on analyzing the pros and cons of directly self-starting algorithms are needed in the near future.

To overcome these issues, the authors will construct and develop a novel implicit family of high-order algorithms, corresponding to the explicit singly diagonally implicit Runge-Kutta methods (ESDIRK), which achieves the following desirable numerical characteristics.

• •

  The self-starting property, avoiding the imposition of an additional starting procedure, has been emphasized in the work of Hilber and Hughes [49]. They underscored that a non-self-starting algorithm often necessitates more in-depth analysis and consideration and the incorporation of an additional starting procedure may result in heightened coding and computational costs. Thus, the algorithm should be crafted with the intention of being (directly) self-starting. However, as mentioned earlier, the directly self-starting algorithms come with certain computational disadvantages due to the loss of $\ddot{\mathbf{U}}_{n}$ in all sub-steps. Consequently, the novel algorithms are intentionally designed to be self-starting but not directly self-starting.

• •

  Identical effective stiffness matrices across all sub-steps, incorporating optimal spectral properties, play a crucial role. The identity of effective stiffness matrices not only decreases computational expenses in solving linear elastic problems within multi-sub-step implicit algorithms but also integrates optimal spectral properties. When traditional implicit algorithms are considered as composite single-sub-step schemes, this property is automatically realized in all traditional algorithms. Remarkably, this property has been conceptually extended to the advancement of multi-sub-step explicit methods [50, 51].

• •

  The ability to control numerical dissipation in the high-frequency range, aiding in the damping out of spurious modes and stabilizing iterative procedures, is recognized as a desirable property in time integration algorithms. To accommodate complicated solutions, the novel methods should be meticulously designed to achieve comprehensive dissipation control, ranging from the non-dissipative scenario to the asymptotically annihilating case.

• •

  Higher-order accuracy in displacement, velocity, and acceleration is crucial for solving general structures. In this paper, novel methods are formulated with real-valued parameters to attain identical high-order accuracy. The sustained high-order accuracy achieved by these novel methods is demonstrated in the resolution of a comprehensive benchmark problem.

• •

  Unconditional stability, rendering the integration step constrained by accuracy rather than stability. Implicit algorithms must unquestionably be devised to possess unconditional stability, as lacking this feature would diminish their competitive edge in practical applications. Note that, achieving identical effective stiffness matrices, controllable numerical dissipation, and high-order accuracy, the proposed methods lack additional parameter design flexibility for further enhancement of BN-stability [52].

• •

  Achieving a favorable equilibrium between computational cost and high-order accuracy is essential. Similar to the existing high-order implicit methods [26, 44, 45, 46, 33, 47, 36], the novel methods necessitate solving linearized equations with the same dimension as the second-order equations of motion per (sub-)step. Consequently, the resulting effective stiffness matrix inherits all desired properties of the global mass, damping, and stiffness matrices, including sparsity, symmetry, and positive definiteness.

The remaining sections of this work are structured as follows. Section 2 introduces the construction of a novel composite $s$-sub-step implicit scheme, where the conditions for achieving $p$th-order accuracy and controllable numerical dissipation are derived. In accordance with these conditions, Section 3 details the development of novel high-order implicit algorithms. Section 4 provides spectral analysis and comparisons with existing high-order algorithms, confirming the well-established unconditional stability and dissipation control of the proposed high-order algorithms. Section 5 presents numerical examples to validate the numerical performance and superiority of the four novel high-order schemes. Lastly, Section 6 summarizes essential conclusions drawn from this study.

For the sake of simplify and clarity, the constructed $s$-sub-step implicit algorithm is described to solve the linear second-order equations of motion (1). Considering the current integration internal $t\in\left[t_{n},~{}t_{n}+\varDelta\!t\right]$, the composite sub-step technique firstly divides the whole time increment into $s$ sub-intervals $\cup_{i=1}^{s}[t_{n}+\gamma_{i-1}\varDelta\!t,~{}t_{n}+\gamma_{i}\varDelta\!t]$, where $\gamma_{i}$ denotes the splitting ratio of the $i$th sub-step. In this paper, $\gamma_{0}=0\ \text{and}\ \gamma_{s}=1$ are always used to provide acceleration responses. Then, an integration scheme in the $i$th sub-step $t\in[t_{n}+\gamma_{i-1}\varDelta\!t,~{}t_{n}+\gamma_{i}\varDelta\!t]~{}(i=1,~{}\dots,~{}s)$ is constructed as

Notice that numerical solutions ($\mathbf{U}_{n+1},~{}\mathbf{\dot{U}}_{n+1}$, and $\mathbf{\ddot{U}}_{n+1}$) at time $t_{n+1}$ are obtained in the last sub-step $t\in[t_{n}+\gamma_{s-1}\varDelta\!t,~{}t_{n}+\gamma_{s}\varDelta\!t]$ due to $\gamma_{s}=1$. That is,

As stressed previously, the composite sub-step methods can also be viewed as a special case of the Runge-Kutta family, and thus the proposed multi-sub-step method (2) can be described in the Butcher tableau [53] as

This paper endeavors to identify algorithmic parameters conducive to achieving identical high-order accuracy and controllable numerical dissipation within the framework of computational structural dynamics. As a result, the developed high-order algorithms given by Eq. (4) are also applicable to the first-order transient problems.

In the initial time $t_{0}$, the initial acceleration vector $\mathbf{\ddot{U}}_{0}$ is given by solving the equilibrium equation at $t_{0}$ with given initial conditions $\mathbf{U}_{0}$ and $\mathbf{\dot{U}}_{0}$.

Hence, the proposed $s$-sub-step implicit algorithm (2) is obviously self-starting, like the known WBZ-$\alpha$ [9], HHT-$\alpha$ [8], and TPO/G-$\alpha$ [10, 11] algorithms.

It is cumbersome and difficult to analyze numerical properties of an integration algorithm by directly manipulating the multi-degree-of-freedom system (1). However, the modal decomposition technique [1] can be employed to uncouple the system (1) by using the orthogonality properties of the free-vibration mode shapes of the undamped system. As a result, the modal damping is often used. It is therefore convenient to analyze numerical properties of an integration method by solving the standard SDOF system [55, 5, 56], which is written as

where $\xi$, $\omega$, and $f(t)$ are the damping ratio, the undamped natural frequency of the system, and the modal force, respectively.

For the subsequent use, considering the above SDOF system with given initial displacement $u_{n}$ and velocity $\dot{u}_{n}$ and external load $f(t)=\exp(t)$, its exact displacement and velocity are analytically obtained as

where $\tau=t-t_{n}$ and $\omega_{d}=\sqrt{1-\xi^{2}}\cdot\omega$. Then, the exact solutions at time $t_{n+1}$ can be collected as

where $\mathbf{D}_{exa}$ is the exact amplification matrix [36, 57], which is

Note that the external load $f(t)=\exp(t)$ is not the only choice in accuracy analysis, and other appropriate external loads, such as $f(t)=\sin(t)+\cos(t)$, can also be used to derive the same conditions as $f(t)=\exp(t)$. Appendix A provides the theoretical explanation for the rationality of using $f(t)=\exp(t)$.

On the other hand, applying the proposed implicit algorithm (2) to solve the SDOF system (8) can yield

Inserting Eq. (12b) into Eq. (12a) to eliminate $u_{n+\gamma_{i}}$ gives

Define the notations $\dot{\mathbf{Y}}$ and $\ddot{\mathbf{Y}}$ as $\dot{\mathbf{Y}}=\begin{bmatrix}\dot{u}_{n+\gamma_{0}}&\dot{u}_{n+\gamma_{1}}&\dots&\dot{u}_{n+\gamma_{s}}\end{bmatrix}^{\mathsf{T}}$ and $\ddot{\mathbf{Y}}=\begin{bmatrix}\ddot{u}_{n+\gamma_{0}}&\ddot{u}_{n+\gamma_{1}}&\dots&\ddot{u}_{n+\gamma_{s}}\end{bmatrix}^{\mathsf{T}}$, respectively, and also define $\mathbf{f}$ as $\mathbf{f}=\begin{bmatrix}f(t_{n}+\gamma_{0}\varDelta\!t)&f(t_{n}+\gamma_{1}\varDelta\!t)&\dots&f(t_{n}+\gamma_{s}\varDelta\!t)\end{bmatrix}^{\mathsf{T}}$. Then, Eq. (13) with varying $i=\{0,~{}1,~{}\cdots,~{}s\}$ can be written as

where $\mathbf{A}$ is defined in Eq. (4); $\mathbf{1}$ is the column vector whose all elements are unity, and $\mathbf{I}$ is the unity matrix with dimension $s+1$. Similarly, Eq. (12c) can be rewritten as

Solving Eqs. (14) and (15) yields

where $\mathbf{0}$ denotes a column vector whose all elements are zero.

Finally, numerical displacement $u_{n+1}$ and velocity $\dot{{u}}_{n+1}$ are updated in the last sub-step due to $\gamma_{s}=1$, that is

The above equations can be further written as

where the column vector $\mathbf{b}$ is defined in Eq. (4). Substituting Eq. (16) into Eq. (18) yields

where $\mathbf{D}_{num}$ and $\mathbf{L}_{num}$ are called the numerical amplification matrix and load operator, respectively.

With the exact and numerical iterative relationships in hand, the conditions for achieving the novel method’ $p$th-order accuracy can be defined as follows.

Notice that Proposition 1 actually requires the numerical scheme in the last sub-step to be $p$th-order accurate due to $\gamma_{s}=1$. To eliminate the redundant algorithmic parameters, an additional requirement is imposed in this paper. In the $i$th sub-step $t\in\left[t_{n}+\gamma_{i-1}\varDelta\!t,~{}t_{n}+\gamma_{i}\varDelta\!t\right]$, the displacement predicted by Eq. (12b) is defined to be $q$th-order consistency if $u(t_{n}+\gamma_{i}\varDelta\!t)-u_{n+\gamma_{i}}=\mathcal{O}(\varDelta\!t^{q})$ is fulfilled, namely,

Calculating Taylor series expansions of $u(t_{n}+\gamma_{i}\varDelta\!t)$ and $\dot{u}(t_{n}+\gamma_{j}\varDelta\!t)$ at time $t_{n}$ gives, respectively,

Substituting Eq. (23) into Eq. (22) and then simplifying the result yield

Then, the following conditions are satisfied to eliminate low-order terms, reaching a third-order consistency [21].

Similarly, the velocity predicted by Eq. (12c) is defined to be $q$th-order consistency if $\dot{u}(t_{n}+\gamma_{i}\varDelta\!t)-\dot{u}_{n+\gamma_{i}}=\mathcal{O}(\varDelta\!t^{q})$ is fulfilled. Adopting the same analysis as Eq. (24), one can easily get

Eq. (26) demonstrates that the conditions (25) also ensure a third-order consistency in velocity (12c). As stated previously, the published directly self-starting implicit algorithms [36] only require a second-order consistency in each sub-step.

The previous subsection has derived the numerical amplification matrix $\mathbf{D}_{num}$, thus it is very easy to design controllable numerical dissipation in the high-frequency range. Firstly, the characteristic polynomial of $\mathbf{D}_{num}$ is computed as

where $\mathbf{I}_{2}$ denotes the unity matrix with dimension $2$, and coefficients $A_{1}$ and $A_{2}$ are two principal invariants of $\mathbf{D}_{num}$. When considering the high-frequency limit ($\omega\to\infty$), the characteristic polynomial (28) can be expressed as

where $A_{1\infty}=\lim\limits_{\omega\to\infty}A_{1}$ and $A_{2\infty}=\lim\limits_{\omega\to\infty}A_{2}$. On the other hand, the optimal high-frequency dissipation [11] is achieved if and only if eigenvalues of the amplification matrix $\mathbf{D}_{num}$ become the same real roots, denoted by $\rho_{\infty}$, in the high-frequency limit. Then,

Comparing Eqs. (29) and (30) yields the conditions to achieve controllable high-frequency dissipation:

It is necessary to point out that two invariants $A_{1}$ and $A_{2}$ given by Eq. (28) are often employed to analyze the stability of an integrator. When the following conditions [1] are fulfilled, the proposed algorithms are said to be unconditionally stable.

After embedding high-order accuracy and controllable numerical dissipation, the developed implicit algorithm (2) will achieve unconditional stability via selecting proper $\gamma_{1}$. It should be emphasized that, after achieving identical high-order accuracy and controllable numerical dissipation, the developed $s$-sub-step method can only achieve fundamental spectral stability instead of the BN-stability [52].

For the three-sub-step scheme, the detailed analysis will be carried out to clearly demonstrate how to determine the unknown algorithmic parameters $\alpha_{ij}$. The proposed $s$-sub-step method (4) in the case of $s=3$ can be rewritten as ($\gamma_{3}=1$)