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

Authors: Heitor Miranda Bottura; Antonio Carlos Rigitano
Journal: Mathematical Problems in Engineering
 (2009)

---

## 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.


---

## Body

## 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)mx••+cx•+kx=f(t),
where mis mass, c is damping, and k is stiffness, while f(t)  is the external force applied as a function of time, which is represented by t. Each dot above the variable represents, as usual, a derivation in t. A comparison with the definition presented in [7] for the linear second-order differential equation demands that coefficients m,c,  and k  are independent from x  or x•  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)mx••+r(x,x•)=f(t).
in which the term r(x,x•) 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)mx••+k(x)x=f(t),
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 ti instant, truncated for the representation of displacement and derived in the ti+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:(4.1)Hm,n=∑j=0n(Δt)jajxi(j)+∑k=0m(Δt)kbkxi+1(k)+R(Δtr+2)=0,(4.2)dHm,ndt=Hm+1,n+1=∑j=1n+1(Δt)jcjxi(j)+∑k=1m+1(Δt)kdkxi+1(k)+R(Δtr+3),
where Δt=ti+1-ti  is the step, ti  is the last instant of known movement, ti+1   is the next unknown instant, R  is the remainder resulting from the truncation of the Hermitian expression, and r  is the order of the local error resulting from efforts (and accelerations). One remarks also that cj=aj-1  and dk=bk-1. The values for m  and n  may lead to guarantee the desired properties, and for the expressions chosen in setting up the HdHr group, resulted for the local error order: r=m+n-1. The coefficients aj  and bk  are obtained by applying Taylor series centered on the i  instant to represent the displacement and its derivatives in the i+1  instant (4.2), as shown in [2]. Specifically, for m=n=1, (4.1) writes(4.3)H1,1=a0xi+b0xi+1+Δta1x•i+Δtb1x•i+1+R(Δt3)=0
or, performing the mentioned expansion in Taylor series,(4.4)xi+1=xi+Δtx•i+Δt22x••i+Δt36x(3)i+⋯,x•i+1=x•i+Δtx••i+Δt22xi(3)+Δt36xi(4)+⋯,
leading to
(4.5)H1,1=a0xi+b0(xi+Δtx•i+Δt22x••i+Δt36xi(3)+⋯)+Δta1x•i+Δtb1(x•i+Δtx••i+Δt22xi(3)+Δt36xi(4)+⋯)+R(Δt3)=0,
and placing a0=1, canceling terms or orders 1 and 2 in Δt  requires(4.6)a0+b0=0,b0+a1+b1=0,b02+b1=0,
and therefore this Hermitian operator in the form of (4.1) has the coefficients a0=1,b0=-1,a1=b1=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(4.7)x••i=-ω2xi,xi(3)=-ω2x•i,…,x••i+1=-ω2xi+1,xi+1(3)=-ω2x•i+1,…
for the free vibration and nondamped mode, taking ω2=k/m. Adequate substitution of (4.7) in (4.1) and (4.2), after determining the necessary aj  and bk,  leads to a pair of relations in the form:(4.8)Hm,n(Δt,xi,x•i,xi+1,x•i+1)≈0,dHm+1,n+1dt(Δt,xi,x•i,xi+1,x•i+1)≈0,
and the solution of this 2×2  linear system provides the values of displacement and its first derivative in the instant i+1.  Another form of (4.8) can be used, that is(4.9)n11xi+n12x•i=m11xi+1+m12x•i+1,n21xi+n22x•i=m21xi+1+m22x•i+1,
or yet, in the traditional form or single-step algorithms:(4.10){xi+1}=[M]-1[N]{xi}=[A]{xi},
where [A]  is the so-called amplification matrix of the method. Several properties of the method may be studied by means of a spectral analysis of [A], what has been done in [2] for various combinations of m  and n  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  r-2.Table 1
Coefficients ofHm,n type operators versus internal forces error order.
Order12345678a01-62460360840672015120b0-16-24-60-360-840-6720-15120a1—-262412036025206720b11-4183624048042008400a2———312603601260b2-1/21-6-9-72-120-1200-2100a3—————420120b3——111216200300a4———————5b4————-1-1-20-25a5————————b5——————11The 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 γ=c/2mω:(4.11)x••i=fi-ω2xi-2γωx•i,xi(3)=f•i-ω2x•i-2γωx••i,…,x••i+1=fi+1-ω2xi+1-2γωx•i+1,xi+1(3)=f•i+1-ω2x•i+1-2γωx••i+1,….The resulting expressions carry now the influence of external forces and take the form:(4.12)-n11xi-n12x•i+m11xi+1+m12x•i+1=l11fi+1+l12f•i+1+li11fi+li12f•i,-n21xi-n22x•i+m21xi+1+m22x•i+1=l21fi+1+l22f•i+1+li21fi+li22f•i,
where the f and its derivative terms are representing the external forces present, so that the shape of the function describing the variation  of  f  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 f functions, may lead to more detailed conclusions in each situation:(4.13){xi+1}=[M]-1([N]{xi}+[L]{fi+1}+[Li]{fi}).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  θ=ω·Δt, for the second-order algorithm (HdH2) the coefficient in (4.12) is written as(4.14)m11=6-θ2,m12=-2(γθ+2),m21=-θ2m12,m22=[6+8γθ-θ2+4(γθ)2],n11=6,n12=2,n21=-θ2n12,n22=6-4γθ,l11=-1,l12=0,l21=4+2γθ,l22=-1,li11=li12=li22=0,li21=2,
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 x(0)=x0,x•(0)=x•0. However, it remains the need of determining(4.15)k(xi+1),fext(x,x•,t)|i+1
in order to attend to equilibrium in the ti+1  instant, and the solution to the problem consists of finding displacement and derivatives for ti+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(4.16)kest=k(xi+1,est),fext(x,x•,t)|i+1,est.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:(4.17)fext(xi+1,ap,x•i+1,ap,ti+1)-(mx••i+1,ap+kestxi+1,ap)=ε(xi+1,ap,x•i+1,ap).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 theti+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 theti+1  instant, constant acceleration was considered throughout the step, that is,(4.18)xi+1,est=xi+Δtx•i+Δt22x••i,x•i+1,est=x•i+Δtx••i,x••i+1,est=x••i=cte.nointervalo
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:(4.19){xi+1,ap}=[Mest]-1([N]{xi}+[L]{fi+1,est}+[Li]{fi})
and for the first-order Hermitian algorithm (HdH1), the expressions for the matrix coefficients [Mest], analogous to those already presented in [3] for the linear case, for free vibration and unitary mass, that is,(4.20)m11=-1+Δt22k(xi+1,est),m12=Δt,m21=-k(xi+1,est)Δt22,m22=[k(xi+1,est)Δt22-1]Δt,
while the [N]  ones continue the same, already presented, taking the values shown below:(4.21)n11=-1,n12=0,n21=0,n22=-Δt.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.Algorithm 1:Steps for the HdHr algorithm in non-linear problems.
1. StartRead x0,x•0,Δt,tf,ξlimCalculate x••0=[fext(x0,x•0,t0)-kx0](1/m)2. Advance in movementn=tfΔt,i=0,n2.1. Estimate of the (restorative) and external forces:2.1.1. Estimate of  displacement in the i+1 instantxi+1,est=xi+Δtx•i+Δt22x••i,x•i+1,est=x•i+Δtx••i.2.1.2. Estimate of the matrixes for the internal forceskest=k(xi+1,est)2.1.3. Estimate of the external forcesCalculate fext(xi+1,ext,x•i+1,est,ti+1)2.2. Evaluation of approximate displacements and velocities{xi+1,ap}=[Mest]-1([N]{xi}+[L]{fi+1,est}+[Li]{fi})2.3. Evaluation of approximate acceleration{x••i+1,ap}=[fext(xi+1,ap,x•i+1,ap,ti+1)-kestxi+1,ap]*(1/m)2.4. Updating of the internal force vectorkest=k(xi+1,ap)2.5. Evaluation of the remainderfext(xi+1,ap,x•i+1,ap,ti+1)-kestxi+1,ap-mx••i+1,ap=ε(xi+1,ap,x•i+1,ap)2.6. Remainder testifε(xi+1,ap,x•i+1,ap)>εlim:returns to 2.2.ifε(xi+1,ap,x•i+1,ap)≤εlim:xi+1⟵xi+1,apx•i+1⟵x•i+1,apx••i+1⟵x••i+1,apifi<nreturns to 2.1.ifi=nEnd


## 4.1. HdHr Algorithms in Linear Form

This group of aforementioned algorithms uses a pair of Hermitian expressions:(4.1)Hm,n=∑j=0n(Δt)jajxi(j)+∑k=0m(Δt)kbkxi+1(k)+R(Δtr+2)=0,(4.2)dHm,ndt=Hm+1,n+1=∑j=1n+1(Δt)jcjxi(j)+∑k=1m+1(Δt)kdkxi+1(k)+R(Δtr+3),
where Δt=ti+1-ti  is the step, ti  is the last instant of known movement, ti+1   is the next unknown instant, R  is the remainder resulting from the truncation of the Hermitian expression, and r  is the order of the local error resulting from efforts (and accelerations). One remarks also that cj=aj-1  and dk=bk-1. The values for m  and n  may lead to guarantee the desired properties, and for the expressions chosen in setting up the HdHr group, resulted for the local error order: r=m+n-1. The coefficients aj  and bk  are obtained by applying Taylor series centered on the i  instant to represent the displacement and its derivatives in the i+1  instant (4.2), as shown in [2]. Specifically, for m=n=1, (4.1) writes(4.3)H1,1=a0xi+b0xi+1+Δta1x•i+Δtb1x•i+1+R(Δt3)=0
or, performing the mentioned expansion in Taylor series,(4.4)xi+1=xi+Δtx•i+Δt22x••i+Δt36x(3)i+⋯,x•i+1=x•i+Δtx••i+Δt22xi(3)+Δt36xi(4)+⋯,
leading to
(4.5)H1,1=a0xi+b0(xi+Δtx•i+Δt22x••i+Δt36xi(3)+⋯)+Δta1x•i+Δtb1(x•i+Δtx••i+Δt22xi(3)+Δt36xi(4)+⋯)+R(Δt3)=0,
and placing a0=1, canceling terms or orders 1 and 2 in Δt  requires(4.6)a0+b0=0,b0+a1+b1=0,b02+b1=0,
and therefore this Hermitian operator in the form of (4.1) has the coefficients a0=1,b0=-1,a1=b1=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(4.7)x••i=-ω2xi,xi(3)=-ω2x•i,…,x••i+1=-ω2xi+1,xi+1(3)=-ω2x•i+1,…
for the free vibration and nondamped mode, taking ω2=k/m. Adequate substitution of (4.7) in (4.1) and (4.2), after determining the necessary aj  and bk,  leads to a pair of relations in the form:(4.8)Hm,n(Δt,xi,x•i,xi+1,x•i+1)≈0,dHm+1,n+1dt(Δt,xi,x•i,xi+1,x•i+1)≈0,
and the solution of this 2×2  linear system provides the values of displacement and its first derivative in the instant i+1.  Another form of (4.8) can be used, that is(4.9)n11xi+n12x•i=m11xi+1+m12x•i+1,n21xi+n22x•i=m21xi+1+m22x•i+1,
or yet, in the traditional form or single-step algorithms:(4.10){xi+1}=[M]-1[N]{xi}=[A]{xi},
where [A]  is the so-called amplification matrix of the method. Several properties of the method may be studied by means of a spectral analysis of [A], what has been done in [2] for various combinations of m  and n  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  r-2.Table 1
Coefficients ofHm,n type operators versus internal forces error order.
Order12345678a01-62460360840672015120b0-16-24-60-360-840-6720-15120a1—-262412036025206720b11-4183624048042008400a2———312603601260b2-1/21-6-9-72-120-1200-2100a3—————420120b3——111216200300a4———————5b4————-1-1-20-25a5————————b5——————11The 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 γ=c/2mω:(4.11)x••i=fi-ω2xi-2γωx•i,xi(3)=f•i-ω2x•i-2γωx••i,…,x••i+1=fi+1-ω2xi+1-2γωx•i+1,xi+1(3)=f•i+1-ω2x•i+1-2γωx••i+1,….The resulting expressions carry now the influence of external forces and take the form:(4.12)-n11xi-n12x•i+m11xi+1+m12x•i+1=l11fi+1+l12f•i+1+li11fi+li12f•i,-n21xi-n22x•i+m21xi+1+m22x•i+1=l21fi+1+l22f•i+1+li21fi+li22f•i,
where the f and its derivative terms are representing the external forces present, so that the shape of the function describing the variation  of  f  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 f functions, may lead to more detailed conclusions in each situation:(4.13){xi+1}=[M]-1([N]{xi}+[L]{fi+1}+[Li]{fi}).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  θ=ω·Δt, for the second-order algorithm (HdH2) the coefficient in (4.12) is written as(4.14)m11=6-θ2,m12=-2(γθ+2),m21=-θ2m12,m22=[6+8γθ-θ2+4(γθ)2],n11=6,n12=2,n21=-θ2n12,n22=6-4γθ,l11=-1,l12=0,l21=4+2γθ,l22=-1,li11=li12=li22=0,li21=2,
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 x(0)=x0,x•(0)=x•0. However, it remains the need of determining(4.15)k(xi+1),fext(x,x•,t)|i+1
in order to attend to equilibrium in the ti+1  instant, and the solution to the problem consists of finding displacement and derivatives for ti+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(4.16)kest=k(xi+1,est),fext(x,x•,t)|i+1,est.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:(4.17)fext(xi+1,ap,x•i+1,ap,ti+1)-(mx••i+1,ap+kestxi+1,ap)=ε(xi+1,ap,x•i+1,ap).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 theti+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 theti+1  instant, constant acceleration was considered throughout the step, that is,(4.18)xi+1,est=xi+Δtx•i+Δt22x••i,x•i+1,est=x•i+Δtx••i,x••i+1,est=x••i=cte.nointervalo
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:(4.19){xi+1,ap}=[Mest]-1([N]{xi}+[L]{fi+1,est}+[Li]{fi})
and for the first-order Hermitian algorithm (HdH1), the expressions for the matrix coefficients [Mest], analogous to those already presented in [3] for the linear case, for free vibration and unitary mass, that is,(4.20)m11=-1+Δt22k(xi+1,est),m12=Δt,m21=-k(xi+1,est)Δt22,m22=[k(xi+1,est)Δt22-1]Δt,
while the [N]  ones continue the same, already presented, taking the values shown below:(4.21)n11=-1,n12=0,n21=0,n22=-Δt.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.Algorithm 1:Steps for the HdHr algorithm in non-linear problems.
1. StartRead x0,x•0,Δt,tf,ξlimCalculate x••0=[fext(x0,x•0,t0)-kx0](1/m)2. Advance in movementn=tfΔt,i=0,n2.1. Estimate of the (restorative) and external forces:2.1.1. Estimate of  displacement in the i+1 instantxi+1,est=xi+Δtx•i+Δt22x••i,x•i+1,est=x•i+Δtx••i.2.1.2. Estimate of the matrixes for the internal forceskest=k(xi+1,est)2.1.3. Estimate of the external forcesCalculate fext(xi+1,ext,x•i+1,est,ti+1)2.2. Evaluation of approximate displacements and velocities{xi+1,ap}=[Mest]-1([N]{xi}+[L]{fi+1,est}+[Li]{fi})2.3. Evaluation of approximate acceleration{x••i+1,ap}=[fext(xi+1,ap,x•i+1,ap,ti+1)-kestxi+1,ap]*(1/m)2.4. Updating of the internal force vectorkest=k(xi+1,ap)2.5. Evaluation of the remainderfext(xi+1,ap,x•i+1,ap,ti+1)-kestxi+1,ap-mx••i+1,ap=ε(xi+1,ap,x•i+1,ap)2.6. Remainder testifε(xi+1,ap,x•i+1,ap)>εlim:returns to 2.2.ifε(xi+1,ap,x•i+1,ap)≤εlim:xi+1⟵xi+1,apx•i+1⟵x•i+1,apx••i+1⟵x••i+1,apifi<nreturns to 2.1.ifi=nEnd


## 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 m  at the end of a rigid, weightless bar of length L, 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(5.1)θ••+K2sinθ=0,K=gL,
and the exact solution for the period  τ for a given initial condition θ0 is written as(5.2)τ=2Lg∫0θ0dθsin2(θ0/2)-sin2(θ/2).Its value can be put in the form:(5.3)τ=2π·K1·Tap,Tap=2πLg,
with Tap  being a good approximation for small values of θ0, that is, when θ0≈sinθ0. Values for 2K1/π, resulting from the solution of (5.2) for a given initial condition θ0, may be found in literature, like in Beer and Johnston [21], where there is a table for some values of θ0. These values were used as the correct ones in the present work.The problem has been solved forK2=g/L=1, for the four methods of the HdHr group, using a truncated Taylor series 0-centered to represent sin(θ), that is:(5.4)sin(θ)=θ-θ36+θ5120-θ75040+⋯.
Regarding (5.1), one may then note, since we truncate the Taylor series, that(5.5)k(θ)=(1-θ26+θ4120-θ65040+θ8362880),θ••+k(θ)θ=0.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: θ0=π/6 and θ0=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-8, 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 rangeTap and total number of iterations (N. Iter.) θ0=π/6 (Texact=6.392568s)—HdH1 and HdH2 methods.
Δt(s)HdH1εr%N. Iter.HdH2εr%N. Iter.1,000005.8373428.6976.3676950,3960,100006.3772510.24286.3909322.6e–2290,010006.3919589.5e–31606.3924132.4e–31600,001006.3924781.4e–315996.3925532.3e–415990,000106.3925217.4e–4159826.3925663.1e–5159820,000016.3925255.7e–41598156.3925681e–5159815Table 3
Approximate rangeTap and total number of iterations (N. Iter.) θ0=π/6 (Texact=6.392568s)—HdH3 and HdH4 methods.
Δt(s)HdH3εr%N. Iter.HdH4εr%N. Iter.1,000006,3178391.276.3453480.7460,100006.3901373.8e–2316.3915631.6e–2310,010006.3923363.6e–33026.3924751.5e–32600,001006.3925453.6e–415996.3925591.4e–415990,000106.3925663.1e–5159826.3925671.5e–5159820,000016.3925681.e–51598156.3925681.e–5159815Table 4
Approximate rangeTap and total number of iterations (N. Iter.) θ0=5π/6 (Texact=11.07226s)—HdH1 and HdH2 methods.
Δt(s)HdH1εr%N. Iter.HdH2εr%N. Iter.1.500005.6940448.668.8075620.4301,000007.25024134.54810.049319.2260,1000010.888121.67511.012780.54720,0100011.055250.1545911.066605.1e–23730,0010011.070561.5e–2276811.071695.1e–327680,0001011.072081.6e–32768111.072196.3e–4276810,0000111.072232.7e–427680611.072241.8e–4276807Table 5
Approximate rangeTap and total number of iterations (N. Iter.) θ0=5π/6 (Texact=11.07226s)—HdH3 and HdH4 methods.
Δt(s)HdH3εr%N. Iter.HdH4εr%N. Iter.1.50008.0603827.2379.5970113.3231,000009.7350312.13210.390866.1220,1000010.984540.798811.035973.3e–1770,0100011.063797.6e–254711.068853.1e–25350,0010011.071407.7e–3487711.071913.2e–338710,0001011.072169.0e–42768111.072214.5e–4276810,0000111.072241.8e–427680611.072241.8e–4276807Figures2 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 θ0=π/6  and Figure 3 to θ0=5π/6.Figure 2
Logarithm for the relative error within rangeX  step estimate (θ0=π/6).Figure 3
Logarithm for the relative error within the rangeX  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 Figure2; 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.



---

*Source: 368309-2009-11-16.xml*