{'fraction_non_alphanumeric': 0.0420802847538066, 'fraction_numerical': 0.020526003559422582, 'mean_word_length': 4.33910472972973, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The TR-BDF2 method for second order problems in structural mechanics

The application of the TR-BDF2 method to second order problems typical of structural mechanics and seismic engineering is discussed. A reformulation of this method is presented, that only requires the solution of algebraic systems of size equal to the number of displacement degrees of freedom. A linear analysis and numerical experiments on relevant benchmarks show that the TR-BDF2 method is superior in terms of accuracy and efficiency to the classical Newmark method and to its generalizations.


Introduction
The Newmark method [21] is an implicit time discretization technique for second order ordinary differential equations that is very widely used in structural dynamics applications, when the stiffness of the system causes standard explicit methods such as the Störmer-Verlet method [11] to be inefficient. The so called Generalized α− methods [7] are extensions of the Newmark method using different time averaging parameters for different forcing terms. A full reference list and discussion of the properties of these methods are reported in [9], along with an analysis of their behaviour in the limit of arbitrarily large frequencies. Generalized α−methods (shortly denoted as G(α) methods in the following) introduce numerical dissipation in the approximation of the highest frequency modes of the solution, thus allowing to achieve unconditional stability and, in some regimes, to avoid overshoots, as discussed in [9]. However, in fully nonlinear regimes spurious oscillations may still be present. Furthermore, the stability and accuracy are dependent on the values of the numerical parameters that define the method, which have to be tuned for each specific application. For these reasons, alternatives to the G(α) methods have been sought in [22], [23] in the class of Singly Diagonally Implicit Runge Kutta (SDIRK) and Rosenbrock methods, respectively. These investigations, which date back more than 25 years by now, do not seem to have changed the attitude of practitioners, who apparently have kept on using G(α) methods for this kind of problems, apart from a limited number of exceptions, see e.g. [5], [6], [12], [13], [18], [20]. One possible reason of this preference is that, in the straightforward application of Runge-Kutta type methods to structural dynamics, the size of the system to be solved is twice that required by G(α) methods.
In this work, we extend the analysis and comparison of [22] to the TR-BDF2 method [1], a second order accurate method whose remarkable accuracy, stability and efficiency properties have been fully analyzed in [16]. In particular, this L-stable method is endowed with an embedded third order method that permits effective time step adaptation. Unconditionally monotonic extensions of the TR-BDF2 method have been recently derived in [4] and its multirate version has been shown in [3] to be quite effective in reducing the computational cost of the time discretization for multirate hyperbolic problems. Furthermore, the combination of this method with high order, adaptive discontinuous finite element discretizations for model problems represented by first order equations in time has been shown in [25] to be extremely effective in reducing the computational cost of these high order methods.
For these reasons, it is of interest to derive a form of this method that allows its application to structural mechanics problems. More specifically, we present a reformulation of the TR-BDF2 method for structural mechanics problems that only implies the solution of nonlinear systems of the same size as the number of displacement degrees of freedom. We also show how the velocity degrees of freedom do not have to be stored explicitly and can be recomputed whenever needed, thus avoiding excessive memory requirements for large scale applications. A similar reformulation was presented in [23] for a two stage Rosenbrock method. An analysis of the accuracy and dissipation properties of the TR-BDF2 method is carried out in the linear regime, extending the classical analyses to consider also the approximation of the damping terms. Both this analysis and a number of numerical experiments on relevant benchmarks show that TR-BDF2 is superior in terms of accuracy and efficiency to the classical G(α) methods.
The paper is organized as follows. In section 2 we show how the TR-BDF2 method can be applied to second order problems and reformulated so as to avoid solving systems of dimension larger than the number of displacement degrees of freedom, while at the same time not requiring extra storage of the velocity degrees of freedom. In section 3, the dissipation properties of the method are compared with those of the G(α) methods. In section 4 results of several numerical simulations are presented, which highlight the accuracy and efficiency of the method for structural mechanics applications. Some conclusions and perspectives for future work are presented in section 5.

The TR-BDFmethod for second order problems
We first consider the generic Cauchy problem y = f (y, t) on the time interval t ∈ [0, T ]. Even though all the results presented in the following also hold for variable time steps, we present for simplicity the TR-BDF2 method as employing a constant time step h = T /N and we follow in this presentation the notation and conventions in [4]. In its original formulation [1], the TR-BDF2 method is defined by the two stages: Here, u n denotes the numerical approximation of the solution at time level t n = nh, γ ∈ (0, 1) is an implicitness parameter and γ 2 = (1 − γ)/(2 − γ), γ 3 = 1/γ(2 − γ). The first stage of (1) is nothing but the application of the trapezoidal rule (or Crank-Nicolson method) over the interval [t n , t n + γh]. The outcome of this stage is then used as input for the BDF2 implicit method. The resulting combination yields a method with several interesting accuracy and stability properties. A detailed analysis of these properties is given in [16], where it is also shown that TR-BDF2 is equivalent to a three stage Diagonally Implicit Runge Kutta (DIRK) method defined by the stages Notice that this method is not Singly Diagonally Implicit Runge Kutta (SDIRK), due to the fully explicit first stage. As shown in [16], the TR-BDF2 method is second order accurate and A-stable for any value of γ. Furthermore, for γ = 2 − √ 2 it is also L-stable. Therefore, with this coefficient value it can be safely applied to problems with eigenvalues whose imaginary part is large, such as typically arise from the discretization of highly oscillatory second order problems. This is not the case for the standard trapezoidal rule (or Crank-Nicolson) implicit method, whose linear stability region is exactly bounded by the imaginary axis. Recently, the method has been applied in [25], [24] to high order accuracy spatial discretization of wave propagation phenomena written as systems of first order in time. Its monotonicity properties have also been studied in [4], where it was shown that an unconditionally monotonic extension of the method can be derived and that the method does not suffer from order reduction when applied to stiff systems. Several specific variants of the method for nonlinear problems have been proposed and analyzed in [8], while an analysis of a multirate implicit method based on TR-BDF2 has been presented in [3].
We now focus on systems of second order ordinary differential equations. In particular, as common in structural mechanics applications, we will consider systems of the form with the initial conditions y(0) = y 0 y (0) = v 0 . Here, M denotes a symmetric and positive definite mass matrix, K denotes the stiffness matrix, which is also assumed to be symmetric and positive definite, while the matrix C represents friction terms. In order to introduce the application of the TR-BDF2 method to this kind of problems, system (3) is rewritten as a first order system letting v = y and setting Γ = M −1 C, B = M −1 K, and f (y) = M −1 g(y), s(t) = M −1 z(t), so as to obtain We now apply method (1) to problem (4). Denoting by u n , w n , respectively, the numerical approximations of y, v at time level n, the trapezoidal rule stage can be written as where we have sets = s(t n+γ ) + s(t n ). The BDF2 stage yields instead Along the lines of what is done for the Newmark and G(α) methods, each of the two stages is now rewritten in terms of a single implicit step for u n+γ , u n+1 , respectively. This amounts to w n+γ = 2 u n+γ − u n − γhw n /2 /γh (9) where we have now set Fixing from now on the value γ = 2 − √ 2, so as to achieve L-stability, one obtains γ 2 = γ/2, so that A 1 = A 2 = A. Notice that this formulation has a number of advantages. First, it allows in practice to avoid doubling the degrees of freedom of the discrete problem, even though a method for first order problems is employed. For the first time step, the value of w n would be recovered from the initial datum, while formula (11) allows the reconstruction of the same term at subsequent time levels. Furthermore, both stages of the method are defined in terms of the same matrix. In order to improve the efficiency of the algebraic solvers, both equations can also be rewritten in terms of the increments δu n+γ = u n+γ − u n δu n+1 = u n+1 − u n+γ , so that the scheme can be rewritten as Finally, it is easy to notice that the mass matrix inversion is not necessary in practice. Indeed, each of the previous equations can be multiplied by M, thus yielding the two nonlinear systems where now A has been redefined as A = M + γhC/2 + γ 2 h 2 K/4 and the corresponding right hand sides as For completeness, we also present the classical Newmark method and its G(α) generalizations. Following [9], we introduce discrete approximations u n , w n , a n of the continuous displacement, velocity and acceleration values, respectively: where for a generic variable φ one has φ n+1−α = (1 − α)φ n+1 + αφ n . where γ N , β, α m , α f are specific averaging parameters. Notice that different definitions of the α averaging can also be employed and that several relationships between the values of the parameters γ N , β, α m , α f must be respected for the method to be second order consistent in the displacement and unconditionally stable. We refer again to [9] for a complete discussion of these issues and of a simple strategy to link the values of these parameters to that of the effective dissipation rate of the highest frequencies, denoted by ρ ∞ . Here, we just summarize the key results, which are reported in the following expression for the parameters β, γ N and for the parameters that define the G(α) method introduced in [7], which we consider in the following as example of the broader G(α) class.

Analysis of the TR-BDF2 method
In [22], several SDIRK methods were compared to some G(α) methods with respect to their dissipation properties. The behaviour of the TR-BDF2 method is very similar to the two stage SDIRK method analyzed in Section 2.2 of [22]. However, the analysis will be repeated here for convenience and extended to the case in which also dissipative terms are taken into account. More specifically, we will consider the methods introduced in section 2 as applied to the linear scalar problem where k, m > 0, c ≥ 0. We will only consider the case in which c 2 < 4km, for which it is well known that the analytic solution has the form where one has As a consequence, one can write for this case the exact evolution operator where we have set y n = [y(t n ), hy (t n )] T , ζ = hζ, ω = hω and Notice that the velocity variable has been rescaled (see e.g. [2]), so as to obtain an evolution operator which only includes dimensionally homogeneous quantities. Defining then x n = [u n , z n ] T , where we set z n = hw n , the discrete evolution operator corresponding to the TR-BDF2 method can be reconstructed from the specific forms of the TR and BDF2 stages, which can be written as where one defines By straightforward algebraic manipulations, one obtains which can then be compared to (27) in order to assess the properties of the TR-BDF2 method. A similar comparison can be carried out for the G(α) methods, taking into account that these methods also employ an approximation of the second derivative. Therefore, one has to consider the exact evolution operator derived from the relationship where ξ = [ω 2 + ζ 2 , ζ] T and equation (26) has been used to introduce the relationship between acceleration, velocity and displacement values. Analogously, for the generic G(α) method one has so that the comparison is between the exact operator The ratios of the spectral norms of the discrete evolution operators to that of the exact evolution operator are reported for different methods in Figure 1 as a function of (ζ, ω), on the domain 10]. Notice that these quantities are symmetric with respect to the ω = 0 axis. The corresponding relative errors in the spectral norm with respect to the exact evolution operators are reported in Figure 2. In particular, the TR-BDF2 method, the Newmark method with three different values of the damping parameters and the G(α) method introduced in [7] (denoted by CH-α) with two different values of the damping parameters are compared in this way.
It can be observed that the TR-BDF2 method introduces significantly less numerical damping in the ranges corresponding to frequencies that need to be accurately resolved, while damping is more effective on the highest frequencies. Furthermore, the TR-BDF2 method also introduces a much smaller error than its counterparts in the same frequency ranges. These findings, which could have been obtained in the undamped case with an analysis analogous to that presented in [22], are also true if damping is added, thus confirming the advantages of the TR-BDF2 method for realistic applications including possibly stiff dissipative terms.

Numerical experiments
A number of numerical experiments have been carried out in order to assess the accuracy and efficiency of the TR-BDF2 method for applications to structural mechanics and to compare it with that of the G(α) methods. We first consider one of the nonlinear numerical benchmarks discussed in [9]. Then, time discretizations of the wave equation are considered, in cases which are representative of possible applications to structural mechanics and seismic engineering. For all the problems concerning the discretization of a wave equation, we have used a finite element spatial discretization and our computations have been performed with the free software FreeFem++ [14], computing the errors in the space-time norms where t n are the time levels used by the time discretization on [0, T ].

Nonlinear system with 2 degrees of freedom
In a first numerical experiment, we consider the strongly nonlinear system with two degrees of freedom y 1 y 2 = − 10 4 y 1 (1 + 10 4 y 2 1 ) − tanh (y 2 − y 1 ) tanh (y 2 − y 1 ) with the initial conditions y(0) = [1, 1.5] T y (0) = 0. The same system was used in [9] to compare the performance of different G(α) methods. We compute a reference solution using the MATLAB ode15s solver and we compare the performance of the TR-BDF2 method, of its parent methods, i.e. the off-centered Trapezoidal Rule or θ− method (setting θ = 0.51) and the BDF2 method, and of 4 different G(α) methods, more specifically the Newmark method with ρ ∞ = 0, 0.25 and the CH-α method with the same values of the dissipation parameter.
As it can be seen in figure 3, while the performance of the G(α) methods and of the TR-BDF2 method is essentially analogous on the the fast and strongly damped variable, the latter is in general the most accurate on the slow variables.

Wave equation in 1 dimension with strongly varying coefficients
We then consider the one dimensional wave propagation problem described by These equations model the displacement u of a clamped-free elastic rod with a unit constant cross-sectional area, length L, constant density ρ and Young modulus given by E(x). The same problem has also been considered in [17], [23], [9] to study the performance of G(α) methods in a stiff case in which strong oscillations and numerical overshoots can arise. In particular, we assume L = 10.5, ρ = 0.01 and a Young modulus such that E(x) =    10 7 0 ≤ x < 0.5 10 2 0.5 ≤ x < 10 10 7 10 ≤ x ≤ 10.5, which yields a rather stiff problem. The system is discretized using P 1 finite elements on the spatial domain [0, L] and either the TR-BDF2 method, the BDF2 method or the Newmark method with β = 1/4, γ N = 1/2 for the time discretization. In all cases we have used a mesh with 21 nodes and ∆t = 0.025s. As a reference solution, we consider the one obtained by computing the exponential of the matrix of the spatial semi-discretization. In table 1 we report the absolute errors between the considered methods and the reference solution, computed at two different final times. It can be observed that the TR-BDF2 method yields consistently the smallest errors among the three second order methods under comparison. The time series of the displacements and velocities corresponding to the end of the rod (x = L) are reported instead in Figures 4, 5, while the corresponding absolute differences between the computed solutions and the reference one are displayed in Figures 6 and 7. Again, the smaller error produced by the TR-BDF2 method is apparent.

Two-dimensional test with analytic solution
In order to verify the implementation used in the more advanced tests presented in Section 4.4, we have considered the the problem u(x, y, 0) = 0 in Ω, ∂u ∂t (x, y, 0) = 2πsin(πx)sin(πy) in Ω, where the domain Ω = [0, 1] × [0, 1] and T = 1. The solution for this problem is given by u(x, y, t) = sin(2πt)sin(πx)sin(πy). This solution has been used to perform a convergence study of the TR-BDF2 time discretization, coupled again to a P 1 finite element discretization in space. To estimate empirically the convergence order, we choose the maximum diameter h of a mesh element equal in size to the time step ∆t and we compute the following expression for two different values of ∆t: r emp = log(err(∆t 2 )/err(∆t 1 )) log(∆t 2 /∆t 1 ) .
The results of the convergence test are displayed in table 2.
We have also repeated the computation using the BDF2 method and the Newmark method with β = 1/4, γ N = 1/2. The results are h = ∆t   Table 3: Convergence behaviour of the BDF2, Newmark and TR-BDF2 method in test (36), absolute errors in L ∞ (L 2 ) and L 2 (H 1 ) norms.
Step Time

Two dimensional elasticity
We consider the plane-strain equations of two dimensional elasticity, presented in [19] as the first order system: ∂σ x,x ∂t = (λ + 2µ) ∂u ∂x + λ ∂v ∂y ∂σ x,y ∂t = µ ∂v ∂x + µ ∂u ∂y ∂σ y,y ∂t = λ ∂u ∂x + (λ + 2µ) ∂v ∂y Here, σ i,j denote the components of the (symmetric) stress tensor Σ, u = (u, v) is the velocity vector with components u, v in the x, y directions, respectively, ρ is the medium density and λ, µ are the medium Lamé coefficients. For the purpose of our work, we reformulate these equations as second order equations in terms of the displacements δ = (δ x , δ y ). As a result, equations (38) are equivalent to which, after substitution of the expression for the stresses into the momentum equation, can be rewritten in vector form as As discussed in [19], this system has P and S-wave solutions with propagation speeds We consider the above equations on the domain Ω = [0, 3] × [0, 3] and on the time interval [0, T ], with T = 10 −2 s. We assume that the density is constant and equal one in all the domain, while the elastic constants are such that a more rigid inclusion is embedded as shown in Figure 9. More specifically, in the darker region, denoted as Zone 2 in Figure 9b), whose boundary is described by the equation is (x − x 0 ) 2 + (y − y 0 ) 2 = 1/100, with (x 0 , y 0 ) = (1.65, 1.65), we assume that λ = 200, µ = 100, while in the region represented in red one has λ = 2, µ = 1. We then proceed to study the wave propagation resulting from imposing at the initial time a displacement in the region denoted as Zone 1 in Figure 9b), whose boundary is described by the equation (x − x 0 ) 2 + (y − y 0 ) 2 = 9/10 4 . More specifically, we consider the function defined as d(x, y) = exp(−(x − x 0 ) 2 /2 − (y − y 0 ) 2 /2) in Zone 1 and d(x, y) = 0 in the rest of the domain and we take as initial conditions interval of size ∆t starting from zero. For the numerical discretization we consider P 1 continuous finite elements on the locally refined unstructured mesh depicted in Figure 9 a), which is composed of 15036 vertices and 7559 elements with an average diameter equal to h = 2.39 × 10 −2 in Zone 2 and equal to h = 0.22 outside. For the time discretization we use a time step of ∆t = 1.25 × 10 −4 , which entails a Courant number based on the P-wave speed of c P ∆t/h ≈ 1.04 in Zone 1 and c P ∆t/h ≈ 0.01 in the rest of the domain, thus mimicking a situation in which an unconditionally stable method is typically applied.
Since an analytic solution is not available, we consider as a reference solution that computed by the same space discretization coupled to a fully implicit Gauss method of order four, see e.g. [11]. Therefore, the errors with respect to the reference solution are an estimate of the time discretization error only. We compare the implicit Euler method, the Crank-Nicolson method, the Newmark method with β = 1/3 and γ N = 1/2 and the TR-BDF2 method. We include also the implicit Euler method as a representative of a robust, unconditionally stable, albeit first order method. A pictorial view of the resulting wave propagation on the portion of the domain around Zone 2 at different instants as computed by the TR-BDF2 method is reported in Figures  10,11.
The global relative errors with respect to the reference solution are then reported in Tables 4 and 5 Figure 12 at the three control points A,B and C shown in Figure 9 b). Point A point is located at the center of Zone 1, where the initial displacement occurs, point B is in Zone 2 and C is the outer region but close to Zone 2. It can be observed that, apart from the solutions computed by the implicit Euler method, the others are all in good agreement. The

Method
Error L ∞ (L 2 ) Error L 2 (H 1 ) Euler 7.07e-2 1.71e-2 Newmark 8.35e-3 5.29e-3 Crank-Nicolson 5.89e-3 3.68e-4 TR-BDF2 2.93e-3 1.82e-4 Table 5: Relative errors on the displacement in the y−direction corresponding absolute errors with respect to the reference solution are shown in Figure 13, while the error norms at the three control points are reported in Tables 6, 7 and 8, respectively. These results highlight the fact that the TR-BDF2 method yields more accurate solutions both in the higher and lower Courant number regions, with errors that are consistently at least 50% smaller than those of the other methods. This is especially important in view to the combination of these methods with higher order discretizations in space.  Method L 2 error on δ x L ∞ error on δ x L 2 error on δ y L ∞ error on δ y Implicit Euler 1.74e-2 0.39 6.31e-2 1.41 Newmark 4.32e-3 0.137 9.39e-3 0.28 Crank-Nicolson 2.99e-3 9.16e-2 6.50e-3 0.19 TR-BDF2 1.46e-3 4.52e-2 3.21e-3 9.43e-2

Conclusions and future work
We have proposed a reformulation of the TR-BDF2 method that allows to apply it without overheads in terms of computational cost or memory to large scale problems in structural mechanics. Our work extends the analysis and the comparison of [22] and a similar reformulation proposed for two-stage Rosenbrock methods in [23]. In particular, we have presented a reformulation of this method that only implies the solution of nonlinear systems of the same size as the number of degrees of freedom necessary to describe the displacement variables. Furthermore, the velocity degrees of freedom do not have to be stored explicitly and can be recomputed whenever needed, thus avoiding the shortcomings of naive implementations of solvers for first order ODE systems. An analysis of the dissipative behaviour of the method, which was carried out considering also the damping term, shows that TR-BDF2 is superior in terms of accuracy and efficiency to the classical G(α) methods. This finding is confirmed by a number of numerical experiments on significant benchmarks of increasing complexity, in which the TR-BDF2 method consistently yields errors least 50% smaller than those of the other methods, both in high and low Courant number regions. This is especially important in view to the combination of robust, unconditionally stable time discretization methods with higher order discretizations in space.
The future developments of this work include the development of a multirate version of the proposed method, based on the self-adjusting multirate technique described in [3], for application to structural mechanics problems with multiple time scales [10] and to monolithic treatment of fluid structure interaction problems [15], as well as the coupling of the proposed time discretization to high order, adaptive discontinuous finite element spatial discretizations, so as to achieve the performance improvements demonstrated in [25] also in the wave propagation problems typical of structural and seismic engineering.