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

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.


Introduction
For decades, direct time integration schemes have been widely used for the effective analysis of linear and nonlinear structural dynamics [1,2,3,4,5,6,7,8,9,10,11]. 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 [8,12,13,14,15,16,17,18,19]. 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 [14,20,21,22,23]. 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 [4,7,14,20,22,24,25,26,27,28,29,30], 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 [2,31,32,33] 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-approximations 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.

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 ≤ ≤ + , where and are the beginning of the time interval and the size of the time interval, respectively. 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 = + , the velocity and displacement vectors are approximated by using the weighted sums of the velocity and acceleration vectors, respectively as where of ̇+ and + denotes that the approximations belong to the th stage, is the collocation parameter that determines the th time point within the time interval, and and are algorithmic weighting parameters for the velocity and acceleration vectors, respectively. It is noted that 0 and 3 should always be chosen as In the current approach, 1 and 2 can be chosen as 1 of the proposed approximations plays the role of the splitting ratio of the existing composite schemes. In the existing composite schemes [12,14,15,34], 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 and should also satisfy and With the restriction on 33 given in Eq. (6), the unified set of approximations does not include the unknown acceleration vector in the third stage (i.e., for = 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 th stage can be stated in terms of other weighting parameters of the th stage. Since 33 is always zero, 3̈+ 3 is not included in the approximations given in Eqs. (1)-(2). 1̈+ 1 and 2̈+ 2 are the two unknown acceleration vectors to be determined by solving two sets of algebraic equations. If 1̈+ 1 and 2̈+ 2 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 = + is expressed as where is the mass matrix, and is the total force vector. In linear structural dynamics, is often expressed as where is the vector of externally applied force, is the damping matrix, 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. This composite scheme can be summarized as follows. In the first sub-step, 1̈+ 1 is computed by using In the second sub-step, 2̈+ 2 is computed by using By using the 1̈+ 1 and 2̈+ 2 , 3̇+ 3 and 3 are updated as where 3 = 1, 33 = 0, and the computation of 3̈+ 3 is not required. To advance another step, ̈, ̇, and of the next time step are updated by using 2̈+ 2 , 3̇+ 3 , and 3 + 3 , 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., , and ) included in Eqs. (9)-(16).

Implicit schemes
To develop improved implicit composite time schemes with preferable numerical characteristics by using Eqs. (9)-(16), the algorithmic parameters , and should be determined carefully. To optimize these algorithmic parameters, the single-degree-of-freedom problem is used [4,5,25]. The single-degree-of-freedom problem is expressed as where is the displacement, is the externally applied force, is the natural frequency, and is the damping ratio. The initial conditions are given by By setting = 1, = 2 , = 2 , = , = = 0, and = 0, the newly proposed approximations given in Eqs.
, and 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 [5] of a scheme is defined by where 0 = 1, 1 is the trace of , 2 is the sum of the principal minors of , and 3 is the determinant of . An algorithm is th-order accurate if is satisfied. In time integration schemes, the spectral radius of is defined by where is the th eigenvalue of . The spectral radius defined in Eq. (22) is a direct measurement of the stability and the numerical dissipation. It is noted that ( ) is a function of = . A scheme is stable if is satisfied. If Eq. (23) is satisfied for all ∈ [0, ∞), 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 By stating one of the algorithmic parameters of a scheme in terms of ∞ according to Eq. (24), desired level of numerical dissipation in the high-frequency limit can be prescribed through ∞ . The most dissipative case is obtained when ∞ is set as zero, and less dissipative cases are obtained as ∞ approaches to one. When ∞ = 1 is used, the nondissipative case is obtained. The non-dissipative case is usually suitable for most of general analyses. Thus, the case of ∞ = 1 is considered standard. The most dissipative case obtained with ∞ = 0 is also called the asymptotic annihilating case. The case of ∞ = 0 is useful for the elimination of the spurious high-frequency modes.

The first family of implicit schemes
In the first family of implicit schemes, 3̇+ given in Eqs. (12)- (13) to eliminate the spurious root of the characteristic polynomial equation of the amplification matrix given in Eq. (19). This condition can be imposed as With the relations given in Eq. (25), the amplification matrix have two eigenvalues. To make two eigenvalues of the amplification matrix complex conjugate, and can be related as = for = 1,2,3 and = 0,1,2,3 To obtain ̈+ , ̇+ , and + from 2̈+ 2 , 3̇+ 3 , and 3 + 3 , 2 and 3 should be chosen as From Eq. (5), 10 and 20 can be determined as When the parameters are determined according to Eqs. (25)- (28), the new approximations given in Eqs. (11)-(16) can be simplified as To guarantee second-order accuracy, 21 should be chosen to satisfy the condition given in Eq. 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 ∞ according to Eq. (31). When 22 ≠ 11 , 22 can be determined from Eq. (31) as For the dissipation control, ∞ can be chosen from zero to one. If ∞ is chosen in the range of 0 ≤ ∞ ≤ 1, the implicit scheme is unconditionally stable for linear problems. The choice of ∞ = 0 gives the most dissipative case, and the choice of ∞ = 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, 11 can be chosen as When 11 is chosen according to Eq. (33), 22 can be determined from Eq. (31) as 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 [23]. 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., = + ) for the first and second dynamic equilibriums. The case 1-3 is equivalent to the scheme proposed by Kim and Reddy [35]. The case 1-4 is equivalent to the scheme proposed by Kim and Choi [15]. It should be noted that the case 1-3 can include the standard Bathe method [8,36] and the ∞ -Bathe method [34]. 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.  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 ( ) to the complete time step ( ), and 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, 11 may be adjusted to have the effect which is the same as adjusting the splitting ratio in the existing schemes. 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 ∞ , and the level of numerical dissipation in the low-frequency range can also be adjusted through 11 . For the case 1-1, 11 can be chosen in the range of is considered standard in the case 1-1. As 11 approaches to 1 2 , weaker numerical damping is introduced into the low-frequency range. As 11 approaches to 1 2 , the period and damping errors of the low-frequency range become almost identical to those of the trapezoidal rule with while a specified level of numerical dissipation can be imposed in the high-frequency limit through ∞ . On the other hand, stronger numerical damping is introduced into the low-frequency range when 11 is chosen in the range of 1 2 < 11 < 1. As 11 approaches to one, strong numerical damping is introduced into the low-frequency range. However, both period and damping errors increase as 11 approaches to one. Useful discussions regarding the role of 1 are presented in Refs. [12,34,35].

The second family of implicit schemes
In the second family of implicit schemes, 3̇+ To make two eigenvalues of the amplification matrix become complex conjugate, and can be related according to Eq. (26), and 33 is chosen as From Eqs. (5) and (35), 11 , 22 , and 32 of the second family can be determined as Since ̈ 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 ̈0 = −1 ( 0 ,̇0, 0) 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 ̈0 = −1 ( 0 ,̇0, 0) 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 To guarantee second-order accuracy, 31 can be determined to satisfy the condition given in Eq. (21) with = 2. Form this condition 31 is computed as For the second family of implicit schemes, the ultimate spectral radius ∞ defined in Eq. (24) is computed as 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 ∞ according to Eq. (40). From Eq. (40), 21 can be determined as For the dissipation control, ∞ can be chosen from zero to one. To obtained an identical effective system matrix for the first and second sub-steps in linear problems, 1 can be chosen as Latin American Journal of Solids and Structures, 2021, 18(4), e367 10/27 In this case, 21 can be determined from Eq. (40) as 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 [12]. The second general form given in Table 2 can also be used to develop various new schemes in this family.

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 [17], because its exact solution is computable at discrete time points based on the conservation of total energy [11,37]. 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 where, = � / , ( ) is the angle between the rod and the vertical line, is the gravitational constant, and is the length of the rod. Dimensionless values = 1.0 and = 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 To enhance the total energy conserving capability in the implicit schemes, we impose a conditions of where is the targeted order of the convergence rate of the total energy error, and 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 = . To satisfy the condition given in Eq. (47), 11 and 1 of the general form given in Table 1 can be chosen as 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 where is given by From Eq. (50), it can be shown that the choice of ∞ = 1 (the non-dissipative case) gives | − | = ( 5 ). Since the condition given in Eq. (50) is satisfied between the exact initial conditions and the solutions at = , the convergence rate of the total energy error at arbitrary time = + becomes one order lower than that of Eq. (50). Thus, the convergence rate of the total energy error with the choice of ∞ = 1 becomes fourth-order. Similarly, 1 and 2 of the general form of the second family given in Table 2 can be chosen as 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 ∞ = 1 is computed as where is given by 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 ≤ ∞ ≤ 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. 0 = 0 is used, and two different values are considered for ̇0 . ̇0 = 1.999999238456499 and ̇0 = 2.000000761543501 are used to synthesize two highly nonlinear cases [17]. The minimum total energy required for the pendulum to make a complete rotation about the pivot point is 1.0. With 0 = 0 and ̇0 = 1.999999238456499, the total energy of the first case becomes about = 0.999998476913288, and the pendulum should not make a complete rotation. For this case, the period is about = 33.7210.
With 0 = 0 and ̇0 = 2.000000761543501 , the total energy of the second case becomes about = 1.000001523087292, and the pendulum passes the peak points slowly. For this case, the period becomes = 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 ( ∞ = 1) is considered to satisfy the condition of | − | = ( 5 ).
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 ( = /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.   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 [38,39] 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 [22]. The softening problem is described by where > 0.0. In this case, = 100.0, 0 = 4.0, and 0 = 0.0 are also used [38]. The period of the problem is computed as = 1.141876 [11,22]. The motion of the hardening elastic spring is stated as where 1 = 100.0, 2 = 10.0, 0 = 1.5, and 0 = 0.0 [38]. The period of the problem is computed as = 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.

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 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, 0 , 1 , 2 , and 3 can be chosen as which is the same as the cases 1-1 and 1-2 as shown in Table 1. The values of 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 become zero for the undamped case (i.e., = 0), 30 and 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 where 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 according to Eq. (62), the desired level of numerical dissipation at the bifurcation point is prescribed through . The most dissipative case is obtained when is zero, and less dissipative cases are obtained as approaches to one. When is one, the non-dissipative case is obtained. In this particular family of explicit schemes, 20 is determined by where 0 ≤ ≤ 1. When the case of = 1 is used to determine 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 changes from one to zero, which is similar to the role of the ultimate spectral radius ∞ 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 . Like cases of the implicit schemes, the choice of = 1 gives the non-dissipative cases, and more dissipative cases are obtained as approaches to zero. Here, is the value of the spectral radius at the bifurcation point. denotes the critical time step. The critical time step is defined as the largest time step that satisfies the condition ( ) ≤ 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 [23], the case 3-2 is spectrally identical Latin American Journal of Solids and Structures, 2021, 18(4), e367 15/27 to the scheme presented in Ref. [13]. The case 3-3 is spectrally identical to the scheme developed by the author [18]. The explicit scheme proposed by Soares [40] 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., = 0 ). However, it becomes only a conditionally stable implicit scheme for damped cases (i.e., ≠ 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.  In these schemes, the level of numerical dissipation at the bifurcation point is determined by specifying the value of . 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 ̈0 = −1 ( 0 ,̇0, 0) 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 Latin American Journal of Solids and Structures, 2021, 18(4), e367 16/27 case 3-1 is the most accurate.
The key values of , / and / of the cases 3-1, 3-2, and 3-3 are presented in Tables 4 and 5, where is the shortest natural period of a given dynamic system.  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 ≤ /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 = 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. [17], and the cases 4-2 and 4-3 are equivalent to the schemes developed in Ref. [23].
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 . 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.

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.

General advantages of composite schemes
To demonstrate general advantages of using the composite schemes, 0 =̇0 = 0.0, = 0.05, = 2 , ( ) = 10.0 sin( ) + 5.0 cos � 3 � 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-methods [5,26]. Then, numerical results are compared with each other. To guarantee enough accuracy in each scheme, = 0.01 is used for all composite schemes, and = 0.005 (i.e., a half time step) is used for the implicit and explicit generalizedmethods.
In general, the computation cost per time step required in implicit composite time schemes is twice higher than that of the implicit generalized-method. As observed in the results presented in Fig. 8 (a), the maximum accuracy level of the non-dissipative cases of implicit composite schemes ( ∞ = 1.0) is almost the same as the cases of using the trapezoidal rule and the generalized-method ( ∞ = 1.0) twice with a half time step. On the other hand, the asymptotic annihilating cases ( ∞ = 0.0) of the implicit composite schemes are presenting much smaller errors than the asymptotic annihilating cases ( ∞ = 0.0) of the generalized-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-method gives less accurate solutions when compared with its nondissipative 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 ∞ 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.

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 [14,41].
In the literature [12,14,40,41], 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. [41], the spring constant of the stiff spring is chosen as 2 = 10.0 , where = 7. 1 = 0.0, 2 = 1.0, 3 = 1.0, = 6/5 and zero initial conditions are used. The spring constant of the soft spring is 3 = 1.0. With = 7, two natural frequencies of the system are given by 1 = 3162.27781828 and 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 ( ∞ = 0) of the implicit composite schemes can effectively filter out the high frequency of the system as presented in Refs. [12,41].
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 is chosen as 2. In the explicit case 3-1 with = 0.0, = 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-method with ∞ = 0.0, three different values are used for to verify that the implicit composite scheme is not suitable for this particular case [42].
With = 2 , two natural frequencies of the system are computed as 1 = 10.0503730777 and 2 = 0.994987939525. With 1 and 2 , two natural periods are computed as 1 = 2 / 1 = 0.625169360244 and 2 = 2 / 2 = 6.31483564532. The bifurcation point of the explicit scheme of the case 3-1 with = 0.0 is computed as = 0.52542780235992, where is the natural period. When is chosen as × 1 , the maximum amount of numerical damping is introduced into the mode associated with 1 . With = 0.52542780235992 × 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 . As shown in Fig. 12, the case 3-1 (explicit) can introduce the maximum numerical dissipation into the mode associated with 1 , while providing enough accuracy for the mode associated with 2 . On the other hand, the case 1-3 (implicit) can introduce only a moderate amount of numerical dissipation in the highfrequency mode. This explains the reason that the explicit case 3-1 can give better high-frequency filtering as shown in Fig. 13 when = 0.328481363056 is used in both cases. Like the explicit case 3-1, the explicit generalized-method also gives good high-frequency filtering capability when = 0.140712 is used. To improve the results of the implicit case 1-3 and the implicit generalized-method for this case, larger time steps are employed. When 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 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 [12]. 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 For the spatial discretization of the problem given in Fig. 15, one thousand uniform linear elements are used. After the spatial discretization, , , , and of this problems are given by   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.

Splitting ratio and alternative parameters
As shown in Table 1, 1 of the case 1-1 and 1-2 is always one. The setting of 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 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 [34,35] with the setting of algorithmic parameters given in the case 1-3. In the existing implicit composite methods, on the other hand, 1 is adjusted to synthesize preferable numerical characteristics, such as the level of numerical dissipation in the low-frequency range [34,35]. Then, the time point of the first dynamic equilibrium given in Eq. (11) also changes accordingly.
In Ref. [34], it has been stated that the choice of 1 = 1.99 and ∞ = 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 1 = 0.5 and ∞ = 0.0. However, the choice of 1 = 1.99 can cause a serious problems in some situations. Basically, the choice of 1 = 1.99 worsens the accuracy. To remedy this, very small time step (due to the use of CFL=0.1) should be employed [34]. 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 11 = 0.995 and ∞ = 0.0 can give the same high-frequency filtering capability of the case 1-2 with 1 = 1.99 and ∞ = 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 11 is safer than varying 1 beyond the range of the current time interval in the case 1-3.  To show a possible misuse of the splitting ratio in the existing implicit composite schemes more clearly, 0 =̇0 = 0.0, = 0.05, and = 2 are used for the single-degree-of-freedom problem given in Eqs. (17)- (18). Then, ( ) = 10.0 is considered for 0.0 ≤ ≤ 1.0, and ( ) = 0.0 is considered for > 1.0 as shown in Fig. 19. For the case 1-3, 1 = 1.01 and ∞ = 0.0 are used. In the case 1-1, however, 11 = 0.505 and ∞ = 0.0 are used to introduce equivalent amounts of numerical damping in the low-frequency range when compares with the case 1-3 with 1 = 1.01 and ∞ = 0.0. To guarantee enough accuracy, = 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 = 1.0 as shown in Fig. 20. To compute the variables at = 1.0, two dynamic equilibrium equations at = 1.0005 and = 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 (1.0005) = 0.0 and (1.0) = 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, ( ) = 0.0 for > 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 = 1.0, because (1.0) = 1.0 is used for the first and second sub-steps.

Computational cost
Among the implicit schemes, the second general form, the case 2-1, and the case 2-2 do not require the calculation ̈0 = −1 ( 0 ,̇0, 0). 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 and ̇+ 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 [32,43], 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 [31,35]. 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. [12,18,23].

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-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.