A Time Finite Element Method Based on the Differential Quadrature Rule and Hamilton ’ s Variational Principle

An accurate and efficient Differential Quadrature Time Finite Element Method (DQTFEM) was proposed in this paper to solve structural dynamic ordinary differential equations. This DQTFEM was developed based on the differential quadrature rule, the Gauss–Lobatto quadrature rule, and the Hamilton variational principle. The proposed DQTFEM has significant benefits including the high accuracy of differential quadrature method and the generality of standard finite element formulation, and it is also a highly accurate symplectic method. Theoretical studies demonstrate the DQTFEM has higher-order accuracy, adequate stability, and symplectic characteristics. Moreover, the initial conditions in DQTFEM can be readily imposed by a method similar to the standard finite element method. Numerical comparisons for accuracy and efficiency among the explicit Runge–Kutta method, the Newmark method, and the proposed DQTFEM show that the results from DQTFEM, even with a small number of sampling points, agree better with the exact solutions and validate the theoretical conclusions.


Introduction
Various algorithms are available for numerical integration of Ordinary Differential Equations (ODEs).The most commonly used algorithms, derived from Finite Difference Method (FDM), are the Direct Integration Methods (DIMs), such as the Central Difference Method (CDM), the Newmark method, and the Runge-Kutta Method (RKM).Feng presented structure-preserving symplectic schemes for the Hamilton system and illustrated the fundamental cause of energy dissipation of classical DIM [1].
As is widely accepted, Finite Element Method (FEM) is superior to FDM and has advantages in space, thus it should also exhibit advantages in time in the mathematical sense.However, Zienkiewicz doubted if it is possible to use Finite Elements (FEs) in the time domain [2].Indeed, when the involved functions are smooth enough, it is not essential to define massive time steps to integrate ODEs, and the time steps should be as many as the FEs to produce a solution of comparable accuracy.In view of the increased storage required, the use of Time Finite Elements (TFEs) to solve ODEs is questionable.However, this problem can be gradually solved due to the rapid development of computer capacity and advanced solution methods.
In many cases, conventional step-by-step algorithms like CDM may call for a very large number of time steps.This is especially true when dealing with excitation and material properties changing rapidly in time.Comparing with the Finite Difference (FD) based method, a variation-based method, with its inherent stability and physical origin, may lower the computational effort considerably.Thus, despite some reservations [2], the extension of FE concept to the time domain is logical.
Many researchers have attempted to employ Time Finite Elements Methods (TFEMs) to discretize ODEs.Similar to FEMs in spatial domain, one can formulate TFEM based on Hamilton's variational principles [3][4][5], Hamilton's law of varying action [6,7], Gurtin's variational principle [8], and weighted residual methods [9,10].Almost all of the reported works dealt with one FE interval and repeated the calculations for subsequent time intervals with new initial conditions.Hence, these TFEMs lead to a step-by-step or recurrence calculation and are conditionally stable, unlike standard FEMs, which produce simultaneous solutions of the whole range of interest.The present work introduces a high-order method, and, generally, the dynamic responses of the whole time domain can be solved in one shot rather than step by step.
The commonly used trial functions in TFEMs based on variational principles are Hermite interpolations [3,4,6,11], rarely used linear functions [12], and cubic spline functions [13][14][15].The advantage of using Hermite interpolations is that the initial conditions can be imposed in the same manner as the boundary conditions in standard FEMs, and the resulting time element is C 1 continuous [16,17].Moreover, the solution accuracy increases as the number of elements increases or the order of trial functions increases.The present TFEM is also based on a variational principle, but usually the nodes or sampling points in an element are not uniformly located, and nodal velocities are only for two end nodes rather than for the interior nodes.
As an alternative to variational principles, weighted residual approaches are also widely used to construct TFEs.The commonly used trial functions in weighted residual approaches are Lagrange interpolations [18,19], linear/quadratic polynomials [20,21], Taylor series [22,23], and Taylor formula [9,24].Zienkiewicz proposed a TFEM based on a weighted residual approach, and gave recurrence schemes of a single step or multiple steps [2].Remarkably, as they reported, the mainly used DIMs can also be derived using weighted residual approaches.Later on, they employed Taylor's expansion to form a unified set of single-step algorithms [9], discussed the accuracy and the stability, and applied the method of backward error analysis to obtain a few general rules to choose parameters.By separately approximating acceleration, velocity and displacement as linear, quadratic, and cubic Taylor series, Hoff and Pahl developed a single-step TFEM based on a weighted residual method, and applied it to linear [25] and nonlinear [22] structural dynamics.Numerical results showed that the method presents some advantages, especially in accuracy and overshooting behavior.Additionally, Krishna and Manjeet presented three and four steps of TFEM using the least squares method as well as Lagrange interpolations for first-order initial-value problems [18].The drawback of the weighted residual approaches based on TFEM is the determination of the optimum weighting functions.
TFEMs were based on the assumption that state variables are continuous in the time domain.However, in the Time-Discontinuous Galerkin Method (TDGM) which was originated by Reed and Hill [26] and Lesaint and Raviart [27] based on Galerkin weighted residual approach, and in the TDGM, the unknown variables are continuous within the time element but discontinuous at time nodes.The TDGM has been employed to solve periodic differential equations and Davey-Stewartson equations [28,29].It was first shown by Lesaint and Raviart [27] and Delfour et al. [30] that TDGM leads to A-stable, higher-order accurate TFEMs.Subsequently, Hughes and Hulbert first applied TDGM to elastodynamics and structural dynamics to develop space-time finite element methods [31][32][33][34] and TFEMs [5,35], respectively.According to their results, TDGM stability is easily proved, but the convergence had been proved only for linear elements, therefore the Time-Discontinuous Galerkin/Least Squares Method (TDGLSM) [36] was developed to solve this problem [31,32,35].Compared with the continuous Galerkin FEM, the drawback of TDGMs is the usage of a higher number of degrees, which has been overcome by designing an improved predictor/multi-corrector algorithm [37] or a suitable adaptive strategy for iterative mesh refinement [38,39].According to the above review, the TFEMs in existence certainly have limitations.As is well known, researchers are always seeking new numerical methods with higher accuracy, efficiency, and flexibility to promote modeling capabilities.Bellman and Casti proposed the Differential Quadrature Method (DQM) as a simply and highly efficient numerical technique for initial-value problems [40], which was regarded as a global and higher-order FD scheme [17,[41][42][43] and could yield accurate results using a considerably smaller number of sampling points.
The general versatility of DQM has been established in a variety of physical problems.A recent survey refers to publications [44,45], and an exhaustive list of the literature on DQM up to 1996 may be found in the survey paper [46].Although DQM's numerical accuracy and computational efficiency has been well demonstrated, its limitations are also significant [46].In 2009, based on the Differential Quadrature (DQ) rule, the Gauss-Lobatto quadrature rule, and the total potential energy variational principle, a DQ Finite Element Method (DQFEM) was proposed by the author and co-workers for free vibration analysis of thin plates [47], and the authors formulated explicit one-to three-dimensional (1-D to 3-D) element matrices [48].The DQFEM uses Lagrange trial functions for C 0 and C 1 FEs and possesses high accuracy of DQM and the generality of standard FEM.
In this context, similar to DQFEM for spatial domain [47,48], Hamilton's variational principle together with the DQ rule and the Gauss-Lobatto quadrature rule are used to develop a TFEM in the present work, where the initial conditions can be readily imposed as in standard FEM.This novel method is called the DQ Time Finite Element Method (DQTFEM).The novelties of the present work compared with the existing TFEMs can be briefly summarized as follows: (1) The general form of Hamilton's variational principle, the DQ rule, and the Gauss-Lobatto quadrature rule are used to formulate a TFEM for structural dynamic analysis.(2) The Gauss-Lobatto quadrature points are used as time sampling points or nodes of the time element, which results in more accurate results with the same number of nodes.(3) The nodal velocities are only involved at two end nodes of the time element, implying that the number of Degrees of Freedom (DOFs) is less than that of the existing TFEMs based on the variational principle and Hermite interpolations.(4) Lagrange interpolations are employed instead of Hermite interpolations in the construction of the C 1 time element, and the stability and accuracy of DQTFEM are studied in detail.
The outline of the present paper is as follows.In Section 2, the formulation of DQTFEM is presented for structural dynamic ODEs.In Section 3, the stability and accuracy are studied.Obtained solutions are compared with analytical solutions, the results of four-order explicit RKM, and the Newmark Average acceleration Method (NMA) in Section 4. Finally, conclusions are drawn in Section 5.

Basic Formulations
For structural dynamic systems with damping and subjected to external loads, the ODE in a discrete form can be derived by using the general form of Hamilton's variational principle [49]: where T, V, and δW stand for the kinetic energy, strain energy, and virtual work of damping forces and external forces, respectively; M = [M ij ] n×n , K = [K ij ] n×n , and C = [C ij ] n×n are the mass, stiffness, and damping symmetrical matrices, respectively; n is the number of spatial DOFs of the system; u(t) and f (t) are the displacement vector and the prescribed load vector corresponding to spatial nodes, respectively; and a time derivative is denoted by a dot over the letter.
Using the variational principle, Equation (1) yields the ODE of discretized system as For the clarity of derivations, the time evolutions of a spatial node "i" are written as for the time sampling points 1 to N, and N is the number of time sampling points.Note that u ij and f ij are the values of u i and f i at time sampling point t j .For a Multiple DOF dynamic system (Equation ( 5)), each spatial DOF is the function of time; its derivative with respect to time can be given through DQ rule as .
where A = [A jk ] N×N is the explicit weighting coefficients matrix [46].It is noteworthy that the weighting coefficients are independent of variables, and dependent only on coordinates of time sampling points.
For the development of the present TFEM, all terms of Equation ( 1) should be dealt with by the DQ and Gauss-Lobatto quadrature rules, that is, .u is expressed in terms of u based on the DQ rule and the integrals are fulfilled by using the Gauss-Lobatto quadrature rule; then we have where where 'cs' indicates column expansion of a matrix, '⊗' is the Kronecker product, E is a n × n identity matrix, and C j is the coefficient of Gauss-Lobatto quadrature.Diagonal matrix C G = diag(C j ).Appendix A presents the derivations of Equations ( 8) to (10).Equation (13) shows that the expanded mass matrix M and stiffness K are symmetrical, but C is not, even though C is symmetrical.Substituting Equations ( 8) to (10) into Equation ( 1), and noting that δU is arbitrary, we have the dynamic equilibrium equation in algebraic form as or in a compact and standard form as where Equation ( 16) can also be written in a different form as which can be reduced to a Sylvester equation whose solutions can be found by the method of Bartels and Stewart [51] if we neglect damping.Then it has no row or column expansion of matrix in Equation (18), or it has no additional storage needed except for A and the diagonal matrix C G when compared with a step-by-step method.
After deriving Equation ( 15), a DQTFEM has been formulated, which is different from classical TFEMs where Hermite interpolations are used and displacement and velocity are involved for all nodes.From the above derivation, one can see that, in DQTFEM, displacements are the only nodal parameter at all time sampling points since Lagrange interpolations are used; refer to Equations ( 11) and (12).In order to impose initial conditions and satisfy C 1 -compatible conditions in the present method, the nodal parameter vector U should be modified to include velocities for all spatial DOFs at the first time sampling point.To achieve this, an accurate and convenient method is to modify U as which is formulated by substituting the displacements u 1N , where T is formed by substituting the first row of A for the Nth row of an unit matrix, that is From Equations ( 16) and ( 20) one can obtain It should be noted that the solution U of Equation ( 23) does not include displacements and velocities of the Nth sampling point, referred to in Equation (19), but these can be easily obtained through Equation ( 20) and the DQ rule, respectively.Rearranging the terms of Equation ( 23), it is rewritten as where subscript 'i' indicates the initial conditions, and 'd' indicates unknowns.From Equation (25) we have where U i is formulated by initial displacements and velocities, thus the initial conditions are dealt with in a general and convenient way.
It is worth noting that here u 1 , u 2 , • • • , u n , corresponding to all time sampling points, are solved concurrently by using the exact initial conditions; this is different from classical step-by-step DIMs, in which the results of the present time step are the initial conditions of the next time step.That is, only the initial conditions of the first time step are exact; all others are approximate.However, DQTFEM can also serve as a step-by-step method; in this situation, a time element is a time step in which there are N time sampling points.Therefore, the initial conditions except for the first time step are also approximate, but the approximate initial conditions are highly accurate since the present method is a high-order method. Remarks: (1) Usually for conservative systems, the actual path will be the one for which the integral of Lagrangian function (T − V) is an extremum, and for the non-conservative system considered here, the general form of Hamilton's variational principle is used and the integral still exists; however, there exists no functional that must be an extremum [49].(2) The method of imposing initial conditions in the present method is standard and similar to the one used in the FEM; see Equation ( 25). ( 3) DQTFEM is a new TFEM, since Gauss-Lobatto quadrature and DQ rules are used for the first time in the construction of a TFEM.The high accuracy and efficiency due to the employment of Gauss-Lobatto quadrature and DQ rules are validated in the following numerical comparisons; see Section 4.

Stability Analysis
Stability analysis is significant for a new TFEM, and for this analysis a spectral approach is generally employed because a linear dynamic system can be decoupled into a set of single DOF systems based on the orthogonality of mode functions or mode vectors.Therefore, in the stability analysis of the present method, consider a single DOF system as where ω denotes the angular frequency.Here, the two dots over u denote the second derivative with respect to time τ, not to t, where τ = t/h, and h and t stand for time step size and the real time, respectively.For Equation ( 27), the displacement vector u in Equation ( 6) becomes where u j (j = 1, 2, . . .N) denotes the displacements at time sampling points, the subscript 'i' in Equation ( 6) denotes the spatial DOF number and is omitted here since i = 1, and The nodal parameter vector in Equation ( 19) has the form of According to Equation (26), one can obtain where column vectors α and β are extracted from G, and u 0 and .u 0 are the initial displacement and velocity, respectively.Moreover, using DQ rule yields .
From Equations ( 31) and ( 32), we have which is the recurrence relation between (u N , . u N ) and (u 0 , .u 0 ), and the elements of the Jacobi matrix According to the relationship between ωh and spectral radius ρ(J), see Figure 1, the first stable intervals for different N can be obtained as given in Table 1.From Table 1 and Figure 1 we can make the following observations: (1) ρ(J) ≥ 1 for any N, and when N becomes larger, the first stable interval also becomes larger.
Although there are several stable intervals with ρ(J) = 1 except for N = 3, the first stable interval is significant and practical; therefore, only the first intervals are listed in Table 1 whose elements are obtained by using the criterion (ρ − 1) < 10 −6 .(2) ρ(J) < 1 in the TFEM by Hulbert [35], but ρ(J) = 1 in DQTFEM for all stable intervals, hence DQTFEM is more accurate than that of Hulbert [35] since DQTFEM is a high-order method that has a smaller phase error than lower order methods.Note that the stability of CDM is ωh ≤ 2, and ωh ≤ 2 √ 2 for four-order explicit RKM, and the stability of the present method is ωh ≤ 2 √ 2.
(3) DQTFEM is highly symplectic.If the Jacobi matrix J of a method satisfies Equation (35) as follows then the method is symplectic.For the present single DOF system, all variables in Equation ( 35) are numbers and Table 2 shows that the precision of Equation ( 35) (especially J 3 ~J6 ) is satisfactory, which indicates that DQTFEM is a symplectic method.In fact, we have numerically verified that DQTFEM is symplectic even if it is not stable.For a stable symplectic method, the magnitudes of all eigenvalues of J are equal to 1, therefore all eigenvalues of J of the present method are equal to 1 when it is stable.
DQTFEM is symplectic even if it is not stable.For a stable symplectic method, the magnitudes of all eigenvalues of J are equal to 1, therefore all eigenvalues of J of the present method are equal to 1 when it is stable.Note: max(ωh) indicates the upper limit of ωh for the first stable interval in Table 1.

Numerical Results and Discussion
This section aims at demonstrating the high accuracy and efficiency of DQTFEM by comparing them with RKM and NMA; note that NMA is equivalent to the Euler mid-point symplectic method.To evaluate the relative computational accuracy of the present method, the relative error for all methods is defined as where 'Exact' denotes the exact solution at each time sampling point, the infinity norm of which denotes the maximum absolute value of a vector, and the elements of this vector are the exact values at all time sampling points; 'Approximate' is the result of DQTFEM at each time sampling point.Note: max(ωh) indicates the upper limit of ωh for the first stable interval in Table 1.

Numerical Results and Discussion
This section aims at demonstrating the high accuracy and efficiency of DQTFEM by comparing them with RKM and NMA; note that NMA is equivalent to the Euler mid-point symplectic method.To evaluate the relative computational accuracy of the present method, the relative error for all methods is defined as where 'Exact' denotes the exact solution at each time sampling point, the infinity norm of which denotes the maximum absolute value of a vector, and the elements of this vector are the exact values at all time sampling points; 'Approximate' is the result of DQTFEM at each time sampling point.

Free Vibration of a Two-DOF System
Consider the following two-DOF system, its governing equilibrium equation is It follows from Figures 2 and 3 that DQTFEM is much more accurate than RKM and NMA, RelErr (DQTFEM) is about 10 −6 , while RelErr (RKM) and RelErr (NMA) are 10 −2 and 10 −1 , respectively; and it is noticeable from Figure 3 that the errors of RKM and NMA accumulate over time.
whose angular frequencies It follows from Figures 2 and 3 that DQTFEM is much more accurate than RKM and NMA, RelErr (DQTFEM) is about 10 −6 , while RelErr (RKM) and RelErr (NMA) are 10 −2 and 10 −1 , respectively; and it is noticeable from Figure 3 that the errors of RKM and NMA accumulate over time.whose angular frequencies It follows from Figures 2 and 3 that DQTFEM is much more accurate than RKM and NMA, RelErr (DQTFEM) is about 10 −6 , while RelErr (RKM) and RelErr (NMA) are 10 −2 and 10 −1 , respectively; and it is noticeable from Figure 3 that the errors of RKM and NMA accumulate over time.

Efficiency Comparison
Computational cost is also a crucial issue for a method.Assume the time domain is [0, 2 × 10 3 ] in this case.In order to achieve the same accuracy of RelErr = 10 −6 , h = 5 × 10 −2 and h = 1 × 10 −3 Appl.Sci.2017, 7, 138 10 of 16 must be used for RKM and NMA, respectively; that is, RKM needs 4 × 10 5 time steps and NMA needs 2 × 10 6 time steps.However, we can achieve the same accuracy by using 1588 elements and each element of size h = 1.26 ≈ 2.828/ω 2 has only three Gauss-Lobatto-Legendre sampling points (N = 3 is the least number of time sampling points in an element; see Table 2).If N > 3, the accuracy of DQTFEM is more than 10 −6 .For example, if N = 20 in an element of size h = 11 ≈ 25.11/ω 2 , the accuracy will be 10 −9 .In other words, the lowest accuracy of DQTFEM is 10 −6 .This is the reason why the data corresponding to RelErr = 10 −2 are not provided in Table 3, which indicates that DQTFEM is much more efficient than the other two methods, especially if high accuracy is required.Another observation is that RKM is more efficient than NMA.
NMA is a symplectic method of second order, in which the amplitude can be rigorously preserved for linear systems, but its phase error is larger than that of RKM.It is well known that dynamic responses must be represented by amplitude and phase together; therefore, RKM is more efficient than NMA for achieving the same accuracy.However, the weakness of RKM is its numerical dissipation and the weakness of NMA is its larger phase error.For long-time kinematics simulations in which only amplitude is taken into account, NMA is a better choice.For further validation of the present method, Rayleigh damping given in Equation ( 38) is added to the system as Equation (37).
N = 33 and h = 0.59375 are also used here.The displacement x 1 and RelErr are given in Figures 4  and 5, respectively, from which one can see that DQTFEM is also the most accurate: its RelErr reaches 10 −6 .Therefore, the present method is also applicable to the analysis of damped problems.
size h = 1.26 ≈ 2.828/ω2 has only three Gauss-Lobatto-Legendre sampling points (N = 3 is the least number of time sampling points in an element; see Table 2).If N > 3, the accuracy of DQTFEM is more than 10 −6 .For example, if N = 20 in an element of size h = 11 ≈ 25.11/ω2, the accuracy will be 10 −9 .In other words, the lowest accuracy of DQTFEM is 10 −6 .This is the reason why the data corresponding to RelErr = 10 −2 are not provided in Table 3, which indicates that DQTFEM is much more efficient than the other two methods, especially if high accuracy is required.Another observation is that RKM is more efficient than NMA.
NMA is a symplectic method of second order, in which the amplitude can be rigorously preserved for linear systems, but its phase error is larger than that of RKM.It is well known that dynamic responses must be represented by amplitude and phase together; therefore, RKM is more efficient than NMA for achieving the same accuracy.However, the weakness of RKM is its numerical dissipation and the weakness of NMA is its larger phase error.For long-time kinematics simulations in which only amplitude is taken into account, NMA is a better choice.For further validation of the present method, Rayleigh damping given in Equation ( 38) is added to the system as Equation (37).and 5, respectively, from which one can see that DQTFEM is also the most accurate: its RelErr reaches 10 −6 .Therefore, the present method is also applicable to the analysis of damped problems.

Free Vibration of a Simply Supported Euler Beam
In order to further explore the accuracy and efficiency of DQTFEM, the free vibration of a simply supported Euler beam of length L = 3 m is studied in this section; flexural rigidity EI = 1 × 10 6 N/m and mass density per unit length ρA = 420 kg/m.The beam is meshed into 16 uniform cubic elements All initial velocities and displacements are equal to zero except the initial velocity at mid-span i 1 m/s; apparently the problem is analogous to the impact problem.
Figure 6 shows an excellent agreement of DQTFEM with the exact method.For reaching the same order of accuracy, NMA needs 2.110 × 10 5 time steps and much more CPU time than DQTFEM see Table 4. Figure 7 presents the RelErr comparisons between DQTFEM and NMA; it follows tha the error of NMA accumulates while the error of DQTFEM does not.
For validating the effectiveness of DQTFEM for long-term simulation, the computational time changes from 3.6 × 10 −3 s to three times the maximum period Tmax = 2π/ω1 = 0.1174 s, where ω1 = 53.51rad/s, the number of elements is 3Tmax/h ≈ 978. Figure 8 gives a comparison of displacements by DQTFEM with the exact solution, which shows that DQTFEM is also highly accurate for this case and the needed CPU time is only 75.53 s.Nevertheless, the NMA results are not presented since NMA requires a prohibitively long CPU time.

Free Vibration of a Simply Supported Euler Beam
In order to further explore the accuracy and efficiency of DQTFEM, the free vibration of a simply supported Euler beam of length L = 3 m is studied in this section; flexural rigidity EI = 1 × 10 6 N/m and mass density per unit length ρA = 420 kg/m.The beam is meshed into 16 uniform cubic elements.All initial velocities and displacements are equal to zero except the initial velocity at mid-span is 1 m/s; apparently the problem is analogous to the impact problem.
Figure 6 shows an excellent agreement of DQTFEM with the exact method.For reaching the same order of accuracy, NMA needs 2.110 × 10 5 time steps and much more CPU time than DQTFEM; see Table 4. Figure 7 presents the RelErr comparisons between DQTFEM and NMA; it follows that the error of NMA accumulates while the error of DQTFEM does not.For validating the effectiveness of DQTFEM for long-term simulation, the computational time changes from 3.6 × 10 −3 s to three times the maximum period T max = 2π/ω 1 = 0.1174 s, where ω 1 = 53.51rad/s, the number of elements is 3T max /h ≈ 978. Figure 8 gives a comparison of displacements by DQTFEM with the exact solution, which shows that DQTFEM is also highly accurate for this case and the needed CPU time is only 75.53 s.Nevertheless, the NMA results are not presented since NMA requires a prohibitively long CPU time.

Conclusions
A highly accurate and highly efficient DQTFEM was developed in the present study by using

Conclusions
A highly accurate and highly efficient DQTFEM was developed in the present study by using the DQ rule, the Gauss-Lobatto integration rule, and the general form of Hamilton's variational principle for solving structural dynamic problems involving damping and external forces.
DQTFEM incorporates the high accuracy and efficiency of DQM into the development of high-order time elements, in which initial conditions can be introduced in the same simple and general form as in standard FEM.Furthermore, DQTFEM can serve as a step-by-step approach; one time element of DQTFEM can be seen as one time step of the conventional step-by-step method.Moreover, it has been shown that DQTFEM is a highly accurate symplectic method.
Extensive numerical comparisons validated the present method, whose solutions can be seen as benchmarks to validate lower-order DIMs and TFEMs.It is expected that DQTFEM will be practical in the future.

) whose angular frequencies ω 1 = √ 2 and ω 2 = √ 5 .
The given initial conditions are u T 0 = [ 1 1 ] and .Let time domain be [0, 19], which is taken as an element with N = 33 non-uniform points in DQTFEM, and 32 uniform time steps are used for both RKM and NMA.The first stable interval for N = 33 is [0, 43.97], hence the stability is satisfied since 19ω 2 = 42.49< 43.97, and the stability 2 √ 2 of RKM is also satisfied by h = 19/32 = 0.59375 since ω 2 h = 1.3277 < 2 √ 2. Additionally, N = 33 and h = 0.59375 implies that the three methods have the same number of interpolation points in the domain [0, 19].The exact solutions used for comparisons are obtained by the mode superposition method.
Let time domain be [0, 19], which is taken as an element with N = 33 non-uniform points in DQTFEM, and 32 uniform time steps are used for both RKM and NMA.The first stable interval for N = 33 is [0, 43.97], hence the stability is satisfied since 19ω2 = 42.49< 43.97, and the stability 2√2 of RKM is also satisfied by h = 19/32 = 0.59375 since ω2h = 1.3277 < 2√2.Additionally, N = 33 and h = 0.59375 implies that the three methods have the same number of interpolation points in the domain [0, 19].The exact solutions used for comparisons are obtained by the mode superposition method.
Let time domain be [0, 19], which is taken as an element with N = 33 non-uniform points in DQTFEM, and 32 uniform time steps are used for both RKM and NMA.The first stable interval for N = 33 is [0, 43.97], hence the stability is satisfied since 19ω2 = 42.49< 43.97, and the stability 2√2 of RKM is also satisfied by h = 19/32 = 0.59375 since ω2h = 1.3277 < 2√2.Additionally, N = 33 and h = 0.59375 implies that the three methods have the same number of interpolation points in the domain [0, 19].The exact solutions used for comparisons are obtained by the mode superposition method.

N
= 33 and h = 0.59375 are also used here.The displacement x1 and RelErr are given in Figures 4