SciELO - Scientific Electronic Library Online

 
vol.16 número6Crack and mechanical behavior of CFRP plate-reinforced bridge roofs under high temperature with different anchoring measuresAnalytical solution of circular tunnel in elastic-viscoplastic rock mass índice de autoresíndice de assuntospesquisa de artigos
Home Pagelista alfabética de periódicos  

Serviços Personalizados

Journal

Artigo

Indicadores

Links relacionados

Compartilhar


Latin American Journal of Solids and Structures

versão impressa ISSN 1679-7817versão On-line ISSN 1679-7825

Lat. Am. j. solids struct. vol.16 no.6 Rio de Janeiro  2019  Epub 22-Jul-2019

https://doi.org/10.1590/1679-78255609 

ORIGINAL ARTICLE

Higher-order explicit time integration methods for numerical analyses of structural dynamics

aDepartment of Mechanical Engineering Korea Army Academy at Yeongcheon, Yeongcheon-si, Gyeongsangbuk-do 38900, Republic of Korea. E-mail: kim.wooram@yahoo.com


Abstract

In this article, third- and fourth-order accurate explicit time integration methods are developed for effective analyses of various linear and nonlinear dynamic problems stated by second-order ordinary differential equations in time. Two sets of the new methods are developed by employing the collocation approach in the time domain. To remedy some shortcomings of using the explicit Runge-Kutta methods for second-order ordinary differential equations in time, the new methods are designed to introduce small period and damping errors in the important low-frequency range. For linear cases, the explicitness of the new methods is not affected by the presence of non-diagonal damping matrix. For nonlinear cases, the new methods can handle velocity dependent problems explicitly without decreasing order of accuracy. The new methods do not have any undetermined algorithmic parameters. Improved numerical solutions are obtained when they are applied to various linear and nonlinear problems.

Keywords: Linear and nonlinear dynamics; explicit time integration method; higher-order method

1 Introduction

Step-by-step direct time integrations are dominantly used for transient analyses of linear and nonlinear dynamic problems described by second-order ordinary differential equations in time. As more sophisticated spatial finite element models [1-4] are developed constantly, demands for more accurate time integration methods are also increasing to take full advantage of improved spatial models in transient analyses. Recently, numerous implicit [5-9] and explicit [10-16] time integration methods have been introduced to effectively analyze challenging dynamic problems.

Direct time integration methods are often categorized into implicit and explicit methods. Usually, implicit methods are unconditionally stable for linear analyses, but they require factorization of non-diagonal system matrices. For large systems, factorizing non-diagonal system matrices requires huge computational efforts. On the other hand, explicit methods are conditionally stable, but matrix factorization can be avoided if the mass matrix is diagonal. Accordingly, computational costs per time step of explicit methods are much lower than implicit methods. Due to this reason, explicit methods are more frequently used in analyses of large systems, such as the wave propagation and impact problems, where sizes of optimal time steps are slightly smaller than critical time steps of explicit methods.

According to the past studies [10,17], preferable attributes of explicit time integration methods can be summarized as (a) explicit methods should not require iterative solution finding procedures for velocity dependent nonlinear problems; (b) at least second-order accuracy should be ensured; (c) number of unspecified algorithmic parameters should be minimized; (d) amplification matrices obtained by applying methods to the single-degree-of-freedom problem should not present spurious eigenvalue, or it should be minimized; (e) explicit methods should not involve matrix factorization when the mass matrix is diagonal; (f) explicit methods should not involve differential or integral type calculus; (g) explicit methods should be applicable to nonlinear problems without any modifications.

The central difference (CD) method is one of the most widely used single-stage explicit time integration methods. The CD method is second-order accurate and non-dissipative. Due to its simplicity and good accuracy, the CD method is the standard explicit method of the well-known software packages [18,19]. However, matrix factorization is required in the CD method if the damping matrix is not diagonal. The velocity can be treated explicitly to remedy this, but additional computation is required to retain second-order accuracy [15]. The CD method cannot satisfy (a) and (b).

Chung and Lee developed a family of second-order accurate explicit method [10] with dissipation control capability. The Chung and Lee (CL) method can maintain explicitness in the presence of the non-diagonal damping matrix. However, the amplification matrix of the CL method has a spurious eigenvalue that influences the quality of solutions when large time steps are used. Hulbert and Chung also developed a family of second-order accurate explicit method [11]. However, the amplification matrix of the Hulbert and Chung (HC) method also has a spurious eigenvalue. The HC method can satisfy all the attributes except (d). The HC method can include a full range of dissipative cases, but the non-dissipative case becomes unconditionally unstable for any choices of time steps in the presence of the damping matrix.

Soares developed an explicit time integration method based on the weighted residual method [13]. The Soares method can be used for the analysis of wave propagation problems, but it becomes an only first-order accurate implicit method in the presence of non-diagonal damping matrix. In addition, it is not easy to use the Soares method for nonlinear analyses, because it has been developed by directly integrating the equations of linear structural dynamics in a weighted residual sense. The Soares method cannot satisfy (a), (b), (d), (e), (f) and (g).

Kim and Lee also proposed a family of second-order accurate explicit methods [14]. The Kim and Lee (KL) method can include a full range of dissipative cases and remain as a second-order explicit method in the presence of a non-diagonal damping matrix. The KL method is very effective, but (c) is not satisfied. Recently, a very accurate two-stage explicit method was introduced by Kim [16] to more effectively tackle challenging nonlinear problems.

For more precise and reliable analyses of challenging dynamic problems [20], higher-order accurate methods can be considered as a good option. Many of the standard time integration methods used in structural dynamics, such as the trapezoidal rule and the central difference method, are second-order accurate, and methods of third- or higher-order accuracy are often called higher-order methods [7,21-24]. Computational cost per time step of higher-order methods is higher than second-order methods because higher-order methods have more stages, while many of second-order methods use only one stage. However, more accurate predictions can be obtained with large time steps in higher-order methods, and they may become more efficient in getting the same prediction when compared with second-methods. In addition, some challenging nonlinear problems can only be effectively solved by using proper higher-order methods [7,8].

One of the most broadly used explicit higher-order methods is the Runge-Kutta (RK) methods. The RK methods have been used in numerous engineering areas. Traditionally, the RK methods have been used to solve first-order ordinary differential equations numerically. Recently, on the other hand, the RK methods have also been used to obtain higher-order accurate transient solutions in structural problems [25].

For second-order initial value problems, the RK methods can provide more accurate numerical solutions than second-order methods if considerably small time steps are used. For example, the convergence rate of the four-stage RK method is fourth-order, and numerical solutions of the four-stage RK method converge to the exact solution with the rate of t4 where t is the size of the time step. In second-order methods, on the other hand, their numerical solutions converge to the exact solution with the rate of t2 as t decreases. However, numerical solutions of the RK methods may become inaccurate by the excessive numerical damping if large time steps are used. For this reason, the RK methods were not recommended for the analysis of structural dynamics, while many second-order explicit methods gained popularity. Adeli, Gere, and Weaver stated that the CD method was more effective than the RK method based on their numerical results [26], and many researchers also agreed with that. In the numerical tests of Adeli et al., the fourth-order RK method introduced too excessive numerical damping in the important low-frequency range when large time steps were used.

In this work, simple and effective third- and fourth-order accurate explicit time integration methods are presented to remedy the shortcomings of the RK methods. For a fair comparison, the three-stage third-order accurate RK method and the four-stage fourth-order accurate RK method are considered. Higher-order accuracies of the new methods are achieved without increasing computational efforts when compared with the equivalent third- and fourth-order RK methods. In the development, the displacement vector is approximated over a time domain by using the known interpolation functions in time and the weighted sums of acceleration vectors. Then, the semi-discrete equations of structural dynamics are discretized in time by adopting the collocation method. After determining optimal collocation points in time, simple difference relations are obtained.

2 New explicit method

The governing equations of many dynamic problems can be expressed as

Mu¨t=fut,u˙t,t, (1a)

and the initial conditions are given by

u0=u0, (1b)

u˙0=u˙0, (1c)

where M is the mass matrix, ut is the displacement vector, u˙t is the velocity vector, u¨t is the acceleration vector, u0 is the initial displacement vector, u˙0 is the initial velocity vector, and fut,u˙t,t is the total force vector. Eq.(1) is suitable for the description of various linear and nonlinear dynamic problems. In linear structural dynamics, fut,u˙t,t in Eq.(1) is frequently expressed as

fut,u˙t,t=qt-Cu˙t-Kut, (2)

where C and K are the damping and stiffness matrices, respectively, and qt is the external force vector. By substituting Eq.(2) into Eq.(1a) and rearranging the equation, the equation of linear structural dynamics is expressed as

Mu¨t+Cu˙t+Kut=qt. (3)

Occasionally, Eq.(3) can also be obtained by spatially discretizing original governing partial differential equations in space and time. Then, Eq.(3) is also called the semi-discrete equation.

In this section, two sets of new higher-order explicit methods are developed to discretize Eqs.(1) and (3) in time. The single-degree-of-freedom case of the linear structural dynamics equation given in Eq.(3) is used in the procedures of the development.

2.1 New third-order accurate explicit method

For the development of the new third-order accurate explicit method, three sets of interpolation functions in time are used to approximate the displacement and velocity vectors over the time domain 0t- t , where t-=t-ts is the local time, and ts is the beginning of the s+1 -th time step or the end of the s -th time step. In the new third-order method, the displacement and velocity vectors at ts+τ1t are approximate as

uts+τ1tu-ts+τ1t=ϕ11τ1tuts+ϕ21τ1tu˙ts+12τ1t2α11u¨ts, (4)

u˙ts+τ1tu-˙ts+τ1t=ϕ11τ1tu˙ts+ϕ21τ1tu¨ts, (5)

where τ1 determines the time point that the first approximated variables are associated with, u¨ts is computed explicitly by using Eq.(1a) at ts , ϕi1t- s are the known linear functions, and α11 is the weighting parameter of the acceleration vector in the first approximation of the displacement vector. ϕi1t- s used in Eqs.(4) and (5) are given by

ϕ11t-=1,ϕ21t-=t-. (6)

The displacement and velocity vectors at ts+τ2t are approximate as

uts+τ2tu-ts+τ2t=ϕ12τ2tuts+ϕ22τ2tu˙ts+ϕ32τ2tu˙ts+τ1t+12τ2t2α21u¨ts+α22u¨ts+τ1t, (7)

u˙ts+τ2tu-˙ts+τ2t=ϕ12τ2tu˙ts+ϕ22τ2tu¨ts+ϕ32τ2tu˙ts+τ1t, (8)

where τ2 determines the time point that the first approximated variables are associated with, u¨ts+τ1t is computed explicitly by using Eq.(1a) at ts+τ1t , ϕi2t- s are the known quadratic interpolation functions, and α21 and α22 are the weighting parameters of the acceleration vectors in the second approximation of the displacement vector. ϕi2t- s used in Eqs.(7) and (8) are given by

ϕ12t-=1,ϕ22t-=-t-t--2τ1t2τ1t,ϕ32t-=t-22τ1t. (9)

Finally, the displacement and velocity vectors at ts+t are updated as

uts+tuts+t=ϕ13tuts+ϕ23tu˙ts+ϕ33tu˙ts+τ1t+ϕ43tu˙ts+τ2t+12t2α31u¨ts+α32u¨ts+τ1t+α33u¨ts+τ2t, (10)

u˙ts+tu˙ts+t=ϕ13tu˙ts+ϕ23tu¨ts+ϕ33tu¨ts+τ1t+ϕ43tu¨ts+τ2t, (11)

where u¨ts+τ2t is computed explicitly by using Eq.(1a) at ts+τ2t , ϕi3t- s are the known cubic interpolation functions, and α31 , α32 , and α33 are the weighting parameters of the acceleration vectors in the third approximation of the displacement vector. ϕi3t- s used in Eqs.(10) and (11) are given by

ϕ31(t¯)=1,ϕ32(t¯)=t¯(2t¯23τ1Δtt¯3τ2Δtt¯+6τ1τ2Δt2)6τ1τ2Δt2,ϕ33(t¯)=t¯2(2t¯+3τ2Δt)6τ1(τ2+τ1)Δt2,ϕ34(t¯)=t¯2(2t¯+3τ1Δt)6τ2(τ1τ2)Δt2. (12)

To achieve third-order accuracy, extended stability, and preferable spectral properties, permissible τi and αij are chosen as

τ1=13,τ2=23,α11=1,α21=23,α22=23,α31=13,α32=23,α33=13. (13)

By using the parameters given in Eq.(13) and rearranging the equations, the new third-order method are summarized as follow: To start the new third-order explicit method, compute the acceleration vector at ts as

u¨ts=M-1futs,u˙ts,ts, (14)

where uts and u˙ts are known properties from s -th time step. By using Eq.(14), the displacement and velocity vectors at ts+t/3 are approximated as

u-ts+t/3=uts+13tu˙ts+118t2u¨ts, (15)

u-˙ts+t/3=u˙ts+13tu¨ts. (16)

By using u-ts+t/3 and u-˙ts+t/3 given in Eqs.(15) and (16), the acceleration vector at ts+t/3 is computed as

u¨ts+t/3=M-1fu-ts+t/3,u-˙ts+t/3,ts+t/3. (17)

By using Eq.(17), the displacement and velocity vectors at ts+2t/3 are approximated as

u-ts+2t/3=uts+23tu˙ts+127t22u¨ts+4u¨ts+t/3, (18)

u-˙ts+2t/3=u˙ts+23tu¨ts+t/3. (19)

By using u-ts+2t/3 and u-˙ts+2t/3 given in Eqs.(18) and (19), the acceleration vector at ts+2t/3 is computed as

u¨ts+2t/3=M-1fu-ts+2t/3,u-˙ts+2t/3,ts+2t/3. (20)

Finally, the displacement and velocity vectors at ts+t are computed as

uts+t=uts+tu˙ts+16t2u¨ts+u¨ts+t/3+u¨ts+2t/3, (21)

u˙ts+t=u˙ts+14tu¨ts+3u¨ts+2t/3. (22)

As shown in Eqs.(14)-(22), M-1 is the only matrix factorization required through out the entire procedure. This trait is commonly found in the existing explicit methods, and total computational effort of analyses can be reduced dramatically when M is diagonal. Also, some recent numerical integration techniques [3] can produce diagonal mass matrix automatically if low order spatial elements are considered. After computing uts+t and u˙ts+t , uts and u˙ts are updated as uts+t and u˙ts+t to advance another step, and procedures can be repeated from Eq.(14) to Eq.(22).

2.2 New fourth-order accurate explicit method

The new fourth-order method is developed by considering an additional stage. In the new fourth-order accurate method, the first and second approximations of the displacement and velocity vectors take the same forms as the approximations used in the third-order method given in Eqs.(4)-(8). The third approximations of the fourth-order method at ts+τ3t take the forms of

uts+τ3tu-ts+τ3t=ϕ13τ3tuts+ϕ23τ3tu˙ts+ϕ33τ3tu˙ts+τ1t+ϕ43τ3tu˙ts+τ2t+12τ3t2α31u¨ts+α32u¨ts+τ1t+α33u¨ts+τ2t, (23)

u˙ts+τ3tu-˙ts+τ3t=ϕ13τ3tu˙ts+ϕ23τ3tu¨ts+ϕ33τ3tu¨ts+τ1t+ϕ43τ3tu¨ts+τ2t, (24)

where τ3 determines the time point that the third approximated variables are associated with, and u¨ts+τ3t is computed explicitly by using Eqs.(7)-(8) and Eq.(1a) at ts+τ3t . The displacement and velocity vectors at ts+t are updated as

uts+tuts+t=ϕ14tuts+ϕ24tu˙ts+ϕ34tu˙ts+τ1t+ϕ44tu˙ts+τ2t+ϕ54tu˙ts+τ3t+12t2α41u¨ts+α42u¨ts+τ1t+α43u¨ts+τ2t+α44u¨ts+τ3t, (25)

u˙ts+tu˙ts+t=ϕ14tu˙ts+ϕ24tu¨ts+ϕ34tu¨ts+τ1t+ϕ44tu¨ts+τ2t+ϕ54tu¨ts+τ3t, (26)

where u¨ts+τ3t is computed explicitly by using Eqs.(23)-(24) and Eq.(1a) at ts+τ3t , ϕi4t- s are the known quartic interpolation functions, and α41 , α42 , α43 , and α44 are the weighting parameters of the acceleration vectors in the fourth approximation of the displacement vector. ϕi4t- s used in Eqs.(25) and (26) are given by

ϕ41(t¯)=1,ϕ42(t¯)=t¯(3t¯34Δtt¯2(τ1+τ2+τ3)+6Δt2t¯(τ1τ2+τ2τ3+τ3τ1))12τ1τ2τ3Δt3+t¯,ϕ43(t¯)=t¯2(3t¯24τ2Δtt¯4τ3Δtt¯+6τ2τ3Δt2)12τ1(τ1τ2+τ2τ3τ3τ1+τ12)Δt3,ϕ44(t¯)=t¯2(3t¯24τ1Δtt¯4τ3Δtt¯+6τ1τ3Δt2)12τ2(τ1τ2τ2τ3+τ3τ1+τ22)Δt3,ϕ45(t¯)=t¯2(3t¯24τ1Δtt¯4τ2Δtt¯+6τ1τ2Δt2)12τ3(τ1τ2τ2τ3τ3τ1+τ32)Δt3. (27)

To achieve fourth-order accuracy, extended stability, and preferable spectral properties, permissible τi and αij are chosen as

τ1=13,τ2=12,τ3=1,α11=1,α21=35,α22=35,α31=35,α32=35,α33=0,α41=α42=α43=α44=0. (28)

With the parameters given in Eq.(28), the new fourth-order method is summarized as follows: To start the new fourth-order explicit method, compute the acceleration vector at ts as

u¨ts=M-1futs,u˙ts,ts, (29)

where uts and u˙ts are known properties from s -th time step. By using Eq.(14), the displacement and velocity vectors at ts+t/3 are approximated as

u-ts+t/3=uts+13tu˙ts+118t2u¨ts, (30)

u-˙ts+t/3=u˙ts+13tu¨ts. (31)

By using u-ts+t/3 and u-˙ts+t/3 given in Eqs.(30) and (31), the acceleration vector at ts+t/3 is computed as

u¨ts+t/3=M-1fu-ts+t/3,u-˙ts+t/3,ts+t/3. (32)

By using Eq.(32), the displacement and velocity vectors at ts+t/2 are approximated as

u-ts+t/2=uts+12tu˙ts+140t22u¨ts+3u¨ts+t/3, (33)

u-˙ts+t/2=u˙ts+18tu¨ts+3u¨ts+t/3. (34)

By using u-ts+t/2 and u-˙ts+t/3 given in Eqs.(33) and (34), the acceleration vector at ts+t/2 is computed as

u¨ts+t/2=M-1fu-ts+t/2,u-˙ts+t/2,ts+t/2. (35)

The interim displacement and velocity vectors at ts+t are approximated as

u-ts+t=uts+tu˙ts+120t2u¨ts+9u¨ts+t/3, (36)

u-˙ts+t=u˙ts+12tu¨ts-3u¨ts+t/3+4u¨ts+t/2. (37)

By using u-ts+t and u-˙ts+t given in Eqs.(36) and (37), the acceleration vector at ts+t is computed as

u¨ts+t=M-1fu-ts+t,u-˙ts+t,ts+t. (38)

Finally, the displacement and velocity vectors at ts+t are computed as

uts+t=uts+tu˙ts+16t2u¨ts+2u¨ts+t/2, (39)

u˙ts+t=u˙ts+16tu¨ts+4u¨ts+t/2+u¨ts+t. (40)

After computing uts+t and u˙ts+t , uts and u˙ts are updated as uts+t and u˙ts+t to advance another step, and procedures can be repeated from Eq.(29) to Eq.(40).

3 Review of the Runge-Kutta methods

The RK Methods are frequently used to analyze various linear and nonlinear dynamic systems described by first-order differential equations. Traditionally, the RK Methods have been considered unsuitable for the analysis of dynamic problems stated by second-order differential equations due to excessive numerical damping. Recently, however, the RK Methods have also been employed to develop a more accurate numerical procedure which can be used for the analysis of structural dynamics [25]. In this section, the third- and fourth-order Runge-Kutta methods are reviewed briefly.

3.1 Third-order Runge-Kutta Method

In the third-order RK method, u¨ts is computed as

u¨ts=M-1futs,u˙ts,ts. (41)

The displacement and velocity vectors at ts+t/2 are approximated as

u-ts+t/2=uts+12tu˙ts, (42)

u-˙ts+t/2=u˙ts+12tu¨ts. (43)

u¨ts+t/2 is computed as

u¨ts+t/2=M-1fu-ts+t/2,u-˙ts+t/2,ts+t/2. (44)

The interim displacement and velocity vectors at ts+t are approximated as

u¯ts+Δt=uts+Δt(u˙ts+2u¯˙ts+Δt)=uts+Δtu˙ts+Δt2u¨ts, (45)

u-˙ts+t=u˙ts+t-u¨ts+2u¨ts+t/2. (46)

u¨ts+t is computed as

u¨ts+t=M-1fu-ts+t,u-˙ts+t,ts+t. (47)

Finally, the displacement and velocity vectors at ts+t are computed as

uts+Δt=uts+16Δt(u˙ts+4u¯˙ts+Δt/2+u¯˙ts+Δt)=uts+Δtu˙ts+16Δt2(u¨ts+2u¨ts+Δt/2), (48)

u˙ts+t=u˙ts+16tu¨ts+4u¨ts+t/2+u¨ts+t. (49)

As shown in Eqs.(41)-(49), computational structures of the new third-order method and the third-order RK method are very similar. It can also be observed that each method contains three major computations (i.e., M-1f ). Due to this reason, computational cost per a step of each method is almost identical.

3.2 Fourth-order Runge-Kutta Method

In the fourth-order RK method, u¨ts is computed as

u¨ts=M-1futs,u˙ts,ts. (50)

The displacement and velocity vectors at ts+t/2 are approximated as

u-ts+t/2=uts+12tu˙ts, (51)

u-˙ts+t/2=u˙ts+12tu¨ts. (52)

u-¨ts+t/2 is computed as

u-¨ts+t/2=M-1fu-ts+t/2,u-˙ts+t/2,ts+t/2. (53)

The displacement and velocity vectors at ts+t/2 are approximated as

u¯¯ts+Δt/2=uts+12Δtu¯˙ts+Δt/2=uts+12Δtu˙ts+Δt/2+14Δt2u¨ts, (54)

u--˙ts+t/2=u˙ts+12tu-¨ts+t/2. (55)

u¨ts+t/2 is computed as

u¨ts+t/2=M-1fu--ts+t/2,u--˙ts+t/2,ts+t/2. (56)

The interim displacement and velocity vectors at ts+t are approximated as

u¯ts+Δt=uts+Δtu¯¯˙ts+Δt/2=uts+Δtu˙ts+12Δt2u¯¨ts+Δt/2, (57)

u-˙ts+t=u˙ts+tu¨ts+t/2. (58)

u¨ts+t is computed as

u¨ts+t=M-1fu-ts+t,u-˙ts+t,ts+t. (59)

Finally, the displacement and velocity vectors at ts+t are computed as

uts+Δt=uts+16Δt(u˙ts+2u¯¯˙ts+Δt/2+2u¯˙ts+Δt/2+u¯˙ts+Δt)=uts+Δtu˙ts+16Δt2(u¨ts+u¯¨ts+Δt/2+u¨ts+Δt/2), (60)

u˙ts+t=u˙ts+16tu¨ts+2u-¨ts+t/2+2u¨ts+t/2+u¨ts+t. (61)

As shown in Eqs.(50)-(61), computational structures of the new fourth-order method and the fourth-order RK method are very similar. It can also be observed that each method contains four major computations (i.e., M-1f ). Due to this reason, computational cost per a step of each method is almost identical.

4 Analysis of the new methods

In this section, numerical characteristics of the new methods are analyzed. To investigate numerical characteristics of time integration methods, the linear single-degree-of-freedom problem [27-29] is frequently used. The linear single-degree-of-freedom problem is given as

u¨t+2ξωu˙t+ω2ut=qt, (62)

with the initial conditions

u0=u0,u˙0=v0. (63)

where ut is the displacement, u0 is the initial displacement, v0 is the initial velocity, ξ is the damping ratio, ω is the natural frequency, and qt is the external force. By setting M=1 , C=2ξω , K=ω2 , u=u , and q=q and applying Eqs.(14)-(22) to Eqs.(62)-(63), fully discrete relation of the new third-order method at Δt is given by

xΔt=A3rdx0+L3rdq3rd, (64)

where A3rd is the amplification matrix of the new third-order method, L3rd is the load operator matrix of the new third-order method, xΔt=uΔt,vΔtT , x0=u0,v0T , and q3rd=q0,qΔt/3,q2Δt/3T . By applying Eqs.(29)-(40) to Eqs.(62)-(63), fully discrete relation of the new fourth-order method at Δt is given by

xΔt=A4thx0+L4thq4th, (65)

where A4th is the amplification matrix of the fourth-order method, L4th is the load operator matrix of the fourth-order method, and q4th=q0,qΔt/3,qΔt/2,qΔtT .

4.1 Order of accuracy

To investigate the order of accuracy of the new methods for linear problems, the characteristic polynomial of the amplification matrices is used. For a time integration method, the characteristic polynomial of the amplification matrix A is expressed as

detA-λI=λ2-2A1λ+A2, (66)

where λ is the eigenvalue of A , and A1 and A2 are the invariants of A . By using A1 and A2 , the local truncation error [28,29] is defined as

τts=1Δt2uts+Δt-2A1uts+A2uts-Δt, (67)

where ut is the exact solution of Eq.(62). For the case of qt=0 , the exact solution is given by

ut=exp-ξωtcosωdt+ξωωdsinωdtu0+1ωdsinωdtv0, (68)

where ωd=1-ξ2 is the damped natural frequency. Then a method is k -th order accurate if

τets=OΔtk. (69)

The new third-order method gives

τets3rd=1108ξω436ξ2ω-5ωu0+72ξ3-28ξv0Δt3+OΔt4. (70)

According to Eq.(69), the new three-stage method is third-order accurate for the damped case ( ξ0 ), and it becomes forth-order accurate for the undamped case ( ξ=0 ). Thus, the new third-order method is expected to give more accurate solutions for the undamped cases when compared with the equivalent third-order RK method. The new four-stage method gives

τets4th=-1720ω5-52ξ2ω+96ξ4ω+ωu0+192ξ5-152ξ3+16ξv0Δt4+OΔt5, (71)

and it is fourth-order accurate for the damped and undamped cases according to Eq.(71).

4.2 Spectral radius and stability limit

The spectral radius is frequently used to investigate linear stability of time integration method. The spectral radius is defined as

ρA=maxλ1,λ2, (72)

where λ1 and λ2 are the eigenvalues of A . A method is stable if 0<ρA1 is provided. The spectral radii of the methods are presented in Fig.1. For the undamped case, the critical time step of the new third-order method is 0.574976T , where T=2π/ω is the natural period. Thus, the new third-order method is stable if 0 Δt<0.574976T is provided. The critical time step of the new fourth-order method is 0.474023T . The stability limit slightly decreases as ξ increases from 0.0 to 1.0 in the new methods. However, similar phenomena have been observed in many well-accepted explicit methods [10,11,30].

Figure 1: Comparison of spectral radius ρA for varying values of Δt/T

4.3 Period and damping error

Figure 2: Comparison of period error T-T-/T for varying Δt/T

Figure 3: Comparison of damping error ξ for varying Δt/T

The order of accuracy is one of the most important measurements for accuracy. However, the period elongation and the damping ratio are more practical measurements, because they can directly describe important characteristics of time integration algorithms. The relative period error is defined as T--T/T where, the exact period is T=2π/ω , and the numerically obtained period is T-=2π/ω- . The algorithmic damping ratio is defined as ξ-=-lnA2/2Ω- , where Ω-=acrtanA2/A12-1 /Δt1-ξ2 and Ω-=ωΔt- . For more details, please see Ref.[28].

The new third- and fourth-order explicit methods have much smaller damping and period errors when compared with the third- and fourth-order RK methods as shown in Figs.2 and 3. Based on the results presented in Fig.2, it can be expected that the new methods will give very small period error in long-term analyses. This will be verified by using various numerical examples in the next section. Fig.3 is showing that the new explicit methods do not introduce excessive numerical (or algorithmic) damping into the important low-frequency range, while considerable numerical damping is introduced in the RK methods.

5 Numerical examples

In this section, six illustrative examples are considered to demonstrate improved performances of the new methods. To keep the simplicity of the numerical analysis, only dimensionless problems are considered. Among the examples, two-degree-of-freedom nonlinear spring-pendulum and double pendulum problems are used to verify the performance of the proposed method for multi-degree-of-freedom problems. It should also be noted that the precision of 16 significant digits is used for evaluations of variables in the computer code.

As shown in Table 1, the computation times of the four different higher-order explicit methods are almost identical when the number of degree-of-freedoms is less than hundred. As degree-of-freedoms increase, on the other hand, the computation times of the third-order methods become approximately 75% of the fourth-order methods. As expected, computational times of the new explicit methods and the RK methods are almost the same when the same numbers of stages are assumed.

Table 1: Computation times (in seconds) taken to obtain nonlinear numerical solutions of 10,000-th time step for varying sizes of system (i.e., degree-of-freedoms). 

Degree-of -freedoms New(3rd) Runge-Kutta(3rd) New(4th) Runge-Kutta(4th)
1 9.5 9.4 9.7 9.7
10 9.8 9.8 10.1 10.0
100 10.6 10.5 10.9 10.9
1,000 107.7 107.6 136.2 135.5
10,000 8186.2 8185.7 10866.2 10865.5

5.1 Single-degree-of-freedom problem

Eqs.(62) and (63) are solved by using the new methods and the RK methods. For undamped cases, u0=1.0 , v0=0.0 , ξ=0.0 , ft=0.0 , and ω=2π are used. For damped cases, u0=1.0 , v0=0.0 , ξ=0.1 , ft=0.0 , and ω=2π are used. Δt=0.1 and 0.2 are used for both cases. All properties are dimensionless. As shown in Figs.4- 7, the new third- and fourth-order methods are presenting more accurate solutions when compared with the RK methods. The third-order RK method gives very inaccurate predictions when Δt=0.2 is used as shown in Figs.4 and 6. The reference solution is obtained by using the exact solution.

Figure 4: Comparison of displacements the single degree of freedom system u¨t+2ξωu˙t+ω2ut=qt . u0=1.0 , v0=0.0 , ξ=0.0 , ft=0.0 , ω=2π and Δt=0.2 are used. 

Figure 5: Comparison of displacements the single degree of freedom system u¨t+2ξωu˙t+ω2ut=qt . u0=1.0 , v0=0.0 , ξ=0.0 , ft=0.0 , ω=2π and Δt=0.1 are used. 

Figure 6: Comparison of displacements the single degree of freedom system u¨t+2ξωu˙t+ω2ut=qt . u0=1.0 , v0=0.0 , ξ=0.1 , ft=0.0 , ω=2π and Δt=0.2 are used. 

Figure 7: Comparison of displacements the single degree of freedom system u¨t+2ξωu˙t+ω2ut=qt . u0=1.0 , v0=0.0 , ξ=0.1 , ft=0.0 , ω=2π and Δt=0.1 are used. 

5.2 Elastic hardening spring problems

The hardening elastic spring problem is considered. The motion of a hardening elastic spring [31-33] is expressed as

u¨+s11+s2u2u=0, (73a)

u0=u0, (73b)

u˙0=v0, (73c)

where s1>0 and s2>0 . The elastic spring problem given in Eq.(73) is a conservative system. The total energy of the spring can be computed exactly by using the initial displacement and velocity. For the numerical test, u0=1.5 , v0=0.0 , s1=100.0 , and s2=10.0 are used [32]. All properties are dimensionless. The nonlinear period T of the problem is computed by using 10th-order implicit Kim and Reddy method [7,20] with time step t=0.000001 . The nonlinear period of the problem is computed as T=0.15153 . All methods used t=T/32 .

Figure 8: Phase portrait of the hardening spring problem u¨+s11+s2u2u=0 , u0=1.5 , v0=0.0 , s1=100.0 , and s2=10.0 are used with t=T/32 to complete fifty cycles (i.e., until t=50T ). (a) New fourth-order method; (b) New third-order method; (c) Fourth-order RK method; (d) Third-order RK method. 

Fifty cycles of displacement and velocity solutions are observed. Since this problem is conservative, the displacement-velocity portrait should form a closed cycle as depicted in Fig.8(a) and (b). The results presented in Fig.8 show that the phase portrait obtained by using the new methods is very accurate. Numerical solutions are forming thirty-two groups on the exact portrait when t=T/32 is used. However, the phase portraits obtained from the RK methods are scatted along the exact phase portrait as shown Fig.8(c) and (d).

Figure 9: Displacement solution of the hardening spring problem u¨+s11+s2u2u=0 , u0=1.5 , v0=0.0 , s1=100.0 , and s2=10.0 are used with t=T/32 . (a) short-term solutions ( 0t2T ); (b) long-term solutions ( 48Tt50T ). 

The results presented in Fig.9 support the fact that the small period error of the new methods can improve the quality of solutions in long-term analyses.

5.3 Elastic softening spring problems

The softening equation [31-33] is also solved by using the new methods. The motion of the softening elastic spring is expressed as

u¨+stanhu=0, (74a)

u0=u0, (74b)

u˙0=v0, (74c)

where s>0.0 . Here we use s=100.0 , u0=4.0 , and v0=0.0 [32]. All properties are dimensionless. Fifty cycles of displacement and velocity solutions are observed. This problem is also conservative, and exact displacement-velocity portrait should form a closed cycle as depicted in Fig.10(a), (b) and (c). The results presented in Fig.10 show that the phase portrait obtained by using the new methods is accurately forming 32 groups on the exact portrait of the displacement-velocity when Δt=T/32 is used. The phase portrait obtained from the fourth-order RK method is less accurate than the new methods. On the other hand, the phase portrait obtained from the third-order RK method is approaching toward the origin. The third-order RK method exhibits the most inaccurate result as shown in Fig.10(d). As shown in Fig.10(d), the third-order RK method becomes excessively dissipative when a large time step is used.

The results presented in Fig.11 also support the fact that the small period and damping errors of the new methods can improve the quality of numerical solutions in long-term analyses.

Figure 10: Phase portrait of the hardening spring problem u¨+stanhu=0 , s=100.0 , u0=4.0 , and v0=0.0 are used with t=T/32 to complete fifty cycles (i.e., until t=50T ). (a) New fourth-order method; (b) New third-order method; (c) Fourth-order RK method; (d) Third-order RK method. 

Figure 11: Displacement solution of the hardening spring problem u¨+stanhu=0 , s=100.0 , u0=4.0 , and v0=0.0 are used with t=T/32 . (a) short-term solutions ( 0t2T ); (b) long-term solutions ( 48Tt50T ). 

5.4 Nonlinear single pendulum problem

Figure 12: Description of nonlinear single pendulum [5

The nonlinear oscillation of single pendulum described in Fig.12 is numerically solved. The motion of the single pendulum is described

θ¨+ω2sinθ=0, (75)

with initial conditions

θ0=θ0, (76b)

θ˙0=θ˙0, (76c)

where ω is g/L , θt is the angle between the rigid rod and the vertical line at time t , g is the gravitational constant, and L is the length of the massless rigid rod. Dimensionless values g=1.0 and L=1.0 are used in the example.

This simple nonlinear pendulum problem is suitable for the test of time integration methods, because important information regarding motions of the pendulum (such as the period and the maximum angle) can be exactly computed by using the initial conditions given in Eqs.(76a) and (76b). For details, please see Refs.[5], [34], and [35].

For the test of the new methods, two sets of initial conditions are used. First, the initial conditions that have been used in Refs.[5,7,14] are also used in this study to synthesize a highly nonlinear case where the pendulum oscillates continuously between two peak points not making complete turns. Second, the initial conditions are chosen as θ0=0.0 and θ˙0=2.000000761543501 to synthesize a highly nonlinear case where pendulum makes complete turns about the point O . These initial conditions can synthesize even more challenging problem than the case used in Ref.[14]. Mathematically computed period of the first set of the initial conditions is 33.72102056501721, and numerically computed period of the second set of the initial conditions is 33.7210. Interestingly, periods of the two cases are almost identical. t=T/400 is used for each set of initial conditions.

Figure 13: Angles of the oscillating nonlinear single pendulum problem described in Fig.12(a) for varying values of t . The new third-order explicit method and the third-order Runge-Kutta method are used, and t is selected as T/400 for both methods, where T=33.72102056485366

Figure 14: Angles of the oscillating nonlinear single pendulum problem described in Fig.12(a) for varying values of t . The new fourth-order explicit method and the third-order Runge-Kutta method are used, and t is selected as T/400 for both methods, where T=33.72102056485366

Figure 15: Angles of the rotating nonlinear single pendulum problem described in Fig.12(b) for varying values of t . The new third-order explicit method and the third-order Runge-Kutta method are used, and t is selected as T/400 for both methods, where T=33.7210

Figure 16: Angles of the rotating nonlinear single pendulum problem described in Fig.12(b) for varying values of t . The new fourth-order explicit method and the third-order Runge-Kutta method are used, and t is selected as T/400 for both methods, where T=33.7210

Figure 17: Comparison of relative error θnumerical-θexact/θexact at t=T/4 . θexact at t=T/4 is analytically computed as 3.13984732433779890888572. 

For the first set of the initial conditions, the numerical solution of the new third-order method is almost superposing the exact solution (the dotted red line) as shown in Fig.13. However, the third-order RK method is presenting noticeable period error as shown in Fig.13. The numerical solutions of the new fourth-order methods are completely superposing the exact solution (the dotted red line) as shown in Fig.14. On the other hand, the fourth-order RK method is showing better results when compared with the third-order RK method, but the period error is still noticeable as shown in Fig.14.

For the second set of the initial conditions, the numerical solution of the new third-order method is almost superposing the exact solution (the dotted red line) as shown in Fig.15. However, the third-order RK method is presenting huge period error and giving a completely incorrect prediction as shown in Fig.15. The result of the third-order RK method indicates that the pendulum is oscillating between the two peak points instead of making complete turns. The numerical solution of the new fourth-order method is completely superposing the exact solution (the dotted red line) as shown in Fig.16, while the fourth-order RK method is showing noticeable period error.

To investigate convergence rate of the nonlinear solutions obtained from various methods, the first set of the initial conditions is used. Errors θnumerical-θexact/θexact at t=T/4 are computed and plotted as shown in Fig.17. The order of accuracies expected in the linear single-degree-of-freedom problem as given in Eqs.(70) and (71) are also achieved in this nonlinear problem.

5.5 Spring-pendulum problem

The two-degree-of-freedom spring-pendulum problem [10] is used for the test of the new methods. The configuration of the two-degree-of-freedom spring-pendulum problem is described in Fig.18. The governing equations are given by

mr¨-mL0+rθ˙2-mgcosθ+kr=0, (76b)

mθ¨+m2r˙θ˙+gsinθL0+r=0, (76c)

and the initial conditions are

r0=r0,θ0=θ0,r˙0=r˙0,θ˙0=θ˙0 (76b)

where r and θ are the displacements in the radial and circumferential directions, respectively, m is the mass of the pendulum, g is the gravitational constant, L0 is the length of the undeformed spring, and k is the spring constant. In this numerical test, m=1.0 , g=9.81 , L=0.5 , k=98.1 , r0=0.25 , θ0=π/2=1.570796326794896 , r˙0=0.0 and θ˙0=0.0 are used. All properties are dimensionless.

Figure 18: Description of spring pendulum [5]. 

Figure 19: Comparison of angles of spring-pendulum. 0.1 is used for t in all methods. 

Figure 20: Comparison of angles of spring-pendulum. 0.05 is used for t in all methods. 

Figure 21: Comparison of relative error rnumerical-rreference/rreference at t=0.1 . rreference is computed by using the 10th-order accurate method of Kim and Reddy [7,20]. 

In this particular problem, the new fourth-order method is presenting accurate solutions while the new third-order method and the third- and fourth order RK methods are not as presented in Figs.19-20. Interestingly, the convergence rate of the new third-order explicit method for the spring-pendulum is third-order as shown in Fig.21. This is due to the velocity dependent nonlinearity included in the governing equations. The results presented in Fig.21 are in good agreement with the order of accuracy computed by using the linear single-degree-of-freedom problem given in Eqs.(70) and (71).

5.6 Double pendulum problem

A double pendulum problem described in Fig.22 is solve to test the new methods. This double pendulum problem can also be used to investigate chaos dynamics [36].

Figure 22: Description of double pendulum. 

Figure 23: Comparison of x1 for 90t95

Figure 24: Comparison of y1 for 90t95

Figure 25: Comparison of x2 for 90t95

Figure 26: Comparison of y2 for 90t95

The governing equations of the double pendulum [37] is given by

θ¨1-m2l2θ˙22sinθ2-θ1-m1+m2gsinθ1l1m1+m2-m2l1cos2θ2-θ1=0, (79a)

θ¨2+m1+m2l1θ˙12sinθ2-θ1+m1+m2gsinθ2l2m1+m2-m2l2cos2θ2-θ1=0, (79b)

θ10,θ20,θ˙10,θ˙20T=ϑ1,ϑ2,ϑ˙1,ϑ˙2T, (79c)

where θ1 and θ2 are angle of the first and second pendulums, respectively, m1 and m2 are masses of the first and second pendulums, respectively, l1 and l2 are weightless rigid rods connected to pendulums as shown in Fig.22, and ϑ1 and ϑ2 are the initial angles. By using θ1 and θ2 , positions of the pendulum are computed as follows:

x1=l1sinθ1, (80a)

y1=-l1cosθ1, (80b)

x2=x1+l2sinθ2, (80c)

y2=y1-l2cosθ2, (80d)

Dimensionless properties m1=1.0 , m2=1.0 , l1=1.0 , l2=1.0 , and g=1.0 are used, and the initial conditions are chosen as ϑ1,ϑ2,ϑ˙1,ϑ˙2T=0.0,π/2,0.0,0.0T . All properties are dimensionless. In this particular problem, the new fourth-order method and the fourth-order RK method are presenting accurate solutions while the new third-order method and the third-order RK method are not as presented in Figs.23- 26. The new fourth-order method is presenting slightly more accurate result when compared with the fourth-order RK method. The convergence rates of the new explicit methods for the double-pendulum problem are presented in Fig.27.

Figure 27: Comparison of relative error y2numerical-y2reference/y2reference at t=1.0 . y2reference is computed by using the 10th-order accurate method of Kim and Reddy [7,20]. 

6 Concluding remarks

The new third- and fourth-order explicit methods presented in this paper could provide accurate numerical solutions when applied to various linear and nonlinear test problems. Simple but meaningful examples were selected from the existing studies where various methods have been developed and tested. The selected illustrative examples were also used to verify improved performances of the new methods. Due to the computational structures of the new methods which are similar to those of the equivalent RK methods, various nonlinear problems of structural dynamics were successfully solved. The new explicit methods could provide more accurate numerical solutions when compared with the well-known third- and fourth-order RK methods.

The advantages of the new explicit methods can be summarized as:

  • (a) The new methods could be applied to linear and nonlinear problems in a consistent and unified manner.

  • (b) The final forms of the new methods did not have any undetermined parameters, which reduced the efforts of a user.

  • (c) For velocity independent nonlinear problems and undamped linear problems, the new third-order method could provide improved solutions with less computational efforts when compared with the fourth-order accurate RK method.

Acknowledgments

The author truly appreciates the support and love from Donghee Son. Thank you.

References

[1] J. N. Reddy. An Introduction to the Finite Element Method. McGraw-Hill New York, 3rd edition, 2006. [ Links ]

[2] O. C. Zienkiewicz and R. L. Taylor. The Finite Element Method for Solid and Structural Mechanics. Butterworth-Heinemann, 2005. [ Links ]

[3] D. Schillinger, J. A. Evans, F. Frischmann, R. R. Hiemstra, M. Hsu, and T. J. R. Hughes. A collocated c0 finite element method: Reduced quadrature perspective, cost comparison with standard finite elements, and explicit structural dynamics. International Journal for Numerical Methods in Engineering, 102(3-4):576-631, 2015. [ Links ]

[4] W. Kim and J. N. Reddy. Novel mixed finite element models for nonlinear analysis of plates. Latin American Journal of Solids and Structures, 7(2):201-226, 2010. [ Links ]

[5] W. Kim and S. Y. Choi. An improved implicit time integration algorithm: The generalized composite time integration algorithm. Computers and Structures, 196:341-354, 2018. [ Links ]

[6] W. Kim and J. N. Reddy. An improved time integration algorithm: A collocation time finite element approach. International Journal of Structural Stability and Dynamics, 17(02):1750024, 2017. [ Links ]

[7] W. Kim and J. N. Reddy. A new family of higher-order time integration algorithms for the analysis of structural dynamics. Journal of Applied Mechanics, 84(7):071008, 2017. [ Links ]

[8] W. Kim and J. N. Reddy. Effective higher-order time integration algorithms for the analysis of linear structural dynamics. Journal of Applied Mechanics, 84(7):071009, 2017. [ Links ]

[9] W. Kim, S. Park, and J. N. Reddy. A cross weighted-residual time integration scheme for structural dynamics. International Journal of Structural Stability and Dynamics, 14(06):1450023, 2014. [ Links ]

[10] J. Chung and J. M. Lee. A new family of explicit time integration methods for linear and non-linear structural dynamics. International Journal for Numerical Methods in Engineering, 37(23):3961-3976, 1994. [ Links ]

[11] G. M. Hulbert and J. Chung. Explicit time integration algorithms for structural dynamics with optimal numerical dissipation. Computer Methods in Applied Mechanics and Engineering, 137(2):175-188, 1996. [ Links ]

[12] D. Soares. An explicit family of time marching procedures with adaptive dissipation control. International Journal for Numerical Methods in Engineering, 100(3):165-182, 2014. [ Links ]

[13] D. Soares. A novel family of explicit time marching techniques for structural dynamics and wave propagation models. Computer Methods in Applied Mechanics and Engineering, 311:838-855, 2016. [ Links ]

[14] W. Kim and J. H. Lee. An improved explicit time integration method for linear and nonlinear structural dynamics. Computers and Structures, 206:42-53, 2018. [ Links ]

[15] W. Kim. A simple explicit single step time integration algorithm for structural dynamics. International Journal for Numerical Methods in Engineering, pages 00-00, 2019. [ Links ]

[16] W. Kim. An accurate two-stage explicit time integration scheme for structural dynamics and various dynamic problems. International Journal for Numerical Methods in Engineering, pages 00-00, 2019. [ Links ]

[17] M. A. Dokainish and K. Subbaraj. A survey of direct time-integration methods in computational structural dynamics i. explicit methods. Computers and Structures, 32(6):1371-1386, 1989. [ Links ]

[18] T. J. R. Hughes. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover, 2012. [ Links ]

[19] K. J. Bathe. Finite Element Procedures. Klaus-Jurgen Bathe, 2006. [ Links ]

[20] W. Kim and J. Lee. A comparative study of two families of higher-order accurate time integration algorithms. International Journal of Computational Mechanics, accepted. [ Links ]

[21] Gregory M Hulbert. A unified set of single-step asymptotic annihilation algorithms for structural dynamics. Computer Methods in Applied Mechanics and Engineering, 113(1-2):1-9, 1994. [ Links ]

[22] T.C. Fung. Unconditionally stable higher-order accurate hermitian time finite elements. International journal for numerical methods in engineering, 39(20):3475-3495, 1996. [ Links ]

[23] T.C. Fung. Weighting parameters for unconditionally stable higher-order accurate time step integration algorithms. part 2second-order equations. International journal for numerical methods in engineering, 45(8):971-1006, 1999. [ Links ]

[24] T.C Fung. Unconditionally stable higher-order accurate collocation time-step integration algorithms for first-order equations. Computer methods in applied mechanics and engineering, 190(13-14):1651-1662, 2000. [ Links ]

[25] J. A. Evans, R. R. Hiemstra, T. J. R. Hughes, and A. Reali. Explicit higher-order accurate isogeometric collocation methods for structural dynamics. Computer Methods in Applied Mechanics and Engineering, 338:208-240, 2018. [ Links ]

[26] H. Adeli, J. M. Gere, and W. Weaver. Algorithms for nonlinear structural dynamics. Journal of the Structural Division, 104(2):263-280, 1978. [ Links ]

[27] T. J. R. Hughes. The finite element method: linear static and dynamic finite element analysis. Courier Corporation, 2012. [ Links ]

[28] T. J. R. Hughes. Analysis of transient algorithms with particular reference to stability behavior. Computational methods for transient analysis(A 84-29160 12-64). Amsterdam, North-Holland, 1983, pages 67-155, 1983. [ Links ]

[29] H. M. Hilber. Analysis and design of numerical integration methods in structural dynamics. PhD thesis, University of California Berkeley, 1976. [ Links ]

[30] G. Noh and K. J. Bathe. An explicit time integration scheme for the analysis of wave propagations. Computers and structures, 129:178-193, 2013. [ Links ]

[31] W. L. Wood and M. E. Oduor. Stability properties of some algorithms for the solution of nonlinear dynamic vibration equations. International Journal for Numerical Methods in Biomedical Engineering, 4(2):205-212, 1988. [ Links ]

[32] Y. M. Xie. An assessment of time integration schemes for non-linear dynamic equations. Journal of Sound and Vibration, 192(1):321-331, 1996. [ Links ]

[33] J. Liu and X. Wang. An assessment of the differential quadrature time integration scheme for nonlinear dynamic equations. Journal of Sound and Vibration, 314(1):246-253, 2008. [ Links ]

[34] T. C. Fung. Solving initial value problems by differential quadrature method-part 2: second and higher-order equations. International Journal for Numerical Methods in Engineering, 50(6):1429-1454, 2001. [ Links ]

[35] T. C. Fung. On the equivalence of the time domain differential quadrature method and the dissipative runge-kutta collocation method. International Journal for Numerical Methods in Engineering, 53(2):409-431, 2002. [ Links ]

[36] S. Maiti, J. Roy, A. K. Mallik, and J. K. Bhattacharjee. Nonlinear dynamics of a rotating double pendulum. Physics Letters A, 380(3):408-412, 2016. [ Links ]

[37] T. Stachowiak and T. Okada. A numerical analysis of chaos in the double pendulum. Chaos, Solitons & Fractals, 29(2):417-422, 2006. [ Links ]

Received: April 29, 2019; Revised: May 18, 2019; Accepted: May 22, 2019; Published: June 03, 2019

* Corresponding author

Editor:

Marco L. Bittencourt

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License