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

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 [¹-⁴] 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 [⁵-⁹] and explicit [¹⁰-¹⁶] 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 [¹⁰,¹⁷], 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 [¹⁸,¹⁹]. 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 [¹⁵]. The CD method cannot satisfy (a) and (b).

Chung and Lee developed a family of second-order accurate explicit method [¹⁰] 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 [¹¹]. 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 [¹³]. 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 [¹⁴]. 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 [¹⁶] to more effectively tackle challenging nonlinear problems.

For more precise and reliable analyses of challenging dynamic problems [²⁰], 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 [⁷,²¹-²⁴]. 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 [⁷,⁸].

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 [²⁵].

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 $∆t^4$ 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 $∆t^2$ 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 [²⁶], 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

$\mathbf{M}\ddot{\mathbf{u}}\left(\mathrm{t}\right)=\mathbf{f}\left(\mathbf{u}\left(\mathrm{t}\right),\dot{\mathbf{u}}\left(\mathrm{t}\right),\mathrm{t}\right),$ (1a)

and the initial conditions are given by

$\mathbf{u}\left(0\right)={\mathbf{u}}_0,$ (1b)

$\dot{\mathbf{u}}\left(0\right)={\dot{\mathbf{u}}}_0,$ (1c)

where $\mathbf{M}$ is the mass matrix, $\mathbf{u}\left(\mathrm{t}\right)$ is the displacement vector, $\dot{\mathbf{u}}\left(\mathrm{t}\right)$ is the velocity vector, $\ddot{\mathbf{u}}\left(\mathrm{t}\right)$ is the acceleration vector, ${\mathbf{u}}_0$ is the initial displacement vector, ${\dot{\mathbf{u}}}_0$ is the initial velocity vector, and $\mathbf{f}\left(\mathbf{u}\left(\mathrm{t}\right),\dot{\mathbf{u}}\left(\mathrm{t}\right),\mathrm{t}\right)$ is the total force vector. Eq.(1) is suitable for the description of various linear and nonlinear dynamic problems. In linear structural dynamics, $\mathbf{f}\left(\mathbf{u}\left(\mathrm{t}\right),\dot{\mathbf{u}}\left(\mathrm{t}\right),\mathrm{t}\right)$ in Eq.(1) is frequently expressed as

$\mathbf{f}\left(\mathbf{u}\left(\mathrm{t}\right),\dot{\mathbf{u}}\left(\mathrm{t}\right),\mathrm{t}\right)=\mathbf{q}\left(\mathrm{t}\right)-\mathbf{C}\dot{\mathbf{u}}\left(\mathrm{t}\right)-\mathbf{K}\mathbf{u}\left(\mathrm{t}\right),$ (2)

where $\mathbf{C}$ and $\mathbf{K}$ are the damping and stiffness matrices, respectively, and $\mathbf{q}\left(\mathrm{t}\right)$ 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

$\mathbf{M}\ddot{\mathbf{u}}\left(\mathrm{t}\right)+\mathbf{C}\dot{\mathbf{u}}\left(\mathrm{t}\right)+\mathbf{K}\mathbf{u}\left(\mathrm{t}\right)=\mathbf{q}\left(\mathrm{t}\right).$ (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 $0\le \stackrel{-}{t}\le ∆t$, where $\stackrel{-}{t}=t-t_s$ is the local time, and $t_s$ is the beginning of the $\left(s+1\right)$-th time step or the end of the $s$-th time step. In the new third-order method, the displacement and velocity vectors at $t_s+{\tau }_1∆t$ are approximate as

$\mathbf{u}\left(t_s+{\tau }_1∆t\right)\approx {\stackrel{-}{\mathbf{u}}}_{t_s+{\tau }_1∆t}={}_{}{}^1{\mathrm{\varphi }}_1\left({\tau }_1∆t\right){\mathbf{u}}_{t_s}+{}_{}{}^1{\mathrm{\varphi }}_2\left({\tau }_1∆t\right){\dot{\mathbf{u}}}_{t_s}+\frac{1}{2}{\left({\tau }_1∆t\right)}^2{\mathrm{\alpha }}_{11}{\ddot{\mathbf{u}}}_{t_s},$ (4)

$\dot{\mathbf{u}}\left(t_s+{\tau }_1∆t\right)\approx {\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+{\tau }_1∆t}={}_{}{}^1{\mathrm{\varphi }}_1\left({\tau }_1∆t\right){\dot{\mathbf{u}}}_{t_s}+{}_{}{}^1{\mathrm{\varphi }}_2\left({\tau }_1∆t\right){\ddot{\mathbf{u}}}_{t_s},$ (5)

where ${\tau }_1$ determines the time point that the first approximated variables are associated with, ${\ddot{\mathbf{u}}}_{t_s}$ is computed explicitly by using Eq.(1a) at $t_s$, ${}_{}{}^1{\mathrm{\varphi }}_{\mathrm{i}}\left(\stackrel{-}{t}\right)$s are the known linear functions, and ${\mathrm{\alpha }}_{11}$ is the weighting parameter of the acceleration vector in the first approximation of the displacement vector. ${}_{}{}^1{\mathrm{\varphi }}_{\mathrm{i}}\left(\stackrel{-}{t}\right)$s used in Eqs.(4) and (5) are given by

${}_{}{}^1{\mathrm{\varphi }}_1\left(\stackrel{-}{t}\right)=1,{}_{}{}^1{\mathrm{\varphi }}_2\left(\stackrel{-}{t}\right)=\stackrel{-}{t}.$ (6)

The displacement and velocity vectors at $t_s+{\tau }_2∆t$ are approximate as

$\mathbf{u}\left(t_s+{\tau }_2∆t\right)\approx {\stackrel{-}{\mathbf{u}}}_{t_s+{\tau }_2∆t}={}_{}{}^2{\mathrm{\varphi }}_1\left({\tau }_2∆t\right){\mathbf{u}}_{t_s}+{}_{}{}^2{\mathrm{\varphi }}_2\left({\tau }_2∆t\right){\dot{\mathbf{u}}}_{t_s}+{}_{}{}^2{\mathrm{\varphi }}_3\left({\tau }_2∆t\right){\dot{\mathbf{u}}}_{t_s+{\tau }_1∆t}+\frac{1}{2}{\left({\tau }_2∆t\right)}^2\left({\mathrm{\alpha }}_{21}{\ddot{\mathbf{u}}}_{t_s}+{\mathrm{\alpha }}_{22}{\ddot{\mathbf{u}}}_{t_s+{\tau }_1∆t}\right),$ (7)

$\dot{\mathbf{u}}\left(t_s+{\tau }_2∆t\right)\approx {\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+{\tau }_2∆t}={}_{}{}^2{\mathrm{\varphi }}_1\left({\tau }_2∆t\right){\dot{\mathbf{u}}}_{t_s}+{}_{}{}^2{\mathrm{\varphi }}_2\left({\tau }_2∆t\right){\ddot{\mathbf{u}}}_{t_s}+{}_{}{}^2{\mathrm{\varphi }}_3\left({\tau }_2∆t\right){\dot{\mathbf{u}}}_{t_s+{\tau }_1∆t},$ (8)

where ${\tau }_2$ determines the time point that the first approximated variables are associated with, ${\ddot{\mathbf{u}}}_{t_s+{\tau }_1∆t}$ is computed explicitly by using Eq.(1a) at $t_s+{\tau }_1∆t$, ${}_{}{}^2{\mathrm{\varphi }}_{\mathrm{i}}\left(\stackrel{-}{t}\right)$s are the known quadratic interpolation functions, and ${\mathrm{\alpha }}_{21}$ and ${\mathrm{\alpha }}_{22}$ are the weighting parameters of the acceleration vectors in the second approximation of the displacement vector. ${}_{}{}^2{\mathrm{\varphi }}_{\mathrm{i}}\left(\stackrel{-}{t}\right)$s used in Eqs.(7) and (8) are given by

${}_{}{}^2{\mathrm{\varphi }}_1\left(\stackrel{-}{t}\right)=1,{}_{}{}^2{\mathrm{\varphi }}_2\left(\stackrel{-}{t}\right)=-\frac{\stackrel{-}{t}\left(\stackrel{-}{t}-2{\tau }_1∆t\right)}{2{\tau }_1∆t},{}_{}{}^2{\mathrm{\varphi }}_3\left(\stackrel{-}{t}\right)=\frac{{\stackrel{-}{t}}^2}{2{\tau }_1∆t}.$ (9)

Finally, the displacement and velocity vectors at $t_s+∆t$ are updated as

$\mathbf{u}\left(t_s+∆t\right)\approx {\mathbf{u}}_{t_s+∆t}={}_{}{}^3{\mathrm{\varphi }}_1\left(∆t\right){\mathbf{u}}_{t_s}+{}_{}{}^3{\mathrm{\varphi }}_2\left(∆t\right){\dot{\mathbf{u}}}_{t_s}+{}_{}{}^3{\mathrm{\varphi }}_3\left(∆t\right){\dot{\mathbf{u}}}_{t_s+{\tau }_1∆t}+{}_{}{}^3{\mathrm{\varphi }}_4\left(∆t\right){\dot{\mathbf{u}}}_{t_s+{\tau }_2∆t}+\frac{1}{2}{\left(∆t\right)}^2\left({\mathrm{\alpha }}_{31}{\ddot{\mathbf{u}}}_{t_s}+{\mathrm{\alpha }}_{32}{\ddot{\mathbf{u}}}_{t_s+{\tau }_1∆t}+{\mathrm{\alpha }}_{33}{\ddot{\mathbf{u}}}_{t_s+{\tau }_2∆t}\right),$ (10)

$\dot{\mathbf{u}}\left(t_s+∆t\right)\approx {\dot{\mathbf{u}}}_{t_s+∆t}={}_{}{}^3{\mathrm{\varphi }}_1\left(∆t\right){\dot{\mathbf{u}}}_{t_s}+{}_{}{}^3{\mathrm{\varphi }}_2\left(∆t\right){\ddot{\mathbf{u}}}_{t_s}+{}_{}{}^3{\mathrm{\varphi }}_3\left(∆t\right){\ddot{\mathbf{u}}}_{t_s+{\tau }_1∆t}+{}_{}{}^3{\mathrm{\varphi }}_4\left(∆t\right){\ddot{\mathbf{u}}}_{t_s+{\tau }_2∆t},$ (11)

where ${\ddot{\mathbf{u}}}_{t_s+{\tau }_2∆t}$ is computed explicitly by using Eq.(1a) at $t_s+{\tau }_2∆t$, ${}_{}{}^3{\mathrm{\varphi }}_{\mathrm{i}}\left(\stackrel{-}{t}\right)$s are the known cubic interpolation functions, and ${\mathrm{\alpha }}_{31}$, ${\mathrm{\alpha }}_{32}$, and ${\mathrm{\alpha }}_{33}$ are the weighting parameters of the acceleration vectors in the third approximation of the displacement vector. ${}_{}{}^3{\mathrm{\varphi }}_{\mathrm{i}}\left(\stackrel{-}{t}\right)$s used in Eqs.(10) and (11) are given by

$\begin{array}{l}{}_{}{}^3\varphi {}_1\left(\overline{t}\right)=\mathrm{1,}{}_{}{}^3\varphi {}_2\left(\overline{t}\right)=\frac{\overline{t}\left(2{\overline{t}}^2-3{\tau }_1\Delta t\overline{t}-3{\tau }_2\Delta t\overline{t}+6{\tau }_1{\tau }_2\Delta t^2\right)}{6{\tau }_1{\tau }_2\Delta t^2},\\ {}_{}{}^3\varphi {}_3\left(\overline{t}\right)=-\frac{{\overline{t}}^2\left(-2\overline{t}+3{\tau }_2\Delta t\right)}{6{\tau }_1\left(-{\tau }_2+{\tau }_1\right)\Delta t^2},{}_{}{}^3\varphi {}_4\left(\overline{t}\right)=\frac{{\overline{t}}^2\left(-2\overline{t}+3{\tau }_1\Delta t\right)}{6{\tau }_2\left({\tau }_1-{\tau }_2\right)\Delta t^2}.\end{array}$ (12)

To achieve third-order accuracy, extended stability, and preferable spectral properties, permissible ${\tau }_i$ and ${\mathrm{\alpha }}_{\mathrm{i}\mathrm{j}}$ are chosen as

$\begin{array}{l}{\tau }_1=\frac{1}{3},{\tau }_2=\frac{2}{3},\\ {\text{α}}_{11}=\mathrm{1,}{}_{\text{α}}=-\frac{2}{3},{}_{\text{α}}=\frac{2}{3},\\ {\text{α}}_{31}=\frac{1}{3},{}_{\text{α}}=-\frac{2}{3},{}_{\text{α}}=\frac{1}{3}.\end{array}$ (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 $t_s$ as

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}}={\mathbf{M}}^{-1}\mathbf{f}\left({\mathbf{u}}_{{\mathrm{t}}_{\mathrm{s}}},{\dot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}},{\mathrm{t}}_{\mathrm{s}}\right),$ (14)

where ${\mathbf{u}}_{{\mathrm{t}}_{\mathrm{s}}}$ and ${\dot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}}$ are known properties from $s$-th time step. By using Eq.(14), the displacement and velocity vectors at $t_s+∆t/3$ are approximated as

${\stackrel{-}{\mathbf{u}}}_{t_s+∆t/3}={\mathbf{u}}_{t_s}+\frac{1}{3}∆t{\dot{\mathbf{u}}}_{t_s}+\frac{1}{18}{∆t}^2{\ddot{\mathbf{u}}}_{t_s},$ (15)

${\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+∆t/3}={\dot{\mathbf{u}}}_{t_s}+\frac{1}{3}∆t{\ddot{\mathbf{u}}}_{t_s}.$ (16)

By using ${\stackrel{-}{\mathbf{u}}}_{t_s+∆t/3}$ and ${\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+∆t/3}$ given in Eqs.(15) and (16), the acceleration vector at $t_s+∆t/3$ is computed as

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/3}={\mathbf{M}}^{-1}\mathbf{f}\left({\stackrel{-}{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/3},{\dot{\stackrel{-}{\mathbf{u}}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/3},{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/3\right).$ (17)

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

${\stackrel{-}{\mathbf{u}}}_{t_s+2∆t/3}={\mathbf{u}}_{t_s}+\frac{2}{3}∆t{\dot{\mathbf{u}}}_{t_s}+\frac{1}{27}{∆t}^2\left(2{\ddot{\mathbf{u}}}_{t_s}+4{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆t/3}\right),$ (18)

${\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+2∆t/3}={\dot{\mathbf{u}}}_{t_s}+\frac{2}{3}∆t{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆t/3}.$ (19)

By using ${\stackrel{-}{\mathbf{u}}}_{t_s+2∆t/3}$ and ${\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+2∆t/3}$ given in Eqs.(18) and (19), the acceleration vector at $t_s+2∆t/3$ is computed as

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+2∆\mathrm{t}/3}={\mathbf{M}}^{-1}\mathbf{f}\left({\stackrel{-}{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+2∆\mathrm{t}/3},{\dot{\stackrel{-}{\mathbf{u}}}}_{{\mathrm{t}}_{\mathrm{s}}+2∆\mathrm{t}/3},{\mathrm{t}}_{\mathrm{s}}+2∆\mathrm{t}/3\right).$ (20)

Finally, the displacement and velocity vectors at $t_s+∆t$ are computed as

${\mathbf{u}}_{t_s+∆t}={\mathbf{u}}_{t_s}+∆t{\dot{\mathbf{u}}}_{t_s}+\frac{1}{6}{∆t}^2\left({\ddot{\mathbf{u}}}_{t_s}+{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆t/3}+{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+2∆t/3}\right),$ (21)

${\dot{\mathbf{u}}}_{t_s+∆t}={\dot{\mathbf{u}}}_{t_s}+\frac{1}{4}∆t\left({\ddot{\mathbf{u}}}_{t_s}+3{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+2∆t/3}\right).$ (22)

As shown in Eqs.(14)-(22), ${\mathbf{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 $\mathbf{M}$ is diagonal. Also, some recent numerical integration techniques [³] can produce diagonal mass matrix automatically if low order spatial elements are considered. After computing ${\mathbf{u}}_{t_s+∆t}$ and ${\dot{\mathbf{u}}}_{t_s+∆t}$, ${\mathbf{u}}_{t_s}$ and ${\dot{\mathbf{u}}}_{t_s}$ are updated as ${\mathbf{u}}_{t_s+∆t}$ and ${\dot{\mathbf{u}}}_{t_s+∆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 $t_s+{\tau }_3∆t$ take the forms of

$\mathbf{u}\left(t_s+{\tau }_3∆t\right)\approx {\stackrel{-}{\mathbf{u}}}_{t_s+{\tau }_3∆t}={}_{}{}^3{\mathrm{\varphi }}_1\left({\tau }_3∆t\right){\mathbf{u}}_{t_s}+{}_{}{}^3{\mathrm{\varphi }}_2\left({\tau }_3∆t\right){\dot{\mathbf{u}}}_{t_s}+{}_{}{}^3{\mathrm{\varphi }}_3\left({\tau }_3∆t\right){\dot{\mathbf{u}}}_{t_s+{\tau }_1∆t}+{}_{}{}^3{\mathrm{\varphi }}_4\left({\tau }_3∆t\right){\dot{\mathbf{u}}}_{t_s+{\tau }_2∆t}+\frac{1}{2}{\left({\tau }_3∆t\right)}^2\left({\mathrm{\alpha }}_{31}{\ddot{\mathbf{u}}}_{t_s}+{\mathrm{\alpha }}_{32}{\ddot{\mathbf{u}}}_{t_s+{\tau }_1∆t}+{\mathrm{\alpha }}_{33}{\ddot{\mathbf{u}}}_{t_s+{\tau }_2∆t}\right),$ (23)

$\dot{\mathbf{u}}\left(t_s+{\tau }_3∆t\right)\approx {\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+{\tau }_3∆t}={}_{}{}^3{\mathrm{\varphi }}_1\left({\tau }_3∆t\right){\dot{\mathbf{u}}}_{t_s}+{}_{}{}^3{\mathrm{\varphi }}_2\left({\tau }_3∆t\right){\ddot{\mathbf{u}}}_{t_s}+{}_{}{}^3{\mathrm{\varphi }}_3\left({\tau }_3∆t\right){\ddot{\mathbf{u}}}_{t_s+{\tau }_1∆t}+{}_{}{}^3{\mathrm{\varphi }}_4\left({\tau }_3∆t\right){\ddot{\mathbf{u}}}_{t_s+{\tau }_2∆t},$ (24)

where ${\tau }_3$ determines the time point that the third approximated variables are associated with, and ${\ddot{\mathbf{u}}}_{t_s+{\tau }_3∆t}$ is computed explicitly by using Eqs.(7)-(8) and Eq.(1a) at $t_s+{\tau }_3∆t$. The displacement and velocity vectors at $t_s+∆t$ are updated as

$\mathbf{u}\left(t_s+∆t\right)\approx {\mathbf{u}}_{t_s+∆t}={}_{}{}^4{\mathrm{\varphi }}_1\left(∆t\right){\mathbf{u}}_{t_s}+{}_{}{}^4{\mathrm{\varphi }}_2\left(∆t\right){\dot{\mathbf{u}}}_{t_s}+{}_{}{}^4{\mathrm{\varphi }}_3\left(∆t\right){\dot{\mathbf{u}}}_{t_s+{\tau }_1∆t}+{}_{}{}^4{\mathrm{\varphi }}_4\left(∆t\right){\dot{\mathbf{u}}}_{t_s+{\tau }_2∆t}+{}_{}{}^4{\mathrm{\varphi }}_5\left(∆t\right){\dot{\mathbf{u}}}_{t_s+{\tau }_3∆t}+\frac{1}{2}{\left(∆t\right)}^2\left({\mathrm{\alpha }}_{41}{\ddot{\mathbf{u}}}_{t_s}+{\mathrm{\alpha }}_{42}{\ddot{\mathbf{u}}}_{t_s+{\tau }_1∆t}+{\mathrm{\alpha }}_{43}{\ddot{\mathbf{u}}}_{t_s+{\tau }_2∆t}+{\mathrm{\alpha }}_{44}{\ddot{\mathbf{u}}}_{t_s+{\tau }_3∆t}\right),$ (25)

$\dot{\mathbf{u}}\left(t_s+∆t\right)\approx {\dot{\mathbf{u}}}_{t_s+∆t}={}_{}{}^4{\mathrm{\varphi }}_1\left(∆t\right){\dot{\mathbf{u}}}_{t_s}+{}_{}{}^4{\mathrm{\varphi }}_2\left(∆t\right){\ddot{\mathbf{u}}}_{t_s}+{}_{}{}^4{\mathrm{\varphi }}_3\left(∆t\right){\ddot{\mathbf{u}}}_{t_s+{\tau }_1∆t}+{}_{}{}^4{\mathrm{\varphi }}_4\left(∆t\right){\ddot{\mathbf{u}}}_{t_s+{\tau }_2∆t}+{}_{}{}^4{\mathrm{\varphi }}_5\left(∆t\right){\ddot{\mathbf{u}}}_{t_s+{\tau }_3∆t},$ (26)

where ${\ddot{\mathbf{u}}}_{t_s+{\tau }_3∆t}$ is computed explicitly by using Eqs.(23)-(24) and Eq.(1a) at$t_s+{\tau }_3∆t$, ${}_{}{}^4{\mathrm{\varphi }}_{\mathrm{i}}\left(\stackrel{-}{t}\right)$s are the known quartic interpolation functions, and ${\mathrm{\alpha }}_{41}$, ${\mathrm{\alpha }}_{42}$, ${\mathrm{\alpha }}_{43}$, and ${\mathrm{\alpha }}_{44}$ are the weighting parameters of the acceleration vectors in the fourth approximation of the displacement vector. ${}_{}{}^4{\mathrm{\varphi }}_{\mathrm{i}}\left(\stackrel{-}{t}\right)$s used in Eqs.(25) and (26) are given by

$\begin{array}{l}{}_{}{}^4\varphi {}_1\left(\overline{t}\right)=\mathrm{1,}\\ {}_{}{}^4\varphi {}_2\left(\overline{t}\right)=-\frac{\overline{t}\left(3{\overline{t}}^3-4\Delta t{\overline{t}}^2\left({\tau }_1+{\tau }_2+{\tau }_3\right)+6\Delta t^2\overline{t}\left({\tau }_1{\tau }_2+{\tau }_2{\tau }_3+{\tau }_3{\tau }_1\right)\right)}{12{\tau }_1{\tau }_2{\tau }_3\Delta t^3}+\overline{t},\\ {}_{}{}^4\varphi {}_3\left(\overline{t}\right)=\frac{{\overline{t}}^2\left(3{\overline{t}}^2-4{\tau }_2\Delta t\overline{t}-4{\tau }_3\Delta t\overline{t}+6{\tau }_2{\tau }_3\Delta t^2\right)}{12{\tau }_1\left(-{\tau }_1{\tau }_2+{\tau }_2{\tau }_3-{\tau }_3{\tau }_1+{\tau }_1^2\right)\Delta t^3},\\ {}_{}{}^4\varphi {}_4\left(\overline{t}\right)=\frac{{\overline{t}}^2\left(3{\overline{t}}^2-4{\tau }_1\Delta t\overline{t}-4{\tau }_3\Delta t\overline{t}+6{\tau }_1{\tau }_3\Delta t^2\right)}{12{\tau }_2\left(-{\tau }_1{\tau }_2-{\tau }_2{\tau }_3+{\tau }_3{\tau }_1+{\tau }_2^2\right)\Delta t^3},\\ {}_{}{}^4\varphi {}_5\left(\overline{t}\right)=\frac{{\overline{t}}^2\left(3{\overline{t}}^2-4{\tau }_1\Delta t\overline{t}-4{\tau }_2\Delta t\overline{t}+6{\tau }_1{\tau }_2\Delta t^2\right)}{12{\tau }_3\left({\tau }_1{\tau }_2-{\tau }_2{\tau }_3-{\tau }_3{\tau }_1+{\tau }_3^2\right)\Delta t^3}.\end{array}$ (27)

To achieve fourth-order accuracy, extended stability, and preferable spectral properties, permissible ${\tau }_i$ and ${\mathrm{\alpha }}_{\mathrm{i}\mathrm{j}}$ are chosen as

$\begin{array}{l}{\tau }_1=\frac{1}{3},{\tau }_2=\frac{1}{2},{\tau }_3=\mathrm{1,}\\ {\text{α}}_{11}=\mathrm{1,}{}_{\text{α}}=-\frac{3}{5},{}_{\text{α}}=\frac{3}{5},\\ {\text{α}}_{31}=\frac{3}{5},{}_{\text{α}}=-\frac{3}{5},{}_{\text{α}}=\mathrm{0,}\\ {\text{α}}_{41}={\text{α}}_{42}={\text{α}}_{43}={\text{α}}_{44}=0.\end{array}$ (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 $t_s$ as

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}}={\mathbf{M}}^{-1}\mathbf{f}\left({\mathbf{u}}_{{\mathrm{t}}_{\mathrm{s}}},{\dot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}},{\mathrm{t}}_{\mathrm{s}}\right),$ (29)

where ${\mathbf{u}}_{{\mathrm{t}}_{\mathrm{s}}}$ and ${\dot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}}$ are known properties from $s$-th time step. By using Eq.(14), the displacement and velocity vectors at $t_s+∆t/3$ are approximated as

${\stackrel{-}{\mathbf{u}}}_{t_s+∆t/3}={\mathbf{u}}_{t_s}+\frac{1}{3}∆t{\dot{\mathbf{u}}}_{t_s}+\frac{1}{18}{∆t}^2{\ddot{\mathbf{u}}}_{t_s},$ (30)

${\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+∆t/3}={\dot{\mathbf{u}}}_{t_s}+\frac{1}{3}∆t{\ddot{\mathbf{u}}}_{t_s}.$ (31)

By using ${\stackrel{-}{\mathbf{u}}}_{t_s+∆t/3}$ and ${\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+∆t/3}$ given in Eqs.(30) and (31), the acceleration vector at $t_s+∆t/3$ is computed as

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/3}={\mathbf{M}}^{-1}\mathbf{f}\left({\stackrel{-}{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/3},{\dot{\stackrel{-}{\mathbf{u}}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/3},{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/3\right).$ (32)

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

${\stackrel{-}{\mathbf{u}}}_{t_s+∆t/2}={\mathbf{u}}_{t_s}+\frac{1}{2}∆t{\dot{\mathbf{u}}}_{t_s}+\frac{1}{40}{∆t}^2\left(2{\ddot{\mathbf{u}}}_{t_s}+3{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆t/3}\right),$ (33)

${\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+∆t/2}={\dot{\mathbf{u}}}_{t_s}+\frac{1}{8}∆t\left({\ddot{\mathbf{u}}}_{t_s}+3{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆t/3}\right).$ (34)

By using ${\stackrel{-}{\mathbf{u}}}_{t_s+∆t/2}$ and ${\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+∆t/3}$ given in Eqs.(33) and (34), the acceleration vector at $t_s+∆t/2$ is computed as

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2}={\mathbf{M}}^{-1}\mathbf{f}\left({\stackrel{-}{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2},{\dot{\stackrel{-}{\mathbf{u}}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2},{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2\right).$ (35)

The interim displacement and velocity vectors at ${\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}$ are approximated as

${\stackrel{-}{\mathbf{u}}}_{t_s+∆t}={\mathbf{u}}_{t_s}+∆t{\dot{\mathbf{u}}}_{t_s}+\frac{1}{20}{∆t}^2\left({\ddot{\mathbf{u}}}_{t_s}+9{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆t/3}\right),$ (36)

${\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+∆t}={\dot{\mathbf{u}}}_{t_s}+\frac{1}{2}∆t\left({\ddot{\mathbf{u}}}_{t_s}-3{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆t/3}+4{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆t/2}\right).$ (37)

By using ${\stackrel{-}{\mathbf{u}}}_{t_s+∆t}$ and ${\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+∆t}$ given in Eqs.(36) and (37), the acceleration vector at $t_s+∆t$ is computed as

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}}={\mathbf{M}}^{-1}\mathbf{f}\left({\stackrel{-}{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}},{\dot{\stackrel{-}{\mathbf{u}}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}},{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}\right).$ (38)

Finally, the displacement and velocity vectors at $t_s+∆t$ are computed as

${\mathbf{u}}_{t_s+∆t}={\mathbf{u}}_{t_s}+∆t{\dot{\mathbf{u}}}_{t_s}+\frac{1}{6}{∆t}^2\left({\ddot{\mathbf{u}}}_{t_s}+2{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆t/2}\right),$ (39)

${\dot{\mathbf{u}}}_{t_s+∆t}={\dot{\mathbf{u}}}_{t_s}+\frac{1}{6}∆t\left({\ddot{\mathbf{u}}}_{t_s}+4{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆t/2}+{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}}\right).$ (40)

After computing ${\mathbf{u}}_{t_s+∆t}$ and ${\dot{\mathbf{u}}}_{t_s+∆t}$, ${\mathbf{u}}_{t_s}$ and ${\dot{\mathbf{u}}}_{t_s}$ are updated as ${\mathbf{u}}_{t_s+∆t}$ and ${\dot{\mathbf{u}}}_{t_s+∆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 [²⁵]. 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, ${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}}$ is computed as

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}}={\mathbf{M}}^{-1}\mathbf{f}\left({\mathbf{u}}_{{\mathrm{t}}_{\mathrm{s}}},{\dot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}},{\mathrm{t}}_{\mathrm{s}}\right).$ (41)

The displacement and velocity vectors at $t_s+∆t/2$ are approximated as

${\stackrel{-}{\mathbf{u}}}_{t_s+∆t/2}={\mathbf{u}}_{t_s}+\frac{1}{2}∆t{\dot{\mathbf{u}}}_{t_s},$ (42)

${\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+∆t/2}={\dot{\mathbf{u}}}_{t_s}+\frac{1}{2}∆t{\ddot{\mathbf{u}}}_{t_s}.$ (43)

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2}$ is computed as

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2}={\mathbf{M}}^{-1}\mathbf{f}\left({\stackrel{-}{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2},{\dot{\stackrel{-}{\mathbf{u}}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2},{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2\right).$ (44)

The interim displacement and velocity vectors at $t_s+∆t$ are approximated as

$\begin{array}{l}{\overline{u}}_{t_s+\Delta t}=u_{t_s}+\Delta t\left(-{\dot{u}}_{t_s}+2{\dot{\overline{u}}}_{t_s+\Delta t}\right)\\ =u_{t_s}+\Delta t{\dot{u}}_{t_s}+\Delta t^2{\ddot{u}}_{t_s},\end{array}$ (45)

${\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+∆t}={\dot{\mathbf{u}}}_{t_s}+∆t\left(-{\ddot{\mathbf{u}}}_{t_s}+2{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2}\right).$ (46)

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}}$ is computed as

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}}={\mathbf{M}}^{-1}\mathbf{f}\left({\stackrel{-}{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}},{\dot{\stackrel{-}{\mathbf{u}}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}},{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}\right).$ (47)

Finally, the displacement and velocity vectors at $t_s+∆t$ are computed as

$\begin{array}{l}{u}_{t_s+\Delta t}=u_{t_s}+\frac{1}{6}\Delta t\left({\dot{u}}_{t_s}+4{\dot{\overline{u}}}_{t_s+\Delta t/2}+{\dot{\overline{u}}}_{t_s+\Delta t}\right)\\ =u_{t_s}+\Delta t{\dot{u}}_{t_s}+\frac{1}{6}\Delta t^2\left({\ddot{u}}_{t_s}+2{\ddot{u}}_{t_s+\Delta t/2}\right),\end{array}$ (48)

${\dot{\mathbf{u}}}_{t_s+∆t}={\dot{\mathbf{u}}}_{t_s}+\frac{1}{6}∆t\left({\ddot{\mathbf{u}}}_{t_s}+4{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2}+{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}}\right).$ (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., ${\mathbf{M}}^{-1}\mathbf{f}$). 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, ${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}}$ is computed as

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}}={\mathbf{M}}^{-1}\mathbf{f}\left({\mathbf{u}}_{{\mathrm{t}}_{\mathrm{s}}},{\dot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}},{\mathrm{t}}_{\mathrm{s}}\right).$ (50)

The displacement and velocity vectors at $t_s+∆t/2$ are approximated as

${\stackrel{-}{\mathbf{u}}}_{t_s+∆t/2}={\mathbf{u}}_{t_s}+\frac{1}{2}∆t{\dot{\mathbf{u}}}_{t_s},$ (51)

${\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+∆t/2}={\dot{\mathbf{u}}}_{t_s}+\frac{1}{2}∆t{\ddot{\mathbf{u}}}_{t_s}.$ (52)

${\ddot{\stackrel{-}{\mathbf{u}}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2}$ is computed as

${\ddot{\stackrel{-}{\mathbf{u}}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2}={\mathbf{M}}^{-1}\mathbf{f}\left({\stackrel{-}{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2},{\dot{\stackrel{-}{\mathbf{u}}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2},{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2\right).$ (53)

The displacement and velocity vectors at $t_s+∆t/2$ are approximated as

$\begin{array}{l}{\overline{\overline{u}}}_{t_s+\Delta t/2}=u_{t_s}+\frac{1}{2}\Delta t{\dot{\overline{u}}}_{t_s+\Delta t/2}\\ =u_{t_s}+\frac{1}{2}\Delta t{\dot{u}}_{t_s+\Delta t/2}+\frac{1}{4}\Delta t^2{\ddot{u}}_{t_s},\end{array}$ (54)

${\dot{\stackrel{-}{\stackrel{-}{\mathbf{u}}}}}_{t_s+∆t/2}={\dot{\mathbf{u}}}_{t_s}+\frac{1}{2}∆t{\ddot{\stackrel{-}{\mathbf{u}}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2}.$ (55)

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2}$ is computed as

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2}={\mathbf{M}}^{-1}\mathbf{f}\left({\stackrel{-}{\stackrel{-}{\mathbf{u}}}}_{t_s+∆t/2},{\dot{\stackrel{-}{\stackrel{-}{\mathbf{u}}}}}_{t_s+∆t/2},{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2\right).$ (56)

The interim displacement and velocity vectors at $t_s+∆t$ are approximated as

$\begin{array}{l}{\overline{u}}_{t_s+\Delta t}=u_{t_s}+\Delta t{\dot{\overline{\overline{u}}}}_{t_s+\Delta t/2}\\ =u_{t_s}+\Delta t{\dot{u}}_{t_s}+\frac{1}{2}\Delta t^2{\ddot{\overline{u}}}_{t_s+\Delta t/2},\end{array}$ (57)

${\dot{\stackrel{-}{\mathbf{u}}}}_{t_s+∆t}={\dot{\mathbf{u}}}_{t_s}+∆t{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2}.$ (58)

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}}$ is computed as

${\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}}={\mathbf{M}}^{-1}\mathbf{f}\left({\stackrel{-}{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}},{\dot{\stackrel{-}{\mathbf{u}}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}},{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}\right).$ (59)

Finally, the displacement and velocity vectors at $t_s+∆t$ are computed as

$\begin{array}{l}{u}_{t_s+\Delta t}=u_{t_s}+\frac{1}{6}\Delta t\left({\dot{u}}_{t_s}+2{\dot{\overline{\overline{u}}}}_{t_s+\Delta t/2}+2{\dot{\overline{u}}}_{t_s+\Delta t/2}+{\dot{\overline{u}}}_{t_s+\Delta t}\right)\\ =u_{t_s}+\Delta t{\dot{u}}_{t_s}+\frac{1}{6}\Delta t^2\left({\ddot{u}}_{t_s}+{\ddot{\overline{u}}}_{t_s+\Delta t/2}+{\ddot{u}}_{t_s+\Delta t/2}\right),\end{array}$ (60)

${\dot{\mathbf{u}}}_{t_s+∆t}={\dot{\mathbf{u}}}_{t_s}+\frac{1}{6}∆t\left({\ddot{\mathbf{u}}}_{t_s}+2{\ddot{\stackrel{-}{\mathbf{u}}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2}+2{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}/2}+{\ddot{\mathbf{u}}}_{{\mathrm{t}}_{\mathrm{s}}+∆\mathrm{t}}\right).$ (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., ${\mathbf{M}}^{-1}\mathbf{f}$). 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 [²⁷-²⁹] is frequently used. The linear single-degree-of-freedom problem is given as

$\ddot{\mathrm{u}}\left(\mathrm{t}\right)+2\mathrm{\xi }\mathrm{\omega }\dot{\mathrm{u}}\left(\mathrm{t}\right)+{\mathrm{\omega }}^2\mathrm{u}\left(\mathrm{t}\right)=\mathrm{q}\left(\mathrm{t}\right),$ (62)

with the initial conditions

$\mathrm{u}\left(0\right)={\mathrm{u}}_0,\dot{\mathrm{u}}\left(0\right)={\mathrm{v}}_0.$ (63)

where $\mathrm{u}\left(\mathrm{t}\right)$ is the displacement, ${\mathrm{u}}_0$ is the initial displacement, ${\mathrm{v}}_0$ is the initial velocity, $\mathrm{\xi }$ is the damping ratio, $\mathrm{\omega }$ is the natural frequency, and $\mathrm{q}\left(\mathrm{t}\right)$ is the external force. By setting $\mathbf{M}=1$, $\mathbf{C}=2\mathrm{\xi }\mathrm{\omega }$, $\mathbf{K}={\mathrm{\omega }}^2$, $\mathbf{u}=\mathrm{u}$, and $\mathbf{q}=\mathrm{q}$ and applying Eqs.(14)-(22) to Eqs.(62)-(63), fully discrete relation of the new third-order method at $\mathrm{\Delta }\mathrm{t}$ is given by

${\mathbf{x}}_{\Delta \mathrm{t}}={}_{}{}^{3rd}\mathbf{A}{\mathbf{x}}_0+{}_{}{}^{3rd}\mathbf{L}{}_{}{}^{3rd}\mathbf{q},$ (64)

where ${}_{}{}^{3rd}\mathbf{A}$ is the amplification matrix of the new third-order method, ${}_{}{}^{3rd}\mathbf{L}$ is the load operator matrix of the new third-order method, ${\mathbf{x}}_{\Delta \mathrm{t}}={\left\{{\mathrm{u}}_{\Delta \mathrm{t}},{\mathrm{v}}_{\Delta \mathrm{t}}\right\}}^T$, ${\mathbf{x}}_0={\left\{{\mathrm{u}}_0,{\mathrm{v}}_0\right\}}^T$, and ${}_{}{}^{3rd}\mathbf{q}={\left\{\mathrm{q}\left(0\right),\mathrm{q}\left(\Delta \mathrm{t}/3\right),\mathrm{q}\left(2\Delta \mathrm{t}/3\right)\right\}}^T$. By applying Eqs.(29)-(40) to Eqs.(62)-(63), fully discrete relation of the new fourth-order method at $\Delta \mathrm{t}$ is given by

${\mathbf{x}}_{\Delta \mathrm{t}}={}_{}{}^{4th}\mathbf{A}{\mathbf{x}}_0+{}_{}{}^{4th}\mathbf{L}{}_{}{}^{4th}\mathbf{q},$ (65)

where ${}_{}{}^{4th}\mathbf{A}$ is the amplification matrix of the fourth-order method, ${}_{}{}^{4th}\mathbf{L}$ is the load operator matrix of the fourth-order method, and ${}_{}{}^{4th}\mathbf{q}={\left\{\mathrm{q}\left(0\right),\mathrm{q}\left(\Delta \mathrm{t}/3\right),\mathrm{q}\left(\Delta \mathrm{t}/2\right),\mathrm{q}\left(\Delta \mathrm{t}\right)\right\}}^T$.

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 $\mathbf{A}$ is expressed as

$\mathit{det}\left(\mathbf{A}-\mathrm{\lambda }\mathbf{I}\right)={\mathrm{\lambda }}^2-2{\mathrm{A}}_1\mathrm{\lambda }+{\mathrm{A}}_2,$ (66)

where $\mathrm{\lambda }$ is the eigenvalue of $\mathbf{A}$, and ${\mathrm{A}}_1$ and ${\mathrm{A}}_2$ are the invariants of $\mathbf{A}$. By using ${\mathrm{A}}_1$ and ${\mathrm{A}}_2$, the local truncation error [²⁸,²⁹] is defined as

$\mathrm{\tau }\left({\mathrm{t}}_{\mathrm{s}}\right)=\frac{1}{{\Delta \mathrm{t}}^2}\left(\mathrm{u}\left({\mathrm{t}}_{\mathrm{s}}+\Delta \mathrm{t}\right)-2{\mathrm{A}}_1\mathrm{u}\left({\mathrm{t}}_{\mathrm{s}}\right)+{\mathrm{A}}_2\mathrm{u}\left({\mathrm{t}}_{\mathrm{s}}-\Delta \mathrm{t}\right)\right),$ (67)

where $\mathrm{u}\left(\mathrm{t}\right)$ is the exact solution of Eq.(62). For the case of $\mathrm{q}\left(\mathrm{t}\right)=0$, the exact solution is given by

$\mathrm{u}\left(\mathrm{t}\right)=\mathit{exp}\left(-\mathrm{\xi }\mathrm{\omega }\mathrm{t}\right)\left(\left(\mathit{cos}\left({\mathrm{\omega }}_{\mathrm{d}}\mathrm{t}\right)+\frac{\mathrm{\xi }\mathrm{\omega }}{{\mathrm{\omega }}_{\mathrm{d}}}\mathit{sin}\left({\mathrm{\omega }}_{\mathrm{d}}\mathrm{t}\right)\right){\mathrm{u}}_0+\frac{1}{{\mathrm{\omega }}_{\mathrm{d}}}\mathit{sin}\left({\mathrm{\omega }}_{\mathrm{d}}\mathrm{t}\right){\mathrm{v}}_0\right),$ (68)

where ${\mathrm{\omega }}_{\mathrm{d}}=\sqrt{1-{\mathrm{\xi }}^2}$ is the damped natural frequency. Then a method is $k$-th order accurate if

${\mathrm{\tau }}_{\mathrm{e}}\left({\mathrm{t}}_{\mathrm{s}}\right)=O\left({\Delta \mathrm{t}}^{\mathrm{k}}\right).$ (69)

The new third-order method gives

${}_{}{}^{3rd}{\mathrm{\tau }}_{\mathrm{e}}\left({\mathrm{t}}_{\mathrm{s}}\right)=\frac{1}{108}\mathrm{\xi }{\mathrm{\omega }}^4\left(\left(36{\mathrm{\xi }}^2\mathrm{\omega }-5\mathrm{\omega }\right){\mathrm{u}}_0+\left(72{\mathrm{\xi }}^3-28\mathrm{\xi }\right){\mathrm{v}}_0\right){\Delta \mathrm{t}}^3+O\left({\Delta \mathrm{t}}^4\right).$ (70)

According to Eq.(69), the new three-stage method is third-order accurate for the damped case ($\mathrm{\xi }\ne 0$), and it becomes forth-order accurate for the undamped case ($\mathrm{\xi }=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

${}_{}{}^{4th}{\mathrm{\tau }}_{\mathrm{e}}\left({\mathrm{t}}_{\mathrm{s}}\right)=-\frac{1}{720}{\mathrm{\omega }}^5\left(\left(-52{\mathrm{\xi }}^2\mathrm{\omega }+96{\mathrm{\xi }}^4\mathrm{\omega }+\mathrm{\omega }\right){\mathrm{u}}_0+\left(192{\mathrm{\xi }}^5-152{\mathrm{\xi }}^3+16\mathrm{\xi }\right){\mathrm{v}}_0\right){\Delta \mathrm{t}}^4+O\left({\Delta \mathrm{t}}^5\right),$ (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

$\mathrm{\rho }\left(\mathbf{A}\right)=\mathit{max}\left(\left|{\mathrm{\lambda }}_1\right|,\left|{\mathrm{\lambda }}_2\right|\right),$ (72)

where ${\mathrm{\lambda }}_1$ and ${\mathrm{\lambda }}_2$ are the eigenvalues of $\mathbf{A}$. A method is stable if $0<\mathrm{\rho }\left(\mathbf{A}\right)\le 1$ 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\pi /\omega$ is the natural period. Thus, the new third-order method is stable if $0\le \Delta t<0.574976T$ is provided. The critical time step of the new fourth-order method is $0.474023T$. The stability limit slightly decreases as $\xi$ increases from 0.0 to 1.0 in the new methods. However, similar phenomena have been observed in many well-accepted explicit methods [¹⁰,¹¹,³⁰].

Figure 1:
Comparison of spectral radius $\mathrm{\rho }\left(\mathbf{A}\right)$ for varying values of $\mathrm{\Delta }t/T$.

4.3 Period and damping error

Figure 2:
Comparison of period error $\left(T-\stackrel{-}{T}\right)/T$ for varying $\mathrm{\Delta }t/T$.

Figure 3:
Comparison of damping error $\mathrm{\xi }$ for varying $\mathrm{\Delta }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 $\left(\stackrel{-}{T}-T\right)/T$ where, the exact period is $T=2\pi /\omega$, and the numerically obtained period is $\stackrel{-}{T}=2\pi /\stackrel{-}{\omega }$. The algorithmic damping ratio is defined as$\stackrel{-}{\xi }=-\mathit{ln}\left(A_2\right)/\left(2\stackrel{-}{\Omega }\right)$, where $\stackrel{-}{\Omega }=\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{t}\mathrm{a}\mathrm{n}\left(\sqrt{A_2/A_1^2-1}\right)/\left(\mathrm{\Delta }t\sqrt{1-{\xi }^2}\right)$ and $\stackrel{-}{\Omega }=\stackrel{-}{\omega \mathrm{\Delta }t}$. For more details, please see Ref.[²⁸].

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, $u_0=1.0$, $v_0=0.0$, $\xi =0.0$, $f\left(t\right)=0.0$, and $\omega =2\pi$ are used. For damped cases, $u_0=1.0$, $v_0=0.0$, $\xi =0.1$, $f\left(t\right)=0.0$, and $\omega =2\pi$ are used. $\Delta 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 $\Delta 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 $\ddot{\mathrm{u}}\left(\mathrm{t}\right)+2\mathrm{\xi }\mathrm{\omega }\dot{\mathrm{u}}\left(\mathrm{t}\right)+{\mathrm{\omega }}^2\mathrm{u}\left(\mathrm{t}\right)=\mathrm{q}\left(\mathrm{t}\right)$. $u_0=1.0$, $v_0=0.0$, $\xi =0.0$, $f\left(t\right)=0.0$, $\omega =2\pi$ and $\Delta t=0.2$ are used.

Figure 5:
Comparison of displacements the single degree of freedom system $\ddot{\mathrm{u}}\left(\mathrm{t}\right)+2\mathrm{\xi }\mathrm{\omega }\dot{\mathrm{u}}\left(\mathrm{t}\right)+{\mathrm{\omega }}^2\mathrm{u}\left(\mathrm{t}\right)=\mathrm{q}\left(\mathrm{t}\right)$. $u_0=1.0$, $v_0=0.0$, $\xi =0.0$, $f\left(t\right)=0.0$, $\omega =2\pi$ and $\Delta t=0.1$ are used.

Figure 6:
Comparison of displacements the single degree of freedom system $\ddot{\mathrm{u}}\left(\mathrm{t}\right)+2\mathrm{\xi }\mathrm{\omega }\dot{\mathrm{u}}\left(\mathrm{t}\right)+{\mathrm{\omega }}^2\mathrm{u}\left(\mathrm{t}\right)=\mathrm{q}\left(\mathrm{t}\right)$. $u_0=1.0$, $v_0=0.0$, $\xi =0.1$, $f\left(t\right)=0.0$, $\omega =2\pi$ and $\Delta t=0.2$ are used.

Figure 7:
Comparison of displacements the single degree of freedom system $\ddot{\mathrm{u}}\left(\mathrm{t}\right)+2\mathrm{\xi }\mathrm{\omega }\dot{\mathrm{u}}\left(\mathrm{t}\right)+{\mathrm{\omega }}^2\mathrm{u}\left(\mathrm{t}\right)=\mathrm{q}\left(\mathrm{t}\right)$. $u_0=1.0$, $v_0=0.0$, $\xi =0.1$, $f\left(t\right)=0.0$, $\omega =2\pi$ and $\Delta 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 [³¹-³³] is expressed as

$\ddot{\mathrm{u}}+{\mathrm{s}}_1\left(1+{\mathrm{s}}_2{\mathrm{u}}^2\right)\mathrm{u}=0,$ (73a)

$\mathrm{u}\left(0\right)={\mathrm{u}}_0,$ (73b)

$\dot{\mathrm{u}}\left(0\right)={\mathrm{v}}_0,$ (73c)

where ${\mathrm{s}}_1>0$and ${\mathrm{s}}_2>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, ${\mathrm{u}}_0=1.5$, ${\mathrm{v}}_0=0.0$, ${\mathrm{s}}_1=100.0$, and ${\mathrm{s}}_2=10.0$ are used [³²]. All properties are dimensionless. The nonlinear period $\mathrm{T}$ of the problem is computed by using 10th-order implicit Kim and Reddy method [⁷,²⁰] with time step $∆\mathrm{t}=0.000001$. The nonlinear period of the problem is computed as $\mathrm{T}=0.15153$. All methods used $∆\mathrm{t}=\mathrm{T}/32$.

Figure 8:
Phase portrait of the hardening spring problem $\ddot{\mathrm{u}}+{\mathrm{s}}_1\left(1+{\mathrm{s}}_2{\mathrm{u}}^2\right)\mathrm{u}=0$, ${\mathrm{u}}_0=1.5$, ${\mathrm{v}}_0=0.0$, ${\mathrm{s}}_1=100.0$, and ${\mathrm{s}}_2=10.0$are used with$∆\mathrm{t}=\mathrm{T}/32$ to complete fifty cycles (i.e., until $t=50\mathrm{T}$). (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 $∆\mathrm{t}=\mathrm{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 $\ddot{\mathrm{u}}+{\mathrm{s}}_1\left(1+{\mathrm{s}}_2{\mathrm{u}}^2\right)\mathrm{u}=0$, ${\mathrm{u}}_0=1.5$, ${\mathrm{v}}_0=0.0$, ${\mathrm{s}}_1=100.0$, and ${\mathrm{s}}_2=10.0$are used with$∆\mathrm{t}=\mathrm{T}/32$. (a) short-term solutions ($0\le t\le 2T$); (b) long-term solutions ($48T\le t\le 50T$).

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 [³¹-³³] is also solved by using the new methods. The motion of the softening elastic spring is expressed as

$\ddot{\mathrm{u}}+\mathrm{s}\tanh \left(\mathrm{u}\right)=0,$ (74a)

$\mathrm{u}\left(0\right)={\mathrm{u}}_0,$ (74b)

$\dot{\mathrm{u}}\left(0\right)={\mathrm{v}}_0,$ (74c)

where $\mathrm{s}>0.0$. Here we use$\mathrm{s}=100.0$,${\mathrm{u}}_0=4.0$, and${\mathrm{v}}_0=0.0$ [³²]. 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 $\mathrm{\Delta }\mathrm{t}=\mathrm{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 $\ddot{\mathrm{u}}+\mathrm{s}\tanh \left(\mathrm{u}\right)=0$, $\mathrm{s}=100.0$,${\mathrm{u}}_0=4.0$, and${\mathrm{v}}_0=0.0$are used with$∆\mathrm{t}=\mathrm{T}/32$ to complete fifty cycles (i.e., until $t=50\mathrm{T}$). (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 $\ddot{\mathrm{u}}+\mathrm{s}\tanh \left(\mathrm{u}\right)=0$, $\mathrm{s}=100.0$,${\mathrm{u}}_0=4.0$, and${\mathrm{v}}_0=0.0$are used with$∆\mathrm{t}=\mathrm{T}/32$. (a) short-term solutions ($0\le t\le 2T$); (b) long-term solutions ($48T\le t\le 50T$).

5.4 Nonlinear single pendulum problem

Figure 12:
Description of nonlinear single pendulum [⁵]

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

$\ddot{\mathrm{\theta }}+{\mathrm{\omega }}^2\sin \left(\mathrm{\theta }\right)=0,$ (75)

with initial conditions

$\mathrm{\theta }\left(0\right)={\mathrm{\theta }}_0,$ (76b)

$\dot{\mathrm{\theta }}\left(0\right)={\dot{\mathrm{\theta }}}_0,$ (76c)

where $\mathrm{\omega }$ is$\sqrt{g/L}$, $\mathrm{\theta }\left(\mathrm{t}\right)$ is the angle between the rigid rod and the vertical line at time$\mathrm{t}$, $g$ is the gravitational constant, and $\mathrm{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.[⁵], [³⁴], and [³⁵].

For the test of the new methods, two sets of initial conditions are used. First, the initial conditions that have been used in Refs.[⁵,⁷,¹⁴] 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 ${\theta }_0=0.0$ and ${\dot{\theta }}_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.[¹⁴]. 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 $\left|\left({\theta }_{\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}}-{\theta }_{\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{c}\mathrm{t}}\right)/{\theta }_{\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{c}\mathrm{t}}\right|$ at $t=T/4$. ${\theta }_{\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{c}\mathrm{t}}$ 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 $\left|\left({\theta }_{\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}}-{\theta }_{\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{c}\mathrm{t}}\right)/{\theta }_{\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{c}\mathrm{t}}\right|$ 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 [¹⁰] 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

$\mathrm{m}\ddot{\mathrm{r}}-\mathrm{m}\left({\mathrm{L}}_0+\mathrm{r}\right){\dot{\mathrm{\theta }}}^2-\mathrm{m}\mathrm{g}\cos \mathrm{\theta }+\mathrm{k}\mathrm{r}=0,$ (76b)

$\mathrm{m}\ddot{\mathrm{\theta }}+\frac{\mathrm{m}\left(2\dot{\mathrm{r}}\dot{\mathrm{\theta }}+\mathrm{g}\sin \mathrm{\theta }\right)}{{\mathrm{L}}_0+\mathrm{r}}=0,$ (76c)

and the initial conditions are

$\mathrm{r}\left(0\right)={\mathrm{r}}_0,\mathrm{\theta }\left(0\right)={\mathrm{\theta }}_0,\dot{\mathrm{r}}\left(0\right)={\dot{\mathrm{r}}}_0,\dot{\mathrm{\theta }}\left(0\right)={\dot{\mathrm{\theta }}}_0$ (76b)

where $\mathrm{r}$ and $\mathrm{\theta }$ are the displacements in the radial and circumferential directions, respectively, $\mathrm{m}$ is the mass of the pendulum, $g$ is the gravitational constant, ${\mathrm{L}}_0$ is the length of the undeformed spring, and $\mathrm{k}$ is the spring constant. In this numerical test, $\mathrm{m}=1.0$, $\mathrm{g}=9.81$, $\mathrm{L}=0.5$, $\mathrm{k}=98.1$, $r_0=0.25$, ${\mathrm{\theta }}_0=\mathrm{\pi }/2=1.570796326794896$, ${\dot{\mathrm{r}}}_0=0.0$ and ${\dot{\mathrm{\theta }}}_0=0.0$ are used. All properties are dimensionless.

Figure 18:
Description of spring pendulum [⁵].

Figure 19:
Comparison of angles of spring-pendulum. 0.1 is used for $∆\mathrm{t}$ in all methods.

Figure 20:
Comparison of angles of spring-pendulum. 0.05 is used for $∆\mathrm{t}$ in all methods.

Figure 21:
Comparison of relative error $\left|\left(r_{\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}}-r_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}\right)/r_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}\right|$ at $t=0.1$. $r_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}$ is computed by using the 10th-order accurate method of Kim and Reddy [⁷,²⁰].

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 [³⁶].

Figure 22:
Description of double pendulum.

Figure 23:
Comparison of $x_1$ for $90\le t\le 95$.

Figure 24:
Comparison of $y_1$ for $90\le t\le 95$.

Figure 25:
Comparison of $x_2$ for $90\le t\le 95$.

Figure 26:
Comparison of $y_2$ for $90\le t\le 95$.

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

${\ddot{\mathrm{\theta }}}_1-\frac{{\mathrm{m}}_2{\mathrm{l}}_2{{\dot{\mathrm{\theta }}}_2}^2\sin \left({\mathrm{\theta }}_2-{\mathrm{\theta }}_1\right)-\left({\mathrm{m}}_1+{\mathrm{m}}_2\right)\mathrm{g}\sin \left({\mathrm{\theta }}_1\right)}{{\mathrm{l}}_1\left({\mathrm{m}}_1+{\mathrm{m}}_2\right)-{\mathrm{m}}_2{\mathrm{l}}_1{\cos }^2\left({\mathrm{\theta }}_2-{\mathrm{\theta }}_1\right)}=0,$ (79a)

${\ddot{\mathrm{\theta }}}_2+\frac{\left({\mathrm{m}}_1+{\mathrm{m}}_2\right){\mathrm{l}}_1{{\dot{\mathrm{\theta }}}_1}^2\sin \left({\mathrm{\theta }}_2-{\mathrm{\theta }}_1\right)+\left({\mathrm{m}}_1+{\mathrm{m}}_2\right)\mathrm{g}\sin \left({\mathrm{\theta }}_2\right)}{{\mathrm{l}}_2\left({\mathrm{m}}_1+{\mathrm{m}}_2\right)-{\mathrm{m}}_2{\mathrm{l}}_2{\cos }^2\left({\mathrm{\theta }}_2-{\mathrm{\theta }}_1\right)}=0,$ (79b)

${\left\{{\mathrm{\theta }}_1\left(0\right),{\mathrm{\theta }}_2\left(0\right),{\dot{\mathrm{\theta }}}_1\left(0\right),{\dot{\mathrm{\theta }}}_2\left(0\right)\right\}}^{\mathrm{T}}={\left\{{\mathrm{\vartheta }}_1,{\mathrm{\vartheta }}_2,{\dot{\mathrm{\vartheta }}}_1,{\dot{\mathrm{\vartheta }}}_2\right\}}^{\mathrm{T}},$ (79c)

where ${\mathrm{\theta }}_1$ and ${\mathrm{\theta }}_2$ are angle of the first and second pendulums, respectively, ${\mathrm{m}}_1$ and ${\mathrm{m}}_2$ are masses of the first and second pendulums, respectively, ${\mathrm{l}}_1$ and ${\mathrm{l}}_2$ are weightless rigid rods connected to pendulums as shown in Fig.22, and ${\mathrm{\vartheta }}_1$ and ${\mathrm{\vartheta }}_2$ are the initial angles. By using ${\mathrm{\theta }}_1$ and ${\mathrm{\theta }}_2$, positions of the pendulum are computed as follows:

${\mathrm{x}}_1={\mathrm{l}}_1\sin {\mathrm{\theta }}_1,$ (80a)

${\mathrm{y}}_1=-{\mathrm{l}}_1\cos {\mathrm{\theta }}_1,$ (80b)

${\mathrm{x}}_2={\mathrm{x}}_1+{\mathrm{l}}_2\sin {\mathrm{\theta }}_2,$ (80c)

${\mathrm{y}}_2={\mathrm{y}}_1-{\mathrm{l}}_2\cos {\mathrm{\theta }}_2,$ (80d)

Dimensionless properties ${\mathrm{m}}_1=1.0$, ${\mathrm{m}}_2=1.0$, ${\mathrm{l}}_1=1.0$, ${\mathrm{l}}_2=1.0$, and $\mathrm{g}=1.0$ are used, and the initial conditions are chosen as ${\left\{{\mathrm{\vartheta }}_1,{\mathrm{\vartheta }}_2,{\dot{\mathrm{\vartheta }}}_1,{\dot{\mathrm{\vartheta }}}_2\right\}}^{\mathrm{T}}={\left\{0.0,\mathrm{\pi }/2,0.0,0.0\right\}}^{\mathrm{T}}$. 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 $\left|\left({y2}_{\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{l}}-{y2}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}\right)/{y2}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}\right|$ at $t=1.0$. ${y2}_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}}$ is computed by using the 10th-order accurate method of Kim and Reddy [⁷,²⁰].

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.
• ^[2]
  O. C. Zienkiewicz and R. L. Taylor. The Finite Element Method for Solid and Structural Mechanics. Butterworth-Heinemann, 2005.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[18]
  T. J. R. Hughes. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover, 2012.
• ^[19]
  K. J. Bathe. Finite Element Procedures. Klaus-Jurgen Bathe, 2006.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[26]
  H. Adeli, J. M. Gere, and W. Weaver. Algorithms for nonlinear structural dynamics. Journal of the Structural Division, 104(2):263-280, 1978.
• ^[27]
  T. J. R. Hughes. The finite element method: linear static and dynamic finite element analysis. Courier Corporation, 2012.
• ^[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.
• ^[29]
  H. M. Hilber. Analysis and design of numerical integration methods in structural dynamics. PhD thesis, University of California Berkeley, 1976.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[37]
  T. Stachowiak and T. Okada. A numerical analysis of chaos in the double pendulum. Chaos, Solitons & Fractals, 29(2):417-422, 2006.

Edited by

Editor:

Marco L. Bittencourt

Publication Dates

• Publication in this collection
  22 July 2019
• Date of issue
  2019

History

• Received
  29 Apr 2019
• Reviewed
  18 May 2019
• Accepted
  22 May 2019
• Published
  03 June 2019