On the optimization of n-sub-step composite time integration methods

A family of n-sub-step composite time integration methods, which employs the trapezoidal rule in the first n-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n-1$$\end{document} sub-steps and a general formula in the last one, is discussed in this paper. A universal approach to optimize the parameters is provided for any cases of n≥2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\ge 2$$\end{document}, and two optimal sub-families of the method are given for different purposes. From linear analysis, the first sub-family can achieve nth-order accuracy and unconditional stability with controllable algorithmic dissipation, so it is recommended for high-accuracy purposes. The second sub-family has second-order accuracy, unconditional stability with controllable algorithmic dissipation, and it is designed for heuristic energy-conserving purposes, by preserving as much low-frequency content as possible. Finally, some illustrative examples are solved to check the performance in linear and nonlinear systems.


Introduction
Direct time integration methods are frequently used to predict accurate numerical responses for general dynamic problems after spatial discretization. Driven by the pursuit of desirable properties, including higher accuracy and efficiency, robust stability, and many others, a number of excellent methods were proposed in the past decades.
In terms of the formulations, existing methods are generally classified into explicit and implicit schemes. Explicit methods are mostly used in wave propagation problems, as their conditional stability limits the allowable time step size to the highest system frequency. Implicit methods have fewer restrictions on the problems to be solved due to the unconditional stability, but they require more computational efforts per step.
In another way, the integration methods can also be divided into single-step, multi-sub-step and multistep techniques. The single-step methods only adopt the states of the last step to predict the current one, while the multi-sub-step methods also need the states at the intermediate collocation points, and the multistep methods require the states of more than one previous step. Each of them has specific advantages and disadvantages.
From the literature, representative single-step methods include the Newmark method [25], the HHT-α method (by Hilbert, Hughes, and Taylor) [17], the WBZ-α method (by Wood, Bossak, and Zienkiewicz) [29], the generalized-α method [9], the GSSSS (gener-alized single-step single-solve) method [34], and many others [28]. These single-step methods were proved to be spectrally identical to the linear multi-step methods [34], so they suffer from the Dahlquist's barrier [10], which states that the methods of higher than secondorder accuracy cannot be unconditionally stable. Therefore, the methods mentioned above are all second-order accurate and unconditionally stable; some of them can also provide controllable algorithmic dissipation.
In the multi-step class, the Dahlquist's barrier certainly works, but in terms of accuracy, the linear twostep method [24,33] is superior to most existing singlestep methods under the same degree of algorithmic dissipation. In this class, BDFs (backward differentiation formulas) [11,16] also represent a widely-used branch, particularly useful for stiff problems owing to the strong algorithmic dissipation. These popular multi-step methods are also second-order accurate and unconditionally stable. However, the multi-step methods are not self-starting, so another method has to be also used to solve the initial steps, which makes the multi-step methods not as convenient to use as the single-step ones.
The multi-sub-step methods, also known as multistage methods, allow more possibilities in terms of properties. The most representative method is the famous Runge-Kutta family [6,7,19], which can be designed to be arbitrarily higher-order accurate and unconditionally stable by choosing proper parameters and enough stages. Besides, Fung [12][13][14][15] provided some methods to reproduce the generalized Padé approximation. These methods can reach up to 2nthorder accuracy by employing n sampling grid points per step, but the dimension of the implicit equation to be solved is n times that of the original, resulting in huge computational costs. In the multi-sub-step class, the composite methods [3], which divide each step into several sub-steps and employ different methods in each sub-step, have received a lot of attention in recent years.
Based on Bank et al.'s work [1], Bathe et al. [3] introduced the concept of the n-sub-step composite method by utilizing the trapezoidal rule in the first n − 1 sub-steps and the (n + 1)-point backward difference scheme at the end of the step. The two-substep scheme is known as the Bathe method, which is asymptotically stable with second-order accuracy. Thanks to its strong dissipation and preferable accuracy, the Bathe method has been found to perform well in many fields [2,4,27]. The three-, and four-sub-step composite methods [8,32], which are asymptotically stable with higher accuracy, were also developed adopting the similar idea. Furthermore, to acquire controllable algorithmic dissipation, the two-sub-step methods [20,21,26], and the controllable three-substep methods [18,23], were proposed by replacing the backward difference scheme with a more general formula. However, with the increase in the number of sub-steps, the number of scalar parameters required to be designed also increases, so the basic requirements, including second-order accuracy, unconditional stability, controllable algorithmic dissipation, are not enough to determine these parameters uniquely. Two optimal sub-families of the controllable three-sub-step method were proposed in [23], since different conditions are considered as a supplement.
On this basis, this paper purposes to provide a universal approach to optimize the parameters of generalized n-sub-step composite method, where n can be any integer greater than 2, and the trapezoidal rule is employed in the first n − 1 sub-steps. Two kinds of optimization goals are considered, producing two optimal sub-families for different purposes. The first one intends to achieve higher-order accuracy, under the premises of unconditional stability and controllable algorithmic dissipation. The second one is dedicated to conserving low-frequency behavior, while still providing controllable high-frequency dissipation. From linear analysis, the resulting schemes in the first subfamily can reach up to nth-order accuracy by using n sub-steps, and the schemes in the second sub-family exhibit very small algorithmic dissipation in the lowfrequency domain. Most of these schemes are developed for the first time, and in each sub-family, the accuracy can be improved by using more sub-steps. Finally, the proposed methods are applied to solve several numerical examples to check the performance. This paper is organized as follows. The formulations of the n-sub-step composite method are shown in Sect. 2. The optimization of the parameters is implemented in Sect. 3. The detailed properties of the two sub-families are discussed in Sect. 4. Numerical examples are provided in Sect. 5, and conclusions are drawn in Sect. 6.

Formulation
In the literature, the composite methods were mostly developed to solve the problems in structural dynamics, as where M is the mass matrix, F collects the damping force, internal force and external load, x,ẋ and x are the displacement, velocity and acceleration vectors, respectively, t is the time, and t 0 , x 0 and v 0 are the given initial time, displacement and velocity, respectively. When this method is applied using n sub-steps, it can be formulated as and denotes the numerical solution at collocation points, h is the step size, and γ , q 0 , q 1 , · · · , q n are the control parameters. The current step [t k , t k + h] is divided into n sub-steps: t k , t k + 2γ h , t k + 2γ h, t k + 4γ h ,· · · , t k + 2(n − 2)γ h, t k + 2(n − 1)γ h , and [t k + 2(n − 1)γ h, t k + h . In the first n − 1 sub-steps, the trapezoidal rule is adopted. In the last one, a general formula containing information about all collocation points is utilized. The present formulation can reduce to the ρ ∞ -Bathe method [26] when n = 2 and to the three-sub-step method [18,23] when n = 3.
In this method, because the same form of assumptions is used to solve x k+1 andẋ k+1 , Eqs. (2) and (3) can be reformulated based on the general first-order differential equation f ( y,ẏ, t) = 0, as f y k+2 jγ ,ẏ k+2 jγ , t k + 2 jγ h = 0 (4a) where {x;ẋ} is replaced by y, and the dynamics equations can be equivalently formulated as first-order differential equations by adding the trivial equatioṅ x =ẋ. Equations (4) and (5) present more general formulations for solving first-order and arbitrarily higherorder differential equations. However, for solving the second-order dynamic problems, Eqs. (2) and (3) are still more recommended, since in the equivalent firstorder expressions, the number of implicit equations to be solved doubles. From the formulation, the first n − 1 sub-steps can share the same procedure in a loop, whereas the last sub-step needs to be implemented separately. The assumption q n = γ is introduced here, which imposes that the last sub-step shares the same form of Jacobi matrix as the first n − 1 sub-steps. This assumption is particularly useful when applied to linear problems, since it allows the constant Jacobi matrix to be factorized only once, like in the single-step methods. For applications, Table 1 shows the computational procedures of the n-sub-step composite method for the general first-order differential equation f ( y,ẏ, t) = 0, where the Newton-Raphson iteration is utilized to solve the nonlinear equation per sub-step.
Besides, by reorganizing the formulations, the composite method can be regarded as a special case of the diagonally-implicit Runge-Kutta methods (DIRKs) with the explicit first-stage. The corresponding Butcher's tableau [6] has the form as Table 1 Computational procedure of the n-sub-step composite method for solving f ( y,ẏ, t) = 0, y(t 0 ) = y 0 A. Initial calculations 1. From the function f ( y,ẏ, t) and its derivative functions with respect to y andẏ, as f y and f˙y, respectively; 2. Initialize t 0 , y 0 andẏ 0 ; 3. Select the time step size h, the algorithmic parameters γ, q 0 , q 1 , · · · , q n−1 , the tolerance error , and the maximum number of iterations N ; 4. Calculate the constant: a = 1 γ h . B. For each time step 1. The first n − 1 sub-steps For j = 1, j < n, j + +: a. Predict y k+2 jγ andẏ k+2 jγ : Prepare the matrices: c. Update y k+2 jγ andẏ k+2 jγ : End.

Optimization
In linear spectral analysis, owing to the mode superposition principle, it is common and enough to consider the single degree-of-freedom equation where ξ is the damping ratio, and ω is the natural frequency. To simplify the analysis, the equivalent first-order differential equation is discussed, aṡ Decomposing the coefficient matrix in Eq. (7) yields the simplified first-order equatioṅ When the composite method is applied, the recursive scheme becomes where the amplification factor A is Since q n = γ is assumed in Sect. 2, Eq. (10) is updated as and a n = γ n−1 For example, n = 5 follows Consequently, the parameters under analysis change from q j ( j = 0, 1, · · · , n − 1) and γ , to a p ( p = 1, 2, · · · , n) and γ in the following. When a p and γ are given, the parameters q j can be obtained uniquely by solving Eqs. (12) and (13). For applications, Table 2 shows Table 2 Formulas of q j ( j = 0, 1, · · · , n − 1) the formulas of q j expressed by a p and γ for the cases n = 2, 3, 4, 5.

Higher-order schemes
A numerical method is naturally expected to be as accurate as possible, so the higher-order schemes are considered first. From the scheme of Eq. (9), the composite method uses the amplification factor A, rewritten as 1 + a 1 z + a 2 z 2 + · · · + a n z n (1 − γ z) n (15) to approximate the exact amplification factorÂ Hence the local truncation error σ can be defined as If σ = O z s+1 , the method is said to be sth-order accurate, which requires that up to sth derivatives of A at z = 0 are all equal to 1, that is To satisfy Eq. (18), a p ( p = 1, 2, · · · , n) can be solved as Therefore, if all a p ( p = 1, 2, · · · , n) follow the relationships in Eq. (19), this method can achieve nth-order accuracy, and then γ becomes the only free parameter to control the stability. A time integration method is said to be uncondi- According to Ref. [19], the bounds on γ can be given by considering the stability on the imaginary axis (ξ = 0), which can result in the unconditional stability when the accuracy order s = n in the DIRKs. Therefore, let z = ±iτ where τ = ωh is a real number, and N (z) = 1 + a 1 z + a 2 z 2 + · · · + a n z n (20a) which are the numerator and denominator of A(z) in Eq. (15), respectively, |A(z)| ≤ 1 is equivalent to Then the condition for unconditional stability can be transformed into where the function S(τ ) is introduced, and the coefficients c 2 j ( j = 0, 1, 2, · · · , n) are expressed as in which a 0 is set to 1. By Eq. (22), the bounds on γ of the cases n = 2, 3, 4, 5 are provided in Table 3 .
It follows that, with s = n, the allowable range of γ narrows as n increases and, in some cases, the n-substep method can achieve (n + 1)th-order accuracy with a fixed γ . Besides, algorithmic dissipation is also a desirable property for a time integration method, to filter out the inaccurate high-frequency content. Generally, it is measured by the spectral radius ρ ∞ at high-frequency limit, that is and it gets stronger with a smaller ρ ∞ . With A(z) from Eq. (15), Eq. (24) can be satisfied if which can be used to solve γ for a given ρ ∞ . Table 4 shows the solutions of γ for several specific ρ ∞ in  Table 3, is selected. So far, the unconditionally stable higher-order accurate schemes with controllable algorithmic dissipation have been developed, whose parameter γ can be solved for a given ρ ∞ by Eq. (25), a p ( p = 1, 2, · · · , n) are determined by γ as shown in Eq. (19), and then q j ( j = 0, 1, 2, · · · , n − 1) can be obtained by solving Eqs. (12) and (13). These information for the cases n = 2, 3, 4, 5 are shown in Tables 2, 3 and 4. The special case n = 2 is identical to the ρ ∞ -Bathe method [26], whereas the other cases are presented here for the first time. In addition, the accuracy and algorithmic dissipation are discussed in more detail in Sect. 4.

Conserving schemes
An original intention of the composite methods was to conserve the energy of the system [2], which explains why the trapezoidal rule is utilized in most sub-steps. Existing two-and three-sub-step methods [18,22,26] really show preferable energy-conserving characteristic over other single-and multi-step methods. In this work, a simple and general approach to determine the parameters, which enable the n-sub-step composite method to conserve as much low-frequency content as possible, is proposed.
First of all, to be competitive, the method needs to have some basically useful properties, including at least second-order accuracy, which requires and controllable algorithmic dissipation, achieved by In addition, unconditional stability also needs to be satisfied, which will be checked last.
Then all parameters of the conserving schemes have been given by combining Eqs. (12), (13), (26), (27) and c 4 = c 6 = · · · = c 2n−2 = 0. The resulting scheme of n = 3 is equivalent to the first sub-family of the threesub-step method proposed in [23]; the other cases are presented here for the first time.
In particular, when ρ ∞ = 1, the resulting scheme is a n-sub-step method with the trapezoidal rule in all substeps, which is supposed to be unconditionally stable in the linear analysis. Empirically, the algorithmic dissipation is acquired by reducing the spectral radius ρ, so the dissipative schemes are likely to be also unconditionally stable, and even present more robust stability. For the undamped case (ξ = 0), the stability can be guaranteed since S(τ ) = (1 − ρ 2 ∞ )γ 2n τ 2n ≥ 0; for other cases, the stability conditions of the schemes listed in Table 5 are checked one by one by considering ξ ∈ (0, 1] and τ ∈ [0, 10000] numerically. As expected, ρ ≤ 1 is satisfied at every point in all schemes, so these methods can be said to possess unconditional stability for linear problems. Other properties are discussed in Sect. 4.

Properties
Two sub-families of the n-sub-step composite method have been presented for different purposes. To identify them, the higher-order schemes are referred to as MSSTH(n), and the conserving schemes are MSSTC(n), where MSST means the multi-sub-step composite method which employs the trapezoidal rule in all substeps except the last one, H and C are utilized to distinguish the two sub-families, and n is the number of sub-steps.
In this section, the representative methods in the literature, including the single-step generalized-α method [9] (G-α) and the linear two-step method [33] (LTS) are also considered for comparison. As the employed methods are all implicit, their computational cost is mainly spent on the iterative calculation when used for nonlinear problems, or the matrix factorization for linear problems. The vector operations brought by the recursive scheme of the method itself is generally considered to have little effect on overall efficiency. Therefore, G-α and LTS are recognized as having equivalent efficiency if the same step size is used. As the composite methods implement a single-step or multi-step scheme in each sub-step, they have the equivalent efficiency to G-α and LTS, if their required number of sub-steps is equal to the number of steps required by G-α and LTS. For this Table 5 a p ( p = 3, 4, · · · , n − 1) and γ for controllable algorithmic dissipation in the conserving schemes reason, to compare the properties under the close computational costs, the same h/n, where n is the number of sub-steps in the composite methods, and n = 1 for the G-α and LTS, is used in these methods. As discussed in Sect. 3, MSSTH(n) has nth-order accuracy under the premises of unconditional stability and controllable algorithmic dissipation. Figures 1 and  2 display the percentage amplitude decay (AD(%)) and period elongation (PE(%)) respectively, of which the definition can refer to [34], of MSSTH(2,3,4,5), G-α and LTS, considering the undamped case (ξ = 0). The abscissa is set as τ/n to compare these methods under the close computational costs.
The results illustrate that the amplitude and period accuracy cannot be improved simultaneously as the order of accuracy increases in MSSTH(n). In terms of amplitude, with a small ρ ∞ , the ρ ∞ -Bathe method (the same as MSSTH (2)) is the most accurate, and when 0.4 < ρ ∞ ≤ 1, LTS shows smaller dissipation error, followed by the G-α and the ρ ∞ -Bathe method. From Fig. 2, MSSTH(3,4,5) have smaller period error than the second-order methods, and MSSTH(5) is the best among them.
In the same way, the percentage amplitude decay and period elongation of MSSTC (2,3,4,5), G-α and LTS for the undamped case are shown in Figs. 3 and 4, respectively. It can be observed that under the similar efficiency, MSSTC(n) presents higher amplitude and period accuracy with a larger n. The gap is more obvious as ρ ∞ decreases, and when ρ ∞ = 1, all the schemes have the same properties as the trapezoidal rule. Both G-α and LTS are less accurate than MSSTC (3,4,5) in the low-frequency range.
As expected, Figs. 5, 6 and 7 demonstrate that MSSTC(n) preserves wider low-frequency range, followed by Padé(n), and MSSTH(n). Note that MSSTH(n) with ρ ∞ = 1 exhibits mild algorithmic dissipation in the medium frequency range, so these schemes are not recommended if all frequencies are requested. Figures 8, 9 and 10 also show that MSSTC(n) has the smallest amplitude dissipation in the low-frequency content. The amplitude decay ratio of MSSTC(5) is very close to 0 over τ ∈ [0, 2]. In terms of period accuracy, Figs. 11, 12 and 13 show that Padé(n) is the most  From the comparison, MSSTC(n) performs really good at conserving the low-frequency content, and its overall accuracy can be improved by using more substeps. MSSTH(n) shows higher period accuracy than the second-order methods, whereas its dissipation error is larger in the low-frequency content.

Numerical examples
To validate the performance, several numerical examples are solved in this section. As the spectral analysis has revealed the properties based on the linear model, this section focuses more on the application and discussion for nonlinear systems.
is considered, and the absolute errors of the displacement x k , velocityẋ k , and accelerationẍ k versus h at t = 10 are plotted in Fig. 14.
The results are consistent with the accuracy order. That is, MSSTC(n) and MSSTH(n) respectively present second-order and nth-order convergence rate. As a result, the higher-order MSSTH(n) enjoys significant accuracy advantage over the second-order methods. However, when h decreases from 10 −2 , it seems that MSSTH(5) cannot maintain fifth-order accuracy. This is because when h is small enough, all effective num-  bers stored in the computer are exactly precise, so if h continues to decrease, the accumulated rounding error can greatly spoil the numerical precision [31]. Van der Pol's equation The van der Pol's equation [19] is solved, where is an adjustable parameter. For the cases = 0.01, 0.001, 0.0001, the absolute errors of x 1,k and x 2,k at t = 1 versus h are plotted in Fig. 15, where the reference solution is obtained by the ρ ∞ -Bathe method with h = 10 −7 . From Fig. 15, in most cases, the second-and nthorder convergence rate can be observed from errors of MSSTC(n) and MSSTH(n), respectively, but for the stiffer case of = 0.0001, MSSTH(3) and MSSTH (5) show obvious order reduction in both x 1 and x 2 . It indi-  cates that the accuracy order also depends on the problem to be solved when applied to nonlinear systems. The order reduction also occurs in other higher-order DIRKs when used for nonlinear problems, see Ref. [19]. Nevertheless, MSSTH(n) still shows significant accuracy advantage over the second-order MSSTC(n) with a small step size.

Multiple degrees-of-freedom examples
In this subsection, some illustrative examples are solved by using the ρ ∞ -Bathe method, MSSTC (3,4,5), MSSTH (3,4,5), G-α and LTS. In these methods, the parameter ρ ∞ is set as 0 uniformly, and the same h/n is used for comparison under close computational costs. The reference solutions are obtained by the ρ ∞ -Bathe method with an extremely small time step.    Fig. 16, the spring-pendulum model, where the spring is fixed at one end and with a mass at the free end, is simulated. Its motion equation can be written as where f (r ) denotes the elastic force of the spring and other system parameters are assumed as m = 1 kg, L 0 = 0.5 m, g = 9.81 m/s 2 . Three kinds of constitutive relations, as Let h/n = 0.01 s; the numerical solutions of E − E 0 (E denotes the system energy and E 0 is the initial value), r and θ for the three cases are summarized in Figs. 17, 18 and 19. From the curves of E − E 0 , it can be observed that MSSTC (3,4,5) can almost preserve the numerical energy from decaying in all cases, despite the oscillations. MSSTH(5) can preserve more energy than the Bathe method, while G-α, LTS, and MSSTH (3,4) show obvious energy-decaying. From the numerical results of r and θ , one can see that with the step size, the numerical solutions of these methods have clearly deviated from the reference solution after a period of simulation. Among these methods, MSSTH(5) predicts the closest solutions to the reference ones, and G-α shows the largest errors. In addition, MSSTC(3,4,5) exhibit good amplitude accuracy thanks to their energy-preserving characteristic. These conclusions are all consistent with the results from linear analysis.
Moreover, to check the algorithmic dissipation, the stiff case, where f (r ) = kr (k = 98.1 × 10 10 N/m), is also simulated with h/n = 0.01 s. The numerical results of E − E 0 , r and θ are plotted in Fig. 20. The results of r indicate that all employed schemes with ρ ∞ = 0 can effectively filter out the stiff component in the first few steps. After the initial decaying, MSSTC (3,4,5) can still preserve the remaining energy in the following simulation. Slider-pendulum model The slider-pendulum model, shown in Fig. 21, is considered in this case. The slider is constrained by the spring, and one end of the pendulum is hinged to the center of mass of the slider. The motion is described by the differential-algebraic equations The system parameters are m 1 = m 2 = 1 kg, L = 1 m, J 2 = 1 12 kg · m 2 , g = 9.81 m/s 2 , k = 1 N/m and 10 10 N/m respectively for the compliant and stiff systems. The slider is excited by the initial horizontal velocity 1 m/s. By using h/n = 0.01 s, the numerical solutions of E − E 0 , x 1 and θ for the compliant and stiff cases are shown in Figs. 22 and 23, respectively. From the results of x 1 and θ , these methods all perform well in terms of accuracy and algorithmic dissipation. However, the numerical energies of MSSTH(4) show a slightly upward trend in the stiff case, so this method cannot give stable results for the problem.
As already discussed in several papers [5,30], the unconditional stability of a time integration method derived from linear analysis cannot be guaranteed when they are applied to nonlinear problems. For nonlinear problems, the stability of a method depends not only on its recursive scheme, but also on the problem itself. Therefore, it is hard to give a definite conclusion about the stability of a method for general problems. From the numerical results, all employed methods, except MSSTH(4), provide stable results when solving stiff problems and differential-algebraic equations, so they can be said to have relatively strong stability.  (4) is not recommended for these problems due to its poorer stability. N -degree-of-freedom mass-spring system The N -degreeof-freedom mass-spring system [23], shown in Fig. 24, is considered to check the computational efficiency. The system parameters are set as With zero initial conditions, three cases, N = 500, 1000 and 1500, are simulated by these methods using h/n = 0.01 s. Figure 25 shows the numerical solutions of x N . It follows that with the step size, all methods can provide reliable results. The CPU time and total number of iterations required by these methods in the simulation of [0, 30 s] are summarized in Table 6. With h/n = 0.01 s, these methods need to proceed 3000 steps (sub-steps for the composite methods) in the whole simulation. Table 6 shows that in addition to MSSTH (4,5), other methods only require one iteration per step/sub-step, so their computational costs are almost equal to each other. One can also see that the required CPU time is approximately proportional to the number of iterations. MSSTH (4,5), especially MSSTH(4), take slightly longer time than other methods.
To check the generality of this conclusion, the required total number of iterations in the above spring-pendulum and slider-pendulum examples are also listed in Table 7. In the two examples, h/n = 0.01 s is adopted, and the simulation of [0, 30 s] is also performed. The results indicate that the total numbers of iterations required by the second-order methods are very close in all cases. Although the higher-order methods need more iterations sometimes, the increased numbers, especially in MSSTH (3,5), are not very large in most cases. Therefore, it is reasonable to say that these methods with the same h/n have similar efficiency for nonlinear problems, and the above comparisons in terms of properties are conducted under the close computational costs.
Overall, the numerical examples in this section demonstrate that when applied to nonlinear problems, the proposed methods can still take advantage of their properties, including the energy-conserving characteristic of MSSTC(n), the high-accuracy of MSSTH(n), and the strong dissipation ability of both sub-families. However, MSSTH(n) shows reduced order and energy instability in some examples. From the presented solutions, MSSTH(5) is more recommended in the higherorder sub-family because of its high accuracy and robust stability, whereas MSSTH(4) is not so prefer-  and energy instability, emerged in MSSTH(n). In this sub-family, MSSTH(5) is more recommended thanks to its high-accuracy and robust stability, and MSSTH (4) is not so preferable, since it shows energy instability and lower efficiency in these examples. However, these conclusions about nonlinear problems are obtained from the existing numerical results. The theoretical analysis is still desired in the future.