About the Use of the HdHr Algorithm Group in Integrating the Movement Equation with Nonlinear Terms

Abstract

This work summarizes the HdHr group of Hermitian integration algorithms for dynamic structural analysis applications. It proposes a procedure for their use when nonlinear terms are present in the equilibrium equation. The simple pendulum problem is solved as a first example and the numerical results are discussed. Directions to be pursued in future research are also mentioned.

1. Introduction

The present work resumes other authors’ studies dedicated to the numerical integration of the movement equation with the use of Hermitian operators defined in [1], aiming application on dynamic structural analysis.

The methods developed as well as the properties used for the solution of free vibrations without damping were described in [2]. Their implementation in a free software for structural analysis with discrete models through the Finite Element Method has been presented in [3]; the damping term inclusion and its effects on properties have been studied and presented in [4, 5].

So far, results confirmed the presence of unconditional stability properties, local error with pre-established order, and asymptotic annihilation; this last property is, by definition, a numerical damping that acts in large steps and eliminates the effects of higher modes of vibration, in multidegree of freedom cases. These modes are artificial ones, introduced by the semidiscrete model employed, so that their contribution is frequently undesirable in the analysis performed. In direct time integration, once a step is chosen for the movement history representation, the higher the mode is the larger is the step related to its period, and more intense on it are the asymptotic annihilation effects.

The fact that these are single-step methods indicates that they are possibly suitable for use in non-linear problems and this is the main motivation for the present text. The study of such theme has a considerable history, but it is still a current matter, once it allows different material behaviors and the influence of displacements in the account of internal forces, as it is shown in [6]. However, what happens is that when nonlinearity is considered, a broad field is opened. Therefore, it is important to limit the problem to guarantee good performance to the model analysis, sufficient quality of response, and usage viability of the method to be employed.

In this context, this is a contribution that begins describing the equations of the non-linear problem, followed by a bibliographical review that allows to formulate a numerical solution, ending with an example of the application method specially developed for single degree of freedom cases, which is applied to the classic pendulum problem.

The conclusion reached is that the method provides indications of its viability and demonstrates the possibility of usage in other cases, which we intend to research in future studies.

2. Objective

The movement equation, employed in dynamic structure analysis, has the following formula:

𝑚𝑥^••+𝑐𝑥^•+𝑘𝑥=𝑓(𝑡), (2.1) where m is mass, c is damping, and k is stiffness, while 𝑓(𝑡) is the external force applied as a function of time, which is represented by 𝑡. Each dot above the variable represents, as usual, a derivation in 𝑡. A comparison with the definition presented in [7] for the linear second-order differential equation demands that coefficients 𝑚,𝑐, and 𝑘 are independent from 𝑥 or 𝑥^• to attend to this classification. For semidiscrete structures, with multiple degrees of freedom, the equation system represented in matrix terms is analogous to (2.1). For a representation of a non-linear equation, the following expression is used:

𝑚𝑥^••+𝑟(𝑥,𝑥^•)=𝑓(𝑡). (2.2) in which the term 𝑟(𝑥,𝑥^•) represents the forces that restore movement. The objective of the present study is to propose a procedure for using the algorithms of the HdHr group in integrating (2.2), which considers, among the phenomena mentioned above, those that allow it to be represented in the following formula:

𝑚𝑥^••+𝑘(𝑥)𝑥=𝑓(𝑡), (2.3) that is, in which damping is not present and stiffness depends on the displacement.

3. Literature Review

In a work from 1977, Zienkiewicz [8] shows that the condition for equilibrium could be understood in terms of weighted residuals; since then, Galerkin’s approach, already traditional in finite elements, has been used to develop methods for structural dynamics.

The properties for ideal algorithms listed in [9] and in other studies are unconditional stability, second-order precision, small errors generation in frequency and damping, to have elevated, or controllable numerical damping for higher frequency response, to be computationally efficient and self-starting, that is, single-step.

In 1994, Piché [10] presents a stable method with second-order precision, using single-step, in which each step is necessary to formulate a matrix of tangential rigidity, to perform the LU decomposition, and to solve two linear systems, but without the need of interactions. Also in 1994, a procedure is proposed by Tarnow and Simo [11] to transform algorithms originally of the second order into fourth order, maintaining their stability and energy conservation properties.

Nawrotzki and Eller [12] developed a new concept of unified stability, based on Lyapunov’s exponents, and an interesting point is to verify how HdHr methods face it, that is, an interesting matter for following works.

A single-step algorithms class was proposed in 2001 by Armero and Romero in [13], for non-linear elasto-dynamics, which contains numerical dissipation controllable at higher frequencies similarly to the asymptotic annihilation of the HdHr group (object of the present study). This property is required for the solution of systems involving conditional stability. These authors also propose [14] another class of second-order algorithms for application in this class of problems, also with numerical dissipation controllable at higher frequencies, with unconditional energy dissipation and momentum conservation, with dissipation properties controlled by some parameters introduced.

Galerkin’s explicit predictive-multicorrective methods, developed for linear analysis, are reviewed in [15] for arbitrary non-linear analysis. The formulation inherits precision properties from the implicit methods from which they are derived and have third-order precision, large limits of stability, and numerical dissipation controlled by a parameter.

Also in 2002, Modak and Sotelino [16] present what they call the generalized method for use in structural dynamics, a single-step one and whose principle is to adopt approximations in Taylor series, truncated for displacement, velocity, and acceleration along each step. In problems involving elasto-plasticity, it may be necessary to adapt the step. There are no limitations in applying the process to geometric nonlinearity.

Recently, the concepts discussed in [15] were reviewed by Mancuso and Ubertini in [17], to present an iterative algorithm of low computational cost, originated from another one for linear problems and inheriting stability and dissipation properties from it; therefore, the iterations are used only to gain precision in the result. The method tends to filter the higher-frequency contributions (in the same way that the HdHr group does).

A conservation analysis of the integration scheme presented in Hilber and Hughes [18] has been used as a starting point to Hauret and Le Tallec [19] present a way of introducing controllable dissipation of energy maintaining momentum conservation in existing methods.

4. Methodology

The following approach is the same employed in several related works, notably Modak and Sotelino [16], which was also based on the use of Taylor series centered on the t_i instant, truncated for the representation of displacement and derived in the t_i+1 instant and in the usage, with adaptations, of the algorithm in linear form. A previous estimative of the displacement is used in the final step only for calculating internal and external forces. The integration method is applied to obtain an approximation for the displacement in this instant and the remainder is calculated in the movement equation. If this is unacceptable, the estimative is refined and the procedure is repeated until the remainder is below a limit. The estimated values of displacement, velocity, and acceleration are then accepted as representing the real values, and the procedure is applied to the next step.

4.1. HdHr Algorithms in Linear Form

This group of aforementioned algorithms uses a pair of Hermitian expressions:

𝐻_𝑚,𝑛=^𝑛_𝑗=0(Δ𝑡)^𝑗𝑎_𝑗𝑥_𝑖^(𝑗)+^𝑚_𝑘=0(Δ𝑡)^𝑘𝑏_𝑘𝑥_𝑖+1^(𝑘)+𝑅(Δ𝑡^𝑟+2)=0,𝑑𝐻_𝑚,𝑛𝑑𝑡=𝐻_{𝑚+1,𝑛+1}=^𝑛+1_𝑗=1(Δ𝑡)^𝑗𝑐_𝑗𝑥_𝑖^(𝑗)+^𝑚+1_𝑘=1(Δ𝑡)^𝑘𝑑_𝑘𝑥^(𝑘)_𝑖+1+𝑅(Δ𝑡^𝑟+3), (4.1)(4.2) where Δ𝑡=𝑡_𝑖+1−𝑡_𝑖 is the step, 𝑡_𝑖 is the last instant of known movement, 𝑡_𝑖+1 is the next unknown instant, 𝑅 is the remainder resulting from the truncation of the Hermitian expression, and 𝑟 is the order of the local error resulting from efforts (and accelerations). One remarks also that 𝑐_𝑗=𝑎_𝑗−1 and 𝑑_𝑘=𝑏_𝑘−1. The values for 𝑚 and 𝑛 may lead to guarantee the desired properties, and for the expressions chosen in setting up the HdHr group, resulted for the local error order: 𝑟=𝑚+𝑛−1. The coefficients 𝑎_𝑗 and 𝑏_𝑘 are obtained by applying Taylor series centered on the 𝑖 instant to represent the displacement and its derivatives in the 𝑖+1 instant (4.2), as shown in [2]. Specifically, for 𝑚=𝑛=1, (4.1) writes

𝐻_1,1=𝑎₀𝑥_𝑖+𝑏₀𝑥_𝑖+1+Δ𝑡𝑎𝑥₁^•_𝑖+Δ𝑡𝑏𝑥₁^•_𝑖+1+𝑅(Δ𝑡³)=0 (4.3) or, performing the mentioned expansion in Taylor series,

𝑥_𝑖+1=𝑥_𝑖+Δ𝑡𝑥^•_𝑖+Δ𝑡²2^••𝑥_𝑖+Δ𝑡³6⁽³⁾𝑥_𝑖+⋯,_•𝑥_𝑖+1=𝑥^•_𝑖+Δ𝑡𝑥^••_𝑖+Δ𝑡²2𝑥_𝑖⁽³⁾+Δ𝑡³6𝑥_𝑖⁽⁴⁾+⋯, (4.4) leading to 𝐻_1,1=𝑎₀𝑥_𝑖+𝑏₀(𝑥_𝑖+Δ𝑡𝑥^•_𝑖+Δ𝑡²2^••𝑥_𝑖+Δ𝑡³6𝑥_𝑖⁽³⁾+⋯)+Δ𝑡𝑎𝑥₁^•_𝑖+Δ𝑡𝑏₁(𝑥^•_𝑖+Δ𝑡𝑥^••_𝑖+Δ𝑡²2𝑥_𝑖⁽³⁾+Δ𝑡³6𝑥_𝑖⁽⁴⁾+⋯)+𝑅(Δ𝑡³)=0, (4.5) and placing 𝑎₀=1, canceling terms or orders 1 and 2 in Δ𝑡 requires

𝑎₀+𝑏₀𝑏=0,₀+𝑎₁+𝑏₁𝑏=0,₀2+𝑏₁=0, (4.6) and therefore this Hermitian operator in the form of (4.1) has the coefficients 𝑎₀=1,𝑏₀=−1,𝑎₁=𝑏₁=1/2.

Equilibrium equation (2.1), as well as its derivatives in time, provides additional relationship between displacement, velocity, and acceleration, that is; one may also write

^•𝑥^•_𝑖=−𝜔²𝑥_𝑖,𝑥_𝑖⁽³⁾=−𝜔𝑥^2•_𝑖,…,_••𝑥_𝑖+1=−𝜔²𝑥_𝑖+1,𝑥⁽³⁾_𝑖+1=−𝜔𝑥^2•_𝑖+1,… (4.7) for the free vibration and nondamped mode, taking 𝜔²=𝑘/𝑚. Adequate substitution of (4.7) in (4.1) and (4.2), after determining the necessary 𝑎_𝑗 and 𝑏_𝑘, leads to a pair of relations in the form:

𝐻_𝑚,𝑛(Δ𝑡,𝑥_𝑖,𝑥^•_𝑖,𝑥_𝑖+1,𝑥^•_𝑖+1)≈0,𝑑𝐻_{𝑚+1,𝑛+1}𝑑𝑡(Δ𝑡,𝑥_𝑖,𝑥^•_𝑖,𝑥_𝑖+1,𝑥^•_𝑖+1)≈0, (4.8) and the solution of this 2×2 linear system provides the values of displacement and its first derivative in the instant 𝑖+1. Another form of (4.8) can be used, that is

𝑛₁₁𝑥_𝑖+𝑛^•₁₂𝑥_𝑖=𝑚₁₁𝑥_𝑖+1+𝑚^•₁₂𝑥_𝑖+1,𝑛₂₁𝑥_𝑖+𝑛^•₂₂𝑥_𝑖=𝑚₂₁𝑥_𝑖+1+𝑚^•₂₂𝑥_𝑖+1, (4.9) or yet, in the traditional form or single-step algorithms:

{𝑥_𝑖+1}=[𝑀]⁻¹[𝑁]{𝑥_𝑖}=[𝐴]{𝑥_𝑖}, (4.10) where [𝐴] is the so-called amplification matrix of the method. Several properties of the method may be studied by means of a spectral analysis of [𝐴], what has been done in [2] for various combinations of 𝑚 and 𝑛 in the Hermitian operators. The goal of unconditional stability and asymptotic annihilation has been achieved when using the operators defined by the coefficients represented in the columns of Table 1, where the error order in accelerations and internal forces is given by 𝑟−2.

Table 1 

Coefficients of H_m,n type operators versus internal forces error order.

The authors performed also a spectral analysis when viscous damping terms are present, for methods in the range from 1 to 4 for the internal forces local error order, and the mentioned desired properties are still present, what is shown in [4, 5].

The forced case requires the inclusion of additional terms in (4.9) brought by the relations provided by the movement equation in the form (2.1), similar to those, present in (4.7), now with the aspect shown below, since we state that 𝛾=𝑐/2𝑚𝜔:

^•𝑥^•_𝑖=𝑓_𝑖−𝜔²𝑥_𝑖−2𝛾𝜔𝑥^•_𝑖,𝑥_𝑖⁽³⁾=𝑓^•_𝑖−𝜔𝑥^2•_𝑖−2𝛾𝜔𝑥^••_𝑖,…,_••𝑥_𝑖+1=𝑓_𝑖+1−𝜔²𝑥_𝑖+1−2𝛾𝜔𝑥^•_𝑖+1,𝑥⁽³⁾_𝑖+1=𝑓^•_𝑖+1−𝜔𝑥^2•_𝑖+1−2𝛾𝜔𝑥^••_𝑖+1,…. (4.11)

The resulting expressions carry now the influence of external forces and take the form:

−𝑛₁₁𝑥_𝑖−𝑛^•₁₂𝑥_𝑖+𝑚₁₁𝑥_𝑖+1+𝑚^•₁₂𝑥_𝑖+1=𝑙₁₁𝑓_𝑖+1+𝑙^•₁₂𝑓_𝑖+1+𝑙𝑖₁₁𝑓_𝑖+𝑙𝑖^•₁₂𝑓_𝑖,−𝑛₂₁𝑥_𝑖−𝑛^•₂₂𝑥_𝑖+𝑚₂₁𝑥_𝑖+1+𝑚^•₂₂𝑥_𝑖+1=𝑙₂₁𝑓_𝑖+1+𝑙^•₂₂𝑓_𝑖+1+𝑙𝑖₂₁𝑓_𝑖+𝑙𝑖^•₂₂𝑓_𝑖, (4.12) where the 𝑓 and its derivative terms are representing the external forces present, so that the shape of the function describing the variation of 𝑓 in time influences the performance of the method. This is not an unusual fact and explains why free vibration problems, like in the present work, are the first choice to study new integration methods, once characteristics present are purely theirs. Once the reliability is verified, further studies, for various 𝑓 functions, may lead to more detailed conclusions in each situation:

{𝑥_𝑖+1}=[𝑀]⁻¹([𝑁]{𝑥_𝑖}+[𝐿]{𝑓_𝑖+1}+[𝐿𝑖]{𝑓_𝑖}). (4.13)

The solution of (4.13) allows advances in understanding the movement history. The aforementioned studies include details from methods, coefficients involved, and resulting properties. Noting 𝜃=𝜔⋅Δ𝑡, for the second-order algorithm (HdH2) the coefficient in (4.12) is written as

𝑚₁₁=6−𝜃²,𝑚₁₂𝑚=−2(𝛾𝜃+2),₂₁=−𝜃²𝑚₁₂,𝑚₂₂=[6+8𝛾𝜃−𝜃²+4(𝛾𝜃)²𝑛],₁₁𝑛=6,₁₂𝑛=2,₂₁=−𝜃²𝑛₁₂,𝑛₂₂𝑙=6−4𝛾𝜃,₁₁𝑙=−1,₁₂𝑙=0,₂₁𝑙=4+2𝛾𝜃,₂₂=−1,𝑙𝑖₁₁=𝑙𝑖₁₂=𝑙𝑖₂₂=0,𝑙𝑖₂₁=2, (4.14) and their determination is sufficient for applying the method. Other order HdHr group defining terms may be found in [5].

4.2. The Nonlinear Case

The previously presented expression (2.3) is used for non-linear case, with known initial conditions 𝑥(0)=𝑥₀,𝑥^•(0)=𝑥^•₀. However, it remains the need of determining

𝑘(𝑥_𝑖+1),𝑓_𝑒𝑥𝑡(𝑥,𝑥^•,𝑡)|_𝑖+1 (4.15) in order to attend to equilibrium in the 𝑡_𝑖+1 instant, and the solution to the problem consists of finding displacement and derivatives for 𝑡_𝑖+1, which result in values given in (4.15) that attend to equilibrium to advance in the movement.

The direct use of the equilibrium equation and its derivatives to follow the path shown for the linear case is not possible, regarding (4.15). One may, however, estimate a displacement value in order to evaluate stiffness and external force terms, leading to

𝑘_𝑒𝑠𝑡=𝑘(𝑥_{𝑖+1,𝑒𝑠𝑡}),𝑓_𝑒𝑥𝑡(𝑥,𝑥^•,𝑡)|_{𝑖+1,𝑒𝑠𝑡}. (4.16)

The integration algorithm applied in these conditions leads to an approximation of the displacement and its derivatives, which generally do not attend to equilibrium. A remainder may be then evaluated, using the following expression:

𝑓_𝑒𝑥𝑡(𝑥_{𝑖+1,𝑎𝑝},𝑥^•_{𝑖+1,𝑎𝑝},𝑡_𝑖+1)−(𝑚𝑥^••_{𝑖+1,𝑎𝑝}+𝑘_𝑒𝑠𝑡𝑥_{𝑖+1,𝑎𝑝})=𝜀(𝑥_{𝑖+1,𝑎𝑝},𝑥^•_{𝑖+1,𝑎𝑝}). (4.17)

This deviation from equilibrium condition may be understood as a residual external force, whose value may be acceptable or not, depending on a previous conveniently established limit. If it exceeds it, a closer displacement estimative must be used to perform the described procedure. If the estimates are successively and consistently better, the remainder will be below the previously established limit after a certain number of trials, and an approximation good enough for the displacement and its derivatives in the 𝑡_𝑖+1 instant will be obtained.

Different methods are basically defined by different ways of estimating displacements in (4.16) and integration schemes used.

4.3. The Proposed Procedure

In order to estimate displacement in the 𝑡_𝑖+1 instant, constant acceleration was considered throughout the step, that is,

𝑥_{𝑖+1,𝑒𝑠𝑡}=𝑥_𝑖+Δ𝑡𝑥^•_𝑖+Δ𝑡²2^••𝑥_𝑖,_•𝑥_{𝑖+1,𝑒𝑠𝑡}=𝑥^•_𝑖+Δ𝑡𝑥^••_𝑖,_••𝑥_{𝑖+1,𝑒𝑠𝑡}=𝑥^••_𝑖=𝑐𝑡𝑒.𝑛𝑜𝑖𝑛𝑡𝑒𝑟𝑣𝑎𝑙𝑜 (4.18) in (4.16) and a set of expressions analogous to those in (4.7), also needed to the Hermitian algorithms already presented for the linear case. The difference is that now estimates were used to represent the forces involved in equilibrium equation. The integration scheme looks now as follows:

{𝑥_{𝑖+1,𝑎𝑝}}=[𝑀_𝑒𝑠𝑡]⁻¹([𝑁]{𝑥_𝑖}+[𝐿]{𝑓_{𝑖+1,𝑒𝑠𝑡}}+[𝐿𝑖]{𝑓_𝑖}) (4.19) and for the first-order Hermitian algorithm (HdH1), the expressions for the matrix coefficients [𝑀_𝑒𝑠𝑡], analogous to those already presented in [3] for the linear case, for free vibration and unitary mass, that is,

𝑚₁₁=−1+Δ𝑡²2𝑘(𝑥_{𝑖+1,𝑒𝑠𝑡}𝑚),₁₂𝑚=Δ𝑡,₂₁=−𝑘(𝑥_{𝑖+1,𝑒𝑠𝑡})Δ𝑡²2,𝑚₂₂=[𝑘(𝑥_{𝑖+1,𝑒𝑠𝑡})Δ𝑡²2−1]Δ𝑡, (4.20) while the [𝑁] ones continue the same, already presented, taking the values shown below:

𝑛₁₁𝑛=−1,₁₂𝑛=0,₂₁𝑛=0,₂₂=−Δ𝑡. (4.21)

The coefficients of these matrixes for the other methods in this Hermitian group, with local error precision of second, third, and fourth orders (HdH2, HdH3, and HdH4), are obtained from those for the linear case, found in the already mentioned work. Once the solution of (4.19) is at hand, the remainder can be calculated through (4.17). If it is lesser than a previously established limit, the approximate value for displacement and its derivatives is accepted, and one advances in the movement history representation. If not, a new estimate of the matrixes in (4.16) is done to perform the integration scheme. The approximate values given by the solution of (4.19) are now used as new estimates and the procedure repeated. In other words, the system is solved again using the new estimates. This generates a new approximation for displacement and velocity, a new remainder is calculated, its value is verified, and this procedure is continued until its limit is satisfied. The iterations for the considered step are then finalized, and one advances, initializing the procedure for the next step until the time interval of interest is completely covered. Algorithm 1summarizes the steps constituting the proposed procedure.

1. Start    Read 𝑥0,𝑥•0,Δ𝑡,𝑡𝑓,𝜉lim    Calculate •𝑥•0=[𝑓𝑒𝑥𝑡(𝑥0,𝑥•0,𝑡0)−𝑘𝑥0](1/𝑚)  2. Advance in movement               𝑡𝑛=𝑓,Δ𝑡𝑖=0,𝑛  2.1. Estimate of the (restorative) and external forces:  2.1.1. Estimate of displacement in the 𝑖+1 instant          𝑥𝑖+1,𝑒𝑠𝑡=𝑥𝑖+Δ𝑡𝑥•𝑖+Δ𝑡22••𝑥𝑖,•𝑥𝑖+1,𝑒𝑠𝑡=𝑥•𝑖+Δ𝑡𝑥••𝑖.  2.1.2. Estimate of the matrixes for the internal forces            𝑘𝑒𝑠𝑡=𝑘(𝑥𝑖+1,𝑒𝑠𝑡)  2.1.3. Estimate of the external forces    Calculate 𝑓𝑒𝑥𝑡(𝑥𝑖+1,𝑒𝑥𝑡,𝑥•𝑖+1,𝑒𝑠𝑡,𝑡𝑖+1)  2.2. Evaluation of approximate displacements and velocities       {𝑥𝑖+1,𝑎𝑝}=[𝑀𝑒𝑠𝑡]−1([𝑁]{𝑥𝑖}+[𝐿]{𝑓𝑖+1,𝑒𝑠𝑡}+[𝐿𝑖]{𝑓𝑖})  2.3. Evaluation of approximate acceleration   {𝑥••𝑖+1,𝑎𝑝}=[𝑓𝑒𝑥𝑡(𝑥𝑖+1,𝑎𝑝,𝑥•𝑖+1,𝑎𝑝,𝑡𝑖+1)−𝑘𝑒𝑠𝑡𝑥𝑖+1,𝑎𝑝]∗(1/𝑚)  2.4. Updating of the internal force vector            𝑘𝑒𝑠𝑡=𝑘(𝑥𝑖+1,𝑎𝑝)  2.5. Evaluation of the remainder   𝑓𝑒𝑥𝑡(𝑥𝑖+1,𝑎𝑝,𝑥•𝑖+1,𝑎𝑝,𝑡𝑖+1)−𝑘𝑒𝑠𝑡𝑥𝑖+1,𝑎𝑝−𝑚𝑥••𝑖+1,𝑎𝑝=𝜀(𝑥𝑖+1,𝑎𝑝,𝑥•𝑖+1,𝑎𝑝)  2.6. Remainder test             if𝜀(𝑥𝑖+1,𝑎𝑝,𝑥•𝑖+1,𝑎𝑝)>𝜀lim:              returns to 2.2.             if𝜀(𝑥𝑖+1,𝑎𝑝,𝑥•𝑖+1,𝑎𝑝)≤𝜀lim:             𝑥𝑖+1⟵𝑥•𝑖+1,𝑎𝑝𝑥𝑖+1⟵𝑥•𝑖+1,𝑎𝑝••𝑥𝑖+1⟵𝑥••𝑖+1,𝑎𝑝   if𝑖<𝑛              returns to 2.1.   if𝑖=𝑛 End

Algorithm 1 

Steps for the HdHr algorithm in non-linear problems.

5. Example

The example shown is that of the simple pendulum, such as the addressed in [20]. Figure 1 illustrates its scheme and notation used. It is a punctual mass 𝑚 at the end of a rigid, weightless bar of length 𝐿, under gravity action, free to rotate around its other end. The angle of the bar in a given time is that measured positive to the right from a vertical line is θ.

Figure 1 

Simple Pendulum.

The equilibrium conditions generate the movement differential equation, given by

^•𝜃^•+𝐾²sin𝜃=0,𝐾=𝑔𝐿, (5.1) and the exact solution for the period 𝜏 for a given initial condition 𝜃₀ is written as

𝜏=2𝐿𝑔^𝜃0₀𝑑𝜃sin²(𝜃₀/2)−sin²(𝜃/2). (5.2)

Its value can be put in the form:

2𝜏=𝜋⋅𝐾₁⋅𝑇_𝑎𝑝,𝑇_𝑎𝑝=2𝜋𝐿𝑔, (5.3) with 𝑇_𝑎𝑝 being a good approximation for small values of 𝜃₀, that is, when 𝜃₀≈sin𝜃₀. Values for 2𝐾₁/𝜋, resulting from the solution of (5.2) for a given initial condition 𝜃₀, may be found in literature, like in Beer and Johnston [21], where there is a table for some values of 𝜃₀. These values were used as the correct ones in the present work.

The problem has been solved for 𝐾²=𝑔/𝐿=1, for the four methods of the HdHr group, using a truncated Taylor series 0-centered to represent sin(𝜃), that is:

𝜃sin(𝜃)=𝜃−³6+𝜃⁵−𝜃120504⁷0+⋯. (5.4) Regarding (5.1), one may then note, since we truncate the Taylor series, that

𝜃𝑘(𝜃)=(1−²6+𝜃⁴−𝜃120⁶+𝜃50403628⁸80),_••𝜃+𝑘(𝜃)𝜃=0. (5.5)

The first change in the signal for 𝜃 was controlled to estimate one fourth of the period. This corresponds to the lowest point in the pendulum trajectory. The method gives a discrete solution, that is, at one observed instant the angle is positive (the pendulum is still to the right of the vertical projection of the point of support) and in the next step its sign has changed (the pendulum has moved to the left of the projection); time has been interpolated proportionally to the observed 𝜃 values. Two situations are given as examples: 𝜃₀=𝜋/6 and 𝜃₀=5𝜋/6. Computational codes were created in Fortran Force 2.0—Fortran Compiler and Editor—language, whose software is available at the site http://force.lepsch.com/, to perform calculations corresponding to the procedure described in Algorithm 1 for several step and remainder limit values found in the pre-established movement (𝜀_lim) equations. Expressions for the first four members of the HdHr group were developed, that is, members HdH1, HdH2, HdH3, and HdH4. The answers of the resulting approximation for the period can be found in Tables 2 to 4 as well as the total number of iterations performed and their relative error, for a value of 𝜀_lim=1×10⁻⁸, chosen as small as necessary to not interfere in the approximation for the period taken with six significant figures. When dealing with forced case problems, the authors’ first approach to be considered will be the usual one of taking the remainder as an external additional force and imposing a limitation to its relative size compared to the actual one. Therefore, the content of these tables serves as the starting point for future discussion and comments as well as for the creation of two figures, also presented and commented on below.

Table 2 

Approximate range 𝑇_𝑎𝑝 and total number of iterations (N. Iter.) 𝜃₀=𝜋/6 (𝑇_exact=6.392568s)—HdH1 and HdH2 methods.

Table 3 

Approximate range 𝑇_𝑎𝑝 and total number of iterations (N. Iter.) 𝜃₀=𝜋/6 (𝑇_exact=6.392568s)—HdH3 and HdH4 methods.

Table 4 

Approximate range 𝑇_𝑎𝑝 and total number of iterations (N. Iter.) 𝜃₀=5𝜋/6 (𝑇_exact=11.07226s)—HdH1 and HdH2 methods.

Table 5 

Approximate range 𝑇_𝑎𝑝 and total number of iterations (N. Iter.) 𝜃₀=5𝜋/6 (𝑇_exact=11.07226s)—HdH3 and HdH4 methods.

Figures 2 and 3 are graphical representations of the values in the tables data and allow to visualize result tendencies. Logarithmic scales were used for the abscissas, where the relationship between the size of the step and the exact period is represented as well as in the ordinates, where the corresponding relative error is represented. Figure 2 relates to the initial angle 𝜃₀=𝜋/6 and Figure 3 to 𝜃₀=5𝜋/6.

Figure 2 

Logarithm for the relative error within range 𝑋 step estimate (𝜃0=𝜋/6).

Figure 3 

Logarithm for the relative error within the range 𝑋 step estimate (𝜃0=5𝜋/6).

6. Discussion and Conclusions

First of all, we can see that error grows more intensely for larger steps in the first-order algorithm, confirming what other authors have observed about the convenience of using higher order methods. This is clearer in Figure 2; as it occurs in small steps, along a weak non-linear path, it is reasonable to conclude that it is the effect of numerical dissipation (asymptotic annihilation) of the algorithm itself, present more intensely than in the higher order methods. In Figure 3, this was not observed as clearly, maybe because of the stronger nonlinearity, and the movement representation error begins to add to the numerical error of the algorithms.

Another noticeable point is that higher-order terms in fact lead to increasing precision, for same time steps. However, larger steps lead to greater effects of the annihilation in all methods.

The apparent convergence of the results for smaller steps must be influenced by truncation errors committed in the representation of values involved; however, we cannot conclude that they will all perform similarly at the limit.

The changing of the adopted remainder limit (for reasonable values) has little effect on results; therefore, if exaggeratedly rigorous values are adopted in comparison to the step, the number of iterations performed increases significantly.

Again, as in linear case, we observe that the HdH2 method presents results whose quality approximates those of higher-order ones. Since it has simpler algebraic expressions and demands fewer operations, studies should focus on its use. The following researches on the matter must definitely contemplate performance evaluation in the application of problems with stronger nonlinearity as well as other iteration control procedures. A stability analysis of these methods and the relative amount of numerical work involved in the application of MDOF problems are also steps that will be taken by the authors in the future.

Following works must also focus on comparing HdHr methods performance with those from traditional ones; for the linear case, this has been done for several examples extracted from literature in works treating on such matter for single-step as well as multistep semidiscrete models analysis via MEF; the results may be seen in the already mentioned references. Advantages of combining unconditional stability, high-order error, and asymptotic annihilation are clear, so far.