BDF-α: A Multistep Method with Numerical Damping Control

When solving numerically the stiff second order ODE system obtained after semidiscretizing the wave-type partial differential equation (PDE) with the finite element method (FEM), and similarly to the HHTα method, which allows the numerical damping of the undesirable high frequency modes associated to FEM semidiscretization, we have constructed a modification of the 2-order BDF method (the BDF2 method), which we have called BDF-α. This new method is second-order accurate and with a smaller local truncation error than the BDF2, it is unconditionally stable for some values of α and it permits a parametric control of numerical dissipation.


Introduction
The second order Ordinary Differential Equation (ODE) system obtained after semidiscretizing the wave-type partial differential equation (PDE) with the finite element method (FEM) shows strong numerical stiffness. Its resolution requires the use of numerical methods with good stability properties and controlled numerical dissipation in the high-frequency range.
Some of the methods developed in the linear range are shown in [1]. Time-stepping algorithms such as the Collocation method [2], the Wilson method [3], the HHT-α method [4], the Houbolt method [5], or more recent works as the generalized-alfa method [6] are some of the methods commonly used.
The development of similar methods for nonlinear problems is more recent and it has been originated by the presence of numerical instabilities when solving nonlinear problems by methods which are unconditionally stable in the linear range. Generally, this instablities are due to the increase of energy of the discrete system and this is the reason why the formulation of energy-momentum conserving schemes, [7] and [8], and Energy Dissipative, Momentum Conserving schemes, EDMC [9,10], have been created.
In this paper, we work in the linear range and we formulate a new method called BDF-α, which is second-order accurate, unconditionally stable for some values of the parameter α and permits a parametric control of numerical dissipation.

BDF methods
We will consider the next initial value problem (IVP): y ′ (t) = f (t, y(t)) , y (t 0 ) = y 0 (1) where T = [t 0 , t n ] is a finite interval and y:[t 0 , t n ] → R m and f:[t 0 , t n ]xR m → R m are continuous functions. The Backward Differentiation Formulae (BDF) were proposed by Gear [11]. They are linear multistep methods useful to solve ODEs of order 1, such as (1). Since they were introduced, the Backward Differentiation Formulae have been widely used due to their good stability properties for solving stiff problems. The BDF of order k can be expressed as follows in terms of backward differences: where: ∇y n+k = y n+k − y n+k−1 and ∇ j y n+k = ∇ ( ∇ j−1 y n+k ) . Backward differences verify the next expression: Developing the backward differences of expression (2) we get this equivalent expression for the BDFs: k ∑ j=0 α j y n+j = hf n+k (4) From (4), it is easily deduced that BDFs are linear multistep methods which respond to the general form: Where the constants α j are chosen in order to verify the order condition, that says that the multistep method will be of order p if the next condition is satisfied [12,13]: The BDFs are unconditionally stable for orders 1 and 2. The stability properties of the BDFs and their error constants can be calculated following the formulae given in [11,12,14]. The A(α)-stability of the BDFs [15,16] and their error constants can be seen in Table 1. If we are interested in studying the amplification factor of the BDFs, the method has to be applied to the test equation y ′ = λy and we have to rewrite the result using the form: where X n+k = (y n+1 , y n+2 , ..., y n+k ) t , X n+k−1 = (y n , y n+1 , ..., y n+k−1 ) t and A is a square matrix of dimension kxk. Eigenvalues of matrix A are calculated and the largest one in module is the one we are interested in, which is the spectral radius: The spectral radius is closely connected to the stability of the method and ρ(A) ≤ 1 is required to prevent amplification of A n as n becomes large.
When working with second order ODEs, we have to consider their test equation d ′′ + ω 2 d = 0, which represents an undamped vibrating physical system with natural frequency f = ω/(2π). Transforming this equation in its equivalent system of two first order ODEs we get: System (9) can be uncoupled in two independent equations of the form: and the application of the general theory, (7) and (8), to the second order problem d ′′ + ω 2 d = 0 is equivalent to apply this theory to the first order test equation y ′ = λy, but now taking into account that λ takes the specific value λ = iω.
In these kind of problems, numerical damping is related to the spectral radius (8). For example, it is easy to see that the trapezoidal method is a linear multistep method with frequency-independent spectral radius ρ = 1. It means that this method does not show numerical damping ( Figure 1). However, BDFs have spectral radii that decrease as frequency increases ( Figure 3); that is to say, BDFs show algorithmic damping ( Figure 1).

Example for calculating the amplification matrix of BDFs
As an example, we have applied this procedure to the BDF4, which expression is given by: We will apply the method given by (11) to the test equation y ′ = λy, reaching the following expression: where:ĥ = λh. Expression (12) can be written in matrix form as: where: And matrix A = A −1 1 A 2 is the amplification matrix of the method (11).
The stability region S of a numerical method is the set of valuesĥ ∈ C for which the spectral radius is less than or equal to one: S = {ĥ ∈ C : ρ(A(ĥ)) ≤ 1 } . Applying this condition to the amplification matrix A of the BDF methods, their stability regions are obtained as it is shown in Figure 2. In the particular case of second order ODEs, if we want to study how the values of the spectral radius depend on frequency, ρ(A(ĥ)) has to be evaluated for different values ofĥ = iωh taking into account that for this specific case λ = iω. Following this condition, the dependence of the spectral radius on ωh/(2π) = h/T can be obtained for BDF methods and others as it is shown in Figure 3. It can be seen that the spectral radii of the BDFs tend to zero when h/T → +∞. It can also be observed that in the case of BDFs of orders 3, 4 and 5, there is a zone in which the spectral radius becomes greater than one. This is because these methods are not stable in a certain zone near to and in the imaginary axis.
When solving a second order ODE, measures of accuracy such as the algorithmic damping ratioξ and the relative period error (T − T )/T are relevant, see Figure 1. The first oneξ, says how much decays the amplitude and the second measures how much periods are elongated or shortened. It is difficult to obtain analytical expressions forξ and (T − T )/T but they can be calculated by computer evaluation.
When solving second order ODEs that come from the FEM semidiscretization of the wave-type PDE, the higher modes of these equations are artifacts of the discretization process and they are not representative of the behaviour of the governing PDE. So it is convenient to introduce algorithmic damping to remove the participation of highfrecuency modal components. Even though the minimum value that BDFs show for ρ ∞ = lim h/T →∞ ρ(A) is the most effective in removing the participation of high-frecuency modal components, it is more convenient to have an algorithm which permits a parametric control of numerical dissipation, allowing the participation of the medium frequencies. Figure 3 shows how the parametric control α = −0.3 of the HHT-α method decreases the medium range frequency numerical damping with respect to the stable BDFs, while maintaining a low spectral radius in the high range, which is enough to dissipate quickly this frequency components.
We are interested in creating a method with the same properties as the HHT-α but taking as basis the BDFs. Next we will review direct methods for second order ODEs and in particular the Newmark parametric family in which HHT-α is based.

Linear multistep methods for second order ODEs
Consider the second order initial value problem: where T = [t 0 , t n ] is a finite interval and y:[t 0 , t n ] → R m and f:[t 0 , t n ]xR m → R m are continuous functions.
Writing equation (15) as a first order system and applying a lineark-step method (5): An elimination of (16) results [14]: where k = 2k. Expression (17) corresponds to a linear multistep method for second order differential equations, which is similar to first order equations but with different coefficients and appearing the factor h 2 instead of the factor h of the first order ODEs.
HHT-α method can be reduced to (17) form, and general properties of second order ODEs multistep methods can be used to study its properties. HHT-α method is an extension of the Newmark family which in turn can be introduced based on the trapezoidal method.
In computational mechanics, the second order linear ODE takes the form: which corresponds to a mass-spring-damper system, where m is the mass, k the spring constant, c the damping coefficient, v = d ′ the velocity and a = d ′′ the acceleration. Converting the second order equation (18) into a first order system: where a verifies (18), and applying the trapezoidal method to it, we get: It can be observed that in relations (20), the displacement and the velocity are updated using the average velocity and acceleration of the instants t n and t n+1 , respectively. By substituting the second formula of (20) in the first one, the exact formula of the uniformly accelerated motion is obtained: If instead of the average acceleration, two weighted means are considered, taking in one of them the parameter 2β and in the other one the parameter γ, we obtain the biparametric family of the Newmark method [17]: where β and γ are free parameters which govern the accuracy, stability and numerical dissipation of Newmark's algorithm.
With the purpose of studying the stability and the numerical damping of the method, we apply the method to the test equation d ′′ + ω 2 d = 0, where w = √ k/m. In this case, (22) can be succinctly written in this recursive way: where: , and: The system (23) can be reduced to a difference equation in the displacementes, which takes the form of a linear multistep method for second order differential equations (17): where the coefficients α j , β j are given by: The order of the method can be calculated by applying the order conditions for linear multistep methods of the type (24) [18]. The method results second-order accurate when γ = 1/2, [1,17]. The method is unstable when γ < 1 2 and it is unconditionally stable when 1 2 ≤ γ ≤ 2β. High frequency dissipation is achieved when: In the second-order accurate Newmark method (γ = 1/2), β ≥ 1/4 retains unconditional stability. If in addition, high frequency dissipation is requiered, β = 1/4 has to be verified. In this case, Newmark's method becomes the trapezoidal method, and although it verifies the high frequency dissipation condition (26), high modes are not damped as ρ ∞ = 1, see Figure 4. So, second-order accurate condition does not allow numerical dissipation.
The HHT-α method is a modification made to the Newmark method, with the aim of obtaining numerical dissipation in the high frequencies while retaining the order and stability conditions. This method was proposed by Hilber-Hughes-Taylor [4] and it is also known as HHT or the α-method. The method consists of maintaining the two first expressions of the Newmark method (22) and varying the third expression. That is to say, modifying the expression of the time-discrete equation of motion with a new parameter α as follows: where: In the case that α = 0, the HHT-α method is reduced to Newmark's method. When applied to d ′′ +ω 2 d = 0, the method takes the recursive form (23), where A = A 1 −1 ·A 2 is the amplification matrix and: Similarly to Newmark, HHT-α method can also be reduced to a three-step linear multistep method for second order differential equations (17): where the coefficients α j , β j are given by: The method is second-order accurate when γ = 1−2α the HHT-α method is a two-order method, unconditionally stable and high frecuency dissipation can be obtained.
In Figure 4 we can see the spectral radii of Newmark and HHT-α methods together with the spectral radii of other typical methods used in computational mechanics, such us, Houbolt's method, Collocation method (among them the Wilson method which is a concrete case of the Collocation method) and Park's method [1,2,5]. In this figure we can see that the spectral radii of Houbolt's and Park's method tend to zero as ∆t/T → ∞ which is typical in backward-differences schemes. It can also be seen, how the A-stable HHT-α method when α ∈ , dampens the high modes more strongly as the value of α decreases, but in a much smoother way than the BDFs for the medium range of frequencies (see Figure 3).

BDF-α
We have seen in Figure 3 that the spectral radii of the BDF2 tends to zero when h/T → +∞. This method has strong dissipation of medium and high frequencies and our aim is to construct a parametrized method wich allows numerical damping control, in the same way that the HHT-α does. The BDF2 is a two-order and A-stable method given by the next expression: We will make some considerations about the modification that have arisen: 1. It has to be based on the BDF2 and it must maintain the order and the stability characteristics of the BDF2; that is to say, the new method needs to be of order 2 and A-stable (unconditionally stable) if it is possible. In addition to this, the new method should lend a wider range of the spectral radius than the ones obtained with the BDF2. 2. The first modification which we have been considering consisted of adjusting each addend of BDF2. In this way a general expression in which y, y ′ evaluated in t n , t n+1 , t n+2 is reached. That is to say, we have been trying with an expression like: We have dismissed this possibility (32), mainly for two reasons. The first one is that this expression's origin is not equivalent to the couple Newmark&HHT-α. And the second one is that even though it is possible to achieve order 2 using (32), expression (32) has all the addends that appear in methods such as Adams Moulton and BDF2, so we would be creating a linear combination that we are not interested in. { Adams Moulton (k = 2): y n+2 = y n+1 + h 3. The HHT-α method, which was constructed taking the Newmark method as its basis, verifies all the characteristics that we want for our method, and the form in which it was constructed may help us in our way of thinking about the most convenient form that the new method should have. Newmark's method is given by the expressions (22). In the method HHT-α, the first two expressions of (22) are maintained and the third one is adjusted: Newmark's form: HHT-α's form: 4. A possible formula which we have tried consisted of adjusting uniquely the right-hand side of expression (31): The problem that we have found is that the numerical method given by the expression (36) has only order 1 when α ̸ = 0 and when α = 0, expression (36) is reduced to the BDF2.
5. Finally, we have considered a modification with three free parameters α, β and γ which can be expressed as follows, in order to clearly see the calculations required: Reordering terms in (37) we obtain the next expression for the modification of the BDF2 that we require: This expression can be written in the following abreviate way, which corresponds with the standard form of a linear multistep method or linear k-step method (5):

Order and error constant
Theorem. The method given by the expression (38) will be of order 2 if the constants α, β, γ verify the following conditions: α = 3 2 β = 2γ.

Proof.
We have already said that BDF-α's expression (38) can be written as a multistep method (5) where the values of the constants are given by (39). So the BDF-α method is a linear multistep method, specifically a 2-step method. This method will be of order 2 if C 0 = C 1 = C 2 = 0 and C 3 ̸ = 0 [12,13], where the C i constants are given by (6). We substitute the values of (39) in the expressions of (6) and we get: We already have C 0 = 0. The method will be of order 2 if C 1 = C 2 = 0 and this happens when α, β, γ are chosen in this way: α = 3 2 β = 2γ.
We can conclude by saying that the expression of the 2-order BDF-α is written as follows: Observe that when α = 0 we have the BDF2 method and for α = −0.5 the trapezoidal method is obtained. When α = a 1−a where a ∈ R the method (41) is the same as the A-BDF method [19] for k = 2. The local truncation error of the method BDF-α is given by the next expression: And the error constant C of the method is given by: where, after substituting the second order condition α = 3 2 β = 2γ in C 3 , we get: and σ(r) is one of the two characteristic polynomials of the method [11,12,14] given by: By substituting (44) and (45) in (43), the error constant of the BDF-α method is obtained: When α = 0 the value of the error constant is C = − 1 3 which, as expected, is the error constant of the BDF2 method. In the case α = −0.5, we have the trapezoidal method and we get the well-known error constant: C = − 1 12 . The value of the error constant in relation to the value of α is shown in Figure 5. In Section 4.2 we will see that the BDF-α method is A-stable for α ≥ −0.5, so, the BDF-α that corresponds to α = −0.5 is the one that has the smallest error constant, which corresponds to the trapezoidal rule. So, the theorem of Dhalquist (1963) is verified; which says that a multistep A-stable method is of order p ≤ 2, being the trapezoidal rule the second order A-stable method with smallest error constant [15,20].
So, we have that the BDF-α is second-order accurate and it has a smaller error constant than the BDF2 for α ∈ [−0.5, 0).

Stability regions
The region of absolute stability of the overall method BDF-α is found using Schur's theorem [18]. To do this, we will apply the method BDF-α to the test equation y ′ = λy. That is to say, hf j = hλy j is introduced in expression (41). ( We will set y p = r p to obtain the characteristic polynomial and after taking outĥ, we get: By substituting r = e iθ in (48) we get the expression of the boundary of the stability regions: For α ≥ −0.5 the denominator ofĥ(θ) is bounded from below by a positive constant. So, having fixed α ≥ −0.5, for a sufficiently large real number which depends on α and independent of θ, R (α) ∈ R, next expression is verified: The boundary of the stability regionĥ(θ) lies in a compact subset of the right half-plane C + . For anyĥ ∈ C − , absolute stability is achieved and for continuity reasons this implies that C − belongs to the stability region, which means that the method is A-stable when α ∈ [−0.5, +∞). Figure 6 shows some stability regions of the BDF-α for some values of α.

Amplification matrix
After applying the BDF-α method to the test equation, the following expression is obtained: where: Matrix A is the amplification matrix whose eigenvalues are given by: The characterization of the numerical dissipation of the method in function of the parameter α requires the study of the limit of the spectral radius whenĥ → ∞: When α ∈ (−∞, −1), ρ ∞ is greater than 1; we have already seen that in this case the method is not A-stable. In the interval α ∈ (−1, 0] expression (54) takes these values: Again, we see that ρ ∞ > 1 when α ∈ (−1, −0.5).
For the A-stable BDF-α method, and this happens when α ≥ −0.5 as has been proved in Section 4.2, the value of ρ ∞ is given by: and from this last expression it is easy to conclude that in the A-stable BDF-α method, ρ ∞ takes all the values of the interval (0, 1]. This means that having fixed a value for ρ ∞ ∈ (0, 1], it is possible to obtain a secondorder accurate and A-stable BDF-α method with that spectral radius inĥ → ∞. So, our new method permits a parametric control of numerical dissipation. In Figure 7 some spectral radii of BDF-α are shown. We have also calculated the value of the parameter α of the method BDF-α, which has the same spectral radius as the method HHT-α when h/T → +∞ for HHT-α = −0.05 and HHT-α = −0.3. For HHT-α = −0.05 the spectral radius when h/T → +∞ is ρ = 0.90476. The same value is obtained in h/T → +∞ for α 1 = −0.475065 and α 2 = 9.50 of the method BDF-α. The Figures of the spectral radii are similar in these cases although it can be seen in Figure 8 that they are not exactly the same. In the same way, the value of the spectral radius in HHT-α = −0.3 is ρ = 0.53846. The method BDF-α for α 1 = −0.35 and α 2 = 1.17 has the same spectral radius in h/T → +∞ as the HHT-α = −0.3. In both cases, HHT-α = −0.05 and HHT-α = −0.3, the BDF-α couples obtained are second-order accurate and unconditionally stable, but the methods BDF-α 2 = 9.50 and BDFα 2 = 1.17 have a lower spectral radius in the medium range frequency with respect to the HHT-α. So, the methods BDF-α 1 = −0.475065 and BDF-α 1 = −0.35 are more similar to their respective HHT-α methods.

Relative period error and algorithmic damping
In Figures 9, 10 and 11 the algorithmic damping and the relative error in period of the method BDF-α are shown, which have been calculated using [1]. In Figure 9 we can see that the greater α, the greater the algorithmic   damping. An algorithmic damping equal to zero is obtained for BDF-α (α = −0.5), which is the trapezoidal method. The trapezoidal rule is again, the one with the smallest relative period error. It can also be seen that figures of the algorithmic damping of the method HHT for α = −0.05 and α = −0.3 are smaller than the ones corresponding to BDF-α with α = −0.475065 and α = −0.35, which were the couple of methods that in h/T → +∞ have the same spectral radius. The same occurs with the relative errors in period as can be seen in Figures 10 and 11: the ones corresponding to the HHT-α methods are smaller than the ones of the BDF-α in the same couples.

Concluding remarks
The BDF-α method is a parametrized second-order accurate multistep method, A-stable when α ∈ [−0.5, +∞) and which allows controlled numerical dissipation in the medium and high-frequency range when applied to second order ODEs modelling vibratory systems. This is a desirable property when dealing with second order ODE systems associated to the FEM semidiscretization of the wave-type PDEs.
In addition, the BDF-α method improves the constant error of the BDF2 and its spectral radii ρ ∞ sweeps the whole interval [0, 1] improving the numerical damping control offered by the HHT-α method.