Methods of Approximation in Hpk Framework for Odes in Time Resulting from Decoupling of Space and Time in Ivps

The present study considers mathematical classification of the time differential operators and then applies methods of approximation in time such as Galerkin method (GM), Galerkin method with weak form (/ GM WF), Petrov-Galerkin method (PGM), weighted residual method (WRY), and least squares method or process (LSM or LSP) to construct finite element approximations in time. A correspondence is established between these integral forms and the elements of the calculus of variations: 1) to determine which methods of approximation yield unconditionally stable (variationally consistent integral forms, VC) computational processes for which types of operators and, 2) to establish which integral forms do not yield unconditionally stable computations (variationally inconsistent integral forms, VIC). It is shown that varia-tionally consistent time integral forms in hpk framework yield computational processes for ODEs in time that are unconditionally stable, provide a mechanism of higher order global differentiability approximations as well as higher degree local approximations in time, provide control over approximation error when used as a time marching process and can indeed yield time accurate solutions of the evolution. Numerical studies are presented using standard model problems from the literature and the results are compared with Wilson's  method as well as Newmark method to demonstrate highly meritorious features of the proposed methodology.


Introduction, Literature Review and Scope of Work
The quest for better, reliable, stable, accurate and efficient methods of approximation for obtaining numerical solutions of the ordinary differential equations ( ODEs ) in time resulting from either decoupling of space and time or using a lumped parameter approximation in spatial coordinates in initial value problems has been a subject of writing in the numerical and computational methods area for over a century.Thus, the published works are voluminous.In this section we provide a summary of historical perspective of the origins of various methods, discuss more recent advances in the various approaches, discuss their merits and shortcomings with the objective of ultimately narrowing down to a single methodology that provides us an infrastructure to treat all ODEs in time in a consistent but rigorous manner without any regard to specific applications or their origin.

Space-Time Decoupling in IVPs
We make a few remarks regarding a strategy for spacetime decoupling in initial value problems ( IVPs ) that lead to ODEs in time.This is helpful in gaining a better understanding regarding the origin of the ODEs .The IVPs are generally a system of partial differential equations in dependent variables, spatial coordinates and time.Thus, these naturally contain spatial as well as time derivatives of the dependent variables.In decoupling of space and time, the spatial derivatives are converted into algebraic expressions, thereby reducing the original partial differential equations ( PDEs ) to a system of ODEs in time.For this purpose we can use finite difference, finite volume or finite element method in spatial coordinates.In a finite difference method one discretizes the spatial domain into a finite number of points.The spatial derivatives are approximated using Taylor series expansions of various orders of truncation errors to obtain algebraic expressions or the finite difference approximations of the spatial derivatives.One then considers the PDEs at each point of the discretization and substitutes the finite difference expressions for the spatial derivatives.The outcome of this process is that spatial derivatives are replaced by dependent variable values (and/or their derivatives) at the points of the discretization and the time derivatives assume values at the points of the discretization and thus we obtain a system of ODEs in time.In principle, the finite volume method of decoupling space and time is similar.In the finite element method of decoupling space and time, one constructs a spatial discretization of the spatial domain using finite elements.We can construct an integral form over the discretized spatial domain based on the fundamental lemma [1][2][3][4][5].The time derivatives for this purpose are treated as constants.One constructs an approximation over the discretized domain and thereby an approximation over each spatial element of the discretization.This approximation, when substituted in the integral form, replaces the spatial derivatives with algebraic expressions and what remains is a system of ODEs in time.
In the following, we present some details of the finite element method of decoupling space and time.Let ( , ) ( , ) 0, ( , ) (0, ) (0, ) be an initial value problem with some boundary conditions ( BCs ) and some initial conditions ( ICs ).
Let T e x x e     be a finite element discretization of x  .Let ( , ) h x t  be an approximation of ( , ) x t  Then based on the fundamental lemma [1][2][3][4][5]  over an ele ment e with domain e x  .At this point if the operator A has even order derivatives with respect to spatial coordinates, then one could use integration by parts to transfer half of the differentiation from e h  to v (Galerkin method with weak form in space) and then the resulting expression could be arranged as   e P is a vector of secondary variables.Thus, , This paper considers the finite element method of approximation for a system of ODEs in time (like (1.10)) resulting from space-time decoupling in IVPs .

Literature Review
The published work addressing the methods of approximation for ODEs in time is heavily dominated by finite difference methods or methods like finite difference methods that originate from the use of Taylor series expansions at discrete points in the time domain to approximate time derivatives in the ODEs by algebraic expressions containing function values and/or their derivatives at the discrete points in the time domain [0, ] x    .In the following we discuss some of the commonly used approaches.
Reference [6] presents various methods of approximation (primarily based on the Taylor series expansions) for ODEs in time.In reference [7] the authors present an account of the various methodologies for approximating solutions of ODEs in time.These include details for first order symmetric systems, second order symmetric systems, as well as nonlinear symmetric and nonsymmetric systems of ODEs .One step algorithms, linear multistep (  )   LMS algorithms, Houbolt method,  -method and partitioning methods are considered as methods of approximation.Some discussion on convergence and stability is also presented.Perhaps the most significant work amongst the earlier works is due to Newmark [8] based on the assumption of linear acceleration in a time interval.Another parallel work by Wilson [9], Wilson's  method, is also as significant.This is also based on the assumption of linear acceleration in a time interval.These two works were strictly aimed at ODEs from decoupling space and time in structural dynamics (second order ODEs in time).In the Newmark method the equilibrium is considered at time t t   at which the evolution state is not known.In Wilson's  method the linear acceleration assumption is extended to the time interval [ ,  ] The value of  is determined from the stability analysis of the method.It is shown that for 1.37   (generally chosen 1.4), the Wilson's  method is unconditionally stable.Many works using the two approaches have appeared in the literature in which accuracy (amplitude decay, base elongation and phase shift of known and specified waves) and applications of these methods have been discussed.A modification of Newmark method referred to as alpha modification is presented in reference [10].The additional parameter alpha is introduced to achieve simultaneous second order accuracy, unconditional stability and positive artificial damping.A family of "one step unconditionally stable integration schemes with improved numerical dissipation" is presented in reference [11].The "generalized  -method in structural dynamics" is presented in reference [12] with "user controlled numerical dissipation".Extensions of LMS methods to L -stable optimal integration methods with 0 u -0 v overshoot properties in structural dynamics are given in reference [13].This approach leads to level-symmetric ( LS ) integration methods.The LS methods are created as symmetric variants of extended three level LMS methods ( 3L -LMS ) with direct use of dynamic equations to obtain algorithmically simple integration schemes that have maximum damping of high frequency modes and thereby without overshoots.The analysis of generalized  -method for non-linear dynamic problems is presented in reference [14].It is shown that generalized α-methods have stability in the energy sense and guaranteed energy annihilation.But, these methods exhibit overshoots and heavy energy oscillations in intermediate frequency ranges.The energy-conserving and decaying algorithms in structural dynamics based on  -method are considered in reference [15].In a series of papers [16][17][18][19][20][21][22], exhaustive developments and details are presented related to time integration of ODEs .
A detailed discussion of various aspects, their merits and shortcomings are not relevant in view of the approach presented in this paper and is also too voluminous to be included here.However, we do remark that various approaches in these papers utilize weighted residual concepts, Taylor series expansions, Taylor series expansion in time and / GM WF and consider single step and multistep algorithms.Stability, accuracy and other features of the various proposed schemes are presented including numerical studies for model problems.

Scope of Present Work and Methodology
We want to take a very pragmatic view i.e., we pose a question "What is our objective?"The answer of course is, we want to consider a method of approximation for ODEs in time that has the following features.
1) The methodology must be equally applicable to all ODEs , i.e., must be independent of the nature of the time operator.
2) The method must be unconditionally stable regardless of the choice of integration time step.
3) The method must be time marching so that the computations of the evolution can be done for an incre-ment of time.This is essential for computational efficiency in applications requiring the evolution for long periods of time as well as in applications in which a large number of ODEs are involved.
4) The method must permit higher degree approximation of the solution in time for an increment of time.
5) To ensure desired global differentiability of the approximation over t  we must ensure that the time approximations of higher degree for each increment of time are such that their union indeed yields desired global differentiability of the approximation in time.For example if an ODE in time contains up to second order derivatives of the dependent variable and if the theoretical solution is analytic, then the approximation of this solution over t  must at least be of class 2 ( ) t C  .6) There must be a measure (and not estimation) of approximation errors in the computed solution regardless of whether we have a theoretical solution of the ( ) ODE s or not.This feature is essential due to the fact that most problems of practical interest do not permit determination of a theoretical solution, yet it is essential to know the approximation errors in the computed solution to determine if the computed solution satisfies the desired requirements of accuracy.
7) There must be a mechanism to reduce the approximation errors to whatever desired level so that a desired accuracy is achievable in the computed solution.
8) The methodology must permit time accurate evolutions.This aspect becomes intrinsic in the approach if features 6) and 7) are present.
The features described above are indeed intrinsically present in the mathematical and computational infrastructure presented in this paper for ODEs in time.Section 2 describes the mathematical details of the formulation, approximations and approximation spaces as well as details of the computational infrastructure.The requirements 1) -8) are described in the paragraph above are all addressed in Section 2. It is shown that the finite element discretization in time in which the local approximations in time are in hpk framework permitting global differentiability of the approximation of order ( 1 k  ) in time when used with integral forms that are variationally consistent provides a mathematical and computational infrastructure for ODEs in time with all of the desired attributes.Numerical studies are presented in Section 3 for commonly used model problems and the numerical results from the present work are compared with Newmark and Wilson's  method.

Finite Element Approximation of ODEs in Time
For presenting the mathematical details and formulations We note that (2.3) requires differential operator in that we transfer all of the time or time differential time differentiation from u to v (left side of (2.3)) using integration by parts.In doing so, we obtain the following in general, ( , ) ( , ) , is called the concomitant which results as a consequence of integration by parts.Thus, based on the definition of symmetry in (2.3) That is, the adjoint A   of the time operator must be the same as the time operator A  and the concomitant must be zero.When the operator A  has even order derivatives in time then we could show that (2.5) holds.Thus, there are time differential oeprators for which (2.5) is possible.(2.6) contains boundary terms and hence it can be expanded.0 , , ,

Methods of Approximation in Time (Based on Time Integral Forms)
Before we consider finite element method of approximation it is perhaps more prudent to consider classical methods of approximation in time (these consider t  without discretization) for the two classes of operators A  to determine which methods of approximation are worthy of consideration in the finite element processes in time.We group the methods of approximation (based on integral form in time) in two categories: 1) those based on fundamental lemma [1][2][3][4][5] 2) those based on minimization of the residual functional [23,24]

Methods Based on Fundamental Lemma
The methods such as GM , PGM , WRM and / GM WF in time are the methods in this category.In all these methods we begin by using fundamental lemma [1][2][3][4][5] In / GM WF we also begin with (2.8) i.e., ( , ) ( , ) but transfer some differentiation from n d to v to obtain a weak form in time, ( , ) ( ) ( , ) ( ) Thus all these methods of approximation in time yield a time integral form ((2.9) or (2.11)).
( ) are substituted in (2.9) or (2.11), we obtain a system of linear or non-linear algebraic equations in i c from which we can determine i c .Remarks 1) These methods of approximation result in an integral form that is converted to an algebraic system upon substituting the approximation.The algebraic system is used to determine the unknowns i c .Thus in these methods we have a "necessary condition" only.
2) Existence and uniqueness of the solution n d or the coefficients i c is never addressed, hence must be considered for each application.
3) This situation can be corrected by establishing a correspondence between the integral forms (necessary conditions) and the calculus of variations.

Variationally Consistent and Variationally
Inconsistent Time Integral Forms Consider extremum of a functional ( ) 3) The sufficient condition or extremum principle is given by 2 0 minimum of ( ) ( ) 0 saddle point of ( ) 0 maximum of ( )


i.e. the first variation of the integral form must yield a unique extremum principle.
It is straight forward to show (see [24]) that the time integral forms resulting from GM , PGM , WRM and / GM WF are all variationally inconsistent.Variationally consistent integral forms yield algebraic systems in which the coefficient matrices are symmetric and positive definite and uniqueness of the solution n d is ensured.VIC integral forms, on the other hand, yield non-symmetric coefficient matrices that are not always ensured to be positive definite.See references [23,24] for more details.Thus, GM , PGM ,WRM and / GM WF in time are not worthy of consideration for finite element processes in time as a general computational methodology for ODEs in time as these do not ensure unconditionally stable computations due to VIC integral forms in time.( ) 2( , ) 2( , )

Methods of Approximation Based on Minimization of Residual Functional
In the following we consider the two categories of the mathematical classification of A  i.e., when A  is non-self adjoint and when A  is non-linear, and determine variational consistency or variational inconsistency of the time integral form given by (2.18).  in Taylor series about 0 n d and retaining only up to linear terms in

Case (b): When
From (2.24), we can solve for We note that 2 2 ( , ) in which  is a scalar generally between 0 and 2 and assumes the largest value between 0 and 2 for which 0 ( ) ( ) holds.This is referred to as line search.The entire process of solving for n d that satisfies ( ) 0 n g d  is called Newton's linear method with line search.

Remarks on Methods of Approximation in
Time Based on Integral Forms in Time 1) The determination of variational consistency or variational inconsistency of an integral form in time allows us to determine unconditional stability of the resulting computations which helps in deciding which integral forms in time are worthy of consideration for the development of finite element method in time for the ODEs in time.
2) The time integral forms resulting from GM , PGM , WRM and / GM WF in time are all variationally inconsistent regardless of whether the time operator A  is non-self adjoint or non-linear.
3) The time integral form resulting from the LSP in time based on residual functional is variationally consistent when the time differential operator A  is non-self adjoint.
4) The time integral form resulting from the LSP in time based on residual functional for non-linear time operators is also variationally consistent provided a) n d in ( ) 0 n g d  is determined using Newton's linear method, and b) the second variation of the residual functional is approximated by 2 ( ) 2( , ) (1)(2)(3)(4), it is obvious to conclude that only LSP in time based on the residual functional is meritorious of consideration as a general numerical solution methodology for developing finite element processes in time for ODEs in time.This would be applicable to the totality of all ODEs as it yields VC integral forms for both non-self adjoint and non-linear operators in time.

Finite Element Processes in Time Based on Minimization of a Residual Functional: Least Squares Finite Element Process
In this section we present details of finite element processes in time based on least squares method using the residual functional.Consider an ODE in time, with some initial conditions.Let ( ) We consider the residual functional E over the discretization ( )

in which v=
That is, g is a non-linear function of nodal degrees of freedom { }  .
Using, Newton's linear method (as described earlier) we obtain, which yields a unique extremum principle and we have, An improved solution is obtained using

Benefits and Meritorious Features of the
Proposed LS Finite Element Formulation for ODEs in Time in hpk Framework 1) Due to the mathematical classification of all time operators into non-self adjoint and non-linear categories, the task of the development of the methods of approximation for ODEs in time is reduced to these two categories that address the totality of all time operators.
2) The investigation of the VC of the time integral forms resulting from various methods of approximation reveals that only the time integral forms resulting from the LSP in time based on the residual are variationally consistent for both classes of differential operators.The finite element processes proposed here are based on LSP in time, hence they yield unconditionally stable computations for all choices of computational and physical parameters for both categories of differential operators in time.
3) The method can be obviously made time marching.One considers only one element in time and time marches in sequel to compute the entire evolution.This results in efficiency of computations.
4) The degree of local approximation p permits desired polynomial for local approximation over an element in time.The local approximations can also be hierarchical.This adds additional efficiency in the computational processes for linear operators if p -levels are increased progressively.7) In the approach presented here, we choose an element in time (thereby choosing h ) and a minimally conforming k and conduct a p -convergence study to ensure that ( ) is as closeto zero as we desire it to be.This may also require one or more adjustments in h (generally reduction).Thus the proposed methodology has a built in mechanism for error measure and control without the knowledge of the theoretical solution.
8) If we choose an element in time and ensure step vii) before we time march, then time accuracy of the evolution is evident.
In summary, we note that the proposed LS finite element processes in time in hpk framework for a single or a system of ODEs based on the residual has all of the desired features that we have listed in Section 1.3 and hence is highly meritorious of consideration as a general methodology for computing numerical solutions of ODEs in time.

Model Problems and Numerical Studies
In this section we consider three model problems that are commonly used as benchmark problems in the published work.For all three model problems, the evolutions computed using the finite element process based on LSP in time in hpk framework are compared with the Newmark method and Wilson's  method.

Model Problem 1
Consider the following non-homogeneous linear ODE in time [12,18]: ; This model problem is typical in finite element processes in structural dynamics when space and time are decoupled (using / GM WF in space) and when the decoupled ODEs are transformed to modal basis in which the ODEs in time become decoupled [24].This model problem is typical of a single ODE in modal basis.
In the numerical studies we choose ( ) 0 f t  and 0.1    used all of the studies conducted using present formulation).Figure 4(a) shows evolutions computed using the Newmark method, Wilson's  method presents formulation ( / 1.6 t T   , 9 p  to 19) and theoretical solution.Evolutions from the Newmark method and Wilson's  method deviate from the true solution but the deviations are not so significant.Both, the Newmark and Wilson's  methods are comparable, however, Wilson's  method has more pronounced phase shift. is greater than 0.1, evolution for both methods are in significant error.Quantitative assessment of amplitude decay, base elongation and phase shift is not conclusive from these graphs due to the highly diffusive mature of the theoretical solution (addressed in model problem 2).A remarkable thing to note here is that in the proposed formulation, extremely high accuracy of the evolution is achievable even for 1.6 t   whereas the Newmark and Wilson's  methods yield evolution with significant errors beyond 0.1 t   .

Model Problem 2
In this model problem we consider ( [10,11,13,19]) This is obviously a 1-D spring, mass, damper system.In the numerical studies we choose . A theoretical solution is given by ( ) sin( )  is defined by       Accuracy of the evolution: Since in this case the theoretical solution is periodic, this model problem serves as a good test to measure the accuracy of the proposed method.Amplitude decay, base elongation and phase shift are good measures of numerical dispersion in the computational process.For this study we choose 1.6 t T   and 13 p  .We note that for this choice of t  , a single time increment contains evolution that is 1.6 times the period of the wave.The evolution is computed for 100 time steps.Figure 10 shows evolution for the last six time steps ( 150. 4  t  to 160 t  ).The periodic nature of the solution is preserved without any measurable amplitude decay, base elongation and phase shift confirming time accuracy of the evolution.

Model Problem 3
Consider the following non-linear ODE in time [21]:    14 and 15.Findings and observations are identical to those reported for model problems 1 and 2 and hence are not repeated.
Accuracy of the evolution: Since in this case the theoretical solution is periodic, this model problem also serves as a good test to measure the accuracy of the proposed method.Amplitude decay, base elongation and phase shift are good measures of numerical dispersion in the computational process.For this study we choose 1.6 t T   and 13 p  .We note that for this choice of t  , a single time increment contains evolution that is 1.6  times the period of the wave.The evolution is computed for 100 time steps.Figure 17 shows evolution for the last six time steps ( 150.4 t  to 160 t  ).The periodic nature of the solution is preserved with virtually no amplitude decay, base elongation and phase shift confirming the time accuracy of the evolution.

Summary and Conclusions
In this paper we have considered methods of approximation for numerical solutions of ODEs in time resulting from decoupling of space and time.All time differential operators in ODEs are mathematically classified as: In this paper we have considered methods of approximation for numerical solutions of ODEs in time resulting from decoupling of space and time.All time differential operators in ODEs are mathematically classified as: non-self adjoint or non-linear.For these two categories  of time differential operators, the methods of approxima-tion based on time integral forms are considered: those based on (1) fundamental lemma such as GM , PGM , WRM and / GM WF in time and (2) minimization of a functional using the residual such as the least squares method in time.A correspondence is established between the time integral forms resulting from these two classes of approximation methods for the two categories of time differential operators and the elements of the calculus of variations.This results in the definitions of VC and VIC integral forms in time.VC integral forms in time yield unconditionally stable computations during the entire evolution whereas unconditional stability of computations is not always ensured in the case of VIC integral forms in time.The time integral forms resulting from GM , PGM , WRM and / GM WF are all VIC for both classes of time differential operators and hence are not meritorious in the development of a general computational infrastructure for ODEs in time.On the other hand, the integral forms resulting from the least squares process in time based on the residual is VC when the time differential operator is non-self adjoint and can be made VC in the case of non-linear time differential operators: by using 1) Newton's linear method for solving the non-linear algebraic equations and 2) approximating second variation of the residual functional by neglecting second variation of the residual function.
The finite element process in time, using least squares in time based on the residual when considered in hpk framework, provides control over the integration time step (i.e. or h or t  ), degree of local approximation ( p ) over an element, and the global differentiability of the evolution through k , the order of the approximation space yielding global differentiability of order ( 1 k  ) so that with the proper choices of h , p and k , quite complex evolutions can be accommodated in a single time step (as in the case of 1.6 t T   ) with desired accuracy.In addition, the residual functional(s) I (or e I ) are a measure of the solution error without the knowledge of the theoretical solution.As 0 I  , the approximation approaches the true solution.This feature is only possible in hpk framework and is a natural outcome of the least squares process in time.This methodology provides a computational infrastructure that addresses all ODEs in time in a consistent and rigorous manner and the resulting computational processes are unconditionally stable regardless of the nature of the time differential operator.Numerical studies are presented for three model problems.In all three problems we intentionally choose 2    so that the time period 1 T  and, hence t t T    .Numerical studies consider 0.1, 0.2, t T   0.4, 0.8 and 1.6.We note that when 1.6 t T   , a single element in time (i.e. the integration interval t  ) contains an evolution that is 1.6 times the time period.Findings in all three model problems are similar.For smaller t T  , lower p -levels suffice.As t T  is increased, higher p -levels are needed (but never beyond 7 or 9) for good accuracy.In all cases, I of the order of O (10 -6 ) or lower is achieved for p -levels of 7 or higher when 3 k  (minimally conforming space, solutions of class 2 ( ) , the computed evolution from these two methods has significant error.In these methods there is no mechanism to accommodate a more complex evolution in a time step than that corresponding to linear acceleration.Hence, for better accuracy there is no other alternative but to reduce the integration time step.In the present approach, the minimally conforming choice of k with progressively increasing p -levels permits accurate computation of more and more complex evolutions for a fixed increment of time.Numerical studies presented for model problems 2 and 3 for 100 time steps with 1.6 t T   at 13 p  , confirm that the evolution remains free of amplitude decay, base elongation and phase shift.Similar computations using Newmark method or Wilson's  method with reasonable accuracy (but worse than the present approach) would require 0.05 t T   , i.e. 3200 time steps.
In conclusion, the methodology presented here addresses all ODEs in time in a uniform and rigorous manner without any special and application dependent adjustments, yields unconditionally stable computations during the entire evolution, has a built in mechanism of error measure without the knowledge of theoretical solution, permits large time steps while maintaining desired accuracy of the evolution, and is free of amplitude decay, base elongation and phase shift for large t T  (with roper choices of k and p ) and, hence is highly meritorious.


10) are a system of ODEs in time.A solution of these in time would yield    and thereby   e  for each element of the spatial discretization and hence, the approximation ( for an element defined by (1.5) completely defines the evolution.


n d and if the derivatives of all orders are continuous then, based on the theorems of calculus of variations ( ) is a necessary condition for the existence of the extremum of ( ) n I d .
, where t  is the time increment, i.e., the length of the element in time.The minimally conforming space h V , i.e., local approximation ( ) e h u t of ( ) u t over e t  , the time domain of an element e of class 2 ( ) e t C  for which the integrals in the entire finite element LSP in time are Riemann.If we choose 2 k  , then all integral measures are Lebesgue.In the numerical studies for this model problem we only consider 3 k  and choose 0.1, 0.2, 0.4, 0.8, 1.6 t t T     with p -levels of 5, 7, ,19.We consider a three node p -version element in time of length t  for the first increment of time.nodal degrees of freedom (as well as those at 0 t  remaining after improving ICs ) are unknown and are to be computed.We compute the solution for the first time increment and study convergence of the solution and time march only upon convergence.Results are presented and discussed in the following.p-convergence study: For the first increment of time t  we choose 0.1, 0.2, 0.4, 0.8 t   and 1.6 and conduct a p -convergence study for each t  by progressively increasing p -levels from 5 to 19 for 3 k  (solutions of class 2 ( ) e t C  ).

Figure 4 . 4 .
Figure 3. Evolutions for model problem 1: 3 k  , , , , / 0.1 0.2 0.4 0.8 t T   and 1.6 for varying p -levels.(a) Evolution: / 0.1 t T   ; (b) Evolution: / 0.2 t T   ; (c) Evolution: / 0.4 t T   ; (d) Evolution: / 0.8 t T   ; (e) Evolution: / t T  .of the order of O (10 −6 ).These studies are instrumental in deciding the choice of p -level for time marching (we only time march when the solution for the current increment has good accuracy).Computations of the evolution: Evolutions are computed using time marching for / 0.1, 0.2, 0.4, t T   0.8 and 7.6 for p -levels of 5 to 19 using local ap-

Figure 4 (
b) and (c) present evolutions obtained using the Newmark method and Wilson's  method for / t T  ranging from 0.01 to 0.4 and a comparison with present solution (for /

Figure 9 . Evolutions for model problem 2
Figure 9. Evolutions for model problem 2: = 3 k and = 2 k , , , , / 0.1 0.2 0.4 0.8 t T   and 1.6 for varying p -levels.(a) Evolution: / 0.1 t T   ; (b) Evolution: / 0.2 t T   ; (c) Evolution: / 0.4 t T   ; (d) Evolution: / 0.8 t T   ; (e) Evolution: / t T  .those at 0 t  remaining after imposing ICs ) are unknown and are to be computed.We compute solutions for the first time increment and study convergence of the solution and time march only upon convergence.Results are presented and discussed in the following.p-convergence study: For the first increment of time t  , we choose 0.1, 0.2, 0.4, 0.8 t   and 1.6 and conduct a p -convergence study for each t  for progressively increasing p -levels from 3 to 19 for 3 k 

Figure 5
shows plots of the residual functional I versus degrees of freedom for 3 k  .The results of Figure 5 are also reported in Figure 6 as the residual functional I versus t T  for each p -level.Similar results for 2 k  (local approximation of class 1 ( )

C
 and p -levels of 3 to 19 for solutions of class 1 ( ) e t C  .Results are shown in Figures 9(a)-(e).First, the solutions of both classes produce almost identical to 19 produce almost identical results.The converged solutions are in excellent agreement with the theoretical solution.

Figure 11 .
Figure 11.Comparison with the Newmark and Wilson's  method: model problem 2, = 3 k .(a) The Newmark and Wilson's  methods: / 0.1 t T   ; (b) The Newmark method; (c) Wilson's  method.thenumerical studies.If we choose 2    , then the

Figure 12 .
Figure 12. p -convergence study for the first time increment ( = 3 k

C
 produce results that are almost as good as those using local approximation of class 2 ( ) e t C  .This is due to the fact that the theoretical solutions are smooth.When comparing I versus degrees of freedom for solutions of class 2 C and 1 C one finds that for a given number of degrees of freedom, lower values of I are obtained in the case of solutions of class 2 C , confirming better accuracy of the evolution.Convergence rates of I versus degrees of freedom are almost the same for 3 k  and 2 k  .When 3 k  (minimally conforming) all integral measures are Riemann and hence are true measures.When 2 k  , all integral measures are in the Lebesgue sense and hence are approximate, but as the approximation approaches the theoretical solution, Lebesgue measures approach Riemann measures.Comparisons with Newmark and Wilson's  methods clearly show that for 0.1 t   operator A has first and second order derivatives of  in time as well as the function  )

2.1. Mathematical Classification of Time Differential Operators
ODEs regardless of whether A is linear or non-linear, for time differential operators A  , (2.6) can not be satisfied, hence the time differential operators can not be symmetric.This leads to the following two categories of mathematical classification for all time differential operators.
7) In (2.7), 0 t  and   are the boundaries of the time domain at 0 t t  (initial conditions) and at t   an open boundary.Let us assume that based on ICs it is possible to show that  exists in case of all ICs of the ODE .When (2.12) and Variationally inconsistent (VIC ) time integral forms.A time integral form resulting from a method of approximation is termed VIC if it is in violation of one or more elements of the calculus of variations.It is sufficient to show the following.Let there exist a functional( ) n I d In this category of methods we consider least squares method or process ( LSM or LSP ) in time based on residual functional.Let

Comparison with the Newmark and Wilson's  methods:
Comparisons with the Newmark and Wilson's  methods are shown in Figure 11.For the first study