## 1 Introduction  
To solve general structural dynamic problems, the direct integration method has been proven to be effective and has been widely used in engineering. The commonly used integration methods can be categorized into the explicit one and the implicit one. The explicit method is conditionally stable, and a smaller time step should be used in the calculation. Thus, for structural vibration analysis, the implicit method is usually considered to be more effective [1]. As for implicit method, it can be designed to be unconditionally stable and to possess controllable numerical dissipation. Although the calculation needs matrix factorization and is cost-consuming [2], the implicit method can allow a larger time step due to its superior numerical stability, which brings much more efficiency in solving dynamic problems [3]. At present, there are some well-known implicit algorithms, such as the Newmark method [4], the Wilson-\(\theta\) method [5], the generalized-\(\alpha\) method [6], the HHT method [7], and the Bathe method [8], etc. However, some algorithms may lead to instability in the nonlinear range [9] due to non-conservation of energy and momentum [10]. Taking trapezoidal rule for example, in the linear range, it is unconditionally stable, and has second-order accuracy without amplitude decay [11, 12]. However, this method can easily become unstable when extended to the nonlinear problems [12, 13], and the response will diverge [3, 10, 14] when using an inappropriate length of time steps. To solve this problem, numerical damping is introduced into the Newmark method by adjusting the parameters, but it results in first-order accuracy [15]. For most of the structural dynamic problems, the main contribution of the response comes from its low-order frequency modes [16]. Considering the impact of discretization, the high-order frequency modes of the structure will change, which affects the computational accuracy. Therefore, the numerical damping which can control numerical dissipation is necessary for both low-order and high-order frequency ranges. In addition, it can also improve the convergence of the iterative solution and increase the numerical stability. Nevertheless, unreasonable numerical damping will result in a great loss of accuracy, which leads to a large attenuation of structural response. It is usually caused by excessive numerical dissipation in a low-frequency range. To evaluate the numerical dissipation of the algorithm, the spectral radius \(\rho\) is used. The change in dissipation is reflected in the spectral radius \(\rho\); that is, the reduction in spectral radius illustrates the dissipation of the algorithm. When the time step is large, the spectral radius will quickly attenuate for some algorithms, such as part of the Newmark method, Houbolt method, Wilson-\(\theta\) method, and part of the generalized \(\alpha\) method, etc. [17]  
In summary, it is of necessity to find a time integration algorithm that has both reasonable and controllable numerical damping. This algorithm can generate good accuracy in a low-frequency range with controlled high dissipation in a high-frequency range [7; 18; 19], and is expected to perform well in nonlinear field. Many researchers have proposed useful methods to control the accuracy of algorithms in different frequency ranges. Some composite algorithms with several steps perform quite well [20; 21; 22]. A typical one is the Bathe method proposed by [PERSON] et al., which divides a time step into two sub-steps to combine two well-known time integration schemes [3; 8; 12; 23]. The trapezoidal integration method and the three-point Euler backward method are respectively used in the first and second sub-step. The Bathe method has reasonable numerical damping and stability but cannot control the dissipation [24]. Based on this, the \(\rho_{\theta}\)-Bathe method and \(\beta_{1}/\beta_{2}\)-Bathe method were proposed [23; 25].