# A Modified Precise Integration Method for Long-Time Duration Dynamic Analysis

Authors: Ce Huang; Minghui Fu
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

This paper presents a modified Precise Integration Method (PIM) for long-time duration dynamic analysis. The fundamental solution and loading operator matrices in the first time substep are numerically computed employing a known unconditionally stable method (referred to as original method in this paper). By using efficient recursive algorithms to evaluate these matrices in the first time-step, the same numerical results as those using the original method with small time-step are obtained. The proposed method avoids the need of matrix inversion and numerical quadrature formulae, while the particular solution obtained has the same accuracy as that of the homogeneous solution. Through setting a proper value of the time substep, satisfactory accuracy and numerical dissipation can be achieved.


---

## Body

## 1. Introduction

In structural dynamic analysis, direct time integration algorithms are widely used to solve the motion equation ‎[1]. Many researchers‎ [2–4] have focused on the improvement of numerical properties of time integration algorithms, such as numerical stability, accuracy, numerical dissipation, and computational efficiency.An alternative way to evaluate the dynamic response is to use the Precise Integration Method (PIM-L)‎ [5], where a numerical exponential matrix function is accurately computed to approximate the fundamental solution matrix and the applied loading is simulated by using piecewise linear functions. For linear dynamic homogeneous systems, the PIM-L can obtain precise numerical results approaching the exact solution. However, its accuracy is limited by the simulation of the applied loading in nonhomogeneous systems. Lin et al. [6] applied Fourier series expansions to simulate the applied loading and proposed the PIM-F method. Both the PIM-L and PIM-F require matrix inversion to obtain the particular solution. However, if the objective matrix is singular or nearly singular, its inverse matrix either does not exist or is hard to obtain steadily.In order to avoid matrix inversion, Gu et al.‎ [7] converted nonhomogeneous equations into homogeneous ones at the expense of dimensional expansion, which significantly increased storage and computational cost. Afterwards, Wang et al.‎ [8] proposed a PIM based on Lagrange piecewise polynomials and improved the dimensional expansion approach. Another practical option to compute the particular solution without matrix inversion is to use the numerical quadrature formulae. Wang et al.‎ [9, 10] applied Gauss quadrature, Simpson quadrature, and Romberg integration formulae to integrate the particular solution. The accuracy of particular solution using these methods depends crucially on the number of integration points used in each time-step and the size of time-step. Basically, the accuracy of the particular solution obtained by these aforementioned methods may not achieve the same accuracy as that of the homogeneous solution, especially when a large time-step is used.In conventional PIMs, the fundamental solution matrix in the first time substep is approximated by the Taylor series expansion of the exponential matrix function. Fung ‎[11] halved the computational effort required using step-response and impulsive-response matrices and considering the symmetry of these matrices. Later on, the Krylov subspace method ‎[12] and the Ritz vectors ‎[13] were also applied to improve the computational efficiency. Gao et al.‎ [14] utilized the special feature of the dynamic system to reduce the computational effort involved in establishing the exponential matrix function. Wang et al.‎ [9] investigated the numerical stability and computational accuracy of the conventional PIMs in detail. Subsequently, Wang et al.‎ [10] proposed the PIM based on Padé approximation, where the generalized Padé approximation was used to replace the Taylor series expansion to achieve unconditional stability and controllable numerical dissipation.This paper aims to modify and improve the existing PIM for solving long-time duration dynamics problems. The numerical fundamental solution and loading operator matrices are obtained by applying an unconditionally stable time integration method and recursive formulae. Using this proposed approach, the same accuracy of particular solution as that of the homogeneous solution can be achieved even for a large time-step, which is considered as a major improvement of PIMs. Besides, the proposed method can obtain approximately the same numerical results as those using the original method with a small time-step; thus it can also be considered as an accelerator of conventional time integration methods.


## 2. The Precise Integration Methods

The motion equation of a discretized linear structural model with n-degrees-of-freedom is given by(1)Mx¨t+Cx˙t+Kxt=Ft,xtt=0=x0,x˙tt=0=x˙0,where M, C, and K are the mass, damping, and stiffness matrices, respectively; F(t) is the vector of applied external loading; x(t), x˙(t), and x¨(t) are the displacement, velocity, and acceleration vectors of the system, respectively.Assumingv=xMx˙, r(t)=0F(t), (1) can be transformed into the following first-order ordinary differential equation:(2)v˙=Hv+rt,where H=0M-1-K-CM-1. An analytical solution of (2) can be expressed as follows‎ [15]: (3)vt=eHtv0+∫0teHt-srsdsApplying time discretization with a time-step τ and assuming the analysis starts at the instant tm=mτ, vm+1 can then be represented by vm: (4)vm+1=Tτvm+v∗m+1Tτ=eHτ,v∗m+1=∫tmtm+1eHtm+1-srsds,where tm=mτ, tm+1=(m+1)τ,vm=v(tm), and vm+1=v(tm+1). T(τ) can be written as(5)Tτ=eHτ=eHΔt2N=TΔt2N,where Δt=τ/2N, N is a given integer. The exponential matrix eHΔt can be approximated by using the Taylor series expansion or the generalized Padé approximation. Taking the Taylor series expansion as an example, T(Δt) can be written as(6)TΔt≈I+T~a0,where(7)T~a0=HΔt+HΔt22!+HΔt33!+⋯+HΔtpp!,where p is the truncation order of the Taylor series expansion. Substituting (6) into (5) gives(8)T2Δt≈I+T~a1=I+2T~a0+T~a0•T~a0T22Δt≈I+T~a2=I+2T~a1+T~a1•T~a1⋮T2NΔt≈I+T~aN=I+2T~aN-1+T~aN-1•T~aN-1T~a(j)(j=0,1,2...N) can be evaluated by applying the following recursive formula:(9)T~aj+1=2T~aj+T~aj•T~ajj=0,1,2…N-1


## 3. A Modified PIM for Long-Time Duration Dynamic Analysis

### 3.1. The Numerical Fundamental Solution Matrix

In the standard finite-element semidiscrete systems, only a few low-frequency modes can represent the original continuous system. The quality of time integration algorithms depends crucially on their dissipative property. It has been found that the PIM is conditionally stable ‎[9]. In order to satisfy the stability conditions, the time substep has to be smaller than the highest modal period of the particular discretized structure. In this case, the dissipative property would be eliminated. Wang et al.‎ [10] proposed the PIM based on Padé approximations to achieve unconditional stability and controllable dissipation, in which the Taylor series expansion of the matrix exponential function is replaced with the generalized Padé approximations. In this work, the fundamental solution matrix is approximated using numerical Green’s function matrix to achieve desirable numerical properties. An analytical expression of x(mτ+Δt), x˙(mτ+Δt) can be written in terms of x(mτ), x˙(mτ)  [16–19]:(10)xmτ+Δt=GΔtC+G˙ΔtMxmτ+GΔtMx˙mτx˙mτ+Δt=G˙ΔtC+G¨ΔtMxmτ+G˙ΔtMx˙mτ,where G(Δt), G˙(Δt) represent Green’s function matrix of the model and its derivative, respectively. Thus, the analytical expression of T(Δt) is given as follows:(11)TΔt=GΔtC+G˙ΔtMGΔtG˙ΔtC+G¨ΔtMG˙ΔtG(Δt) and G˙(Δt) can be calculated numerically by solving the following equation:(12)MG¨t+CG˙t+KGt=0G0=0,G˙0=M-1In order to satisfy the unconditional stability, a known unconditionally stable time integration method (referred to as original method in the following discussion) with time-stepΔt is used to solve (12). By selecting a proper original method and setting the value of time substep, one can obtain satisfactory accuracy, numerical dissipation, and other numerical properties.SubstitutingG(Δt) into the following equation, Ta(0) can be obtained: (13)Ta0=TΔt-ISubstitutingTa(0) into (8) and (9), T(τ) can be evaluated. On the other hand, there is another method to evaluate T(Δt), which does not need the inverse of the mass matrix. The j column of T(Δt) can be obtained through solving the following equation: (14)v˙=Hv,v0=1j,where 1ij is the i entry of the unity vector 1j. 1ij is zero if i≠j and 1 if i=j (1ij=δij, δij being the Kronecker delta).It can be seen that the method incorporates the PIM with Green’s function matrix based method ‎[16–19]. For this method, the sparsity of the fundamental matrix would be weakened as the number of time substeps increases. As a result, this method requires more storage usage compared to conventional time integration methods. However, much larger time-step can be used, which causes less computational cost. Users can find a balance between storage usage and computation cost by setting proper N. The improvement of the storage efficiency is feasible, which needs further investigation.


### 3.2. The Loading Operator Matrix

For homogeneous systems, the PIMs can give a precise numerical solution. However, for nonhomogeneous dynamic systems, the accuracy of numerical results depends essentially on the accuracy of the particular solution. The accuracy of the PIM-L and PIM-F methods is mainly limited by the error in the matrix inversion and the applied loading simulation. For the PIM proposed by Wang et al. ‎[10], the main error was caused by the process of the numerical quadrature formulae. In order to obtain a high accuracy, a large number of integration points were used, which required extra computational effort. When a larger time-step is adopted, the particular solution obtained may not attain the same accuracy as that of the homogeneous solution. In this study, a novel method is proposed to evaluate the particular solution, which avoids the matrix inversion. It can retain the same accuracy of the particular solution as that of the homogeneous solution even though the time-step is very large.In practical engineering problems, abounding loadings can be expressed in a periodical form, which can be approximated by using the Fourier series expansion(15)Ft≈B0+∑i=1lAisin⁡iωt+Bicos⁡iωt,where l is the truncation order of the Fourier series expansion. According to Euler’s formula, the trigonometric functions in (15) can be transformed as an exponential function. (16)eiωt=cos⁡ωt+isin⁡ωti2=-1An exponential function loading can be expressed as follows in the time interval [tm,tm+1],(17)rt=r~ebt,where r~ is a constant vector and b is a real or complex number. Substituting (17) into (4), the following equation can be obtained:(18)v∗m+1=∫tmtm+1eHtm+1-tr~ebtdtAssuming s=tm+1-t, (18) becomes(19)v∗m+1=∫0τeHsebtm+1-sdsr~=ebtm+1S~r~(20)S~=∫0τeHsds=∑j=02N-1∫jΔtj+1ΔteHse-bsds,where S~r~ is the loading operator matrix. It can be seen that S~r~ is not associated with the integration interval [tm,tm+1], thus it is calculated only once during the whole step-by-step integration procedure. As (19) is the analytical expression of the particular solution, once S~r~ is properly computed, the particular solution would become extremely accurate. In order to ensure the same accuracy of the particular solution as that of the homogeneous solution, eHΔt in (20) can be approximated by T(Δt), which is generated in Section 3.1. S~r~ can be written as (21)S~r~≈∑j=02N-1TΔtj-1e-bj-1Δt∫0ΔteHse-bsdsr~It can be seen that S~r~ is the sum of 2N terms. Let S~(j)r~(j=0,1,2...N) denote the sum of the first 2j terms of (21). The following recursive formula can be used to evaluate S~(N)r~ from S~(0)r~: (22)S~j+1r~=Z2jS~jr~Z2j=e-2jbΔtI+Taj,where Z=T(Δt)e-bΔt. Ta(j)(j=0,1,2...N) is obtained in Section 3.1. The matrices Ta(j) and S~(j)r~ are computed using the following looping algorithm:(23)forj=1:NTa=2∗Ta+Ta∗TaSr=Sr+e-2jbΔt∗I+Ta∗SrendFor conventional PIMs, the algorithm (23) can also be used to acquire S~(N)r~. Thus the particular solution would possess the same accuracy as the homogeneous solution. In this work, S~(0)r~ is obtained numerically by solving the following equation using the original method with step-size Δt: (24)v˙=Hv+r~ebt,vtm=0(25)vtm+Δt=ebΔtS~0r~During the recursive process(22), only matrix-vector multiplications are required. Therefore, it can be summarized that only N times of matrix-matrix multiplication are needed to obtain the numerical fundamental solution and loading operator matrices (as the computational effort of matrix-vector multiplication is negligible compared to that of matrix-matrix multiplication, only matrix-matrix multiplication is considered).


### 3.3. Properties Analysis

As the recurrence relations ((9) and (22)) are analytical expressions, the accuracy, stability, dissipation, and other numerical properties of the proposed method are determined only by the original method and the size of time substep. In this section, the stability of the proposed method is investigated. For the analysis of other properties, one can refer to [16] for detailed information.Considering the following undamped unforced single degree of freedom (SDOF) problem:(26)x¨+ω2x=0,where ω is the natural frequency of the model and the natural period is Tn=2π/ω. Assuming the time-step is τ, time substep is Δt (Δt=τ/2N). The following recursive relationship can be obtained using the original method:(27)xtm+τ2Nx˙tm+τ2N=Aτ2Nxtmx˙tm,where A(τ/2N) stands for the amplification matrix of the original method. According to (5), the amplification matrix of the proposed method is obtained as A(τ)=A(τ/2N)2N. The stability condition is ρ(A)≤1, where ρ(A)=max⁡(λ1,2) is the spectral radius of A and λ1,2 are the eigenvalues of A.In order to investigate the stability of the proposed method, the Newmark (α=3/10, δ=11/20), the Bathe [2] and the fourth-order Runge-Kutta (RK-4) methods are selected as the original method, respectively.Figure1 shows the spectral radius of the proposed method, where MPIM-Newmark, MPIM-Bathe, and MPIM-RK denote the results obtained using the Newmark, Bathe, and RK-4 method as the original method, respectively. For all curves the value of spectral radius is 1 until τ/Tn is equal to a certain value and thereafter it rapidly diminishes. As N increases, this value increases. The MPIM-Newmark has lower limit values than the Newmark method. Compared to RK-4, the stability condition of MPIM-RK4 is relaxed.Figure 1
Spectral radius of the proposed method. (a)MPIM-Newmark, (b)MPIM-Bathe, and (c)MPIM-RK4.
(a) (b) (c)


## 3.1. The Numerical Fundamental Solution Matrix

In the standard finite-element semidiscrete systems, only a few low-frequency modes can represent the original continuous system. The quality of time integration algorithms depends crucially on their dissipative property. It has been found that the PIM is conditionally stable ‎[9]. In order to satisfy the stability conditions, the time substep has to be smaller than the highest modal period of the particular discretized structure. In this case, the dissipative property would be eliminated. Wang et al.‎ [10] proposed the PIM based on Padé approximations to achieve unconditional stability and controllable dissipation, in which the Taylor series expansion of the matrix exponential function is replaced with the generalized Padé approximations. In this work, the fundamental solution matrix is approximated using numerical Green’s function matrix to achieve desirable numerical properties. An analytical expression of x(mτ+Δt), x˙(mτ+Δt) can be written in terms of x(mτ), x˙(mτ)  [16–19]:(10)xmτ+Δt=GΔtC+G˙ΔtMxmτ+GΔtMx˙mτx˙mτ+Δt=G˙ΔtC+G¨ΔtMxmτ+G˙ΔtMx˙mτ,where G(Δt), G˙(Δt) represent Green’s function matrix of the model and its derivative, respectively. Thus, the analytical expression of T(Δt) is given as follows:(11)TΔt=GΔtC+G˙ΔtMGΔtG˙ΔtC+G¨ΔtMG˙ΔtG(Δt) and G˙(Δt) can be calculated numerically by solving the following equation:(12)MG¨t+CG˙t+KGt=0G0=0,G˙0=M-1In order to satisfy the unconditional stability, a known unconditionally stable time integration method (referred to as original method in the following discussion) with time-stepΔt is used to solve (12). By selecting a proper original method and setting the value of time substep, one can obtain satisfactory accuracy, numerical dissipation, and other numerical properties.SubstitutingG(Δt) into the following equation, Ta(0) can be obtained: (13)Ta0=TΔt-ISubstitutingTa(0) into (8) and (9), T(τ) can be evaluated. On the other hand, there is another method to evaluate T(Δt), which does not need the inverse of the mass matrix. The j column of T(Δt) can be obtained through solving the following equation: (14)v˙=Hv,v0=1j,where 1ij is the i entry of the unity vector 1j. 1ij is zero if i≠j and 1 if i=j (1ij=δij, δij being the Kronecker delta).It can be seen that the method incorporates the PIM with Green’s function matrix based method ‎[16–19]. For this method, the sparsity of the fundamental matrix would be weakened as the number of time substeps increases. As a result, this method requires more storage usage compared to conventional time integration methods. However, much larger time-step can be used, which causes less computational cost. Users can find a balance between storage usage and computation cost by setting proper N. The improvement of the storage efficiency is feasible, which needs further investigation.


## 3.2. The Loading Operator Matrix

For homogeneous systems, the PIMs can give a precise numerical solution. However, for nonhomogeneous dynamic systems, the accuracy of numerical results depends essentially on the accuracy of the particular solution. The accuracy of the PIM-L and PIM-F methods is mainly limited by the error in the matrix inversion and the applied loading simulation. For the PIM proposed by Wang et al. ‎[10], the main error was caused by the process of the numerical quadrature formulae. In order to obtain a high accuracy, a large number of integration points were used, which required extra computational effort. When a larger time-step is adopted, the particular solution obtained may not attain the same accuracy as that of the homogeneous solution. In this study, a novel method is proposed to evaluate the particular solution, which avoids the matrix inversion. It can retain the same accuracy of the particular solution as that of the homogeneous solution even though the time-step is very large.In practical engineering problems, abounding loadings can be expressed in a periodical form, which can be approximated by using the Fourier series expansion(15)Ft≈B0+∑i=1lAisin⁡iωt+Bicos⁡iωt,where l is the truncation order of the Fourier series expansion. According to Euler’s formula, the trigonometric functions in (15) can be transformed as an exponential function. (16)eiωt=cos⁡ωt+isin⁡ωti2=-1An exponential function loading can be expressed as follows in the time interval [tm,tm+1],(17)rt=r~ebt,where r~ is a constant vector and b is a real or complex number. Substituting (17) into (4), the following equation can be obtained:(18)v∗m+1=∫tmtm+1eHtm+1-tr~ebtdtAssuming s=tm+1-t, (18) becomes(19)v∗m+1=∫0τeHsebtm+1-sdsr~=ebtm+1S~r~(20)S~=∫0τeHsds=∑j=02N-1∫jΔtj+1ΔteHse-bsds,where S~r~ is the loading operator matrix. It can be seen that S~r~ is not associated with the integration interval [tm,tm+1], thus it is calculated only once during the whole step-by-step integration procedure. As (19) is the analytical expression of the particular solution, once S~r~ is properly computed, the particular solution would become extremely accurate. In order to ensure the same accuracy of the particular solution as that of the homogeneous solution, eHΔt in (20) can be approximated by T(Δt), which is generated in Section 3.1. S~r~ can be written as (21)S~r~≈∑j=02N-1TΔtj-1e-bj-1Δt∫0ΔteHse-bsdsr~It can be seen that S~r~ is the sum of 2N terms. Let S~(j)r~(j=0,1,2...N) denote the sum of the first 2j terms of (21). The following recursive formula can be used to evaluate S~(N)r~ from S~(0)r~: (22)S~j+1r~=Z2jS~jr~Z2j=e-2jbΔtI+Taj,where Z=T(Δt)e-bΔt. Ta(j)(j=0,1,2...N) is obtained in Section 3.1. The matrices Ta(j) and S~(j)r~ are computed using the following looping algorithm:(23)forj=1:NTa=2∗Ta+Ta∗TaSr=Sr+e-2jbΔt∗I+Ta∗SrendFor conventional PIMs, the algorithm (23) can also be used to acquire S~(N)r~. Thus the particular solution would possess the same accuracy as the homogeneous solution. In this work, S~(0)r~ is obtained numerically by solving the following equation using the original method with step-size Δt: (24)v˙=Hv+r~ebt,vtm=0(25)vtm+Δt=ebΔtS~0r~During the recursive process(22), only matrix-vector multiplications are required. Therefore, it can be summarized that only N times of matrix-matrix multiplication are needed to obtain the numerical fundamental solution and loading operator matrices (as the computational effort of matrix-vector multiplication is negligible compared to that of matrix-matrix multiplication, only matrix-matrix multiplication is considered).


## 3.3. Properties Analysis

As the recurrence relations ((9) and (22)) are analytical expressions, the accuracy, stability, dissipation, and other numerical properties of the proposed method are determined only by the original method and the size of time substep. In this section, the stability of the proposed method is investigated. For the analysis of other properties, one can refer to [16] for detailed information.Considering the following undamped unforced single degree of freedom (SDOF) problem:(26)x¨+ω2x=0,where ω is the natural frequency of the model and the natural period is Tn=2π/ω. Assuming the time-step is τ, time substep is Δt (Δt=τ/2N). The following recursive relationship can be obtained using the original method:(27)xtm+τ2Nx˙tm+τ2N=Aτ2Nxtmx˙tm,where A(τ/2N) stands for the amplification matrix of the original method. According to (5), the amplification matrix of the proposed method is obtained as A(τ)=A(τ/2N)2N. The stability condition is ρ(A)≤1, where ρ(A)=max⁡(λ1,2) is the spectral radius of A and λ1,2 are the eigenvalues of A.In order to investigate the stability of the proposed method, the Newmark (α=3/10, δ=11/20), the Bathe [2] and the fourth-order Runge-Kutta (RK-4) methods are selected as the original method, respectively.Figure1 shows the spectral radius of the proposed method, where MPIM-Newmark, MPIM-Bathe, and MPIM-RK denote the results obtained using the Newmark, Bathe, and RK-4 method as the original method, respectively. For all curves the value of spectral radius is 1 until τ/Tn is equal to a certain value and thereafter it rapidly diminishes. As N increases, this value increases. The MPIM-Newmark has lower limit values than the Newmark method. Compared to RK-4, the stability condition of MPIM-RK4 is relaxed.Figure 1
Spectral radius of the proposed method. (a)MPIM-Newmark, (b)MPIM-Bathe, and (c)MPIM-RK4.
(a) (b) (c)


## 4. Numerical Examples

In this section, two numerical examples are investigated to evaluate the validity of the proposed method. The proposed method solution is compared with those obtained by the SPIM (q=2, m=8, μ=0.8), the ExGA-RK method [16], and the original method with a small time-step. The Newmark (α=3/10, δ=11/20), the Bathe [2], and the fourth-order Runge-Kutta (RK-4) methods are selected as the original method, respectively. The RK-4 method is selected just for the sake of the comparison with the ExGA-RK method. An unconditionally stable method can be used when solving particular dynamic problems.In Figures2–11, MPIM-RK, MPIM-Newmark, and MPIM-Bathe denote the solution obtained by the proposed modified PIM using the RK-4, Newmark, and the Bathe method as the original method, respectively.Figure 2
Displacement of DOF1 in Case1 using different methods: (a) displacement; (b) displacement error.
(a) (b)Figure 3
Velocity of DOF1 in Case1 using different methods: (a) velocity; (b) velocity error.
(a) (b)Figure 4
Error of numerical fundamental solution matrix versus N.Figure 5
Displacement of DOF2 in Case2 (N=4): (a) displacement; (b) displacement error.
(a) (b)Figure 6
Velocity of DOF2 in Case2 (N=4): (a) velocity; (b) velocity error.
(a) (b)Figure 7
Displacement of DOF2 in Case2 (N=6): (a) displacement; (b) displacement error.
(a) (b)Figure 8
Velocity of DOF2 in Case2 (N=6): (a) velocity; (b) velocity error.
(a) (b)Figure 9
Displacement of DOF2 in Case2 (N=8): (a) displacement: (b) displacement error.
(a) (b)Figure 10
Velocity of DOF2 in Case2 (N=8): (a) velocity; (b) velocity error.
(a) (b)Figure 11
Log of maxerror versus log of time-step: (a) displacement of DOF2; (b) velocity of DOF2.
(a) (b)


### 4.1. Two-DOF Model

The first example considers the following two-DOF system:(28)2001x¨1x¨2+6-2-24x1x2=f1f2The natural periods of the system areT1=4.443,T2=2.810 (some values are rounded). In the following cases, the solution obtained by the mode superposition method is referred to as the reference solution. In Case 1, a constant loading is used, while the applied loading is a sine function in Case 2.Case 1.
The initial conditions are set asx1(0)=1, x2(0)=2, x˙1(0)=0, and x˙2(0)=0 with f1=100,f2=0. In this case, the time-step is selected as τ=0.5236, the recursive parameter N=6, and the time substep Δt=0.0082.Figures2 and 3 show the simulation results obtained using various methods. The error is defined as follows:(29)error=x^it-xiti=1,2,where x^i(t) is the response obtained using numerical methods and xi(t) is the reference solution. The ExGA-RK method employs the RK-4 to calculate the fundamental solution matrix and Simpson’s 1/3 rule to obtain the particular solution. The SPIM applies the generalized Padé approximation (q=2, μ=0.8) to calculate the fundamental solution matrix and the repeated Simpson integration formulae (using 8 integration points) to obtain the particular solution. In this case, the applied loading is a constant vector and the integrand of the particular solution remains unchanged. As a result, all methods perform very well. The error follows the sequence of MPIM- - Bathe > SPIM >ExGA-RK=MPIM-RK.The accuracy of numerical fundamental solution matrices obtained by different methods is also investigated by fixing the time-stepτ=0.262 and choosing different recursive parameters N. The error of T~(τ) is defined as follows:(30)errorofT~=T~τ-Tτ2,where T~(τ) is the numerical fundamental solution matrices obtained using different methods and T(τ) is the reference fundamental solution matrix. Figure 4 depicts the error versus the recursive parameters N. It can be seen that the error decreases with the increasing of N. With the same N, the error follows the sequence of MPIM- Newmark >MPIM- Bathe >SPIM>MPIM-RK.Case 2.
The initial conditions are set asx1(0)=1, x2(0)=2, x˙1(0)=0, and x˙2(0)=0 and the load f1=100sin⁡(3t),f2=0. In this case, the time substep is set as Δt=0.0082. The recursive parameter is chosen as N=4,6,8. Thus τ=0.131,0.524,2.094 and τ/Tp=1/16,1/4,1, where Tp=2.094 is the period of the applied loading.Figures5–10 show the simulation results using different N. It can be seen that the errors of the ExGA-RK, SPIM solution increase as the time-step increases, while the solutions obtained by the proposed methods (MPIM-Bathe, MPIM-RK) is always the same as the reference solution. For instance, for N=6,τ=Tp/2, the loading is zero at the start and end points of each time-step; thus it would be considered as zero within each time-step. Although the SPIM uses 8 integration points to calculate the particular, the numerical results have large errors when N=6,8.In this case, the numerical convergence property of the proposed method is also investigated by calculating the time order of convergence rate. The time order of convergence rate is calculated by applying the following formula:(31)s=log⁡maxerror1-log⁡maxerror2log⁡τ1-log⁡τ2,where τ1, τ2 are different time-steps, and maxerror1, maxerror2 are the maximum error of corresponding numerical result. Figure 11 shows the log of the maxerror versus the log of time-step for displacement and velocity of node 2. It can be seen that the slopes of SPIM-RK and SPIM-Bathe curves remain unchanged as N increases. The slopes of SPIM-RK and SPIM-Bathe curves are about 3.99 and 1.97, respectively.


### 4.2. Regular Diagonally Positioned Pyramid

The second example considers the regular diagonally positioned pyramid shown in Figure12, which is composed of aluminium tubes. The section dimensions and the material properties are as follows: the outer diameter 150mm, the wall thickness 11mm, the density ρ=2.699×103kg/m3, the Poisson ratio ν=0.32, and Young’s modulus E=6.895×104MPa. The top level of the pyramid is fixed along four sides, while node 1 is subjected to a vertical load F=5×103sin⁡(1000t)kN. By using spatial discretization (spatial bar element), a structural model with 543 degrees of freedom is obtained. A consistent mass matrix is used in the analysis and the material is assumed to be linearly elastic. The damping matrix is assumed to be C=a0K+b0M (Rayleigh damping), where the coefficients a0=0.001 and b0=0.45 [3]. The maximum natural period of the model is Tmax=0.0565s.Figure 12
Regular diagonally positioned pyramid.The time-step is set asτ=0.0029s, the recursive parameter N=5, and the time substep Δt=9.062×10-5s. Thus τ/Tmax=0.0513, τ/Tp=0.462, where Tp=6.283×10-3s is the period of the loading. The solution obtained by the Bathe method with time-step 9.062×10-5s is set as the reference solution. The model was solved for over 10 seconds. Table 1 shows the CPU time occupied by different methods (using a desktop PC with Inter(R) Core(TM) i5-5200U CPU @ 2.2 GHz, Ram@16G). The time-step for the RK-4 and Bathe method is Δt=9.062×10-5s. It can be seen that the proposed method requires much less CPU time than the original method with the small time-step.Table 1
CPU time (s) occupied by different methods.
MPIM-RK4 RK4 MPIM-Bathe Bathe N=5 13.101 373.608 9.969 149.257 N=10 8.961 5.534Figure13 shows the simulation results of the vertical displacement of node 1 and node 2, respectively. Similar conclusions can be obtained for other nodes. It demonstrates that the MPIM-Bathe and MPIM-RK solution exactly follows the reference solution. As the time-step is set as nearly half the period of the load, larger error occurs in the ExGA-RK and SPIM solutions.Figure 13
Simulation results for the vertical displacement: (a) node 1; (b)node 2.
(a) (b)From the foregoing discussion, it can be concluded that the accuracy of the proposed method is determined by the size of time substep and the properties of the original method. Even with a large time-step, the particular solution has the same accuracy as the homogeneous solution. On the premise of obtaining almost the same numerical result, the proposed method can greatly improve the calculative efficiency of the original method.


## 4.1. Two-DOF Model

The first example considers the following two-DOF system:(28)2001x¨1x¨2+6-2-24x1x2=f1f2The natural periods of the system areT1=4.443,T2=2.810 (some values are rounded). In the following cases, the solution obtained by the mode superposition method is referred to as the reference solution. In Case 1, a constant loading is used, while the applied loading is a sine function in Case 2.Case 1.
The initial conditions are set asx1(0)=1, x2(0)=2, x˙1(0)=0, and x˙2(0)=0 with f1=100,f2=0. In this case, the time-step is selected as τ=0.5236, the recursive parameter N=6, and the time substep Δt=0.0082.Figures2 and 3 show the simulation results obtained using various methods. The error is defined as follows:(29)error=x^it-xiti=1,2,where x^i(t) is the response obtained using numerical methods and xi(t) is the reference solution. The ExGA-RK method employs the RK-4 to calculate the fundamental solution matrix and Simpson’s 1/3 rule to obtain the particular solution. The SPIM applies the generalized Padé approximation (q=2, μ=0.8) to calculate the fundamental solution matrix and the repeated Simpson integration formulae (using 8 integration points) to obtain the particular solution. In this case, the applied loading is a constant vector and the integrand of the particular solution remains unchanged. As a result, all methods perform very well. The error follows the sequence of MPIM- - Bathe > SPIM >ExGA-RK=MPIM-RK.The accuracy of numerical fundamental solution matrices obtained by different methods is also investigated by fixing the time-stepτ=0.262 and choosing different recursive parameters N. The error of T~(τ) is defined as follows:(30)errorofT~=T~τ-Tτ2,where T~(τ) is the numerical fundamental solution matrices obtained using different methods and T(τ) is the reference fundamental solution matrix. Figure 4 depicts the error versus the recursive parameters N. It can be seen that the error decreases with the increasing of N. With the same N, the error follows the sequence of MPIM- Newmark >MPIM- Bathe >SPIM>MPIM-RK.Case 2.
The initial conditions are set asx1(0)=1, x2(0)=2, x˙1(0)=0, and x˙2(0)=0 and the load f1=100sin⁡(3t),f2=0. In this case, the time substep is set as Δt=0.0082. The recursive parameter is chosen as N=4,6,8. Thus τ=0.131,0.524,2.094 and τ/Tp=1/16,1/4,1, where Tp=2.094 is the period of the applied loading.Figures5–10 show the simulation results using different N. It can be seen that the errors of the ExGA-RK, SPIM solution increase as the time-step increases, while the solutions obtained by the proposed methods (MPIM-Bathe, MPIM-RK) is always the same as the reference solution. For instance, for N=6,τ=Tp/2, the loading is zero at the start and end points of each time-step; thus it would be considered as zero within each time-step. Although the SPIM uses 8 integration points to calculate the particular, the numerical results have large errors when N=6,8.In this case, the numerical convergence property of the proposed method is also investigated by calculating the time order of convergence rate. The time order of convergence rate is calculated by applying the following formula:(31)s=log⁡maxerror1-log⁡maxerror2log⁡τ1-log⁡τ2,where τ1, τ2 are different time-steps, and maxerror1, maxerror2 are the maximum error of corresponding numerical result. Figure 11 shows the log of the maxerror versus the log of time-step for displacement and velocity of node 2. It can be seen that the slopes of SPIM-RK and SPIM-Bathe curves remain unchanged as N increases. The slopes of SPIM-RK and SPIM-Bathe curves are about 3.99 and 1.97, respectively.


## 4.2. Regular Diagonally Positioned Pyramid

The second example considers the regular diagonally positioned pyramid shown in Figure12, which is composed of aluminium tubes. The section dimensions and the material properties are as follows: the outer diameter 150mm, the wall thickness 11mm, the density ρ=2.699×103kg/m3, the Poisson ratio ν=0.32, and Young’s modulus E=6.895×104MPa. The top level of the pyramid is fixed along four sides, while node 1 is subjected to a vertical load F=5×103sin⁡(1000t)kN. By using spatial discretization (spatial bar element), a structural model with 543 degrees of freedom is obtained. A consistent mass matrix is used in the analysis and the material is assumed to be linearly elastic. The damping matrix is assumed to be C=a0K+b0M (Rayleigh damping), where the coefficients a0=0.001 and b0=0.45 [3]. The maximum natural period of the model is Tmax=0.0565s.Figure 12
Regular diagonally positioned pyramid.The time-step is set asτ=0.0029s, the recursive parameter N=5, and the time substep Δt=9.062×10-5s. Thus τ/Tmax=0.0513, τ/Tp=0.462, where Tp=6.283×10-3s is the period of the loading. The solution obtained by the Bathe method with time-step 9.062×10-5s is set as the reference solution. The model was solved for over 10 seconds. Table 1 shows the CPU time occupied by different methods (using a desktop PC with Inter(R) Core(TM) i5-5200U CPU @ 2.2 GHz, Ram@16G). The time-step for the RK-4 and Bathe method is Δt=9.062×10-5s. It can be seen that the proposed method requires much less CPU time than the original method with the small time-step.Table 1
CPU time (s) occupied by different methods.
MPIM-RK4 RK4 MPIM-Bathe Bathe N=5 13.101 373.608 9.969 149.257 N=10 8.961 5.534Figure13 shows the simulation results of the vertical displacement of node 1 and node 2, respectively. Similar conclusions can be obtained for other nodes. It demonstrates that the MPIM-Bathe and MPIM-RK solution exactly follows the reference solution. As the time-step is set as nearly half the period of the load, larger error occurs in the ExGA-RK and SPIM solutions.Figure 13
Simulation results for the vertical displacement: (a) node 1; (b)node 2.
(a) (b)From the foregoing discussion, it can be concluded that the accuracy of the proposed method is determined by the size of time substep and the properties of the original method. Even with a large time-step, the particular solution has the same accuracy as the homogeneous solution. On the premise of obtaining almost the same numerical result, the proposed method can greatly improve the calculative efficiency of the original method.


## 5. Discussion

In recent years, many integration methods have been applied to calculate Green’s function matrix, resulting in different versions of time marching method, such as the Newmark family methods [17], the generalized-α method [18], the Runge-Kutta, and the central difference schemes [16]. Using a substepping procedure, these aforementioned approaches can obtain accurate homogeneous solution. However, the technologies used to evaluate the particular solution are basically based on some approximations that are only suitable for a small time-step. For instance, in [16] the external loading vector is assumed to be linear functions in each time-step. In ‎[19] the integrand of the particular solution is assumed to be a linear or parabolic function within each time-step. When a larger time-step is considered, these approximations are oversimplified and results in unexpected errors. For instance, the time-step is set as an integral multiple of the period of the periodical loading; the loading would be considered to be unchanged within each integral interval. As a result, the accuracy of the particular solution obtained would be inadequate compared to that of the homogeneous solution.It is widely known that the computational effort required is inversely proportional to the size of time-step in the structural dynamic analysis. Compared to these aforementioned approaches, the proposed method can use a large time-step and can preserve the same accuracy of particular solution as that of the homogeneous solution; thus it is more suitable for the long-time duration dynamic analysis.


## 6. Conclusions

An improved version of PIM is presented in this paper. This approach replaces the Taylor series expansion or the generalized Padé approximation of the matrix exponential function in the conventional PIMs with numerical fundamental solution matrix and employs loading operator matrix to obtain particular solutions. Efficient recurrence expressions are developed to evaluate the loading operator matrices in the first time-step. The highlights of the approach can be summarized as follows,(1)
The proposed method avoids matrix inversion operations and numerical quadrature formulae; thus it expands the scope of application. Meanwhile, the accuracy of particular solution obtained can be maintained the same as that of the homogeneous solution even for a large time-step.(2)
As the intermediate variables required for evaluation of the loading operator matrix have already been generated in calculating the fundamental solution matrix, the proposed method improves the computational efficiency of conventional PIMs.(3)
By selecting an appropriate original method and time substep, the satisfactory numerical properties, such as accuracy, numerical dissipation, and unconditional stability, can be obtained.(4)
Although the proposed method uses a large time-step, it can obtain almost the same numerical results as the original method using a small time-step. Thus, the proposed method is actually an accelerator of conventional time integration methods.



---

*Source: 9646017-2018-06-14.xml*