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
^{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

and the initial conditions are given by

where

where

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

where

The displacement and velocity vectors at

where

Finally, the displacement and velocity vectors at

where

To achieve third-order accuracy, extended stability, and preferable spectral properties, permissible

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

where

By using

By using Eq.(17), the displacement and velocity vectors at

By using

Finally, the displacement and velocity vectors at

As shown in Eqs.(14)-(22),
^{3}] can produce diagonal mass matrix automatically if low order spatial elements are considered. After computing

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

where

where

To achieve fourth-order accuracy, extended stability, and preferable spectral properties, permissible

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

where

By using

By using Eq.(32), the displacement and velocity vectors at

By using

The interim displacement and velocity vectors at

By using

Finally, the displacement and velocity vectors at

After computing

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,

The displacement and velocity vectors at

The interim displacement and velocity vectors at

Finally, the displacement and velocity vectors at

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

3.2 Fourth-order Runge-Kutta Method

In the fourth-order RK method,

The displacement and velocity vectors at

The displacement and velocity vectors at

The interim displacement and velocity vectors at

Finally, the displacement and velocity vectors at

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

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

with the initial conditions

where

where

where

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

where
^{28},^{29}] is defined as

where

where

The new third-order method gives

According to Eq.(69), the new three-stage method is third-order accurate for the damped case (

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

where
^{10},^{11},^{30}].

4.3 Period and damping error

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
^{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.

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,

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

where
^{32}]. All properties are dimensionless. The nonlinear period
^{7},^{20}] with time step

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

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

where
*. Here we use*
*,*
*, and*
^{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
*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.

5.4 Nonlinear single pendulum problem

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

with initial conditions

where

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
^{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.

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

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

and the initial conditions are

where

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}].

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

where

Dimensionless properties

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.