A re-evaluation of overshooting in time integration schemes: The neglected effect of physical damping in the starting procedure

An important property of implicit time integration algorithms for structural dynamics is their tendency to “overshoot” the exact solution in the first few steps of the computed response due to high-frequency components in the initial excitations. The typical analysis technique for overshooting involves the study of asymptotic response of the algorithm’s first step in the limiting high frequency case. This article finds that the prior analysis of overshooting in much of the engineering literature is incomplete in that it neglects the effect of physical damping. With physical damping included, first-order overshooting components enter into several well-known time integration algorithms which were previously thought to exhibit zero-order overshooting in displacement. The Newmark method, Wilson- 𝜃 method, Bazzi- 𝜌 method, HHT- 𝛼 method, WBZ- 𝛼 method, and three parameter optimal/generalized- 𝛼 method are analyzed, as well as the generalized single-step single-solve (GSSSS) framework which encompasses all of the prior schemes and other new and optimal algorithms and designs based upon the issues under consideration. The additional overshooting component is eliminated in the novel amended GSSSS V0 family (which is noteworthy and cannot be derived by conventional means), while the numerically dissipative schemes in the GSSSS U0 family (encompassing traditional methods such as the HHT- 𝛼 method, WBZ- 𝛼 method, and three parameter optimal/generalized- 𝛼 method among other new algorithm designs) are shown to be irremediable as the additional overshooting component from


INTRODUCTION
Many time integration or time-stepping schemes have been proposed for the numerical solution of second-order dynamical systems such as arising from the spatial semi-discretization of the equations of structural dynamics (the reader can refer to a series of papers by Zhou, Tamma et al. [1][2][3][4][5] detailing a unified theory regarding the development and analysis of algorithms for time-dependent problems). Of these, the most popular and widely used are the single-step (SS) 3-value methods which use the displacement, velocity, and acceleration values at one time point to compute the values at the next time point. These SS methods admit equivalent linear multistep (LMS) representations by way of the Cayley-Hamilton theorem, and as such any SS method actually pertains to the class of LMS methods. Thus, the Lax equivalence theorem can be straightforwardly employed to prove the convergence of any SS method to the exact solution as well as its order of accuracy. This involves a proof of consistency, which requires that the exact solution approximately satisfies the discrete time integration scheme to some order, as well as a proof of stability, which requires that the errors introduced by the discrete time integration scheme do not grow without bound. A consistent LMS method is convergent if and only if it is stable. The consistency of an LMS method can be analyzed using the local truncation error, the order of which can be used to determine the order of accuracy if the method converges. The stability can be analyzed using the method's characteristic polynomial, and a method is conditionally stable if it is stable for a finite range of the time step size, and unconditionally stable if it is stable for all time step sizes. However, consistency order and stability are not the only properties that impinge on the design of a successful time integration scheme. The computational complexity of a scheme is important, particularly the number and size of implicit systems of equations that must be solved in each time step, as well as the starting procedure that is used to begin the time-stepping process. Dahlquist 6 proved that there can be no explicit and unconditionally stable LMS methods for 1st-order differential equation systems, and that the highest consistency order attainable by an unconditionally stable LMS method is 2. This "Dahlquist barrier theorem" can be straightforwardly extended to LMS methods (and hence SS methods due to their equivalent LMS representations) for second-order differential equation systems as well. 7 In addition, several analysis techniques have been proposed for the study and verification of other numerical properties of SS methods for second-order systems, such as numerical dissipation (amplitude distortion), numerical dispersion (period error), and a special overshooting property of SS methods which is of particular note here. Much of the development of SS methods for second-order systems has been centered on methods with controllable numerical dissipation of high frequency modes which is capable of attenuating spurious numerical oscillations due to high-frequency errors that naturally occur due to spatial semi-discretization of the equations of structural dynamics. One early example was the Wilsonmethod, 8 which is implicit and unconditionally stable with second-order accuracy and featured controllable numerical dissipation of the high-frequency modes. However, Goodreau and Taylor 9 discovered that the Wilson-method tended to significantly overshoot the amplitude of the exact solution to a free vibration problem as the time step size was increased despite its unconditional stability property, which should maintain bounded solutions for all time step sizes. Hilber and Hughes 10 further elaborated on the nature of this overshooting, revealing that it occurs principally due to the first step of the time-stepping procedure and proposing an analysis technique to detect and quantify overshooting in SS methods via an asymptotic representative solution in the high-frequency limit. In the time since, overshooting analysis has become a staple in papers proposing and analyzing unconditionally stable implicit SS methods, 4,[11][12][13][14][15] with zero-order overshooting in displacement and/or velocity seen as a desirable property for a time integration scheme. However, most of these papers neglect physical damping in the overshooting analysis, likely because its influence was omitted in the original paper by Hilber and Hughes, despite appearing in Hilber's PhD thesis dissertation 16 2 years prior. Furthermore, the analysis technique employed by Zhou and Tamma 3 to prove the overshooting properties of several families of time integration schemes was sufficient only for physically undamped problems, in which case it produces representative solutions with the proper order of overshoot but still incorrect asymptotic behavior in the high-frequency limit.
This article re-examines the effect of physical damping on the overshooting behavior of several well-known SS methods, including Newmark's method, 17 the Wilson-method, 8 and Bazzi-method, 11 and in particular demonstrates that the numerically dissipative Hilber-Hughes-method, 18 Wood-Bossak-Zienkiewicz-method, 19 and three parameter optimal/generalized-method, 15,20,21 which have been claimed to exhibit zero-order overshooting in displacement and first-order overshooting in velocity, actually exhibit first-order overshooting in displacement in the presence of damping. It then extends the analysis to the generalized single-step single-solve (GSSSS) framework of algorithms, including the prior three algorithms in the U0 family as well as their counterparts in the V0 family. First, the structure of a SS method and its basic components are examined. A distinction is made between the integration procedure (defined by the LMS representation) which determines the consistency and stability of the time integration scheme, and the starting procedure (defined by the SS representation) which determines the overshooting behavior. Next, the first-order displacement overshooting behavior of the selected algorithms under physically damped situations is proven analytically by the proper extension of the asymptotic representative solution to the physically damped case, including a novel perspective on overshooting that occurs in the second step of the time integration. Modal decomposition is elaborated to extend this analysis to multi-degree of freedom systems, with practical considerations for finite element discretizations. The analytical proofs are verified numerically with a single-degree of freedom spring-mass-dashpot system undergoing free vibration due to initial conditions and a multi-degree of freedom beam vibration problem with Rayleigh damping, with implications for practical problems in structural dynamics. The conclusion of the study recommends additional treatment of the damping matrix in the development of numerically dissipative time integration schemes and suggests optimal algorithms for problems with various combinations of initial conditions and physical damping.

SINGLE-STEP METHODS AND LMS REPRESENTATIONS
After spatial semi-discretization, the initial value problem governing linear structural dynamics can be written as in Equation (1) Mü + Cu + Ku = F(t), The single-step (SS) 3-value methods are commonly employed to numerically approximate the solution to Equation (1). One of the most well-known methods of this type is Newmark's method as given in Equation (2), which is implicit, unconditionally stable, and second-order accurate with the average acceleration parameters = 1 2 and = 1 4 .
The typical linear analyses of SS methods apply modal analysis concepts in which properties such as accuracy, stability, and overshooting are analyzed for the single-degree of freedom problem given by Equation (3), with natural frequency n = √ k∕m and damping ratio = c∕(2 √ km) conforming to the case where Equation (1) only has one degree of freedom. With these properties known for the SDOF case, it is straightforward to extrapolate to the MDOF case via modal superposition (this is described in further detail in Section 5.2).
When applied to the single-degree of freedom case, an SS method such as Newmark's method can be rewritten in the form of Equation (4), where A is a 3 × 3 amplification matrix and L n is a 3 × 1 load vector.
A straightforward application of the Cayley-Hamilton theorem allows for the SS method to be recast into an equivalent LMS representation as in Equation (5) where A 1 , A 2 , and A 3 are the principal invariants defining the characteristic polynomial of A as in Equation (6).
The local truncation error can be derived by substitution of the exact solution into Equation (5) to prove consistency and order of accuracy. The roots of Equation (6) can be shown to have modulus less than or equal to unity to prove stability, and thus ensure that discrete solutions remain bounded. By the Lax equivalence theorem, this is necessary and sufficient to prove that the integration procedure as defined by the LMS representation converges to the exact solution as Δt → 0. Moreover, an unconditionally stable (implicit) scheme will have bounded solutions in the limit Δt → ∞ as well. However, any LMS integration procedure must be equipped with a distinct starting procedure to provide the initial backstep values (if A 1 , A 2 , and A 3 are all nonzero, then at least two backstep values must be given). Hence, the SS representation as in Equation (4) defines the method's starting procedure, while the LMS representation as in Equation (5) defines the integration procedure. Note that only the spectral eigenvalue properties of A play a role in the homogeneous integration procedure, while the full matrix structure of A plays a role in the homogeneous starting procedure. Special care must be taken in the design of SS methods to ensure not only that the integration procedure is consistent and stable, but also that the starting procedure does not introduce unacceptable errors which propagate into the rest of the solution-the starting values must possess the same order of accuracy as the integration procedure to maintain the same accuracy globally (see Zhou and Tamma's 4 discussion on the accuracy of the primary variables).

OVERSHOOTING ANALYSIS
We have discussed how both the starting procedure and integration procedure of an algorithm affect its order of convergence in the limit as Δt → 0, as well as how the unconditional spectral stability of an implicit algorithm's integration procedure ensures that its solutions remain bounded even for Δt → ∞. However, we have left out the effect of the starting procedure in the limit as Δt → ∞. Presently, the effect of the homogeneous terms of the starting procedure in the high-frequency limit is understood in the literature by the concept of overshooting. The homogeneous integration procedure is determined purely by the characteristic equation of a method, which for an unconditionally stable method has roots of modulus less than or equal to unity. This guarantees that initial integration errors will not grow without bound. However, it is possible for the starting procedure to introduce errors larger than O(1) into the solution. Goodreau and Taylor 9 discovered that the second-order time accurate and unconditionally stable Wilson-method tended to significantly overshoot the amplitude of the exact solution to a free vibration problem. Hilber and Hughes 10 further elaborated on the nature of this overshooting, revealing that it occurs principally due to the first step of the analysis, and can be measured by the full matrix norm of the amplification matrix, which can generally be larger than its spectral radius. Equivalently, high-frequency overshooting can simply be detected by applying the entries of the amplification matrix to the homogeneous initial conditions u 0 , v 0 , and a 0 = m −1 (−cv 0 − ku 0 ) = −2 n v 0 − 2 n u 0 for modal frequency tending to infinity, Ω = n Δt → ∞ to obtain a representative solution. For a 3-value SS method, the representative solution for the limiting high-frequency case is obtained as in Equation (7), and the maximum order in Ω and/or Δt of coefficients of the solution determines the order of overshoot for displacement and velocity.
where y n = [ u n v n Δt a n Δt 2 ] T .
In order to assess the "order" of overshooting associated with each dependent variable, the representative solution for n = 1 or n = 2 can be written as in Equation (8), where p is the order of displacement overshoot due to initial displacement, q is the order of displacement overshoot due to initial velocity, r is the order of velocity overshoot due to initial displacement, and s is the order of velocity overshoot due to initial velocity.
This representative solution in Equation (7) differs from the one advocated by all previous papers on overshooting, as it includes the second step of the starting procedure in addition to the first, which is generally required for 3-value SS methods with a nonzero spurious root. In Section 3.5, we will demonstrate how this can be critical in assessing the overshooting behavior of an algorithm in the presence of physical damping, as the second step can introduce overshooting components that were not present in the first step. To be precise, if the first step exhibits zero-order overshooting in displacement and first-order overshooting in velocity, it is possible for the first-order overshoot in velocity to be propagated into the second step displacement due to physical damping, as the O( Ω −1 ) term in the (1, 2) position of the 3 × 3 amplification matrix will act on v ∞ 1 Δt ∼ O(Ω 2 ) to produce u ∞ 2 ∼ O( Ω) uniquely in the second step. In the interest of brevity, we will only present results for the first step analysis for the schemes that follow except in Section 3.5 on the GSSSS framework of algorithms, which is capable of encapsulating all of the other cases and for which the second step overshooting is a salient design feature.
It is worth noting that more than zero-order overshooting is only possible for displacement and velocity in 3-value SS methods (the maximum orders of overshooting are (1)). Hence, acceleration overshooting is not considered here, but for general k-value methods with unconditional stability (k > 3 is rare in the literature as Dahlquist's barrier theorem limits the convergence order of an implicit unconditionally stable algorithm to 2), higher order derivatives should also be assessed for overshooting that could occur due to the starting procedure.

Newmark's method
When the technique prescribed in Section 3 is employed to study the overshooting properties of Newmark's method in Equation (2), the representative solution for the first step is obtained as in Equation (9).
In general, Newmark's method exhibits first-order overshoot in displacement and first-order overshoot in velocity in the presence of damping, reducing to zero-order overshoot in displacement and first-order overshoot in velocity for physically undamped problems. Additionally, if average acceleration parameters = 1 2 and = 1 4 are adopted, then Newmark's method exhibits zero-order overshoot in both displacement and velocity with or without physical damping. This makes Newmark's average acceleration method very competitive, as it is unconditionally stable, second-order accurate, and features zero-order overshoot in all dependent variables for all problems regardless of physical damping. However, it does not possess controllable numerical dissipation, the introduction of which imposes a tradeoff on the order of overshooting as proven by Zhou and Tamma. 4 Hence, in the sections that follow we will revisit several numerically dissipative time integration schemes appearing in the literature whose overshoot properties have been incompletely understood.

The Wilson-method
As an example, we consider the Wilson-method, 8 which was derived by extrapolating the integrator equation to a collocation point, satisfying the equation of motion at n + instead of at n + 1, where ≥ 1. This endows the method with controllable numerical dissipation of high frequency content (i.e., the algorithm's spectral radius at Ω → ∞ is less than unity). The Wilson-method can be written as in Equation (10).
When the analysis technique prescribed in Section 3 is employed to study the overshooting properties of the Wilsonmethod, the representative solution for the first step is obtained as in Equation (11) Hence, the Wilson-method exhibits second-order overshoot in displacement and first-order overshoot in velocity. If there is no physical damping component, this overshoot is due only to the initial displacement, reducing to first-order overshoot in displacement and zero-order overshoot in velocity for pure initial velocity problems. However, in the presence of physical damping, the Wilson-method exhibits second-order overshoot in displacement and first-order overshoot regardless of the initial conditions. In general, the Wilson-method introduces large errors into the first step of the analysis where high frequency modes are concerned, which is undesirable as it defeats the purpose of high-frequency numerical dissipation. Hilber and Hughes 10 reached the same conclusion as the above in the physically undamped case, proving that a generalized collocation scheme with controllable numerical dissipation would be unable to escape second-order overshoot in displacement and first-order overshoot in velocity, thus ruling out collocation as an effective procedure to derive time integration schemes with desirable overshooting properties. When the analysis of this generalized collocation scheme is extended to the physically damped case, the overshooting properties remain identical to those given above for the Wilson-method. Thus, we look to other numerically dissipative algorithm designs to assess their overshooting properties, particularly noting the influence of the physical damping coefficient.

The Bazzi-method
Bazzi and Anderheggen 11 developed the unconditionally stable Bazzi-method with controllable numerical dissipation as given by Equation (12) (M + The Bazzi-method was originally claimed to be second-order accurate for physically undamped problems by way of its homogeneous local truncation error, but Zhou and Tamma 4 pointed out that it actually suffers from degraded global accuracy due to its starting procedure, which is only first-order accurate. Here, we further point out that the original overshooting analysis was also performed incorrectly. When the analysis technique prescribed in Section 3 is employed to study the overshooting properties of the Bazzi-method, the representative solution for the first step is obtained as in Equation (13).
The Bazzi-method actually exhibits zero-order overshoot in displacement and first-order overshoot in velocity in general, where it was originally claimed to exhibit first-order overshoot in velocity only in the damped case (note that both displacement and velocity exhibit zero-order overshooting if = 1, which leads to the numerically non-dissipative classical midpoint rule method with endpoint acceleration. These results are maintained in the second step analysis, which is omitted here). In place of Equation (7), Bazzi et al. considered only the limiting behavior of the amplification matrix as applied to arbitrary values of displacement, velocity, and acceleration, ignoring the initial acceleration term a 0 Δt 2 , in which u 0 introduces a potential O(Ω 2 ) component and v 0 introduces a potential O( Ω 2 ) component. This omission illustrates the importance of correctly applying the initial conditions during overshooting analysis, as they are what introduce high-frequency components into the starting procedure. While the Bazzi-method is an improvement on the Wilsonmethod in terms of overshooting behavior, its first-order accuracy is a major drawback to its competitiveness. In the following sections, we discuss families of numerically dissipative methods which preserve second-order accuracy while managing to attain similarly improved overshooting properties.

The HHT-, WBZ-, and three-parameter optimal/generalized-method
Hilber et al. 18 developed the HHT-method via a modification of Newmark's integrator equation that introduces a convex combination on the stiffness term (in subsequent iterations, 22,23 the damping and force terms as well) between the time points t n and t n+1 . Wood et al. 19 developed the WBZ-method by similarly introducing a convex combination on the inertial term in Newmark's integrator equation. In both methods, the Newmark parameters , , and a new parameter are re-evaluated to form one-parameter families of algorithms with second-order accuracy, unconditional stability, and controllable numerical dissipation. Later, Shao and Cai 20,21 developed the three-parameter optimal (TPO) scheme, which combines the two approaches under a single family of algorithms with optimal numerical dissipation. The same approach is taken by Chung and Hulbert 15 to develop the generalized-method, which is actually identical to the original and first developed TPO method and thus we use TPO/G-to refer to the family of algorithms henceforth. Since HHT-and WBZare contained in its structure, an overshooting analysis of the TPO/G-method with arbitrary parameters (Newmark's and , HHT-'s f , and WBZ-'s m ) suffices for all three. The TPO/G-method can be written as in Equation (14) M When the analysis technique prescribed in Section 3 is employed to study the overshooting properties of the TPO/Gmethod, the representative solution for the first step is obtained as in Equation (15).
In the special subcase where Newmark average acceleration parameters = 1 2 , = 1 4 , m = f = 0 are adopted, Equation (15) obtains the same zero-order overshoot in displacement and velocity as described in Section 3. For more general algorithm parameters, the TPO/G-method (along with its subcases in HHT-and WBZ-) exhibits zero-order overshoot in displacement and first-order overshoot in velocity for physically undamped problems. This is line with analysis presented by Hilber et al. 18 and Zhou and Tamma. 4 However, in the presence of physical damping the TPO/Gmethod exhibits first-order overshoot in both displacement and velocity. While Hilber made note of this fact for the HHTmethod in his PhD thesis dissertation 16 (which Bazzi and Anderheggen 11 also briefly reference in their own overshooting analysis), physical damping is neglected in the majority of the overshooting analyses of the TPO/G-methods and its subcases-Hilber et al. 18 did not include physical damping in their paper introducing overshooting analysis; Chung and Hulbert 15 simply assert that the TPO/G-method does not exhibit overshooting in displacement (which can be assumed to mean that it exhibits only zero-order overshooting); and Zhou and Tamma 4 applied an incorrect modification of the analysis technique presented in Section 3 which will be discussed in the next section.

The GSSSS family of algorithms
The GSSSS framework of algorithms 4,24 consists of second-order accurate algorithms which are implicit and unconditionally stable and are designed for optimal high-frequency dissipation and overshooting properties. The GSSSS "U0" family encompasses many of the prior traditional methods, such as the Newmark average acceleration method, classical midpoint rule method, HHT-method, WBZ-method, and TPO/generalized-method, among many others which were new and novel, which all exhibit zero-order overshooting in displacement in the physically undamped case and no more than first-order overshooting in velocity. In addition, the GSSSS "V0" family includes counterparts (or competitors based on initial conditions) of every U0 algorithm, possessing the same spectral characteristics but exhibiting zero-order overshooting in velocity and no more than first-order overshooting in displacement. In this section, we will study the design of second-order accurate schemes in the GSSSS framework for overshooting and discuss the resulting numerical behavior in the high-frequency limit. The GSSSS framework can be written as in Equation (16) Integrator Updates u n+1 = u n + 1 v n Δt + 2 a n Δt 2 + 3 ΔaΔt 2 , v n+1 = v n + 4 a n Δt + 5 ΔaΔt, a n+1 = a n + Δa, where W i Λ j and k are scalar coefficients. The choice of these coefficients determines the numerical properties of the resulting algorithm, such as accuracy, stability, overshooting, and numerical dissipation. The parametric constraints given in Equation (17) are imposed to guarantee second-order accuracy of the starting and integration procedure.
When the analysis technique prescribed in Section 3 is employed to study the overshooting properties of the GSSSS framework, the representative solution for the first step is obtained as in Equation (18).
It should first be noted that this representative solution differs from the one given by Zhou and Tamma. 4 Rather than evaluating the limiting behavior of the entire first step as in Equation (7), Zhou and Tamma obtained a representative solution for the GSSSS framework by evaluating the limiting behavior of the amplification matrix before applying it to the initial conditions, which has the consequence of eliminating O( Ω −1 ) terms in the amplification matrix which can be magnified by the O(Ω 2 ) and O( Ω 2 ) components in the initial acceleration term and thus produce additional overshooting with physical damping. Based on this simplification, Zhou and Tamma proposed the following constraints to control overshooting: The U0 family employs conditions 1 and 2, resulting in algorithms with zero-order overshoot in displacement and at most first-order overshoot in velocity in the physically undamped case. The V0 family employs conditions 1 and 3, resulting in algorithms with zero-order overshoot in velocity and at most first-order overshoot in displacement in the physically undamped case. There also exist algorithms satisfying conditions 1-3, dubbed the U0V0 algorithms, which can feature no numerical dissipation of the high frequencies, while the so-called "U0V1" algorithms and their "V0U1" counterparts have controllable numerical dissipation of high frequency components. The optimal U0 and V0 families of algorithms are given as in Equation (19) according to the dissipation parameters 0 ≤ s ∞ ≤ min ∞ ≤ max ∞ ≤ 1 which determine the absolute value of the three spectral roots of the amplification matrix at Ω = n Δt → ∞. Each algorithm can thus be compactly designated as GSSSS U0( min ∞ , max ∞ , s ∞ ) or GSSSS V0( min ∞ , max ∞ , s ∞ ).

U0 family of algorithms
, , , V0 family of algorithms As is clear from Equation (18), conditions 1-3 are not sufficient to control overshooting in the physically damped case, where additional components enter into the analysis. The U0 family of algorithms evidently suffers from first-order overshoot in displacement in the physically damped case unless W 2 Λ 5 = 2W 1 3 , and the V0 family of algorithms will similarly exhibit first-order overshoot for initial displacement problems that would otherwise be avoided in the physically undamped case. Since the W 2 Λ 5 parameter was constrained by Zhou and Tamma 4 via the arbitrary condition that the approximation functions equate to the update functions (which does not lead to any particular desirable numerical property), it would seem possible to apply the new condition as a simple fix, stated below.
However, an unforeseen problem arises when the U0 family is amended in this manner-while the first step has its displacement overshoot under control, the second step of the analysis still produces first-order overshoot in displacement (recall from Section 2 that the starting procedure for a 3-value SS method consists of the first two steps) when physical damping is present. The representative solution for the second step with the original U0 conditions and W 2 Λ 5 = 2W 1 3 condition applied is given by Equation (20) u Evidently, even though first-order overshoot in displacement can be eliminated from the first step of the amended U0 family in the presence of physical damping by way of an additional constraint on W 2 Λ 5 , the second step still exhibits first-order overshoot in displacement, owing directly to the first-order overshoot in velocity that was present in the first step. The only two possible constraints that could eliminate the first-order overshoot in displacement with damping are given as follows: It turns out that condition A is only possible for the one scheme in the U0 family, namely the numerically non-dissipative classical midpoint rule, while condition B is only possible for algorithms in the U0 family which satisfy not just conditions 1 and 2 for zero-order displacement in the undamped case, but also condition 3 for zero-order overshoot in velocity (the numerically non-dissipative U0V0 algorithms). Hence, there is no numerically dissipative algorithm in the U0 family which exhibits zero-order overshoot in displacement and first-order overshoot in velocity for physically damped problems. The only possibilities in the U0 family in the presence of damping are zero-order overshoot in both displacement and velocity (the numerically non-dissipative U0V0 algorithms), or first-order overshoot in both displacement and velocity (the numerically dissipative U0 algorithms).
Remark 1. The traditional first step overshooting analysis of a 3-value SS method is sufficient to establish the overshooting behavior in the physically undamped case, as well as to reveal overshooting components that arise due to the first step in the presence of physical damping. However, an algorithm that produces zero-order overshoot in displacement but first-order overshooting in velocity in the first step may still be capable of producing first-order overshoot in displacement in the second step. Hence, the second step must also be analyzed in the physically damped case to determine the leading overshooting behavior in displacement for such algorithms.
The V0 family is capable of preserving a special property in the physically damped case when condition 4 is additionally satisfied. Namely, for problems involving only initial displacements, this amended V0 family exhibits zero-order overshoot in both displacement and velocity. Thus, the amended V0 family possesses a distinct advantage for initial displacement problems, whether or not physical damping is present (note that while the U0 family can similarly exhibit zero-order overshooting in both displacement and velocity for physically undamped initial velocity problems, this advantage disappears in the presence of physical damping). The amended V0 family is given as in Equation (21), and will be henceforth abbreviated as the V0* family.

Amended V0* family of algorithms
, The overshooting behavior of algorithms in the GSSSS framework for damped and undamped problems with initial displacement and/or velocity is summarized in Table 1. The optimal choices in terms of zero-order overshooting in both displacement and velocity are highlighted in bold.

EXTENSION TO MULTI-DEGREE OF FREEDOM SYSTEMS
The preceding analysis has centered on single-degree of freedom problems. This greatly simplifies the analysis, but is not a realistic depiction of overshooting behavior-the time step will always be chosen to resolve the frequency at hand in a practical SDOF problem. Hence, it is important to use modal analysis techniques to extend the SDOF analysis to MDOF systems. While usually considered trivial, it is elaborated in depth here, as the implications of overshooting on MDOF structures (such as arising from finite element discretizations of continuum bodies) has been only discussed obliquely in papers that mention it. Hence, we seek to offer a rigorous perspective on the consequences of an algorithm's starting procedure for practical MDOF problems.

TA B L E 1 GSSSS framework overshooting behavior
Physical damping? Initial conditions Algorithm family Overshooting behavior First, for simplicity we assume that the initial value problem in Equation (1) uses Rayleigh (i.e., proportional) damping such that C = M + K, as this allows the problem to be decomposed into N orthogonal modes (where N is the number of DOFs) as in Equation (22) By mode orthogonality, we have that T M = I and T K = , enabling decomposition into N independent modal coordinates with associated mode shapes. The initial conditions can be mapped onto the modal basis by (0) = T u(0) anḋ(0) = Tu (0), leading to a separate initial value problem for each mode as in Equation (23).
For free vibration problems, the system response u(t) = (t) can be viewed as a collection of the responses of each mode to initial excitation-that is, each mode participates according to the excitation of each mode shape. If the mode shapes corresponding to low frequency modes are activated, then only a low-frequency response will be observed. However, if there is high-frequency content in the initial conditions, there will be high-frequency content in the solution.
Implicit methods are typically employed for problems in inertial dynamics in which only a few low-frequency modes dominate the solution, as the unconditional stability property of such methods allow for the selection of a large time step size which is capable of resolving those lower frequencies but which exceeds the period of the maximum system frequency. By contrast, explicit methods are typically employed for wave propagation problems in which high frequencies significantly participate in the solution, as a smaller time step size below the critical stability limit is needed to resolve the dynamic response anyway, and it is much more efficient to perform explicit calculations at each time step. It is well understood that the spatial discretization of a continuum body by, for example, the finite element method introduces errors in the frequency spectrum, particularly for the largest frequencies. This leads to spurious high-frequency participation even for inertial dynamics problems, resulting in noisy oscillations in the numerical solution. Hence, controllable numerical dissipation to attenuate high frequencies is seen as a desirable property for a time integration algorithm. However, as we have seen in Section 3, numerically dissipative single-step algorithms with second-order accuracy suffer from high-frequency overshooting (depending on the initial conditions). In particular, this overshooting comes from high frequency content in the initial condition that the chosen time step size is incapable of resolving, resulting in the Ω → ∞ type behavior for high frequency modes that can exceed the magnitude of the initial excitation. Hence, while the chosen step size may be sufficient to resolve the dominant low-frequency response, high-frequency errors from spatial discretization can induce large errors in the computed solution. This does not pose a difficulty for algorithmic convergence; as a consequence of algorithmic consistency, an infinitesimal time step size tending to Δt → 0 is capable of resolving any finite frequency content, no matter how large (and while driving the spatial mesh size to zero will tend toward an infinite frequency spectrum, consistency of the spatial discretization method ensures that the error in the associated frequency content approaches zero as well). Rather, these large errors are encountered in the practical application of implicit time integration schemes to semi-discrete structural dynamics problems for finite time step sizes.
This leads to a conundrum in which spurious high-frequency content entails the use of numerically dissipative algorithms, but these in turn induce large errors proportional to the highest frequency. The fix is appropriate algorithm selection-in many cases, the algorithm can be matched to the initial conditions via Table 1 so as to eliminate overshooting. Assuming no physical damping is present, the GSSSS U0 algorithms (including the dissipative HHT-, WBZ-, and TPO/G-families) should be matched with initial velocity problems; the GSSSS V0* algorithms (containing counterparts to the listed dissipative algorithms) should be matched with initial displacement problems; and a tradeoff between numerical dissipation and overshooting is required for problems with both initial displacement and velocity (the GSSSS U0V0 algorithms do not overshoot but cannot attenuate spurious high frequencies). If physical damping is present, the GSSSS U0 algorithms lose this comparative advantage and either GSSSS V0* or GSSSS U0V0 methods may be generally preferred.
Rayleigh damping is usually designed to achieve a desired damping ratio over a given frequency range. As the damping ratio varies for each mode by i = 2 i + 1 2 i , it is important to select the coefficients and so that some frequencies are not excessively damped. Moreover, the stiffness proportional component of damping can lead to large O( Ω) and O( 2 Ω) overshooting components if the damping ratio is not properly controlled. A typical recommendation is given by Equation (24), which ensures that the damping ratio never exceeds desired (other schema may be necessary in the case that the minimum and maximum frequencies cannot be easily estimated). In general, the user must take care when using numerically dissipative algorithms alongside Rayleigh damping to avoid the excessive overshooting behavior that this combination can produce. Similar conclusions hold for other damping models in which the modes of vibration cannot be fully decoupled, but this analysis is out of scope for the present work.

NUMERICAL EXAMPLES
In the sections that follow, we numerically verify the overshooting properties for various selected algorithms in the dynamic response of SDOF and MDOF systems.

Overshoot response of spring-mass-dashpot SDOF system
First we examine a simple spring-mass-dashpot system with a single degree of freedom to numerically verify the overshooting behavior of selected algorithms with and without physical damping. The material parameters are m = 1 kg, k = 10 4 N/m leading to a natural frequency of n = 100 rad/s, and a damping ratio of = 0.5 is applied in the physically damped case. The time step size is chosen as Δt = 0.5 s, leading to Ω = n Δt = 50, which is sufficient to approach limiting high-frequency behavior as the time step size is nearly 8 times larger than the natural period of vibration. Initial displacements of u 0 = 1.0 m and initial velocities of v 0 = (1.0 m) n are combined to produce three test cases for both damping situations.
To clearly illustrate the effect of the second step on the damped overshooting response, the response of the SDOF spring-mass-dashpot system is first depicted in Figure 1 for the GSSSS U0V0(1.0, 1.0, 1.0) algorithm (a.k.a. classical midpoint rule method with endpoint acceleration), the GSSSS U0(0.5, 0.5, 0.5) algorithm (a.k.a. TPO/G-with = 0.5), its GSSSS V0 counterpart, and the amended GSSSS U0* and GSSSS V0* versions of these algorithms in which W 2 Λ 5 = 2W 1 3 . It can be seen that the amendment to the U0 family is incapable of removing the first-order displacement overshooting that occurs in the second step for physically damped problems, while the amendment to the V0 family results in zero-order displacement overshooting for initial displacement problems as expected.
Next, we numerically verify the overshooting behavior for a selection of algorithms in Section 3; Newmark's average acceleration method (i.e., GSSSS U0(1, 1, 0) which belongs to the U0V0 family), the GSSSS U0(0.5, 0.5, 0.5) method (which is equivalent to the TPO/G-method with = 0.5), its counterpart GSSSS V0*(0.5, 0.5, 0.5) in the amended V0 family, the Bazzi-method with = 0.5, and the Wilson-method with = 1.4. The physically undamped results are shown in Figure 2, while the physically damped results are shown in Figure 3. The velocity results are scaled by the natural frequency to enable straightforward comparisons. Due to O(Ω 2 ) overshooting, the Wilson-method's displacement response greatly exceeds all others, so it is left out of the figure views. The Newmark average acceleration method exhibits zero-order overshooting in displacement or velocity in both cases, regardless of initial conditions. The Wilsonmethod exhibits second-order overshooting in displacement in general; and first-order overshooting in velocity for problems with initial displacements. The Bazzi-method exhibits zero-order overshooting in displacement in general; and zero-order overshooting in velocity for undamped initial velocity problems, otherwise exhibiting first-order overshooting in velocity (but is still uncompetitive due to its first-order time accuracy). The GSSSS U0(0.5, 0.5, 0.5) method (i.e., TPO/G-method with = 0.5) exhibits zero-order overshooting in displacement for undamped problems but first-order overshooting in the presence of damping; and zero-order overshooting in velocity for undamped initial velocity problems with initial velocities, but otherwise first-order overshooting in velocity. The GSSSS V0*(0.5, 0.5, 0.5) method exhibits zero-order overshooting in displacement for initial displacement problems but first-order overshooting otherwise; and zero-order overshooting in velocity, regardless of physical damping in either case. This confirms the behavior stipulated by the analysis in Section 3.

Notched beam with high-frequency excitation
As a demonstration of overshooting for MDOF problems in structural dynamics, we consider the notched beam with a finite element bilinear quadrilateral mesh shown in Figure 4. The beam is made of structural steel (E = 210 GPa, = 8050 kg/m 3 , = 0.3) with a thickness of 0.1 m. The lowest frequency of the system is min = 1.2886 × 10 3 rad/s, and the highest frequency of the system is max = 2.1416 × 10 6 rad/s. The time step size is chosen as Δt = 1 10 2 min so as to resolve the lowest frequency while being over 166 times larger than the period of the highest frequency mode, thus producing asymptotic high-frequency behavior in the presence of the associated excitations.
The system response is computed by with and without physical damping. In the latter case, Rayleigh proportional damping is applied so as to produce a damping ratio not exceeding = 0.5 across the range of system frequencies as in Equation (24). The damping ratio is plotted between min and max in Figure 5.
The initial conditions for displacement and velocity are constructed as a sum of the mode shapes associated with the smallest and largest frequencies of the system so as to produce a low-frequency and high-frequency excitation simultaneously. Initial values of i (0) anḋi(0) are applied to the minimum and maximum normal modes of vibration such that the excitation of the highest mode is 5% of the excitation of lowest mode, where i = 1 corresponds to min and i = N corresponds to max . The final initial conditions are then given by Equation (25) Exact solution bounds are represented with dashed lines. The Newmark average acceleration method exhibits zero-order overshooting in general, while the other schemes exhibit at least first-order overshooting in displacement or velocity except when matched with special initial conditions. The GSSSS U0 algorithm (i.e., TPO/G-) exhibits first-order overshooting in both displacement and velocity and in general, while its amended GSSSS V0* counterpart shows advantages for initial displacement problems in terms of zero-order overshooting in both variables  Exact solution bounds are represented with dashed lines. The Newmark average acceleration method exhibits zero-order overshooting in general, while the other schemes exhibit at least first-order overshooting in displacement or velocity except when matched with special initial conditions. The GSSSS U0 algorithm (i.e., TPO/G-) exhibits first-order overshooting in both displacement and velocity and in general, while its amended GSSSS V0* counterpart shows advantages for initial displacement problems in terms of zero-order overshooting in both variables (x-direction) and transverse (y-direction) initial conditions given in Figure 6. Just as in Section 5.1, we consider both initial displacement problems, initial velocity problems, and combined problems to verify the overshooting components that arise in each case for the same selected algorithms.
The physically undamped results are shown in Figure 7, while the physically damped results are shown in Figure 8. The L 2 norm is used to depict the magnitude of the displacement and velocity response over the entire system, and the velocity results are scaled by the natural frequency to enable straightforward comparisons. Due to O(Ω 2 ) overshooting, the Wilson-method's displacement response greatly exceeds all others, so it is left out of the figure views. The same overshooting behavior described in Section 5.1 is displayed here, confirming the analysis in Section 3 applies equally to the MDOF case due to high-frequency excitations, particularly when the physical damping is of a Rayleigh type that preserves orthogonal modal decomposition. The Newmark average acceleration method (falling into the GSSSS U0V0 family) maintains zero-order overshooting in both displacement and velocity for all cases. The Bazzi-method maintains zero-order overshooting in displacement and first-order overshooting in velocity, except for physically undamped initial velocity problems in which it has zero-order overshooting in both variables (but is still uncompetitive due to its first-order time accuracy). The GSSSS U0 algorithm displays advantages only for physically undamped initial velocity problems in terms of zero-order overshooting in both displacement and velocity, while the GSSSS V0 algorithm maintains its advantage for physically damped and undamped initial displacement problems with zero-order overshooting in both displacement and velocity.

CONCLUSIONS
In this article, it is shown that the analysis of overshooting for many implicit time integration algorithms in the engineering literature is incomplete, particularly with respect to the influence of physical damping on overshooting behavior. New clarity is brought to the proper analysis technique for detecting overshoot in the presence of physical damping, which generally must include both the first and second step of the starting procedure in the common case of 3-value methods with displacement, velocity, and acceleration as dependent variables. The article examines the role of physical damping in the overshooting behavior of the Wilson-method, Bazzi-method, Hilber-Hughes-Taylormethod, Wood-Bossak-Zienkiewicz-method and three parameter optimal (a.k.a. generalized-) method, as well as the U0 and V0 families of the GSSSS framework of algorithms, proving that previously undisclosed overshooting components as a consequence of the full analysis. The GSSSS V0 family of methods is amended to the novel V0* family in order to remove the additional overshooting component, while the GSSSS U0 family (encompassing the HHT-, WBZ-, and generalized-methods) is unable to be fully remedied due to overshooting in the second step of the starting procedure, which has not been seen in any prior design and analysis. This is particularly noteworthy, as the algorithm framework is fully general, and as such there is no further improvement possible for the entire space of 3-value implicit single-step methods with second-order accuracy. Practical considerations for multi-degree of freedom structures, that is, those arising from finite element discretizations, are discussed, and the overshooting behavior of selected algorithms for SDOF and MDOF structures under various initial conditions are numerically demonstrated. Based on the results of the study, we recommend that new attention be paid to the analysis of overshooting due to physical damping and the second step of the starting procedure in the design of schemes, and those schemes with optimal overshooting behavior be chosen for problems based on their physical damping content and initial conditions.

DATA AVAILABILITY STATEMENT
Data sharing is not applicable to this article as no new data were created or analyzed in this study.