A novel explicit three-sub-step time integration method for wave propagation problems

A novel explicit three-sub-step time integration method is proposed. From linear analysis, it is designed to have at least second-order accuracy, tunable stability interval, tunable algorithmic dissipation and no overshooting behaviour. A distinctive feature is that the size of its stability interval can be adjusted to control the properties of the method. With the largest stability interval, the new method has better amplitude accuracy and smaller dispersion error for wave propagation problems, compared with some existing second-order explicit methods, and as the stability interval narrows, it shows improved period accuracy and stronger algorithmic dissipation. By selecting an appropriate stability interval, the proposed method can achieve properties better than or close to existing second-order methods, and by increasing or reducing the stability interval, it can be used with higher efficiency or stronger dissipation. The new method is applied to solve some illustrative wave propagation examples, and its numerical performance is compared with those of several widely used explicit methods.

NB method and single-step explicit methods. As τ b decreases, the proposed method can offer better period accuracy and stronger algorithmic dissipation. A comparative study of the proposed method with NB, EG-α and TW is presented.
Since the explicit methods are commonly used to solve wave propagation problems, the dispersion analysis for these problems is also performed, to support the selection of the parameters, including τ b , ρ b , and the corresponding CFL number. For these problems, high dispersion accuracy is expected to provide accurate solutions, while strong algorithmic dissipation is also important to filter out the inaccurate high-frequency dynamics, which can greatly spoil the overall accuracy as the errors accumulate. For a given ρ b , the proposed method with a selected τ b shows better or close properties compared to several widely used explicit methods. As τ b increases, it has higher efficiency and smaller dispersion error, and as τ b gets smaller, it can offer stronger high-frequency dissipation. Some benchmark problems are simulated to assess the performance of the proposed method.
This paper is organized as follows. The formulations of the proposed method are presented in Sect. 2. Accuracy, stability, overshoot and dispersion analysis are discussed in Sect. 3. Numerical examples are presented in Sect. 4, and conclusions are drawn in Sect. 5.

Formulations
Linear structural dynamics problems have the general form Mẍ + Cẋ + K x = R(t) (1) where M, C and K are the constant mass, damping and stiffness matrices, respectively,ẍ,ẋ and x are the acceleration, velocity and displacement vectors, respectively, and R is the external load vector, usually a function of the time t. The initial displacement x 0 and velocityẋ 0 are the initial conditions of the corresponding Cauchy problem. The proposed method divides each time step, spanning the interval [t, t + t], into three sub-steps, as [t, t + γ 1 t], [t + γ 1 t, t + γ 2 t] and [t + γ 2 t, t + t] (0 ≤ γ 1 ≤ γ 2 ≤ 1). In the first sub-step, the displacement, velocity and acceleration vectors are updated according to Mẍ t+γ 1 t + Cẋ t+γ 1 t + K x t+γ 1 t = R t+γ 1 t (2a) x t+γ 1 t = x t + γ 1 tẋ t + 1 2 γ 2 1 t 2ẍ t (2b) x t+γ 1 t =ẋ t + γ 1 tẍ t (2c) In the second sub-step, we use and in the last one Mẍ t+ t + Cẋ t+ t + K x t+ t = R t+ t (4a) Finally, the velocityẋ t+ t is updated bẏ Here γ 1 , γ 2 , γ 3 , γ 4 , γ 5 , γ 6 , γ 7 , γ 8 , β 1 , β 2 , β 3 are control parameters to be determined. If the external load R is known as a function of time, t, its exact values at the internal points can be used. Otherwise, if it is only defined and sampled at discrete time points, the loads at t + γ 1 t and t + γ 2 t can be obtained for example through linear interpolation, as An approach for selecting the optimal loads at the internal points was given in [22]. From Eqs.
(2)-(5), the proposed method employs a form similar to Taylor expansion at time t to predict the displacement and velocity of each sub-step, and replaces the acceleration term with a linear combination of accelerations at known time points. The method is essentially explicit, with a diagonal mass matrix. By evaluating the internal force and the damping force at the element level prior to assembling, this method can completely avoid matrix operations.
In terms of computational effort, the amount of calculation required for each sub-step of the proposed method is equivalent to that required for each step of a single-step explicit method. Consequently, if the step size of the proposed method is set as three times that of a single-step explicit method, or 3/2 times that of a two-sub-step explicit method, their computational costs can be considered equivalent.

Accuracy and stability
When applied to the single-degree-of-freedom test modelẍ + 2ξ 0 ω 0ẋ + ω 2 0 x = 0, where ξ 0 is the damping ratio and ω 0 is the natural frequency, the proposed method can be formulated in compact form as where A is the amplification matrix, whose formulation, not presented here for the sake of conciseness, can be obtained from Eqs. (2)- (5). The characteristic polynomial of the proposed method takes the form where λ denotes the characteristic root, and A j , j = 1, 2, 3, for the undamped case (ξ 0 = 0) are The local truncation error is defined as where x(t) denotes the exact solution at time t. The method is said to be th-order accurate if σ = O( t ) [10]. Expanding Eq. (10) yields Thus, the method is second-order accurate if and third-order accurate if, in addition to Eq. (12), When physical damping is considered (ξ 0 = 0), the method is naturally second-order accurate with Eq. (12), whereas in addition to Eq. (13), the conditions for third-order accuracy also require The proposed method is expected to be at least second-order accurate; hence, Eq. (12) is considered in the following. Taking it into account, the coefficients A j , j = 1, 2, 3 in Eq. (9) can be rewritten as used in the following for the sake of conciseness. After obtaining p 1 , p 2 , q 1 and q 2 , the original parameters β i , γ i can be determined, although not uniquely, according to Eqs. (16). From Eq. (9c), A 3 = 0, so the two nonzero characteristic roots can be expressed as Since an oscillatory solution is expected for the undamped oscillator problem considered herein, λ 1,2 should be a pair of conjugate complex numbers for small values of τ and bifurcate into two real numbers when approaching the stability limit. At the bifurcation point, denoted as τ b , λ 1,2 should be two equal real numbers, and the spectral radius at this point, denoted as ρ b = |λ 1,2 | ∈ [0, 1], is used to denote the degree of algorithmic dissipation of the explicit method. According to |λ 1,2 | = ρ b at τ = τ b , which imposes p 1 and q 1 can be expressed as To simplify the expression, P b ∈ [−1, 1] is used to replace ±ρ b in the following analysis. Consequently, p 1 and q 1 can be rewritten as 0 needs to be satisfied to ensure that λ 1,2 are conjugate complex numbers, to avoid internal bifurcation. With the conditions in Eqs. (20), A 2 1 − 4 A 2 can be expressed as From Eq. (21), the condition for A 2 1 − 4 A 2 ≤ 0 for τ < τ b changes to f (τ ) ≥ 0. Using x instead of τ 2 , f (τ ) can be rewritten as an equivalent quartic function g(x), as Since a 4 ≥ 0, the condition g(x) ≥ 0 imposes that the plot of g(x) cannot cross the real axis for x ≥ 0. In other words, the unary quartic equation g(x) = 0 can only have double real roots or complex roots in pairs when x ≥ 0. Considering the critical situation that g(x) = 0 has exactly two double roots, which means that g(x) can be expressed as the form (ax 2 + bx + c) 2 , we can obtain the following conditions From Eq. (24), two sets of critical values of p 2 and q 2 are solved, as and If p 2 and q 2 satisfy Eqs. (25) or (26), g(x) ≥ 0 is naturally satisfied for x ≥ 0.
Once the relations p 1 , p 2 , q 1 , q 2 with respect to ρ b and τ b are determined, Eq. (27) can be used to give the allowable range of τ b for a given P b . If they are set by Eqs. (20) and (25), the range of τ b for several P b can be solved as (28a) If p 1 , p 2 , q 1 , q 2 are set by Eqs. (20) and (26), the range of τ b for several P b can be solved as (29b) As can be seen, the coefficients in Eqs. (20) and (26) with P b = −ρ b offer a considerable range of τ b . When ρ b = 1, the maximum τ b can reach 6. Since the CDM has τ b = 2, which was proved to be the largest stability domain among single-step explicit methods [20], and the two-sub-step NB method can achieve τ b = 4 at the most [28], it is reasonable to say that the coefficients in Eqs. (20) and (26) with P b = −ρ b can provide the largest stability interval for the proposed three-sub-step method. Other possible values of p 2 and q 2 that satisfy A 1 − 4 A 2 ≤ 0 always make the maximum allowable value of τ b smaller. Therefore, the coefficients in Eqs. (20) and (26), where P b is replaced by −ρ b , are finally selected, expressed as From the analysis, the coefficients in Eqs. (30) can ensure that λ 1,2 are conjugate complex numbers in τ ∈ [0, τ b ] and bifurcate into two real numbers at τ = τ b . Moreover, they can maximize the allowable range of τ b , which needs to satisfy the condition The allowable range of τ b can be solved from As ρ b changes from 0 to 1, the maximum τ b of the proposed method increases from 5.5425 to 6, while for ρ b ∈ [0, 1], the NB method [28] has τ b ∈ [3.4142, 4], the EG-α method [12] has τ b ∈ [1.4142, 2], and the TW method [34] [25] has τ b ∈ [1, 2]. Since 5.5425/3 > 3.4142/2 > 1.4142 > 1, or 1.8475 > 1.7071 > 1.4142 > 1, the proposed method can offer a broader stability interval for a given 0 ≤ ρ b < 1. Besides, from Eqs. (13) and (16), the method has third-order accuracy for undamped systems if p 1 − q 1 = 1/12, which requires For example, if ρ b = 0, the method has third-order accuracy with τ b ≈ 5.1451; if ρ b = 0.5, it has third-order accuracy with τ b ≈ 5.4495. For simplicity, the value of τ b that makes the method achieve third-order accuracy is denoted as τ b3 , and the maximum value of τ b is denoted as τ bm . Generally, τ b3 < τ bm for a given ρ b ∈ [0, 1], so the second-order schemes can offer a broader stability interval. Considering the undamped case, Figs. 1 and 2, respectively, plot the spectral radii of the proposed threesub-step method (τ b = τ bm ) and of the two-sub-step NB method for different values of ρ b . It can be clearly seen that under the same computational cost, that is τ New /3 = τ NB /2, the proposed method allows a broader stability interval for a given ρ b and can better preserve low-frequency dynamics with τ b = τ bm .
Figures 3-4 present the amplitude decay ratio, which is the damping ratio ξ in the numerical solution for the undamped model, of the proposed method (τ b = τ bm ) and the NB method for different values of ρ b , respectively. Figures 5-6 present their period elongation ratio, which is the relative error of period, (T −T 0 )/T 0 , of the numerical solution applied to the undamped model. As can be seen, as ρ b decreases from 1 to 0, the two methods both exhibit stronger amplitude decay and better period accuracy. The results also illustrate that for a given ρ b < 1, the proposed method with τ b = τ bm exhibits much smaller amplitude decay than the NB method, but it has larger period error.
In addition to ρ b , the new method has another adjustable parameter, τ b . For comparison, Fig. 7 plots the spectral radii, and Fig. 8 shows the amplitude decay ratios and period elongation ratios of the proposed method with  [25] to avoid excessive loss of accuracy. Since in the NB method p = 0.54, which corresponds to ρ b = 0.45, is recommended [28], Figs. 9 and 10 also plot the spectral characteristics of these methods, where the second-order methods use ρ b = 0.45 and the TW method still uses ϕ = 1.033. From Eqs. (31) and (32), with ρ b = 0.0, the new method has τ bm ≈ 5.5425 and τ b3 ≈ 5.1451, and with ρ b = 0.45, the new method has τ bm ≈ 5.7728 and τ b3 ≈ 5.4241. In these figures, the abscissa is set to τ/n, where n = 3 for the proposed method, n = 2 for the NB method, and n = 1 for the EG-α and TW methods, to compare under equivalent computational costs. The results illustrate that as τ b decreases from τ bm to τ b3 , the proposed method shows stronger algorithmic dissipation and better period accuracy. With τ b = τ bm , its amplitude decay ratios are very close to 0 for τ/n ≤ 0.4, and with τ b = τ b3 , its period elongation ratios are quite small for τ/n ≤ 0.4. When τ b is less than τ b3 , the period accuracy no longer improves as τ b decreases, so τ b is best selected in the interval [τ b3 , τ bm ].
From Figs. 7-8, with ρ b = 0, the NB method shows period accuracy close to that of the proposed method with τ b = τ b3 , but its amplitude accuracy and stability interval are not as good as those of the new method with τ b = 5.3. The EG-α method shows a narrow stability interval and poor accuracy in this case. From Figs. 9-10, with ρ b = 0.45, the NB and EG-α methods show spectral characteristics close to those of the new method with τ b = 5.6 or τ b = 5.7. The TW method shows larger amplitude and period errors, especially for τ/n < 0.5, in both cases.  From the comparisons, for a given 0 ≤ ρ b < 1, the new method can offer higher-amplitude accuracy and a broader stability interval than the EG-α and NB methods with τ b close to τ bm and higher period accuracy with τ b close to τ b3 . With a suitable τ b ∈ (τ b3 , τ b3 ), it can reproduce the spectral characteristics of these secondorder methods. Therefore, the schemes with τ b = τ bm are recommended for high-efficiency purposes because of their maximum stability interval. These schemes are more conservative in the low-frequency domain, but their phase accuracy is slightly reduced compared with the NB method. The schemes with τ b = τ b3 are more suitable for high-accuracy or strong-dissipation purposes, as they possess favourable period accuracy and large-amplitude decay.

Overshoot
So far, the parameters in the formulation cannot be uniquely determined yet. For linear analysis, as long as they satisfy Eqs. (16) and (30), the proposed method can achieve the accuracy and stability characteristics described in Sect. 3.1. However, the spectral characteristics only control the long-term behaviour of a method. For the short-term behaviour at the initial steps or after sudden switches, the norm of the amplification matrix A also needs to be considered to avoid overshoot [10]. If A has a large norm for the high-frequency dynamics, the solutions may produce pathological growth at the initial steps, the so-called overshoot.
Considering the test problemẍ + ω 2 0 x = 0, the recursive scheme can also be expressed as with As can be seen, a 11 , a 12 , a 21 , a 22 contain the terms τ 4 , and τ 6 , so for a large τ close to τ b , these elements may have large absolute values, even though the spectral radius of A is small. The large elements may cause the initial values to be enlarged abnormally during the first few steps, resulting in overshoot. To avoid overshoot as much as possible, |a 11 |, |a 12 |, |a 21 |, and |a 22 | need to be as small as possible. From the spectral analysis, a 11 , a 12 , a 21 , a 22 satisfy the following conditions where A 1 and A 2 are determined in Sect. 3.1. Thus, to make |a 11 | and |a 22 | as small as possible, an artificial assumption a 11 = a 22 is introduced. From Eq. (34), it yields two conditions Unfortunately, no condition can be used to constrain |a 12 | and |a 21 |. After combining Eqs. (12), (16) and (36), one condition is still missing to determine all parameters, so another assumption γ 2 = 2γ 1 , which imposes that the first two sub-steps have equal size, is used. Then, the parameters can be obtained as To evaluate the overshoot characteristic of the proposed scheme with the parameters in Eq. (37), the 2-norm of A, expressed as is employed, and A 2 is expected to be small to avoid large values in |a 11 |, |a 12 |, |a 21 | and |a 22 |. Figure 11 plots A 2 versus τ/τ b of the proposed scheme with the parameters in Eq. (37) and τ b = τ bm , as well as the NB method with p = 0.54. It can be seen that although the peak values of A 2 of the proposed schemes are higher than that of the NB method, they are not greater than 3, which is still too small to cause observable overshoot. When τ > 0.8τ b , A 2 in all schemes reduces to a small value less than 2, which further prevents the occurrence of high-frequency overshoot. Therefore, with the set of parameters in Eq. (37), the proposed method possesses satisfactory overshoot characteristics. Generally, overshoot occurs in implicit, unconditionally stable methods, but owing to the broad stability region, the proposed explicit method also risks suffering from overshoot. For example, if a set of parameters is chosen that only satisfies the spectral characteristics, as Eqs. (12) and (16), such as the corresponding A 2 with τ b = τ bm is shown in Fig. 12. The peak values can reach about 25 at τ ≈ 0.9τ b , which is large enough to cause noticeable overshoot. Therefore, the analysis of the overshoot characteristic is also necessary for explicit methods with large stability limits. For damped systems (ξ 0 = 0), the remaining three parameters γ 4 , γ 7 and γ 8 can be set as where γ 7 and γ 8 satisfy the third-order accuracy condition in Eq. (14) to gain as high accuracy as possible.
Using the parameters in Eqs. (37) and (40), Fig. 13 plots the spectral radii of the proposed method with τ b = τ bm for the damped case ξ 0 = 0.05. As can be seen, the internal bifurcation appears in a small interval, and the stability limits are reduced compared with the undamped case. No optimal values of γ 4 , γ 7 and γ 8 were found for arbitrary values of ξ , so Eq. (40) shows a selection considering high accuracy. According to the analysis in Sects. 3.1 and 3.2, the parameters in Eqs. (37) and (40) are finally the recommended ones for the proposed method considering accuracy, stability, algorithmic dissipation and overshoot characteristics.

Dispersion analysis of wave propagating problems
In terms of wave propagation problems, Refs. [21,28] present a method to measure the dispersion error using the displacement-based spatial discretizations. Following the procedure in Ref. [28], the dispersion analysis of the proposed method is discussed in this subsection, to support the selection of suitable parameters, including ρ b , τ b and the CFL number, for these problems.
The scalar wave propagation equation is employed as the model problem, where u is the field variable and c 0 is the exact wave velocity. Considering the 2-dimensional case, the exact solution has the form u(x, y, t) = A 0 e i(k 0 x cos θ +k 0 y sin θ −ω 0 t) (42) where (x, y) is the spatial coordinate of the solution, ω 0 is the exact frequency, k 0 = ω 0 /c 0 is the exact wave number, and θ is the propagation angle measured from the x-axis.
In terms of the numerical solution, the governing equation after spatial discretization can be written as where U summarizes the numerical solutions of all nodes, and M and K are assembled by the corresponding element matrices M e and K e , respectively. Using 4-node elements, and assuming the element length x = y = h, M e and K e can be written as where the lumped mass matrix is used to gain the efficiency advantage of explicit methods. From Sect. 3.1, the modal degree of freedom x satisfies the recursive scheme where p 1 , p 2 , q 1 and q 2 are obtained in Eq. (30) and τ = ω t. Considering all modal degrees of freedom of the problem in Eq. (43), we have where X is the vector collecting all modal degrees of freedom, Λ is the corresponding diagonal matrix composed of all ω 2 i , and I is the unit matrix. The finite element degrees of freedom can be expressed as where Φ is the matrix composed of all corresponding eigenvectors, and it satisfies Multiply Eq. (46) at the left by Φ follows According to the element matrices M e and K e in Eq. (44), the row of global mass matrix M for nonboundary nodes is So the term M −1 in Eq. (49) can be replaced by 1/ h 2 . Using CFL = c 0 t/ h, we can obtain Considering the middle node of a 4 4-node elements patch, the row of global stiffness matrix K is which means that the term K U t in Eq. (51) has the form where u(n x h, n y h, n t t) denotes the numerical solution for the middle node (n x h, n y h) at time n t t. The expressions of the terms +Â 2 u(n x h, n y h, (n t − 1) t) = 0 (58) whereÂ 1 =2 − CFL 2 κ + p 1 CFL 4 κ 2 + p 2 CFL 6 κ 3 (59a) It follows that Consequently, when CFL, kh, θ and the algorithmic parameters are specified, the dispersion error (c − c 0 )/c 0 can be solved by Eq. (62a). The dispersion error is expected to be close to 0 for all participating modes. However, since the higher modes of about kh/π > 0.6 [21,30] cannot be resolved spatially or cannot be accurately described by the time step size, their participation can greatly pollute the overall solution and cause loss of accuracy. Therefore, the algorithmic dissipation is expected to increase rapidly when kh/π > 0.6 to filter out the higher modes. Considering different CFL numbers, algorithmic parameters and θ , the dispersion error (c − c 0 )/c 0 and the degree of algorithmic dissipation, represented by the spectral radius ρ = |λ| = Â 2 , versus kh/π are discussed in the following. Using θ = 0 and different CFL numbers, Fig. 14 shows the dispersion errors and spectral radii of the proposed method with ρ b = 0.0 and τ b = τ bm . The results illustrate that as the CFL number increases, the dispersion error gets smaller, and the algorithmic dissipation gets stronger. With CFL = τ b /2, the bifurcation point of the spectral radius is exactly at kh/π = 1.0, so the algorithmic dissipation can play its role to a great extent. The conclusion also applies to other explicit schemes. As shown in Eq. (59), if kh/π = 1, θ = 0 and   16 Dispersion errors and spectral radii versus kh/π of the proposed method with ρ b = 0.0 and several existing methods at θ = 0 using CFL= τ b /2 CFL = τ b /2, the spectral radius is exactly at the bifurcation point. With CFL > τ b /2, the spectral radius will bifurcate when kh/π < 1. Consequently, the CFL number is best set as τ b /2.
Using θ = 0 and CFL = τ b /2, Fig. 15 plots the dispersion errors and spectral radii of the proposed method with τ b = τ bm and different values of ρ b . The results show that as ρ b gets smaller, the dispersion error increases for the modes with kh/π ≤ 0.6, and the algorithmic dissipation gets stronger for those with kh/π > 0.6. When ρ b = 1, the scheme with τ b = τ bm has no dispersion error, but also no dissipation, like the CDM. Since in the proposed method also τ b is adjustable, two cases, ρ b = 0.0 and ρ b = 0.45, with different τ b ∈ [τ b3 , τ bm ] are discussed in the following.
Considering θ = 0 and CFL = τ b /2, Fig. 16  As can be seen, the proposed method with τ b = τ bm shows the smallest dispersion error for the modes with kh/π ≤ 0.5, and as τ b decreases, the dispersion accuracy worsens, but the range of discarded modes increases. With ρ b = 0.0, the proposed method with τ b = 5.40 presents a spectral radius very close to that of the NB and EG-α methods, but its dispersion errors are clearly smaller for the participating modes. With ρ b = 0.45, the case τ b = 5.70 shows a spectral radius almost identical to that of the NB method, stronger algorithmic dissipation than the EG-α method for kh/π > 0.5 and also smaller dispersion error than these two methods for kh/π ≤ 0.5. The TW method exhibits good dispersion accuracy, but the smallest algorithmic dissipation for kh/π > 0.6. Despite this fact, the TW method shows powerful filtering ability in the numerical examples in Sect. 4.
From the comparisons, among the second-order methods, with ρ b = 0.0 and ρ b = 0.45, the proposed method with τ b = 5.40 and τ b = 5.70, correspondingly, is superior to the EG-α and NB methods, because they have better dispersion accuracy for the participating modes, and stronger or at least equivalent algorithmic dissipation for the discarded modes. With a larger τ b ≤ τ bm , the proposed method has smaller dispersion error, and with a smaller τ b ≥ τ b3 , it better dissipates the unwanted modes.
Considering θ = π/6 and CFL = τ b /2, Figs. 18 and 19 show the dispersion errors and spectral radii of the employed methods with ρ b = 0.0 and ρ b = 0.45, respectively. Figures 20 and 21 show the analogous results considering θ = π/4 and CFL = τ b /2. These results indicate that the proposed method with τ b = 5.40 for the case ρ b = 0.0, and τ b = 5.70 for ρ b = 0.45, still holds a certain accuracy advantage for different angles over the NB and EG-α methods, but their gap decreases as θ increases. From the plots, with θ > 0 the spectral radii have not dropped to ρ b when kh/π = 1.0, so the effect of increasing θ is similar to that of reducing the CFL number, which ensures the stability of the methods at different angles, and also reduces the differences between these methods.  Table 1 Step-by-step solution using the proposed three-sub-step explicit method for solving linear problems Certainly, if higher efficiency is required, a larger τ b ≤ τ bm can also be used with the proposed method.

Illustrative solutions
In this section, several single degree-of-freedom linear and nonlinear examples are solved to demonstrate the properties of the proposed method, and some of the benchmark wave propagation problems, used in Refs. [14,21,28,29,39], are simulated to assess the performance of the proposed method. The proposed method with ρ b = 0.45 and different τ b , the NB method with p = 0.54 [28], the EG-α method with ρ b = 0.45 [12] and the TW method with ϕ = 1.033 [34] [25] are employed for comparison. The EG-α method with the explicit treatment of physical damping [12] is selected, and its improved formulation illustrated in [1] is used for the computations. The step-by-step computational procedure of the proposed method used for the linear problems is presented in Table 1.

One-dimensional wave propagation in a clamped-free bar
The 1-dimensional wave propagation in a clamped-free bar excited by an external load F at the free end [29] shown in Fig. 22 is considered. The elastic modulus, density, cross-sectional area and length of the bar and the  τ b = τ bm shows the closest solution to the reference one at the switch point, followed by the scheme with τ b = 5.7, the NB method, the EG-α method and the scheme with τ b = τ b3 , but as one would expect, the scheme with τ b = τ bm shows more oscillation. Considering accuracy and dissipation capability, the proposed method with τ b = 5.7 performs best. These conclusions are consistent with the dispersion error analysis in Sect. 3.3.

Two-dimensional wave propagation in a pre-stressed membrane
As shown in Fig. 24, the 2-dimensional membrane [21,28] excited by a lumped load at the centre, is considered. The wave propagation equation has the form The wave velocity c 0 is assumed as 1, and the load f (0,  Figure 25 summarizes the snapshots of displacement at t ≈ 9139/800 resulting from the employed methods. In this example, CDM using CFL = 1 is also employed, but its solution is clearly unsatisfactory due to the high-frequency pollution. Among the remaining methods, the proposed one with τ b = τ b3 and TW shows smoother solutions; the solutions of the other second-order methods are also quite accurate. The displacements and velocities of θ = 0 and π/4 in the interval 8 ≤ x 2 + y 2 ≤ 12 are plotted in Fig. 26. Due to the use of constant time steps, the results given by these methods are the solutions of the step closest to t = 9139/800, but not exactly at t = 9139/800, so the offset between the numerical solution and the reference one cannot be used to account for the dispersion error of the method. For example, the scheme with τ b = τ b3 shows the results at t ≈ 11.1646 actually. From Fig. 26, it can be observed that the second-order methods all exhibit slight oscillations near the position x 2 + y 2 ≈ 10. The NB method, the EG-α method and the proposed method with τ b = 5.70 show very close solutions, while the scheme with τ b = τ b3 shows slighter oscillations because of its stronger dissipation. These conclusions are also consistent with the comparison of algorithmic dissipation in Sect. 3.3.

A lamb problem
The lamb problem used in [21,28] is considered. As shown in Fig. 27, the elastic medium is excited by the line load F at x = 0, y = 0 in plane strain conditions. The mass density is set as ρ = 2200. The velocities of P-wave and S-wave are c P = 3200, and c S = 1847.5. The CFL number is computed using c P . Two receivers are placed at x = 640, y = 0 and x = 1280, y = 0. In the first case, the Ricker wavelet line force, as is applied, where the parameters are assumed asf = 12.5 and t 0 = 0.1. The plane is meshed into 1280 × 640 4-node elements. Figure 28 shows the time history of displacements at two receivers. Figure 29 presents the snapshots of the von Mises stress at the fixed time point t ≈ 0.999. All employed methods predict very accurate solutions, the TW method showing the largest errors at the peaks in Fig. 28. The wave fronts of the P-wave, S-wave and Rayleigh wave can be clearly seen in Fig. 29. In this case, since the contribution of high frequencies to the solution is very small, the proposed method with τ b = τ bm is a better choice, owing to its high efficiency and high dispersion accuracy.
In the second case, the line force is changed to Since more wave modes can be excited by the force, a finer mesh of 3200 × 1600 finite elements is used for the spatial discretization. Figures 30-31 show the displacements at the two receivers and the snapshots of the von Mises stress at t ≈ 0.999, respectively. Oscillations can be observed in this case. From the zoom-in figures in Fig. 30, the first-order TW method can filter out oscillations faster; among the second-order methods, the proposed scheme with τ b = τ b3 shows stronger dissipation ability. In this case, because of the high-frequency pollution, τ b ∈ [τ b3 , 5.70] can be used in the proposed method to make the solution smoother.  To compare the computational efficiency, Table 2 lists the required CPU time and the total number of sub-steps of the employed methods for solving the period t ∈ [0, 0.999]. In single-step methods, each step is counted as 1 sub-step. The simulations were run on an Intel i5-8265U CPU @ 1.60 GHz with 8.00 GB RAM. As can be seen, the effort of each sub-step for these explicit methods is very close, so the proposed method with τ b = τ bm and τ b = 5.70, as well as the TW method, is slightly more efficient than the remaining schemes.

Conclusions
A novel second-order accurate three-sub-step explicit time integration method is proposed. The optimal parameters, controlled by τ b and ρ b , i.e. the position of the bifurcation point, and the spectral radius at that point, are defined according to the requirements of accuracy, stability, algorithmic dissipation and overshooting characteristics.
A distinctive feature of the proposed method is that, in addition to ρ b , also τ b can be adjusted to control the properties. For a given ρ b , when τ b = τ bm , the proposed method has a broader stability interval and smaller amplitude decay than the NB and EG-α methods. As τ b decreases from τ bm to τ b3 , the proposed method shows higher period accuracy and stronger algorithmic dissipation. By selecting an appropriate τ b ∈ (τ b3 , τ bm ), the proposed method can reproduce spectral characteristics close to those of existing second-order methods. The analysis of the spectral characteristics is illustrated by the convergence rates of several single degree-of-freedom examples, where the scheme with τ b = τ b3 exhibits an appreciable accuracy advantage.
For the analysis of wave propagation, the dispersion error is discussed based on the displacement-based spatial discretizations. As ρ b or τ b decreases, the proposed method has larger dispersion error for the participating modes, but the range of discarded modes gets wider. For a good trade-off between accuracy and algorithmic dissipation, one set of parameters, as ρ b = 0.45, τ b = 5.70, and CFL = τ b /2 = 2.85, is selected in the proposed method. The recommended scheme shows better dispersion accuracy than the NB method ( p = 0.54, CFL = 1.85) and the EG-α method (ρ b = 0.45, CFL = 0.90), and almost the same degree of algorithmic dissipation of the NB method ( p = 0.54), and stronger dissipation than the EG-α method (ρ b = 0.45). If higher efficiency is expected, a larger τ b ≤ τ bm can be used in the proposed method, whereas for strong dissipation demand, a smaller τ b ≥ τ b3 is more useful.
These considerations are illustrated by solving some benchmark wave propagation problems, where the proposed method with ρ b = 0.45 and different τ b , the NB method with p = 0.54, the EG-α method with ρ b = 0.45 and the TW method with ϕ = 1.033 are considered. The comparison indicates that if the solutions are mainly contributed by low-frequency dynamics, the scheme with τ b = τ bm is a better choice, since it has higher dispersion accuracy and requires less computational effort than the other schemes. On the other hand, if the solutions contain significant high-frequency pollution, a smaller τ b ∈ [τ b3 , 5.70] can be used to quickly filter out the high-frequency dynamics. The NB and EG-α methods show performances similar to those of the proposed method with τ b = 5.70. The TW method shows powerful filtering ability, but it cannot predict accurate peak values.