An Improved Time Integration Method Based on Galerkin Weak Form with Controllable Numerical Dissipation

: Based on the newly proposed generalized Galerkin weak form (GGW) method, a two-step time integration method with controllable numerical dissipation is presented. In the ﬁrst sub-step, the GGW method is used, and in the second sub-step, a new parameter is introduced by using the idea of a trapezoidal integral. According to the numerical analysis, it can be concluded that this method is unconditionally stable and its numerical damping is controllable with the change in introduced parameters. Compared with the GGW method, this two-step scheme avoids the fast numerical dissipation in a low-frequency range. To highlight the performance of the proposed method, some numerical problems are presented and illustrated which show that this method possesses superior accuracy, stability and efﬁciency compared with conventional trapezoidal rule, the Wilson method, and the Bathe method. High accuracy in a low-frequency range and controllable numerical dissipation in a high-frequency range are both the merits of the method.


Introduction
To solve general structural dynamic problems, the direct integration method has been proven to be effective and has been widely used in engineering. The commonly used integration methods can be categorized into the explicit one and the implicit one. The explicit method is conditionally stable, and a smaller time step should be used in the calculation. Thus, for structural vibration analysis, the implicit method is usually considered to be more effective [1]. As for implicit method, it can be designed to be unconditionally stable and to possess controllable numerical dissipation. Although the calculation needs matrix factorization and is cost-consuming [2], the implicit method can allow a larger time step due to its superior numerical stability, which brings much more efficiency in solving dynamic problems [3]. At present, there are some well-known implicit algorithms, such as the Newmark method [4], the Wilson-θ method [5], the generalized-α method [6], the HHT method [7], and the Bathe method [8], etc.
However, some algorithms may lead to instability in the nonlinear range [9] due to non-conservation of energy and momentum [10]. Taking trapezoidal rule for example, in the linear range, it is unconditionally stable, and has second-order accuracy without amplitude decay [11,12]. However, this method can easily become unstable when extended to the nonlinear problems [12,13], and the response will diverge [3,10,14] when using an inappropriate length of time steps. To solve this problem, numerical damping is introduced into the Newmark method by adjusting the parameters, but it results in first-order accuracy [15].
For most of the structural dynamic problems, the main contribution of the response comes from its low-order frequency modes [16]. Considering the impact of discretization, the high-order frequency modes of the structure will change, which affects the computational accuracy. Therefore, the numerical damping which can control numerical dissipation is necessary for both low-order and high-order frequency ranges. In addition, it can also improve the convergence of the iterative solution and increase the numerical stability. Nevertheless, unreasonable numerical damping will result in a great loss of accuracy, which leads to a large attenuation of structural response. It is usually caused by excessive numerical dissipation in a low-frequency range. To evaluate the numerical dissipation of the algorithm, the spectral radius ρ is used. The change in dissipation is reflected in the spectral radius ρ; that is, the reduction in spectral radius illustrates the dissipation of the algorithm. When the time step is large, the spectral radius will quickly attenuate for some algorithms, such as part of the Newmark method, Houbolt method, Wilson-θ method, and part of the generalized α method, etc. [17] In summary, it is of necessity to find a time integration algorithm that has both reasonable and controllable numerical damping. This algorithm can generate good accuracy in a low-frequency range with controlled high dissipation in a high-frequency range [7,18,19], and is expected to perform well in nonlinear field.
Many researchers have proposed useful methods to control the accuracy of algorithms in different frequency ranges. Some composite algorithms with several steps perform quite well [20][21][22]. A typical one is the Bathe method proposed by Bathe K.J. et al., which divides a time step into two sub-steps to combine two well-known time integration schemes [3,8,12,23]. The trapezoidal integration method and the three-point Euler backward method are respectively used in the first and second sub-step. The Bathe method has reasonable numerical damping and stability but cannot control the dissipation [24]. Based on this, the ρ θ -Bathe method and β 1 /β 2 -Bathe method were proposed [23,25]. According to the evaluation, these methods possess effective numerical dissipation, and the β 1 /β 2 -Bathe method is the special case of the ρ θ -Bathe method [26]. Compared with the generalized-α method, this series of Bathe methods show much better performance [24]. Some multi-step algorithms are further proposed based on the methods mentioned above [27][28][29][30]. They are usually composed of a non-dissipative algorithm, e.g., trapezoidal rule, and a strongly dissipative one [29].
Focusing on the nonlinear problems, a time integration method based on Galerkin weak form, i.e., the Galerkin weak form (GGW) method is proposed [31], considering weighted residuals of equations [32]. The Newton-Raphson iteration method is used first so that the iterative linear procedure can be constructed. To solve the resulting linear equations, a time-domain finite element method (FEM) based on Galerkin weak form is established [33][34][35], in which the Gaussian numerical integration is used to solve the integral term. For different degrees m (m = 1, 2, 3) of the time elements, the implicit time integration schemes are derived. At the same time, the numerical damping is controlled by modifying the numerical integration form. According to the analysis, the GGW method is an unconditionally stable, self-starting and one-step method. Compared with the Newmark method, Wilson-θ method, and Bathe method, the GGW method has better stability and efficiency. However, for some representative nonlinear problems, it is found that in the GGW method, the numerical dissipation may decay too fast and the speed of the dissipation cannot be controlled in some cases. Hence, the GGW method is expected to be improved to possess a more effective numerical dissipation performance.
Based on the deep study of the GGW method and inspired by the idea of the Bathe method, an improved two-step time integration method is proposed hereinafter. The time step of this method is divided into two sub-steps. In the first sub-step, the GGW method is used, and in the second sub-step, another parameter is introduced by the trapezoidal integral, i.e., a non-dissipative algorithm is combined with the GGW method so that the improved algorithm can solve the defects of the GGW method but retaining its excellent characteristics. After that, the corresponding iterative procedure is established and presented, and then the accuracy and numerical stability of the algorithm are discussed.
Typical numerical examples including an undamped single-degree-of-freedom (single-DOF) system, the nonlinear pendulum, and a large-span plane skeletal system are solved, through which it can be seen that the numerical dissipation is controllable in the proposed two-step time integration method and sensitivity of parameters is reduced. It dramatically improves the performance of the GGW method. Compared with other algorithms, such as the trapezoid rule, the Wilson-θ method, and the Bathe method, the results sufficiently show the high accuracy and computational efficiency of the proposed method.

Governing Equations and Iterative Schemes
The nonlinear dynamic problems can be usually attributed to the following nonlinear second-order ordinary differential equations (ODEs), which are named dynamic equations hereinafter: u are the nodal displacement vector, velocity vector, and acceleration vector respectively, M is the mass matrix, f d is the damping force vector, f s is the resilience vector, f ext is the vector of externally applied forces, and u 0 , . u 0 are the known displacement and velocity vectors at the initial time, respectively. The time domain is (0, t f ]. Considering time t in the time interval [t n , t n+1 ], when f d and f s are the nonlinear functions related to u and . u, the Newton-Raphson method can be used to transform the nonlinear second-order ODEs (Equation (1)) into linear equations as follows: where In Equations (2)-(4) u (k) , int are all functions of time t. Therefore, through Newton's method, the process of solving nonlinear Equation (1) is transformed into an iterative solution of linear Equation (2). For convenience, the superscript k standing for the k-th iteration step is omitted in Equation (2), and Equation (2) is rewritten as follows:

Galerkin Finite Element Method
Many methods can be used to solve Equation (5), and Ref. [31] uses the time FEM based on the Galerkin weak form to solve Equation (5). Setting v ∈ H 1 v to be the test function and u h ∈ H 1 E to be the trial function, the weighted residual weak form of Equation (5) in [t n , t n+1 ] is established as follows: where H 1 E and H 1 v are the trial function space and the test function space, respectively. Regarding the time interval [t n , t n+1 ] as a single time element of length ∆t = t n+1 − t n , the test function and trial function can be expressed with Equation (7), in which the time interval [t n , t n+1 ] is divided into m (m = 1, 2, 3, . . . ) intervals so that the m-degree interpolation polynomial is used.
where N = [N 1 (t)I, N 2 (t)I, . . . , N m+1 (t)I] is the shape function matrix and m is the degree of time elements. N i (t) (i = 1, 2, . . . , m + 1) is the shape function of Lagrange polynomial form and I is the unit matrix.
Substituting Equation (7) into Equation (6) and implementing partial integration, the following equation can be derived by coordinate mapping, in which time From Equation (8), for any positive integer m, we can derive and obtain a time integration scheme. In fact, considering both the calculation efficiency and accuracy, letting m = 1 or 2 or 3 is sufficient [31]. At the same time, it can be proved that when reduced Gaussian integration is used, the derived scheme is unconditionally stable for linear problems [33]. When m = 1 and using the Gaussian integration, the time integration scheme from Equation (8) can be presented as where where ω l , ζ l are the Gaussian integral coefficients and locations, respectively, and n l is the number of Gaussian points. Finally, when the reduced Gaussian integration with one Gaussian integration point is used, it can be simplified to the following equations: In summary, the nonlinear dynamic equation, Equation (1), is converted into the linear dynamic equation, Equation (5), with Newton-Raphson iteration, and then the Galerkin FEM is used to solve this ODE. This is the basic idea and the key strategy of the GGW method. However, the numerical damping is not introduced until Equations (14) and (15).

Improvement of Algorithm
In order to introduce numerical dissipation into the GGW method, Ref. [31] proposes a feasible time integration scheme. A parameter α is introduced to control the numerical damping by changing the form of the test function as Equation (16). At the same time, to calculate the integral terms in Equations (12) and (13), a specific numerical integration rule as Equation (17) is constructed using Equation (16).
In this way, the time integration scheme with numerical dissipation for linear time elements is established based on the Galerkin weak form. However, the parameter α in Equations (16) and (17) turns out to be too sensitive under certain circumstances, which will result in too fast decay (this is demonstrated in Section 4.3).
To overcome this shortcoming of the GGW method, here we present an improved method in which the time step is divided into two subsections and another parameter β (0 < β < 1) is introduced in addition to α.
In the first step, we take the step length to be β∆t, and by using the GGW method, we obtain the solutions u n+β , u n+1 at t n+1 (t n+1 = t n + ∆t). The detailed computational strategy of the second sub-step is as follows: Considering the definition of integral as Equation (18) shown as follows: u(t)dt (18) and using the algorithm of trapezoidal integral, Equation (18) can be rewritten into Equation (19).
It is obvious that, when β = 1, the second sub-step vanishes and the proposed method degenerates to the GGW method; when β = 0, the first sub-step vanishes and the proposed method degenerates to the Newmark method.
The iterative algorithm of the proposed method is finally summarized into the following steps (Algorithm 1). Determine parameters in the algorithm:β, α, and calculate c 1 , c 2 , c 3 . Aim Solve the equation at t n + ∆t when knowing the result at t n . 2 //--First sub-step: solve the dynamic response u n+β , Use the step-by-step time integration algorithm based on the GGW method: . u n ; Let k be the number of iteration steps, starting from 1 to meet the requirements.

2.2
Since u n+β are known, the k+1 step of iteration is performed: Determine the parameters: n+β using the following equation: n+β > Tol iter then k = k + 1 and return to 2.2, otherwise go to the next step.
. u n ; k is the number of iterations, starting from 1 to meet the requirements.
n+1 are known, the (k + 1)th step of iteration is performed: Determine the parameters: Solve the response by the following equations: n+1 > Tol iter then k = k+1 and return to step 3.2, otherwise go to the next step. 3.7 Get the answers u n+1 , . u n+1 and .. u n+1 3.8 End

Analysis of Numerical Stability and Accuracy
The following single-DOF system is considered to analyze the stability and accuracy of the proposed method: where ω is the circular frequency of free vibration. The relationship between the response of the (n + 1)-th time step and the n-th time step can be obtained by referring to Equations (11) and (20), that is which leads to As we know, the transfer matrix A in the following equation is generally used to analyze the stability of the method: Using ρ(A) to represent the spectral radius of A, it is well-known that when ρ(A) ≤ 1, the algorithm is stable. Referring to Equation (23), for the proposed two-step algorithm, the transfer matrix A can be written into the following formula: where B is the transfer matrix computed in the first sub-step, i.e., B satisfies As the transfer matrix of the GGW method, B was derived from Ref. [31]: where Λ = ωβ∆t, α is the parameter in the GGW method to control the numerical damping. It has already been concluded that the spectral radius ρ(B) ≤ 1 [33,34] when using reduced Gaussian integration in Equation (8) and the GGW method is unconditionally stable [31]. By using Equations (25) and (27), the spectral radius of the present two-step algorithm is derived: where Ω = ω∆t. With T denoting the period of free vibration, when ∆t/T tends to infinity, the spectral radius ρ(A) turns out to be Appl. Sci. 2021, 11, 1932 8 of 21 Through Equation (29), the stability of the algorithm with different parameters can be predicted. Generally, we take β = 1/2 and then c 1 = c 3 = 1/4, c 2 = 1/2, which results in the following expression: When ρ ∞ both and β are given, the corresponding parameter α can be computed from Equation (29). Generally, when ∆t/T tends to infinity, the suggested ρ ∞ satisfies ρ ∞ ∈ [0.5, 0.8]. Therefore, according to the Equation (30), when β = 1/2, α ∈ [2/19, 2/7]. In this case, the change in ρ(A) relating to ∆t/T is plotted in Figure 1 with ρ ∞ = 0, 0.4, 0.8 and β = 0.5, 0.65, 0.8, respectively, which aims to evaluate the impact of α and β on ρ(A). Furthermore, to evaluate the errors of this method, the period elongation (PE) and amplitude decay (AD) are plotted in Figures 2 and 3 as well.
When both   and β are given, the corresponding parameter α can be computed from Equation (29). Generally, when tT tends to infinity, the suggested   satisfie . Therefore, according to the Equation (30) When both   and β are given, the corresponding parameter α can be computed from Equation (29). Generally, when tT tends to infinity, the suggested   satisfie . Therefore, according to the Equation (30)     From Figure 1, we can find that the value of   is controlled by changing α and β. I the same time, if   is a constant, the decay rate of this algorithm turns out to be large when β is smaller. As shown in Figure 2 and Figure 3, the PE and AD are also related t   and β. When   maintains a constant value, the trends of PE and β are positively co related, while those of AD and β are opposite. Moreover, when   becomes larger, bot AD and PE become smaller.
The comparison of the GGW method and the proposed method is shown in Figure  When tTis small (e. g. , Δ / ≈ 10 −2 in Figure 4), the spectral radius of the GGW metho starts to reduce, while the spectral radius of the proposed method starts to reduce whe tTis larger (e. g. , Δ / ≈ 10 −1 in Figure 4). This illustrates that the proposed metho improved the disadvantage of too fast numerical dissipation in the GGW method.

Numerical Examples
In this section, three numerical examples are considered to demonstrate the accurac and efficiency of the proposed two-step method. Section 4.1 gives an undamped standar single-DOF system, Section 4.2 gives an oscillating nonlinear single pendulum, and Se tion 4.3 gives a long span planar skeletal structure. These examples are solved with th Wilson method (θ = 1.4), Bathe method, GGW method, and the trapezoidal rule as well i From Figure 1, we can find that the value of ρ ∞ is controlled by changing α and β. In the same time, if ρ ∞ is a constant, the decay rate of this algorithm turns out to be larger when β is smaller. As shown in Figures 2 and 3, the PE and AD are also related to ρ ∞ and β. When ρ ∞ maintains a constant value, the trends of PE and β are positively correlated, while those of AD and β are opposite. Moreover, when ρ ∞ becomes larger, both AD and PE become smaller.
The comparison of the GGW method and the proposed method is shown in Figure  4. When ∆t/T is small (e.g., ∆t/T ≈ 10 −2 in Figure 4), the spectral radius of the GGW method starts to reduce, while the spectral radius of the proposed method starts to reduce when ∆t/T is larger (e.g., ∆t/T ≈ 10 −1 in Figure 4). This illustrates that the proposed method improved the disadvantage of too fast numerical dissipation in the GGW method. From Figure 1, we can find that the value of   is controlled by changing α and β. I the same time, if   is a constant, the decay rate of this algorithm turns out to be large when β is smaller. As shown in Figure 2 and Figure 3, the PE and AD are also related t   and β. When   maintains a constant value, the trends of PE and β are positively cor related, while those of AD and β are opposite. Moreover, when   becomes larger, bot AD and PE become smaller.
The comparison of the GGW method and the proposed method is shown in Figure 4 When tTis small (e. g. , Δ / ≈ 10 −2 in Figure 4), the spectral radius of the GGW metho starts to reduce, while the spectral radius of the proposed method starts to reduce whe tTis larger (e. g. , Δ / ≈ 10 −1 in Figure 4). This illustrates that the proposed metho improved the disadvantage of too fast numerical dissipation in the GGW method.

Numerical Examples
In this section, three numerical examples are considered to demonstrate the accurac and efficiency of the proposed two-step method. Section 4.1 gives an undamped standar single-DOF system, Section 4.2 gives an oscillating nonlinear single pendulum, and Sec

Numerical Examples
In this section, three numerical examples are considered to demonstrate the accuracy and efficiency of the proposed two-step method. Section 4.1 gives an undamped standard single-DOF system, Section 4.2 gives an oscillating nonlinear single pendulum, and

A Single-DOF System
In order to evaluate the accuracy of the proposed method, an undamped standard single-DOF system is firstly considered [19]. Its governing equation is: ..
Letting ω = π, the exact solution can be written as: The free-vibration period T of this example is 2 s, the time step ∆t = T/10 = 0.2 s, and the parameters are chosen to be α = 0, β = 0.5, c 1 = c 3 = 1/4, c 2 = 1/2. The vibrational responses within 40-50s of the proposed two-step method, the exact solution, and other time integration methods are shown in Figures 5-7. The computational time of these methods are listed in Table 1. order to give a brief analysis of the performance of the proposed method. All examples are programed in MATLAB 2020b under the same computing environment.

A Single-DOF System
In order to evaluate the accuracy of the proposed method, an undamped standard single-DOF system is firstly considered [19]. Its governing equation is: Letting ω = π, the exact solution can be written as: The free-vibration period T of this example is 2 s, the time step Δt = T/10 = 0.2 s, and the parameters are chosen to be α = 0，β = 0.5, c1 = c3 = 1/4，c2 = 1/2. The vibrational responses within 40-50s of the proposed two-step method, the exact solution, and other time integration methods are shown in Figures 5-7. The computational time of these methods are listed in Table 1.
It can be observed that the Wilson method has a larger amplitude decay when the time step is a little longer, e.g., Δt = 0.2 s in this example, whereas the trapezoid rule, Bathe method, GGW method, and the proposed method can still maintain good accuracy. However, the Bathe method will bring a certain attenuation. Among all these methods, the proposed algorithm shows a lower period elongation. Furthermore, as shown in Table 1, the computational times of these methods are close to each other, thus the numerical efficiency of the proposed method can be proved.

Nonlinear Single Pendulum
In order to verify the effectiveness of the proposed algorithm in nonlinear cases, an oscillating nonlinear single pendulum [36] described in Figure 8 is considered in this subsection.  It can be observed that the Wilson method has a larger amplitude decay when the time step is a little longer, e.g., ∆t = 0.2 s in this example, whereas the trapezoid rule, Bathe method, GGW method, and the proposed method can still maintain good accuracy. However, the Bathe method will bring a certain attenuation. Among all these methods, the proposed algorithm shows a lower period elongation. Furthermore, as shown in Table 1, the computational times of these methods are close to each other, thus the numerical efficiency of the proposed method can be proved.

Nonlinear Single Pendulum
In order to verify the effectiveness of the proposed algorithm in nonlinear cases, an oscillating nonlinear single pendulum [36] described in Figure 8 is considered in this subsection. The governing equation is where = gL, θ(t) is the angle between the rod and the vertical plane at time t, L is the θ 0 (33) where ω = g/L, θ(t) is the angle between the rod and the vertical plane at time t, L is the length of the rod regardless of mass, and g is the acceleration of gravity, assuming a dimensionless case of L = g = 1 so that ω = 1. Large motion causes the nonlinearity, the degree of which can be adjusted by changing the initial conditions. In this case, the energy E total is conserved, which can be written as [37,38] 1 2 Considering the initial conditions θ 0 and . θ 0 , the equation can be expressed as: Letting θ 0 = 0, the relationship of Equation (35) can be rewritten as: where κ = 2ω/ . θ 0 and Ei (z,k) is the elliptic integral of the first kind. With the change in κ, the pendulum has two movement modes: if κ < 1, the pendulum will rotate around the center and the angle of rotation will keep increasing; if κ ≥ 1, the pendulum will oscillate around the center of rotation and the angle will be limited into a certain range and form a periodic change. The maximum angle becomes: Here we consider a highly nonlinear case; that is, when the speed of rotation is zero, the maximum oscillating angle reaches 179.9 • . According to the proposed equation, the initial condition is set to be . θ 0 = 1.999999238456499 such that the movement mode of the pendulum is close to the critical state. This case can be used to verify the validity of the algorithm. For the proposed method, the parameters are c 1 = c 3 = 1/4, α = 0, β = 0.5 and for the GGW method, the parameter is α = 0.
The corresponding results are shown in Figure 9 with ∆t/T = 10 −4 . It can be seen that the trapezoid rule, Wilson method, and Bathe method all give a result that the pendulum tends to rotate around the center, which indeed does not reflect the actual movement. Instead, the proposed method and GGW method generate numerical solutions which are very close to the exact solution. It can be also found that when the time step is ∆t = T/50000, as shown in Figure 10, the solution of the proposed method possesses high accuracy compared with the exact solution.
the trapezoid rule, Wilson method, and Bathe method all give a result that the pendulum tends to rotate around the center, which indeed does not reflect the actual movemen Instead, the proposed method and GGW method generate numerical solutions which ar very close to the exact solution. It can be also found that when the time step is Δt = T/50000 as shown in Figure 10, the solution of the proposed method possesses high accuracy com pared with the exact solution.   Table 2. We can find that th computational time of the proposed method, the GGW method, and the Bathe method ar close to each other, which is a little longer than that of Wilson method, since these thre methods are two-step ones. It takes the shortest time for computation with trapezoida rule, but as shown in Figure 11, no matter how small the Δt is, the trapezoidal rule canno obtain stable results. the trapezoid rule, Wilson method, and Bathe method all give a result that the pendulum tends to rotate around the center, which indeed does not reflect the actual movemen Instead, the proposed method and GGW method generate numerical solutions which ar very close to the exact solution. It can be also found that when the time step is Δt = T/5000 as shown in Figure 10, the solution of the proposed method possesses high accuracy com pared with the exact solution.   Table 2. We can find that th computational time of the proposed method, the GGW method, and the Bathe method ar close to each other, which is a little longer than that of Wilson method, since these thre methods are two-step ones. It takes the shortest time for computation with trapezoid rule, but as shown in Figure 11, no matter how small the Δt is, the trapezoidal rule canno obtain stable results. The computational time of these methods are listed in Table 2. We can find that the computational time of the proposed method, the GGW method, and the Bathe method are close to each other, which is a little longer than that of Wilson method, since these three methods are two-step ones. It takes the shortest time for computation with trapezoidal rule, but as shown in Figure 11, no matter how small the ∆t is, the trapezoidal rule cannot obtain stable results.

Long Span Planar Skeletal Structure
A 25 m (length) × 1 m (height) large span truss model is considered to demonstrate the numerical stability and dissipation of the proposed algorithm. The model shown in Figure 12 is composed of 101 member bars, and each bar is a single element with hinged constraints at both ends. As for the mechanical properties of the model, the elastic modulus is E = 2.0  10 5 MPa, the density is 7800 kg/m 3 , and the cross-section area of each bar is 100 mm 2 . An initial velocity v = 100 m/s is set on the points A and B located near the middle of the span, as shown in Figure 12.

Long Span Planar Skeletal Structure
A 25 m (length) × 1 m (height) large span truss model is considered to demonstrate the numerical stability and dissipation of the proposed algorithm. The model shown in Figure 12 is composed of 101 member bars, and each bar is a single element with hinged constraints at both ends. As for the mechanical properties of the model, the elastic modulus is E = 2.0 × 10 5 MPa, the density is 7800 kg/m 3 , and the cross-section area of each bar is 100 mm 2 . An initial velocity v = 100 m/s is set on the points A and B located near the middle of the span, as shown in Figure 12.

Long Span Planar Skeletal Structure
A 25 m (length) × 1 m (height) large span truss model is considered to demonstrate the numerical stability and dissipation of the proposed algorithm. The model shown in Figure 12 is composed of 101 member bars, and each bar is a single element with hinged constraints at both ends. As for the mechanical properties of the model, the elastic modulus is E = 2.0  10 5 MPa, the density is 7800 kg/m 3 , and the cross-section area of each bar is 100 mm 2 . An initial velocity v = 100 m/s is set on the points A and B located near the middle of the span, as shown in Figure 12.      As shown in Figure 14, the velocity response calculated by the GGW method is larger than that calculated by other methods, i.e., the GGW method does not filter the effects of high-frequency oscillations. In order to demonstrate this effect more clearly, the velocity response from 3 to 8 s in Figure 14 is replotted in Figure 16. The presented method has much improved the performance of the GGW method on high-frequency oscillation problems. Figure 16 also shows that after filtering out high-frequency effects, the velocity responses calculated by the proposed method, Bathe method and Wilson method oscillate between −10 and 10 m/s, while the response of the GGW method oscillates between −50 and 50 m/s.
To study the influence of parameter α on the results of the proposed method, setting the time step Δt = 0.001 s, the response of this example is computed with different α (α = As shown in Figure 14, the velocity response calculated by the GGW method is larger than that calculated by other methods, i.e., the GGW method does not filter the effects of high-frequency oscillations. In order to demonstrate this effect more clearly, the velocity response from 3 to 8 s in Figure 14 is replotted in Figure 16. The presented method has much improved the performance of the GGW method on high-frequency oscillation problems. Figure 16 also shows that after filtering out high-frequency effects, the velocity responses calculated by the proposed method, Bathe method and Wilson method oscillate between −10 and 10 m/s, while the response of the GGW method oscillates between −50 and 50 m/s. To study the influence of parameter α on the results of the proposed method, setting the time step ∆t = 0.001 s, the response of this example is computed with different α (α = 0.005, 0.1, 0.25, 0.5, 0.667) respectively, and the results are plotted in Figures 17-19. According to Figures 17-19, as α becomes larger, the high-frequency oscillation gradually weakens while the low-frequency response is well preserved. It demonstrates the superiority of the proposed method with controllable numerical damping. method, Bathe method, GGW method, and the proposed method.
As shown in Figure 14, the velocity response calculated by the GGW method is larger than that calculated by other methods, i.e., the GGW method does not filter the effects of high-frequency oscillations. In order to demonstrate this effect more clearly, the velocity response from 3 to 8 s in Figure 14 is replotted in Figure 16. The presented method has much improved the performance of the GGW method on high-frequency oscillation problems. Figure 16 also shows that after filtering out high-frequency effects, the velocity responses calculated by the proposed method, Bathe method and Wilson method oscillate between −10 and 10 m/s, while the response of the GGW method oscillates between −50 and 50 m/s.
To study the influence of parameter α on the results of the proposed method, setting the time step Δt = 0.001 s, the response of this example is computed with different α (α = 0.005, 0.1, 0.25, 0.5, 0.667) respectively, and the results are plotted in Figures 17-19. According to Figures 17-19, as α becomes larger, the high-frequency oscillation gradually weakens while the low-frequency response is well preserved. It demonstrates the superiority of the proposed method with controllable numerical damping.  To verify the decaying speed of the numerical solution with the proposed method, the time steps are set to be Δt = 0.001, 0.005, 0.01, and 0.02 s, respectively, and the corresponding results are shown in Figure 20 and Figure 21, where α = 2/7 so that the corresponding final spectral radius is 0.5, and low frequencies are retained while high frequencies are filtered out. To verify the decaying speed of the numerical solution with the proposed method, the time steps are set to be Δt = 0.001, 0.005, 0.01, and 0.02 s, respectively, and the corresponding results are shown in Figure 20 and Figure 21, where α = 2/7 so that the corresponding final spectral radius is 0.5, and low frequencies are retained while high frequen- To verify the decaying speed of the numerical solution with the proposed method, the time steps are set to be ∆t = 0.001, 0.005, 0.01, and 0.02 s, respectively, and the corresponding results are shown in Figures 20 and 21, where α = 2/7 so that the corresponding final spectral radius is 0.5, and low frequencies are retained while high frequencies are filtered out. To verify the decaying speed of the numerical solution with the proposed method the time steps are set to be Δt = 0.001, 0.005, 0.01, and 0.02 s, respectively, and the corre sponding results are shown in Figure 20 and Figure 21, where α = 2/7 so that the corre sponding final spectral radius is 0.5, and low frequencies are retained while high frequen cies are filtered out.  It can be found that when the time step increases, the high-frequency part of the nu merical solution gradually disappears while the low-frequency part is retained. This con clusion shows that the numerical result obtained from the proposed method graduall attenuates, but the accurate and effective low-frequency solution can still be obtained This result can help to improve the computational efficiency.
When filtering out high-frequency numerical solutions, it is also necessary to con sider the influence of numerical damping on low-frequency solutions. Due to the rapi numerical damping and amplitude decay, the low-frequency solution of some algorithm will gradually lose its value during calculation. When using the Wilson method, with th passage of time, the amplitude of numerical solution gradually decays. As for the GGW method, the numerical damping can be controlled by changing the parameter α, but it solution attenuates too fast, which will also cause numerical dissipation in the low-fre quency mode when filtering out the high-order modes.
As a significant improvement to the GGW method, the proposed method can effec tively avoid such a problem. To verify the superior accuracy of the proposed method, th responses are computed with the Wilson method, GGW method, and proposed metho when Δt = 0.01 s. The corresponding results are shown in Figure 22, where α = 0.05 for th GGW method and α = 2/7 for the proposed method. It can be found that when the time step increases, the high-frequency part of the numerical solution gradually disappears while the low-frequency part is retained. This conclusion shows that the numerical result obtained from the proposed method gradually attenuates, but the accurate and effective low-frequency solution can still be obtained. This result can help to improve the computational efficiency.
When filtering out high-frequency numerical solutions, it is also necessary to consider the influence of numerical damping on low-frequency solutions. Due to the rapid numerical damping and amplitude decay, the low-frequency solution of some algorithms will gradually lose its value during calculation. When using the Wilson method, with the passage of time, the amplitude of numerical solution gradually decays. As for the GGW method, the numerical damping can be controlled by changing the parameter α, but its solution attenuates too fast, which will also cause numerical dissipation in the low-frequency mode when filtering out the high-order modes.
As a significant improvement to the GGW method, the proposed method can effectively avoid such a problem. To verify the superior accuracy of the proposed method, the responses are computed with the Wilson method, GGW method, and proposed method when ∆t = 0.01 s. The corresponding results are shown in Figure 22, where α = 0.05 for the GGW method and α = 2/7 for the proposed method. method, the numerical damping can be controlled by changing the parameter α, b solution attenuates too fast, which will also cause numerical dissipation in the low quency mode when filtering out the high-order modes.
As a significant improvement to the GGW method, the proposed method can e tively avoid such a problem. To verify the superior accuracy of the proposed method responses are computed with the Wilson method, GGW method, and proposed me when Δt = 0.01 s. The corresponding results are shown in Figure 22, where α = 0.05 fo GGW method and α = 2/7 for the proposed method. It can be found from Figure 22 that the accuracy of the proposed method is b than both the GGW method and the Wilson method with the same time step. This i cause of the unnecessary numerical dissipation of these two methods. Using a sm It can be found from Figure 22 that the accuracy of the proposed method is better than both the GGW method and the Wilson method with the same time step. This is because of the unnecessary numerical dissipation of these two methods. Using a smaller time step can improve the accuracy of these two methods, but this will increase the cost of computation.
To compare the efficiency of the proposed method and the Bathe method, the displacement responses are computed with different time steps: ∆t = 0.03, 0.02, and 0.01 s. Letting α = 1/7, the results are shown in Figure 23 and the computational time versus different ∆t is listed in Table 3.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 2 time step can improve the accuracy of these two methods, but this will increase the of computation.
To compare the efficiency of the proposed method and the Bathe method, the placement responses are computed with different time steps: Δt = 0.03, 0.02, and 0 Letting α = 1/7, the results are shown in Figure 23 and the computational time versu ferent Δt is listed in Table 3. As shown in Table 3, the computational time of the Bathe method is shorter tha proposed method when the same time step is used. However, we can find from Figu that the accuracy of the proposed method is better than the Bathe method with the time step. That is to say, under the conditions of the same accuracy, the proposed me calculates faster with larger time step than Bathe method with smaller time step. Fo ample, in Figure 23, the Bathe method uses 10.186s with Δt = 0.01 s to obtain results o same accuracy as the proposed method, which takes 5.039s with Δt = 0.03 s. From point of view, the proposed method is sufficiently efficient in both the computational and accuracy.   Table 3, the computational time of the Bathe method is shorter than the proposed method when the same time step is used. However, we can find from Figure 23 that the accuracy of the proposed method is better than the Bathe method with the same time step. That is to say, under the conditions of the same accuracy, the proposed method calculates faster with larger time step than Bathe method with smaller time step. For example, in Figure 23, the Bathe method uses 10.186s with ∆t = 0.01 s to obtain results of the same accuracy as the proposed method, which takes 5.039s with ∆t = 0.03 s. From this point of view, the proposed method is sufficiently efficient in both the computational time and accuracy.

Conclusions
This paper proposes a new implicit time integration method based on the GGW method. To control the numerical damping effectively, the time step is divided into two substeps. In the first step, the GGW method is used, and in the second step, a new parameter is introduced according to the idea of trapezoidal integration. The numerical analysis of numerical stability and accuracy shows that spectral radius, period extension rate, and amplitude attenuation rate of this method can all be controlled by changing the parameters α and β. Generally, we take β = 1/2 and α ∈ [2/19, 2/7], thus ρ ∞ ∈ [0.5, 0.8]. Several numerical examples are given to demonstrate the effectiveness of the proposed method.
The first representative example is an undamped single-DOF problem. The results show that the proposed method has better accuracy, smaller period extension, and smaller amplitude attenuation compared with other well-known time integration methods. The second example is a large displacement nonlinear problem. When using the same time step, the Bathe method, trapezoidal rule, and Wilson method cannot compute reasonable results, while the results of the proposed method are very similar to the exact solutions, which proves that the method is suitable for nonlinear problems. The third example is a large-span planar structure problem. The results show that the proposed method can not only possess high dissipation in a high-frequency range, but also preserve good accuracy in a low-frequency range. Hence, a larger time step can be chosen to solve such problems.
In summary, the proposed method has the following merits: 1.
Reasonable and controllable numerical damping. It possesses high dissipation in a high-frequency range and preserves good accuracy in a low-frequency range as well; 2.
Applicability of both linear and nonlinear problems. This method can maintain high accuracy in large-displacement nonlinear problems; 3.
Due to the advantages mentioned above, a stable and reasonable numerical result can still be obtained when a larger time step is selected, which results in the superior efficiency.
The proposed method is expected to have a bright future in nonlinear vibration problems.

Conflicts of Interest:
The authors declare no conflict of interest.