Improved second-order unconditionally stable schemes of linear multi-step and equivalent single-step integration methods

Second-order unconditionally stable schemes of linear multi-step methods, and their equivalent single-step methods, are developed in this paper. The parameters of the linear two-, three-, and four-step methods are determined for optimal accuracy, unconditional stability and tunable algorithmic dissipation. The linear three- and four-step schemes are presented for the first time. As an alternative, corresponding single-step methods, spectrally equivalent to the multi-step ones, are developed by introducing the required intermediate variables. Their formulations are equivalent to that of the corresponding multi-step methods; their use is more convenient, owing to being self-starting. Compared with existing second-order methods, the proposed ones, especially the linear four-step method and its alternative single-step one, show higher accuracy for a given degree of algorithmic dissipation. The accuracy advantage and other properties of the newly developed schemes are demonstrated by several illustrative examples.


Introduction
Direct time integration methods are powerful numerical tools for solving the ordinary differential and differential-algebraic equations arising in structural dynamics, multi-body dynamics and many other fields. A lot of excellent integration methods have been developed and are available in the literature; they are often classified as explicit and implicit methods. Explicit methods are easy to implement, but their intrinsic conditional stability limits the step size. Implicit methods can be designed to be unconditionally stable, so they are suitable when the accuracy requirements are not very strict.
In addition, time integration methods also can be categorized as multi-step and multi-stage techniques. The multistep methods employ the states of several previous steps to express the current one, whereas the multi-stage utilize the states of selected collocation points. The single-step methods, which only use the information of the last step, may be B Huimin Zhang huimin.zhang@polimi.it 1 without acceleration vectors [22], were subsequently proposed. All these methods include the standard Bathe method as a special case; they present the same properties with an appropriate choice of the parameters [26]. To achieve higher accuracy, the three-and four-sub-step implicit composite methods [8,18,29,30,42] were also developed, adopting similar ideas. In addition, many efforts have been devoted to the development of explicit composite methods, producing several good candidates, such as the Noh and Bathe [34], the Soares [37], and many other methods [21,24,43].
In the multi-step class, including many existing singlestep methods, Dahlquist [11] has revealed that the methods of higher than second-order accuracy cannot achieve unconditionally stable, known as Dahlquist's second barrier. As a result, for implicit methods in this class, second-order unconditionally stable schemes are of most concern.
Owing to the simplicity, the single-step methods have received a lot of attention after the 1950s, generating many widely-used schemes, including the Newmark method [33], the Wilson-θ method [39], the HHT-α method proposed by Hilbert, Hughes, and Taylor [15], the WBZ-α method by Wood, Bossak, and Zienkiewicz [40], the generalizedα (G-α) method [9,17], the GSSSS (generalized single-step single-solve) method [45], and others reviewed in [38]. These single-step methods can display the identical spectral characteristics of the corresponding linear multi-step methods [45]. However, the linear two-step method with optimal parameters [31,44] can provide better accuracy than most existing single-step methods under the same degree of algorithmic dissipation.
Besides, representative multi-step methods also include the explicit Adams-Bashforth methods [2], the implicit Adams-Moulton methods [13] and the backward differentiation formulas (BDFs) [12]. The Adams-Bashforth and Adams-Moulton methods are designed to achieve higher-order accuracy, but most of them have poor stability. In practical applications, the implicit Adams-Moulton methods are usually used in tandem with the Adams-Bashforth methods as predictor-corrector pairs, which are essentially explicit. The BDFs, although unconditionally stable only up to second order, are particularly useful for solving stiff problems and differential-algebraic equations, owing to their strong algorithmic dissipation. One drawback of the multi-step methods is that they are not self-starting, so another method has to be used to solve the first steps, until there are enough to start the multi-step procedure. As a consequence, multi-step methods may not be as convenient to use as single-step ones.
This paper focuses on second-order accurate and unconditionally stable schemes of the linear multi-step methods. To acquire these accuracy and stability properties, the linear two-step method has been originally studied in [31,32,44]. Seeking higher accuracy, the linear three-, and four-step methods are discussed in this paper. From linear analysis, their optimal parameters are determined by minimizing the global error, under the constraints of second-order accuracy, unconditional stability and tunable algorithmic dissipation. The resulting optimal schemes present better accuracy than the linear two-step method, the linear four-step scheme being the best among them. Furthermore, to avoid additional startup procedures, the single-step methods, spectrally equivalent to the linear multi-step ones, are proposed as well. Instead of utilizing the states of previous steps, the single-step methods get additional information by introducing intermediate variables. With proper parameters, they can achieve the same properties as the linear multi-step methods. Finally, the newly developed methods are applied to solve a few illustrative examples, and the solutions are compared with those of some widely-used implicit second-order schemes.
This paper is organised as follows. The formulations of the linear multi-step methods and the corresponding singlestep ones are presented in Sect. 2. The optimization of the parameters is implemented in Sect. 3. The detailed properties of the proposed schemes are discussed in Sect. 4. Numerical examples are provided in Sect. 5, and conclusions are drawn in Sect. 6.

Formulations
The linear multi-step method for solving the initial-value problem of the generic first-order implicit differential equation f (x,ẋ, t) = 0, x(t 0 ) = x 0 can be cast as where x k andẋ k are the numerical approximations of x(t k ) andẋ(t k ), respectively, t = t k − t k−1 is the time step size, assumed as constant for simplicity, α j and β j are scalar parameters. The initial condition x 0 is supposed to be given in advance.
Applied to the linear structural dynamics problem, expressed as where M, C, K are the mass matrix, damping matrix, and stiffness matrix, respectively, R(t) is the external load, q, q andq denote the displacement, velocity and acceleration vectors, respectively, q 0 and v 0 are the known initial displacement and velocity vectors, this method can be formulated as From Eq. (3), q k ,q k andq k can be solved by combining the linear multi-step scheme and the discrete dynamics equation at t k . Table 1 provides the step-by-step solution of the linear r -step method for solving the linear structural dynamics problems. To start the procedure of the linear r -step method, the initial k < r steps need to be prepared by another method. The single-step scheme which can be obtained by setting α 1 = 1, β 1 = 1 − β 0 and α j = 0, β j = 0 ( j = 2, 3, · · · , r ) in Eq. (1b), is adopted in Table 1. Certainly, other methods, such as the trapezoidal rule, and the combination of single-step and multi-step methods, can also be used for solving the first r − 1 steps. The scheme in Eq. (4) is employed since it shares the same effective stiffness matrix as the linear r -step method. When applied to linear systems, no additional matrix factorization is requested. With the start-up procedure, the linear multi-step method almost does not require additional computational costs compared with the existing single-step methods. The employed single-step scheme is first-order accurate and unconditionally stable when β 0 > 1/2, and it becomes the second-order trapezoidal rule when β 0 = 1/2. When applied to nonlinear problems, the Newton-Raphson iteration needs to be used to solve the nonlinear equations per step.
To completely save the start-up procedure of initial steps, the single-step method, which can be used as an alternative of the linear r -step method in Eq. (1b), is also proposed as . . .
Clearly, a series of intermediate variables are introduced by this method. These variables are only useful to proceed the current and next step, but do not need to be stored and output. From Eq. (5), they can be seen as the shift approximations ofẋ k , soẋ 1 0 =ẋ 2 0 = · · · =ẋ r −1 0 =ẋ 0 is assumed at the initial time. By choosing appropriate values for the γ j parameters, the single-step method of Eq. (5) can achieve spectral characteristics identical to those of the linear r -step method of Eq. (1b), which will be demonstrated in Sect. 3. Table 2 presents the step-by-step solution of the singlestep method for solving the linear structural dynamics problems. The intermediate variables of step k − 1 are useful to update the state vectors of the current step, and after obtaining the intermediate variables of step k, the storage occupied by these variables of step k − 1 can be released. Therefore, this single-step method only requires to store 2(r −1) vectors additionally. The additional storage is negligible compared to the overall storage. In terms of costs, although the single-step methods need to update the intermediate variables additionally, the required efforts of vector addition and subtraction operations are also negligible compared to the matrix operations per step, especially for large finite element systems. Therefore, it is reasonable to say that the proposed single-step methods have equivalent computational costs to the existing single-step methods. Due to being free from the need of additional start-up procedures, the proposed single-step method can be regarded as a more convenient alternative of the linear r -step method for applications.

Optimization
In this section, the linear spectral analysis is executed to determine the parameters and optimize the properties. According to Dahlquist's second barrier [11], only the second-order unconditionally stable schemes are concerned here.
Generally, the test equation for linear analysis is q + 2ξωq + ω 2 q = 0 ( 7 ) Table 1 Step-by-step solution using the linear r -step method for solving the linear structural dynamics problems A. Initial calculations 1. Form mass matrix M, stiffness matrix K and damping constant C; 2. Initialize q 0 ,q 0 andq 0 ; 3. Select time step t, the algorithmic parameters α 1 , α 2 , · · · , α r , β 0 , β 1 , β 2 , · · · , β r ; 4. Calculate integration constants: B. For each time step The first r − 1 steps (k < r ) 1. Calculate effective loads at time t k : 2. Solve for displacement at time t k : 3. Calculate velocity and acceleration at time t k : From the r th step (k ≥ r ) 1. Calculate effective loads at time t k : Solve for displacement at time t k : 3. Calculate velocity and acceleration at time t k : where ξ is the damping ratio, and ω is the natural frequency. It can be equivalently reduced to the first-order equatioṅ Decomposing the coefficient matrix in Eq. (8) gives the simplified first-order equatioṅ where i = √ −1. To simplify the analysis, Eq. (9) is used as the test equation here.
Applied to Eq. (9), the recurrence equation of the linear r -step method, shown in Eq. (1b), is The characteristic polynomial can be written as where the generic characteristic root is denoted by μ. According to Eq. (11), this polynomial has r characteristic roots μ 1 , μ 2 , · · · , μ r . One of them, called the principal root, dominates the quality of numerical solutions; the other roots are spurious.
Similarly, applying the single-step method of Eq. (5) to the test equation (9) yields the compact scheme where A ∈ C r ×r is the amplification matrix, which governs the stability and accuracy characteristics of a method. From Eq. (12), the difference scheme only with respect to x can be written as where A j ( j = 1, · · · , r ) is the sum of the jth-order principal minors of A. The characteristic polynomial of A is Equation (14) shows a form similar to that of Eq. (11), and they can be identical by choosing appropriate values for the Table 2 Step-by-step solution using the single-step method equivalent to the linear r -step method for solving the linear structural dynamics problems A. Initial calculations 1. Form mass matrix M, stiffness matrix K and damping constant C; 2. Initialize q 0 ,q 0 ,q 0 ,q 1 0 ,q 1 0 ,q 2 0 ,q 2 0 · · · ,q r −1 0 ,q r −1 0 ; 3. Select time step t, the algorithmic parameters γ 0 , γ 1 , γ 2 , · · · , γ 2r −3 , γ 2r −2 ; 4. Calculate integration constants: B. For each time step 1. Calculate effective loads at time t k : 3. Calculate velocity and acceleration at time t k : respective parameters.Then these two methods can be called spectrally equivalent to each other. As representatives of the proposed formulation, the linear two-, three-, and four-step methods, respectively referred to as LMS2, LMS3 and LMS4 in the following, and the alternative single-step methods, correspondingly named SS2, SS3 and SS4, are discussed here. After determining the parameters of the linear multi-step methods, the parameters of the single-step methods are directly obtained by comparing the respective characteristic polynomials. To our knowledge, apart from LMS2 [31,44], the second-order unconditionally stable schemes of the higher-step methods have not been proposed yet.

LMS2 and SS2
Setting r = 2 in Eq. (1b) yields the formulation of LMS2 as It is second-order accurate and unconditionally stable with the optimal parameters The detailed process of parameter optimization is presented in [44]. In Eq. (16), all parameters are controlled by a single parameter ρ ∞ , which is the asymptotic spectral radius, defined as max{|μ j |, j = 1, 2, · · · , r } when |λ t| → +∞. Generally, ρ ∞ is used to denote the degree of algorithmic dissipation. A smaller ρ ∞ means stronger algorithmic dissipation for the high-frequency content. As shown in Eq. (5), the corresponding single-step scheme SS2 has the form From Eq. (14), SS2's characteristic polynomial is where Meanwhile, the characteristic polynomial of LMS2 can be obtained by Eq. (11) and the parameters in Eq. (16) as Equation (18) is the same as Eq. (20) when the parameters satisfy Therefore, SS2 with the parameters in Eq. (21) is spectrally equivalent to LMS2 with the parameters in Eq. (16), so the single-step method SS2 also has second-order accuracy, unconditional stability and tunable algorithmic dissipation.

LMS3 and SS3
According to Eq. (1b), LMS3 can be formulated as This method is also expected to have second-order accuracy, unconditional stability, and to offer tunable algorithmic dissipation. With respect to Eq. (10), the local truncation error is defined as Expanding σ at t k gives where so this method is second-order accurate if To maximize the low-frequency accuracy for a given degree of algorithmic dissipation, the characteristic roots of Eq. (11) when |λ t| → +∞ need to be equal [9], as −ρ ∞ , which implies that Equations (26) and (27) impose six conditions, so only one parameter, e.g. β 0 , is now adjustable.
To evaluate the stability, μ = (1 + z)/(1 − z) is assumed; the stability condition |μ| ≤ 1 is equivalent to Re(z) ≤ 1. With this assumption, the characteristic polynomial of Eq. (11) can be transformed into where the real and imaginary parts of the coefficients, denoted by p j and q j , respectively, are determined by the scalar parameters and λ t. When r = 3, the generalized Routh array [41] with respect to Eq. (28) is listed in Table 3.
According to the generalized Routh stability criterion [41], the roots of Eq. (28) satisfy Re(z) ≤ 1 only when Then this method is unconditionally stable if the constraints in Eq. (29) hold for any λ in the complex left-half-plane, including the imaginary axis. To simplify the analysis, the undamped case, i.e. ξ = 0, which implies λ = ±ωi, is considered first. The desirable range of β 0 that makes this method unconditionally stable can be obtained as To determine β 0 uniquely, the global error, defined as E k = x(t k ) − x k , is checked. Subtracting the recursive scheme, as in Eq. (10), from the right-hand side of Eq. (23), yields so the global error can be measured by the constant The optimal β 0 that minimizes EC is finally obtained as Then, the other parameters can be obtained from Eqs. (26)- (27).
With the set of parameters obtained above, the characteristic polynomial of Eq. (11) can be written as The unconditional stability for the damped cases is checked by solving Eq. (35) with different values of ξ . As expected, the condition |μ| ≤ 1 is satisfied for any ξ > 0 and ω t > 0, so this method is unconditionally stable. For illustration, the absolute values of the characteristic roots |μ| as a function of t/T , where T = 2π/ω, for the cases ξ = 0.0 and ξ = 0.1 are plotted in Figs. 1 and 2 , respectively.
Consequently, the optimal parameters of LMS3 have been obtained, as shown in Eqs. (26), (27) and (34). The resulting scheme has optimal accuracy, unconditional stability and tunable algorithmic dissipation.
According to Eq. (5), the single-step scheme SS3 is Its characteristic polynomial is where By comparing Eqs. (37) and (35), SS3 is spectrally equivalent to LMS3 when its parameters are where γ 1 and γ 3 can be solved from the quadratic equation Note that if ρ ∞ < 1, γ 1 and γ 3 are complex conjugate numbers. As such, SS3 yields complex intermediate variables, nonetheless resulting in real-valued x k andẋ k , as explained in the following. Asẋ 1 0 =ẋ 2 0 =ẋ 0 is assumed, the recurrence equations of the first few steps after eliminating the intermediate variables can be written as where a 11 = a 22 = a 33 = a 44 (41a) 21 (41d) The parameters in Eq. (41) are all real numbers, so the computed state vectors are also real. From Eq. (40), the schemes of the first two steps reduce to first-order accurate, and from the third step on, SS3 shows the same properties of LMS3. Actually, it can be verified that x 3 and x 4 in Eq. (40) can be reorganized as The same formulation can also be found in the fifth step, the sixth step and so on, which further validates the equivalence of SS3 and LMS3. In another way, SS3 can be seen as LMS3 with the specific start-up procedure shown in Eq. (40).

LMS4 and SS4
Similar to LMS2 and LMS3, LMS4 is cast as Its local truncation error is which can be also expanded in the form of Eq. (24) with the coefficients This method is also second-order accurate if Likewise, the requirement of tunable algorithmic dissipation imposes four conditions, as Then two parameters remain free, selected as β 0 and α 1 .
To investigate the stability, the generalized Routh array for Eq. (28) with r = 4 is listed in Table 4.
By usingẋ 1 0 =ẋ 2 0 =ẋ 3 0 =ẋ 0 , the first several steps without the intermediate variables can be also written as where a 11 = a 22 = a 33 = a 44 = β 0 (57a) One can check that the first three steps only have first-order accuracy, but from the fourth step on, the same formulation as LMS4 can be reproduced. Therefore, SS4 can be also regarded as LMS4 with the specific start-up procedure shown in Eq. (56).
In conclusion, the parameters of LMS2, LMS3, LMS4, SS2, SS3, and SS4 for optimal accuracy, unconditional stability and tunable algorithmic dissipation have been determined in this section. From the parameters of LMS2, LMS3 and LMS4, we can see that when 0 ≤ ρ ∞ < 1, β 0 > 1/2 and when ρ ∞ = 1, β 0 = 1/2 , so the single-step scheme used as the start-up procedure in Table 1 is first-order unconditionally stable as 0 ≤ ρ ∞ < 1, and recovers the trapezoidal rule as ρ ∞ = 1. As opposed to LMS2, other schemes are here presented for the first time. The single-step schemes SS2, SS3 and SS4 can present properties identical to those of the corresponding multi-step method, except for the first few steps.

Accuracy and stability
Due to the equivalence of single-step and linear multi-step methods, the proposed single-step methods are not considered in this section. Firstly, the error constants (ECs), which are employed to measure the global errors, as in Eqs. (33) and (49), of LMS2, LMS3 and LMS4 are plotted in Fig. 5. As can be seen, LMS4 gives the smallest global error, followed by LMS3, and LMS2, which means that the accuracy can be improved by using the information of more previous steps. As ρ ∞ increases, the gap between them gradually decreases; when ρ ∞ = 1, these schemes all become the equivalent multi-step expressions of the trapezoidal rule.
Considering the test equation (9), the exact solution has the form Using the linear multi-step method, the numerical solution can be expressed as where μ p and μ s respectively denote the principal and spurious roots. As shown in Figs. 1, 2, 3 to the first r − 1 steps of the recurrence equation of the linear multi-step method with ρ ∞ = 1, as shown in Eq. (10), yields that It follows that, for this case, the large spurious roots do not contribute to the solution, that is c 1 = c 2 = · · · = c r −1 = 0 in Eq. (59). Therefore, only the contribution of the principal root μ p is considered in the following.
Representing μ p = e −ξ ±i 1−ξ 2 ω t , the numerical solution can be expressed in a form similar to that of the exact solution, as where ξ and ω are the damping ratio and natural frequency of the numerical solution, and they can be obtained from the norm |μ p | and phase μ p as By comparing Eqs. (58) and (62), ξ and (T − T )/T (T = 2π/ω), known as the amplitude decay ratio (AD) and period elongation ratio (PE), can be used to measure the amplitude and period accuracy, respectively. Figures 6, 7 and 8 display the spectral radii (SR=|μ p |), AD(%) and PE(%), respectively, of LMS2, LMS3, LMS4 and two representative implicit methods, the generalized-α method and the ρ ∞ -Bathe method, referred to as G-α and Bathe respectively, for the undamped case. Because of the two-sub-step scheme, the computational cost required by the Bathe method per step is twice that of other methods approximately, so the abscissa is set as t/(nT ) (n = 2 for Bathe and n = 1 for other methods) to compare the methods under equivalent computational costs. In addition to ρ ∞ , the ρ ∞ -Bathe method can also change the spectral characteristics by adjusting the splitting ratio γ , so three cases of γ = γ 0 , γ = γ 0 − 0.2 and γ = γ 0 + 1 are considered in these figures, where the value of γ 0 can refer to [35]. When γ = γ 0 , its two sub-steps share the same effective stiffness matrix for linear problems.
As illustrated in Fig. 6, all methods can provide tunable algorithmic dissipation in the high-frequency domain by adjusting the parameter ρ ∞ . For a specified 0 ≤ ρ ∞ < 1, LMS4 shows the spectral radius closest to 1 in the lowfrequency region, so its dissipation error is the smallest, followed by LMS3, LMS2 and G-α. With a small ρ ∞ , the Bathe method with γ = γ 0 and γ = γ 0 −0.2 have the spectral radii close to that of LMS3 in the low-frequency region, but when 0.6 ≤ ρ ∞ < 1, they have larger dissipation errors than G-α. The Bathe method with γ = γ 0 + 1 shows the largest dissipation error for any 0 ≤ ρ ∞ < 1. When ρ ∞ = 1, all methods, except the Bathe method with γ = γ 0 − 0.2 and γ = γ 0 + 1, present the same properties as the trapezoidal rule. Figures 7 and 8 focus on the low-frequency behaviour; they display results similar to those of Fig. 6. LMS4 is the most accurate among these methods in terms of both amplitude and period, followed by LMS3, LMS2, and G-α. The Bathe method with γ = γ 0 + 1 shows the largest errors, and when ρ ∞ is close to 1, the Bathe method with γ = γ 0 and γ = γ 0 − 0.2 are less accurate than G-α. From [35], for a given ρ ∞ and γ ∈ (0, 1), the Bathe method with γ = γ 0 possesses the largest amplitude decay and the smallest period error, which also can be observed by comparing the cases of γ = γ 0 and γ = γ 0 − 0.2.
Owing to the higher accuracy, the proposed LMS3 and LMS4 are really better candidates than the existing secondorder methods. Since the developed single-step schemes, SS2, SS3 and SS4, are spectrally equivalent to LMS2, LMS3 and LMS4, respectively, SS3 and SS4 also have the same accuracy advantages.

Overshoot
From [14], the methods of unconditional stability are possible to show pathological growth at the initial steps, that is the so-called "overshoot". For a convergent method, there is no overshoot as ω t → 0, so generally, only the case of ω t → +∞ needs to be concerned. With the assumption, the recursive schemes of the linear r -step method, applied to the test modelq + ω 2 q = 0, can be written as Since the overshooting phenomenon only occurs at the first few steps, the numerical results of the r th step are considered the approximate results at step r of LMS2, LMS3 and LMS4 can be obtained as respectively, where the parameters are represented by ρ ∞ .
The results indicate that LMS2, LMS3 and LMS4 have zeroorder overshoot in both displacement and velocity. If the trapezoidal rule is used to solve the first r − 1 steps, that is the approximate results at step r of LMS2, LMS3 and LMS4 have the form as Also, no overshoot can be observed from the results. Therefore, as long as the method used for the start-up procedure does not have overshoot, the proposed multi-step methods also exhibit no overshoot from the r th step.
For the proposed single-step methods, SS2, SS3 and SS4, the results of the first step, applied to the modelq +ω 2 q = 0, are checked. With the assumption ω t → +∞, they can be expressed as respectively. The numerical results indicate that the proposed single-step methods also have zero-order overshoot in both displacement and velocity.
In addition, the numerical example used in [14] to measure the overshooting behaviour is also considered here, which has the form as With t/T = 10, the numerical displacements and velocities given by the proposed multi-step methods, single-step methods and the multi-step methods with the trapezoidal rule (TR) as the start-up procedure, are summarized in Figs. 9 and 10 .
From the results, when ρ ∞ < 1, the displacements show a decaying trend from the beginning; the velocities also begin to decay after the initial disturbance caused by the non-zero initial displacement. When ρ ∞ = 1, both the displacements and velocities present periodic increases and decreases, since no algorithmic dissipation exists. In all the results, no pathological growth can be observed, which further validates that the proposed methods have no overshoot in both displacement and velocity.

Illustrative solutions
To check the performance of the proposed LMS2, LMS3, LMS4, SS2, SS3, and SS4, several illustrative examples are simulated in this section, and the G-α and Bathe method are also employed for comparison. In order to compare under the close computational costs, the step size of the two-substep Bathe method is set as twice that of the other methods. The G-α developed in [1], which improves the acceleration accuracy to second-order, is used in the computations. Unless otherwise specified, the single-step method in Eq. (4) is used as the start-up procedure for the linear multi-step methods. Single degree-of-freedom example. The single-degree-offreedom linear equation [22] q + 2ξωq + ω 2 q = R 1 sin(ω 1 t) + R 2 cos(ω 2 t) where ξ = 0.1, ω = 2π, R 1 = 10, R 2 = 15, is considered to check the convergence rate. The global errors (GE), defined as within [0, 10] versus t/n for these methods using ρ ∞ = 0 are plotted in Fig. 11. Second-order convergence rate can be clearly observed in all methods. The Bathe method with γ = γ 0 + 1 shows significantly larger errors than others. Besides, the absolute errors (AE), calculated by with t = 0.01 for the case ρ ∞ = 0 and ρ ∞ = 0.6 are plotted in Figs. 12 and 13, respectively. In both cases, LMS4 and SS4 predict the most accurate solutions, followed by LMS3 and SS3, LMS2 and SS2, consistent with the discussion in Sec. 4. If the Bathe method with γ = γ 0 +1 is excluded, G-α shows the largest errors when ρ ∞ = 0, and the Bathe method with γ = γ 0 − 0.2 is the most inaccurate when ρ ∞ = 0.6. The Bathe method with γ = γ 0 shows higher accuracy than the method with γ = γ 0 − 0.2, so only this scheme is considered in the following examples. To exhibit the effect of different start-up procedures on the accuracy, Fig. 14 plots the absolute errors given by LMS4 fol-   Table 1, SS4, and LMS4 using the trapezoidal rule for the first three steps with t = 0.01. In addition to the trapezoidal rule, the other two start-up procedures are first-order accurate. From the results, the accuracy difference can be observed at the beginning, but gradually disappears over time. The single-step method SS4 has very close accuracy to LMS4 with the trapezoidal rule as the start-up procedure in the whole time. Therefore, the start-up procedure has a very limited impact on the overall accuracy. For linear systems, the procedure shown in Table 1 and the single-step method are more recommended, since they only need one factorization of the effective stiffness matrix. Stiff-soft mass-spring system. The mass-spring oscillator [5], shown in Fig. 15, is considered as a model problem to test the dissipation ability of a method. The motion of m 0 is prescribed as q 0 = sin(1.2t), and other parameters are k 1 = 10 7 , k 2 = 1, m 0 = 0, m 1 = m 2 = 1. The governing equation of m 1 and m 2 can be written as With zero initial conditions and t/n = 0.1309, the numerical solutions of m 1 by these methods using ρ ∞ = 0 are summarized in Fig. 16. The results of q 1 are overlapping with each other, and fromq 1 (v 1 ) andq 1 (a 1 ), these methods all can filter out the rigid component in the first several steps. Because the first r − 1 steps of the linear r -step method and the corresponding single-step method cannot offer the prescribed degree of algorithmic dissipation, the proposed methods need to experience more steps of high-frequency oscillations. Fig. 17 Acceleration of m 1 given by LMS4 with different start-up procedures Figure 17 plots the acceleration of m 1 given by LMS4 with different start-up procedures. One can see that SS4 shows stronger dissipation in the first three steps, and LMS4 using the trapezoidal rule as the start-up procedure exhibits larger oscillation due to being non-dissipative. At the sixth step, the stiff component has been filtered out in all schemes. Wave propagation in a clamped-free bar. As shown in Fig. 18, the clamped-free bar [18,35], excited by the load F at the free end, is considered. The system parameters are the elastic modulus E = 3 × 10 7 , cross-sectional area A = 1, density ρ = 7.3 × 10 −4 , length L = 200, and the step load F = 10 4 . The bar is meshed by 1000 linear two-node finite elements.
With different CFL numbers (the ratio of the propagation length per step √ E/ρ t to the element length L/1000), the velocity at the midpoint x = 100 given by these methods using ρ ∞ = 0 are plotted in Figs. [19][20][21][22][23]. From the results, the spurious oscillation can be observed in all methods with CFL=1.0 (CFL=2.0 for the Bathe method), and it lasts longer for LMS3, SS3, LMS4 and SS4. However, with a smaller CFL number, which is different for each method, all methods can rapidly filter out the high-frequency oscillation. Among these methods, LMS4 and SS4 provide favorable results with   a larger CFL number, so they require fewer steps than other methods. A two-dimensional wave propagation problem. The two-dimensional wave propagation problem, a pre-stressed membrane excited by the concentrated load at the center [20,28], as shown in Fig. 24, is considered. Its transverse governing equation can be expressed as   where the exact wave velocity c 0 is assumed as 1. The concentrated load has the form as With zero initial conditions, the displacements at t ≈ 13 given by these methods with ρ ∞ = 0 and the CFL number suggested from the last example are shown in Figs. 25-27, where the membrane is meshed by 70 × 70, 105 × 105 and 140 × 140 elements, respectively. By the mesh of 140 × 140 elements, Fig. 28 plots the obtained displacements and velocities at the angle of θ = 0 and θ = π/4, where θ denotes the angle between the wave propagation direction and the x-axis. Except for the trapezoidal rule, other methods all can predict reasonably accurate solutions. The linear multi-step method shows almost identical results to its corresponding single-step method. From the above two examples, the proposed methods perform very well in the wave propagation problems with an appropriate CFL number. A more detailed analysis of the proposed methods for wave propagation problems, such as the optimal CFL numbers for different ρ ∞ , other advantages and disadvantages, is still deserved in the future. Spring-pendulum model. The spring-pendulum model [10] shown in Fig. 29 is solved. Its motion equation is expressed  With t/n = 0.01 s, Figs. 30-31 show the numerical results of r and θ , as well as their absolute errors for the compliant case, where ρ ∞ = 0 is adopted in all methods, and the reference solution is obtained using the Bathe method with t = 10 −5 s. For the nonlinear problem, LMS4 and SS4 still provide the most accurate results, and G-α yields the largest errors. The single-step schemes also show close accuracy to the corresponding multi-step methods. Figure 32 plots the computed r and θ for the stiff case by using t/n = 0.01 s and ρ ∞ = 0. Since the algorithmic dissipation of the proposed methods works from the r th step, LMS4 and SS4 needs more steps to damp out the stiff component. From about the sixth step, all methods can give reasonable results for r . Truss element. In this example, the classic pendulum model [27], as shown in Fig. 33, is considered. It is discretized as 1 truss element considering the Green strain. This example is usually used to check the energy-conserving characteristic of a method for nonlinear structural dynamics. The system parameters are the length L = 3.0443 m, density per unit length ρ A = 6.57 kg/m, axial stiffness E A = 10 4 N, and 10 10 N for the compliant and stiff cases. The pendulum is excited by the horizontal velocity v 0 = 7.72 m/s at the free end. The gravity is not considered, so the pendulum makes continuous rotations in this plane, along with axial straining.  Figure 34 displays the axial deformation L and total energy for the compliant case. Significant amplitude and energy decay can be observed with the G-α, LMS2 and SS2 methods. LMS4 and SS4 exhibit better energy-conserving properties for the nonlinear system owing to their higher accuracy. Figure 35 gives the axial deformation L and rotation angle θ for the stiff case. Plots similar to those of Fig. 32 can be observed. All these methods can effectively filter out the stiff component after a few oscillations.

Conclusions
In this paper, second-order unconditionally stable schemes of linear multi-step methods, represented by LMS2, LMS3, and LMS4, are discussed. From linear analysis, the parameters, controlled by ρ ∞ , of LMS3 and LMS4 are given for optimal accuracy, unconditional stability and controllable algorithmic dissipation. In addition, single-step methods, spectrally equivalent to the linear multi-step ones, are proposed by using intermediate state variables. The parameters of the latter methods, named SS2, SS3, and SS4, are also obtained by enabling them to have the same characteristic polynomials as LMS2, LMS3, and LMS4, respectively. The resulting schemes show equivalent formulations of the corresponding linear multi-step methods, except for the first few steps, in which the linear multi-step methods require a dedicated startup procedure. Being self-starting, the equivalent single-step schemes are more convenient to use.
Compared with some representative second-order methods, including LMS2, G-α and the Bathe method, the newly developed LMS3 and LMS4 show higher amplitude and period accuracy for a given degree of algorithmic dissipation. Among them, LMS4 is the most accurate, followed

Compliance with ethical standards
Conflicts of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.