## 1 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). [PERSON] 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, [PERSON] 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 changingrapidly 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 [PERSON]'s variational principles [3; 4; 5], [PERSON]'s law of varying action [6; 7], [PERSON]'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]. [PERSON] 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.