A Posteriori Error Computations in Finite Element Method for Initial Value Problems

K. S. Surana, J. Abboud

Department of Mechanical Engineering, University of Kansas, Lawrence, KS, USA.

DOI: 10.4236/ajcm.2025.151005                

Abstract

A posteriori error computations in the space-time coupled and space-time decoupled finite element methods for initial value problems are essential: 1) to determine the accuracy of the computed evolution, 2) if the errors in the coupled solutions are higher than an acceptable threshold, then a posteriori error computations provide measures for designing adaptive processes to improve the accuracy of the solution. How well the space-time approximation in each of the two methods satisfies the equations in the mathematical model over the space-time domain in the point wise sense is the absolute measure of the accuracy of the computed solution. When $L_2$ -norm of the space-time residual over the space-time domain of the computations approaches zero, the approximation ${\varphi }_h\left(x,t\right)\to \varphi \left(x,t\right)$ , the theoretical solution. Thus, the proximity of ${\Vert E\Vert }_{L_2}$ , the $L_2$ -norm of the space-time residual function, to zero is a measure of the accuracy or the error in the computed solution. In this paper, we present a methodology and a computational framework for computing ${\Vert E\Vert }_{L_2}$ in the a posteriori error computations for both space-time coupled and space-time decoupled finite element methods. It is shown that the proposed a posteriori computations require $h$ , $p$ , $k$ framework in both space-time coupled as well as space-time decoupled finite element methods to ensure that space-time integrals over space-time discretization are Riemann, hence the proposed a posteriori computations can not be performed in finite difference and finite volume methods of solving initial value problems. High-order global differentiability in time in the integration methods is essential in space-time decoupled method for posterior computations. This restricts the use of methods like Euler’s method, Runge-Kutta methods, etc., in the time integration of ODE’s in time. Mathematical and computational details including model problem studies are presented in the paper. To authors knowledge, it is the first presentation of the proposed a posteriori error computation methodology and computational infrastructure for initial value problems.

Keywords

A Posteriori Error Computation, Space-Time Coupled, Space-Time Decoupled, A Priori Error Estimation, A Posteriori Error Estimation, hpk Scalar Product Spaces, Minimally Conforming Scalar Product Spaces

1. Introduction, Literature Review, and Scope of Work

1.1. Introduction

The mathematical models of initial value problems describe evolutions in which dependent variables naturally exhibit simultaneous dependence on spatial coordinates and time. Thus, the space-time coupled finite element method [1] with space-time variationally consistent integral forms is ideally suited for obtaining numerical solutions to the initial value problems. In space-time coupled finite element methods, one constructs space-time discretization using space-time finite elements with space-time p-version hierarchical local approximations [1] in higher-order scalar product spaces in space and time. In this approach, it is possible to maintain space-time integrals over space-time discretization as Riemann integrals. One considers a space-time strip (in ${ℝ}^1$ ) or a space-time slab (in ${ℝ}^2$ ) and computes solution for the space-time discretization using $h$ , $p$ , and $k$ processes beginning with the first space-time strip $0\le t\le \text{Δ}t$ . Upon obtaining a converged solution for this space strip (requires a posteriori error computations to quantify the errors in the solution), the solution is time marched for subsequent space-time strips, keeping $h$ , $p$ , and $k$ constant. In this approach, an initial value problem in ${ℝ}^1$ becomes a 2D problem with $x_1$ and $t$ and an initial value problem in ${ℝ}^2$ ($x_1$ , $x_2$ coordinates) becomes a 3D problem with spatial coordinates $x_1$ , $x_2$ and time $t$ . Treatment of initial value problems in ${ℝ}^3$ require consideration of $x_1$ , $x_2$ , $x_3$ , and time $t$ , leading to significant complexities in the details of the computations. Thus, for initial value problems in ${ℝ}^3$ , space-time coupled finite element methods are generally avoided. We note that in the space-time coupled finite element methods for initial value problems, we construct a space-time discretization for an increment of time. Upon obtaining a converged solution for the current space-time strip (or slab), we proceed to the next increment of time i.e. next space-time strip using the calculated solution for the current space-time strip to determine initial conditions for the next space-time strip. This process continues till desired time is reached. We shall see in a subsequent section that a posteriori error computations in space-time coupled finite element method are quite straight forward when using $h$ , $p$ , and $k$ scalar product spaces for space-time coupled local approximations for the space-time elements and can be designed to yield precise a posteriori errors in the computed solution.

To overcome the difficulties in space-time coupled methods for initial value problems in ${ℝ}^3$ , we generally consider space-time decoupled methods for initial value problems. In space-time decoupled finite element methods, we construct a discretization in space for the spatial domain of the initial value problem i.e. a 1D discretization in $x_1$ , for the 1D spatial domain $\left[0,L\right]$ of the initial value problem in $x_1$ and $t$ , a 2D discretization in $x_1,x_2$ for the 2D spatial domain of the initial value problem in $x_1,x_2$ and time $t$ and so on. We construct an integral form in space for the spatial discretization using fundamental Lemma [1] [2] of the Calculus of Variations, in particular using Galerkin method with weak form (GM/WF). The local approximation for an element in space consists of approximation functions dependent on spatial coordinates and nodal degrees of freedom dependent on time. Upon substituting local approximations in the integral form and carrying out the integration in space, we obtain element expressions in which space coordinates are eliminated, and only degrees of freedom and their time derivatives remain. Assembly of element equations yields a system of ordinary differential equations in time. Thus, reducing partial differential equations in space and time to ordinary differential equations in time, a consequence of decoupling space and time.

The ordinary differential equations in time resulting from decoupling of space and time are integrated using explicit or implicit time integration methods or finite element methods in time to obtain solutions for the nodal degrees of freedom and their time derivatives. Upon substituting the solutions for nodal degrees of freedom in space-time decoupled local approximation, we obtain space-time approximation of the solution for each element of the spatial discretization. For the finite element discretization in space, the resulting ordinary differential equations in time are integrated for an increment of time. Using the calculated solution for the current increment of time, we proceed to the next increment of time and continue till the desired time is reached. Thus, we see that a posteriori error computations in space-time decoupled methods also need to be performed for an increment of time as well, same as in the case of space-time coupled methods. Designing a mathematical and computational infrastructure for a posteriori error computation for the space-time decoupled finite element solution is not as simple and as straight forward as for space-time coupled finite element processes as shown in subsequent section.

From the details of space-time coupled and space-time decoupled finite element methods for obtaining solutions, we note that in both methods we determine evolution for an increment of time and then time march. Thus, a posteriori computation methodologies and infrastructures for both space-time coupled and space-time decoupled require consideration for an increment of time.

1.2. Literature Review

Before the introduction of k-version of finite element method in which $k$ is the order of the approximation space, in addition to $h$ and $p$ as the independent parameters in all finite element computational processes by Surana et al. [1] [3]-[9], all finite element computations for boundary value problems as well as for initial value problems were performed using approximation of class $C^0$ in space for boundary value problems and of $C^0$ in space and time for initial value problems. Such solutions of global differentiability $C^0$ are not admissible in the partial differential equations describing the boundary value problem or the initial value problem, due to lack of required global differentiability, dictated by the orders of derivatives of the dependent variables in space and time. Inter-element discontinuities of the derivatives at the inter element boundaries lead to the integrals over the discretization in Lebesgue sense that are not reliable measure of residuals. Thus, in the $h$ -$p$ framework in the absence of $k$ (i.e. solution of class $C^0$ only), the order of the approximation space, a posteriori error computation based on space-time residual functionals obtained by substituting space-time approximations in the equations of the mathematical model and then obtaining $L_2$ -norm is not possible. This led to investigations and research into what is called “a posteriori estimation” methods [10]. In this approach, it is accepted that a posteriori computations are not possible, and one resorts to establishing a bound of the error in the computed solution. More compact or shorter interval bounds of a posteriori error estimates are obviously more useful. Unfortunately, the a posteriori error estimation methods were largely developed using weaknesses in the global approximation as the major factor driving the determination of the estimates. As an illustration, consider solutions of class $C^0$ . Such finite element solutions will exhibit discontinuity of the first derivative of the dependent variable(s) normal to the inter-element boundaries. Thus, for solutions of class $C^0$ , one could undertake derivation of a posteriori error estimates based on the jumps of the first derivatives of the solution normal to the inter-element boundaries. In the study of convergence of the finite element solutions, how this measure based on the jumps in the first derivatives approaches zero is some indication of the solution accuracy and convergence. This methodology has obvious shortcomings. If we consider solutions of class $C^1$ over the discretization, then the first derivative of the dependent variable in the tangential and normal direction to the inter-element boundaries is continuous. Hence, a posteriori error estimates derived based on the jumps in the first derivatives normal to the inter-element boundaries are meaningless as these don’t exist. Thus, we note that with the introduction of $k$ , the order of approximation space as an independent parameter in addition to $h$ and $p$ , it is possible to construct local approximations that yield desired global differentiability required by the orders of derivatives in the partial differential equations. We examine the highest order of derivatives of the dependent variable(s) in the partial differential equations in space and time and then choose minimum value $k=\left(k_1,k_2\right)$ , $k_1$ in space and $k_2$ in time, so that the residual function is continuous everywhere over the space-time domain or dicsretization. This allows us to maintain integrals in $L_2$ -norm computations in the Riemann sense over the discretization, essential to ensure accuracy of $L_2$ -norms. To our knowledge, there are no published works on a posteriori error computations except Surana et al. [1] [9] [11]. It is rather obvious that solutions of class $C^0$ cannot be solutions of any partial differential equation, as a partial differential equation at the very least must have first order derivatives of the dependent variables, requiring approximations of class $C^1$ if the integrals are to be Riemann over the discretization, thus it is not surprising that there are no published works on a posteriori error computations before the introduction of $k$ -version of finite element methods by Surana et al. [1] [3]-[9]

1.3. Why Proposed Approach? Proof of Accuracy and Convergence

Consider an initial value problem $A\varphi -f=0$ in ${\Omega }_{xt}$ with ${\overline{\Omega }}_{xt}^T={\displaystyle \underset{e}{\cup }{\overline{\Omega }}_{xt}^e}$ , as space-time discretization of space-time domain ${\overline{\Omega }}_{xt}$ in which ${\varphi }_h={\displaystyle \underset{e}{\cup }{\varphi }_h^e}$ is approximation of $\varphi$ over ${\overline{\Omega }}_{xt}^T$ and ${\varphi }_h^e$ , the approximation of $\varphi$ over ${\overline{\Omega }}_{xt}^e$ , a space-time element $e$ . Then, in the absence of the theoretical solution of the initial value problem, the best measure of the accuracy of ${\varphi }_h$ is how well it satisfies the initial value problem $A\varphi -f=0$ over ${\overline{\Omega }}_{xt}^T$ . This requires that we substitute ${\varphi }_h$ in $A\varphi -f=0$ , and since ${\varphi }_h$ is the approximation of $\varphi$ , there will be a residual function $E=A{\varphi }_h-f$ $\forall x,t\in {\overline{\Omega }}_{xt}^T$. The corresponding residual functional $I\left({\varphi }_h\right)$ i.e. ${\Vert E\Vert }_{L_2}^2$ is obtained by calculating:

$I\left({\varphi }_h\right)={\left(E,E\right)}_{{\overline{\Omega }}_{xt}^T}={\displaystyle \sum }_e{\left(E^e,E^e\right)}_{{\overline{\Omega }}_{xt}^T}={\displaystyle \sum }_e{\left(A{\varphi }_h^e-f,A{\varphi }_h^e-f\right)}_{{\overline{\Omega }}_{xt}^e}$

When ${\varphi }_h$ converges to the theoretical solution $\varphi$ , then $I\left({\varphi }_h\right)=I\left(\varphi \right)=0$ . Thus, proximity of $I\left({\varphi }_h\right)$ to zero is a measure of the accuracy or error in the solution ${\varphi }_h$ . This is undeniably the best measure of solution accuracy or error as $I\left({\varphi }_h\right)$ indicates how well $A\varphi -f=0$ is satisfied in the pointwise sense over the space-time discretization ${\overline{\Omega }}_{xt}^T$ .This is a posteriori error computation, however, for the a posteriori error computations to be accurate, the approximation ${\varphi }_h$ must have the required differentiability over ${\overline{\Omega }}_{xt}^T$ as dictated by the highest order of the derivatives of $\varphi$ in space and time in $A\varphi -f=0$ . With the choice of minimally conforming approximation spaces for ${\varphi }_h^e$ , the integrals in the calculation of $I\left({\varphi }_h\right)$ over ${\overline{\Omega }}_{xt}^T$ are ensured to be Riemann, in which case when $I\left({\varphi }_h\right)\to 0$ , $A\varphi -f=0$ is satisfied over ${\overline{\Omega }}_{xt}^T$ in the pointwise sense i.e. $A{\varphi }_h-f\to 0$ holds everywhere in ${\overline{\Omega }}_{xt}^T$ . With minimally conforming approximation spaces, $h,p$ refinements are carried out to ensure that $I\left({\varphi }_h\right)\to 0$ , generally of the order $\mathcal{O}\left({10}^{-6}\right)$ or lower when using dimensionless form of the equations in the mathematical model describing the initial value problem is found to be adequate but can be reduced as low as $\mathcal{O}\left({10}^{-15}\right)$ with much improved accuracy.

In the currently published works on a posteriori literature, only a posteriori error estimates (bounds on error) are reported, generally derived using jumps in the derivatives of the dependent variable(s) normal to the inter-element boundaries. Obviously, these estimates are incapable of determining quantitatively how well the equations in the mathematical model of the initial value problem are satisfied by the computed solution. There is no reported work on a posteriori error computations because valid, meaningful and accurate a posteriori computations are not possible with solution of class $C^0$ in space and time used currently for initial value problems. Thus, there is nothing in the published works to compare the residual reported in this paper. Convergence and proof of convergence of ${\Vert E\Vert }_{L_2}$ to computed zero based on $h,p$ refinements for minimally conforming approximation spaces follows the details on convergence reported by Surana et al. [1] [6] [9] [11]. For a fixed mesh (fixed $h$ ), progressively increasing $p$ results in progressively increasing convergence rate of ${\Vert E\Vert }_{L_2}$ . Likewise, for a fixed $p$ with $h$ -refinement, convergence rate does not change. Increased in $p$ -level results in increased convergence rate with $h$ -refinement. Details of the convergence of ${\Vert E\Vert }_{L_2}$ and its rate of convergence can be found in references [1] [6] [9] [11].

1.4. Scope of Work

In this paper, we present a posteriori error computation methodology and computational infrastructure for space-time coupled as well as space-time decoupled finite element processes for initial value problems strictly based on $L_2$ -norm of the space-time residual function obtained using the mathematical model for the discretized space-time computational domain. In space-time coupled finite element method, the $p$ -version hierarchical local approximations with required regularity are substituted using computed solutions in the partial differential equations to compute residual functions followed by the computation of $L_2$ -norm of the residual functions over the space-time elements of the discretization yielding ${\Vert E\Vert }_{L_2}$ , the $L_2$ -norm of residual function for the space-time discretization. When the computed solution ${\varphi }_h\left(x,t\right)$ approaches $\varphi \left(x,t\right)$ , the theoretical solution, ${\Vert E\Vert }_{L_2}\to 0$ . Thus, proximity of ${\Vert E\Vert }_{L_2}$ to zero is a measure of error in the solution or is a measure of accuracy of the solution. Thus, in space-time coupled methods, a posteriori error computation involves computing ${\Vert E\Vert }_{L_2}$ over the space-time discretization using the computed solution for an increment of time.

In space-time decoupled methods, we have finite element discretization in space, followed by integral form over the spatial discretization using fundamental lemma of the calculus of variations (generally GM/WF). Substitution of local approximation, in which approximations are function of space coordinates and are $p$ -version hierarchical with higher-order global differentiability and the degrees of freedom are function of time. After carrying out integration in space and summing element equations, we obtain ordinary differential equations in time. These are integrated using implicit time integration methods such as Newmark’s linear method. When the approximation functions in space and the time approximation of solution in the time integration method have desired global differentiability in space and time, the space-time decoupled approximations $\varphi \left(x,t\right)$ can also be substituted in the partial differential equations of the mathematical model to obtain space-time residual functions for the element in space and the time increment constituting space-time domain. This is followed by $L_2$ -norm computation for the space-time domain. Proximity of ${\Vert E\Vert }_{L_2}$ to zero in this case is a measure of error or accuracy of the solution.

Thus, the proposed framework of “a posteriori” error computation requires higher-order global differentiability approximations in space and time in space-time coupled as well as decoupled methods. When we choose approximations in space and time that are minimally conforming, the space-time integrals over the space-time discretization in both methods are ensured to be Riemann, essential for accurate computation of ${\Vert E\Vert }_{L_2}$ . In choosing time integration method for ordinary differential equations in time in space-time decoupled method, we must ensure that the time integration scheme has desired global differentiability in the entire discretized time domains.

Thus, the proposed a posteriori error computations are based on space-time residual functions for the discretized space-time domains in both space-time coupled and space-time decoupled methods. These error measures are true measures of how well the partial differential equations in the mathematical model are satisfied by the computed solutions in the pointwise sense, ensured by the Riemann integrals over the space-time discretization. Clearly, these a posteriori error computations require $h$ , $p$ and $k$ framework in which desired global differentiability in space and time as necessitated by the partial differential equations is possible. The proposed a posteriori error computations are not possible in finite difference or finite volume methods of solving initial value problems.

Mathematical and computational details associated with the proposed a posteriori error computation are presented in this paper. This is followed by model problem studies to demonstrate the applications, effectiveness, and usefulness of a posteriori computation frameworks presented in the paper for space-time coupled as well as space-time decoupled finite element methods.

2. Theoretical Consideration and Computational Infrastructure

In this section, we present some theoretical considerations in space-time coupled and space-time decoupled finite element methods that are necessary to design a possible a posteriori error computation strategy in these methods for initial value problems. We present details in the following.

2.1. A Posteriori Error Computation in Space-Time Coupled Finite Element Method for Initial Value Problems

Let

$A\varphi -f=0\text{in}{\Omega }_{xt}={\Omega }_x\times {\Omega }_t$ (1)

be an initial value problem over the space-time domain ${\Omega }_{xt}$ in which $A$ is the space-time differential operator containing $2m_1$ and $2m_2$ as the highest orders of derivatives in space and time. Let ${\overline{\Omega }}_{xt}={\Omega }_{xt}\cup \Gamma$ , in which $\text{Γ}$ is the closed boundary of ${\Omega }_{xt}$ . Surana et al. [1] have shown that space-time differential operators are non self-adjoint or non-linear but not self-adjoint as in space-time differential equations, $A^{\text{*}}$ adjoint of the operator is not the same as the operator $A$ and there are issues in the concomitant due to open boundary. Thus, for initial value problems, all space-time integral form based on fundamental lemma (Space-time Galerkin Method, Space-time Galerkin Method with Weak Form, Space-time Petrov-Galerkin Method, Space-time Weighted Residual Method) are space-time variationally inconsistent due to lack of unique extremum principle. Thus, for initial value problems, only the integral form based on space-time residual functional are space-time variationally consistent [1], hence are considered in the present work for a posteriori error computations.

(a) Space-time domain ${\overline{\Omega }}_{xt}$ (b) Discretization of ${\overline{\Omega }}_{xt}$ into space-time strips

(c) Discretization of ${\overline{\Omega }}_{xt}^{\left(n\right)}$ into space-time elements.

Figure 1. Space-time domain, space-time strips, and discretization of the $n^{th}$ space-time element.

Let ${\left({\overline{\Omega }}_{xt}^{\left(n\right)}\right)}^T={\displaystyle \underset{e}{\cup }{\overline{\Omega }}_{xt}^{\left(e\right)}}$ be discretization of the $n^{th}$ space-time strip or the $n^{th}$ space-time slab ${\overline{\Omega }}_{xt}^{\left(n\right)}$ in which ${\overline{\Omega }}_{xt}^{\left(e\right)}$ is a $p$ -version hierarchical space-time finite element (see Figure 1). Let ${\varphi }_h^e\left(x,t\right)$ be local approximation of $\varphi$ over ${\overline{\Omega }}_{xt}^{\left(e\right)}$ given by

${\varphi }_h^e\left(x,t\right)={\displaystyle \sum }_{i=1}^n{N}_i\left(x,t\right){\delta }_i^e$ (2)

in which $N_i\left(x,t\right)$ are space-time local approximation functions for an element ${\overline{\Omega }}_{xt}^{\left(e\right)}$ derived using its map in the natural coordinate space $\xi$ , $\eta$ and $\zeta$ in a one unit length, one unit square or a one unit cube in ${ℝ}^1$ , ${ℝ}^2$ and ${ℝ}^3$ respectively.

Let

${\varphi }_h\left(x,t\right)={\displaystyle \underset{e}{\cup }{\varphi }_h^e\left(x,t\right)}$ (3)

be approximation of $\varphi \left(x,t\right)$ over the discretization ${\left({\overline{\Omega }}_{xt}^{\left(n\right)}\right)}^T$ of ${}_{\overline{\text{Ω}}}^{xt}$, the space-time domain of $n$ -th space-time strip.

Since $2m_1$ and $2m_2$ are the highest orders of the derivatives of $\varphi$ in the space and time, then we have:

${\varphi }_h\left(x,t\right)\in V\subset H^{\left(k,p\right)}\left({\left({\overline{\Omega }}_{xt}^{\left(n\right)}\right)}^T\right)$ (4)

in which

$\begin{array}{l}{k}_1\ge 2m_1+1,k_2\ge 2m_2+1,\hfill \\ p_1\ge 2k_1-1,p_2\ge 2k_2-1\hfill \end{array}$ (5)

Thus,

${\varphi }_h^e\left(x,t\right)\in V\subset H^{\left(k,p\right)}\left({\overline{\Omega }}_{xt}^{\left(e\right)}\right)$ (6)

must hold i.e. the union of local approximation ${\varphi }_h^e$ in (3) must yield the desired global differentiability of ${\varphi }_h\left(x,t\right)$ over ${\left({\overline{\Omega }}_{xt}^{\left(n\right)}\right)}^T$ . Equation (6) implies that the local approximation function $N_i\left(.,.\right)$ must belong to the scalar product space $H^{\left(k,p\right)}\left({\Omega }_{xt}^e\right)$ defined by (5), in which

$k_1=2m_1+1,k_2=2m_2+1\text{and}{p}_1=2k_1-1,p_2=2k_2-1$

define minimally conforming scalar product approximation space for which the space-time integrals over ${\left({\overline{\Omega }}_{xt}^{\left(n\right)}\right)}^T$ will be Riemann.

In the space-time finite element method based on space-time residual functional, the space-time integral forms that are space-time variationally consistent for linear as well as non-linear space-time differential operators [1], hence are meritorious for obtaining solutions of initial value problems. Details of space-time coupled finite element methods for linear and nonlinear differential operators can be found in [1] and are omitted here for the sake of brevity.

If ${\varphi }_h\left(x,t\right)$ in (4) is minimally conforming space or is of class higher than minimally conforming, then the space-time integrals over ${\left({\Omega }_{xt}^{\left(n\right)}\right)}^T$ are Riemann and the residual function over ${\left({\Omega }_{xt}^{\left(n\right)}\right)}^T$ and is given by:

$A{\varphi }_h-f=E\left(x,t\right)\forall \left(x,t\right)\in {\left({\Omega }_{xt}^{\left(n\right)}\right)}^T$ (7)

residual functional $I\left({\varphi }_h\right)$ is given by

$I\left({\varphi }_h\right)=\left(E,E\right)$ (8)

or

$I\left({\varphi }_h\right)={\displaystyle \sum }_e{I}^e\left({\varphi }_h^e\right)={\displaystyle \sum }_e{\left(E^e,E^e\right)}_{{\Omega }_{xt}^e}$ (9)

in which

$E^e=A{\varphi }_h^e\left(x,t\right)-f\forall x\in {\overline{\Omega }}_{xt}^e$(10)

and $I^e$ is the residual functional and $E^e$ is the residual function over ${\overline{\Omega }}_{xt}^e$ , domain a space-time element e.

For a space-time strip we calculate $I^e$ and $I$ using (9). When ${\varphi }_h\to \varphi$ , the theoretical solution, then

$I\left({\varphi }_h\right)\to I\left(\varphi \right)=0$ (11)

Thus, $I^e\left({\varphi }_h\right)$ is a posteriori measure of error in space-time element $e$ and $I\left({\varphi }_h\right)$ is the a posteriori measure of error for the space-time discretization of a space-time strip.

2.2. A Posteriori Error Computations Procedure

The procedure for calculating a posteriori error based on space-time residual functional for space-time coupled finite element methods for initial value problem is given in the following:

1) Choose minimally conforming spaces for local approximation based on the highest orders of the derivatives of dependent variable(s) in space and time so that the space-time integral over the space-time finite element discretization is Riemann. If the evolution is sufficiently smooth, the order of the space in space and time can be lower by one, if the space-time integral in Lebesgue sense suffices over the space-time discretization depending on the smoothness of the solution.

2) Calculate the solution for the first space-time strip or slab, i.e., ${\varphi }_h^e$ and ${\varphi }_h={\displaystyle {\cup }_e{\varphi }_h^e}$ . In cases of nonlinear initial value problems, calculate the converged solution using iterative solution methods (such as Newton’s linear method).

3) Using ${\varphi }_h^e\left(x,t\right)$ for each space-time element, calculate the residual functionals $I^e\left({\varphi }_h^e\right)$ for ${\overline{\Omega }}_{xt}^e$ and $I\left({\varphi }_h\right)$ for ${\left({\overline{\Omega }}_{xt}^{\left(1\right)}\right)}^T$ space-time discretization of the first space-time strip.

$I^e\left({\varphi }_h^e\right)={\left(E^e,E^e\right)}_{{\overline{\Omega }}_{xt}^e},E^e=A{\varphi }_h^e-f$(12)

$I\left({\varphi }_h\right)={\displaystyle \sum }_e{I}^e\left({\varphi }_h^e\right)$ (13)

$I\left({\varphi }_h\right)$ is a measure of error in the computed solution for the first space-time strip or slab.

4) When ${\varphi }_h\to \varphi$ , $I\left({\varphi }_h\right)\to I\left(\varphi \right)=0$ , thus proximity of calculated $I\left({\varphi }_h\right)$ to computed zero ($\mathcal{O}\left({10}^{-4}\right)$ or lower) is a measure of error in the solution for the space-time discretization of ${\overline{\Omega }}_{xt}^{\left(1\right)}$ and $I^e\left({\varphi }_h^e\right)$ is a measure of error in the space-time finite element with space-time domain ${\overline{\Omega }}_{xt}^e$ .

5) We use $h$ , $p$ adaptive refinement for ${\overline{\Omega }}_{xt}^{\left(1\right)}$ to obtain converged solutions (i.e., $I\left({\varphi }_h\right)\le \mathcal{O}\left({10}^{-4}\right)$ ), keeping minimally conforming $k=\left(k_1,k_2\right)$ . $I^e\left({\varphi }_h^e\right)$ are helpful in refining mesh or changing $p$ -levels in space and time. Upon obtaining converged solutions for ${\overline{\Omega }}_{xt}^{\left(1\right)}$ , we time march to the next space-time strip $\left(t_0+\text{Δ}t\le t\le t_0+2\text{Δ}t\right)$ , keeping $h$ , $p$ and $k$ fixed. In general, this procedure works well. If a change in $h$ or a change in $p$ -level is necessary for subsequent space-time strips, then, unfortunately, we must again begin with ${\overline{\Omega }}_{xt}^{\left(1\right)}$ with the new mesh and $p$ -levels.

6) We note that for a space-time strip or a space-time slab the residual functional $I^e$ for an element $e$ with domain ${\overline{\Omega }}_{xt}^e$ is a measure of the error in the solution for element e. These residual functionals $I^e,e=1,2,\cdots$ are quite effective in designing an $h,p$ (possibly $k$ ) adaptive strategy. The element with the largest value of $I^e$ is the first candidate for $h$ or $p$ or $hp$ -refinements. This process can be designed such that all elements with a value of $I^e$ larger than $I^{\text{*}}$ (a threshold value) are considered for $h$ or $p$ or $hp$ -refinements. How these refinements are carried out is beyond the scope of this work. Here we only wish to point out that quantities $I^e$ calculated in ‘a posteriori’ error computations indeed provide excellent measures to design adaptive processes.

3. A Posteriori Error Computation in Space-Time Decoupled Finite Element Processes

3.1. Mathematical Model, Integral Form in Space

Let

$A\varphi -f=0\text{in}{\Omega }_{xt}={\Omega }_x\times {\Omega }_t$ (14)

be the initial value problem over the space-time domain ${\Omega }_{xt}$ , in which $A$ is a space-time differential operator of the highest order $2m_1$ in space and $2m_2$ in time. ${\overline{\Omega }}_{xt}={\Omega }_{xt}\cup \Gamma$ , where $\text{Γ}$ is the closed boundary of ${\Omega }_{xt}$ . We discretize the spatial domain ${\overline{\Omega }}_x$ with finite elements in ${ℝ}^1$ , or in ${ℝ}^2$ or ${ℝ}^3$ in space such that ${\overline{\Omega }}_x^T$ is the discretization of ${\overline{\Omega }}_x$ , and consists of ${\overline{\Omega }}_x^T={\displaystyle \underset{e}{\cup }{\overline{\Omega }}_x^e}$ , in which ${}^{\overline{\text{Ω}}}\left(x\right)$ is a finite element in the spatial domain ${\overline{\Omega }}_x$ . Consider an increment of time $\text{Δ}t$ , and let ${\varphi }_h\left(x,t\right)$ be the approximation of $\varphi$ over the spatial domain ${\overline{\Omega }}_x^T={\displaystyle \underset{e}{\cup }{\overline{\Omega }}_x^e}$ for time $t_0\le t\le t_0+\text{Δ}t$ in which $t=t_0$ is the initial time when the evolution described by (14) commences. We consider ${\varphi }_h\left(x,t\right)$ to be the union of ${\varphi }_h^e\left(x,t\right)$

${\varphi }_h\left(x,t\right)={\displaystyle \underset{e}{\cup }{\varphi }_h^e\left(x,t\right)}$ (15)

in which ${\varphi }_h^e\left(x,t\right)$ is the approximation of $\varphi \left(x,t\right)$ over ${\overline{\Omega }}_x^e$ for time $t_0\le t\le t_0+\text{Δ}t$ , given by:

${\varphi }_h^e\left(x,t\right)={\displaystyle \sum }_{i=1}^n{N}_i\left(x\right){\delta }_i^e\left(t\right)\text{and}{\varphi }_h\left(x,t\right)={\displaystyle \underset{e}{\cup }{\varphi }_h^e\left(x,t\right)}$ (16)

In (16), $N_i\left(x\right)$ are the approximation functions that are functions of space coordinates only and ${\delta }_i^e\left(t\right)$ are nodal degrees of freedom that are a functions of time only. The local approximation ${\varphi }_h^e\left(x,t\right)$ (16) decouples space and time for ${\overline{\Omega }}_{xt}^e$ and ${\varphi }_h\left(x,t\right)$ decouples space and time for the discretization ${\left({\overline{\Omega }}_x\right)}^T\times {\overline{\Omega }}_t={\left({\overline{\Omega }}_x\right)}^T\times \left[t_0,t_0+\Delta t\right]$ . Using fundamental lemma of the calculus of variations, we construct the integral form of (14) over the discretization ${\overline{\Omega }}_x^T$ of the spatial domain ${\overline{\Omega }}_x$ using a test function $v\left(x\right)$ and treating time as a constant.

${\left(A{\varphi }_h\left(x,t\right)-f,v\left(x\right)\right)}_{{\left({\overline{\Omega }}_x\right)}^T}=0$ (17)

Since (17) is a functional, we can write this as a sum of the functional over ${\overline{\Omega }}_x^e$ i.e.

${\left(A{\varphi }_h\left(x,t\right)-f,v\right)}_{{\overline{\Omega }}_x^T}={\displaystyle \sum }_e{\left(A{\varphi }_h^e\left(x,t\right)-f,v\right)}_{{\overline{\Omega }}_x^e}=0$ (18)

Consider $\left(A{\varphi }_h^e\left(x,t\right)-f,v\left(x\right)\right)$ in Galerkin methods with weak form in space treating time constant, we perform integration by parts in space if beneficial [1] [9] to obtain the following:

${\left(A{\varphi }_h^e\left(x,t\right)-f,v\left(x\right)\right)}_{{\overline{\Omega }}_x^T}\underset{\text{IBP}}{=}{B}^e\left({\varphi }_h^e,v\right)-l^e\left(v\right)$(19)

We perform some integration by parts in ${\left(A{\varphi }_h^e\left(x,t\right)-f,v\right)}_{{}_{\overline{\text{Ω}}}^x}$ only for those terms that contain even order derivatives in spatial coordinate $x$ to transfer half of the differentiation from ${\varphi }_h^e\left(x,t\right)$ to the test function $v\left(x\right)$ . Integral form (19) for ${\overline{\Omega }}_x^e$ is called the weak form of the initial value problem (14) in space. In the weak form the order of differentiation on ${\varphi }_h^e\left(x,t\right)$ , with respect to $x$ may have been lowered or weakened, hence the name weak form, but increased on the test function. After performing the integration in space for an element $e$ , we obtain a system of ordinary differential equations in time and $\left\{{\delta }^e\right\}$ and its time derivatives. The order of the time derivatives of $\left\{{\delta }^e\right\}$ in the ordinary differential equations depends upon the order of the time derivatives of $\varphi$ in (14), that is (19) can be written as

$\begin{array}{l}{\left(A{\varphi }_h^e\left(x,t\right),v\left(x\right)\right)}_{{\overline{\Omega }}_x^e}=B^e\left({\varphi }_h^e,v\left(x\right)\right)-l^e\left(v\right)\\ =H_1^e\left\{{\delta }^e\right\}+H_2^e\left\{{\dot{\delta }}^e\right\}+\cdots +H_n^e\left\{{\delta }_{\left[n\right]}^e\right\}-\left\{P^e\right\}-\left\{F^e\right\}\end{array}$ (20)

Substituting (20) into (18), we obtain

${\left(A{\varphi }_h\left(x,t\right)-f,v\right)}_{{\overline{\Omega }}_{xt}^T}=\left[H_1\left]\left\{\delta \right\}+\right[H_2\right]\left\{\dot{\delta }\right\}+\cdots +\left[H_n\right]\left\{{\delta }_{\left[n\right]}\right\}-\left\{P\right\}-\left\{F\right\}=0$ (21)

In which

$H_1={\displaystyle \sum }_e\left[H_1^e\right]$ , $H_2={\displaystyle \sum }_e\left[H_2^e\right]$ … etc. and $\left\{\delta \right\}={\displaystyle \underset{e}{\cup }\left\{{\delta }_e\right\}}$ , $\left\{\dot{\delta }\right\}={\displaystyle \underset{e}{\cup }\left\{{\dot{\delta }}_e\right\}}$ … etc.

In (21), we have a system of ordinary differential equations in time for the initial value problem (14) for ${\overline{\Omega }}_{xt}^T={\overline{\Omega }}_x^T\times \left[t_0,t_0+\Delta t\right]$

3.2. Methods for Obtaining Solution of Ordinary Differential Equations in Time

The system of ordinary differential equations in (21) is linear if the space-time differential operator $A$ in (14) is linear, otherwise (21) represents a system of non-linear ordinary differential equations in time.

The ordinary differential equations in time in (21) can be integrated in time using:

1) Finite element methods in time

2) Explicit time integration methods

3) Implicit time integration methods

Explicit time integration methods have poor accuracy and stability issues but are much simpler and are quite fast computationally. We do not consider these methods in this work because of poor accuracy, but also due to other requirements discussed later.

3.3. Finite Element Method in Time

Consider space-time coupled finite element method for initial value problems. The local approximation of ${\varphi }^e\left(x,t\right)$ over a space-time element ${\overline{\Omega }}_x^e$ or ${\overline{\Omega }}_{\xi }^e$ is given by:

${\varphi }_h^e\left(x,t\right)=\left[N\left(x,t\right)\right]\left\{{\delta }^e\right\}$ (22)

If $\left[N\left(x,t\right)\right]$ are generated using tensor product, then we can write (22) as follows:

${\varphi }_h^e\left(x,t\right)=\left[N\left(x\right)\right]\otimes \left[N\left(t\right)\right]\left\{{\delta }^e\right\}$ (23)

The $\otimes$ symbol indicates tensor product [1] [6] as we use it in finite elements to generate 2D approximation functions from 1D approximation functions.

We can also write (23) as:

${\varphi }_h^e\left(x,t\right)=\left[N\left(x\right)\right]\left(\left[N\left(t\right)\right]{\delta }^e\right)$ (24)

or

${\varphi }_h^e\left(x,t\right)=\left[N\left(x\right)\right]\left\{{\delta }^e\left(t\right)\right\}$ (25)

where

${\delta }^e\left(t\right)=\left[N\left(t\right)\right]\left\{{\delta }^e\right\}$ (26)

In (24) and (26) $\left\{{\delta }^e\right\}$ is not a function of time.

Remarks

1) Equation (24) is the space-time coupled finite element method.

2) Equation (25) is the space-time decoupled finite element method.

3) Equation (26) is the finite element method in time.

4) Substitution from (26) into (25) yields (24).

That is in space-time decoupled methods, if we solve ordinary differential equations in time using finite element method in time, then the space-time coupled and the space-time decoupled methods are the same provided the space-time approximation functions are generated using tensor product of 1D space and time approximation functions. This has been investigated by Surana et al. [12] and is shown to be true. For further details on this approach, see reference [12].

3.4. Implicit Time Integration Methods for Ordinary Differential Equations in Time

The a posteriori error computation proposed in the present work considers $L_2$ -norm of the space-time residual function for space-time decoupled finite element methods. This technique requires:

(i) Desired regularity or differentiability of ${\varphi }_h^e\left(x,t\right)=\left[N\left(x\right)\right]{\delta }^e\left(t\right)$ in space and time. Since the spatial discretization uses finite elements in space, the desired regularity in space can be easily satisfied by appropriate choices of $k$ and $p$ in space.

(ii) We also need desired global regularity or differentiability of ${\varphi }_h\left(x,t\right)$ in time as well. That is, we must consider a time integration method in which approximation of $\varphi$ over a time increment is a function or an expression in time (for example, polynomial in time).

(iii) At the inter-element boundaries, desired order of differentiability required by the initial value problem must be satisfied so that the space-time integrals are Riemann.

(iv) When the requirements (i) - (iii) are satisfied, the space-time decoupled approximation ${\varphi }_h^e\left(x,t\right)=\left[N\left(x\right)\right]\left\{{\delta }^e\left(t\right)\right\}$ can be substituted in the initial value problem (14) for any value of time in the interval $\left(n-1\right)\text{Δ}t\le t\le n\text{Δ}t$ . Thus, permitting $L_2$ -norm of the space-time residual functional over ${\overline{\Omega }}_x^e\times \left[\left(n-1\right)\Delta t\le t\le n\Delta t\right]$ for the $n^{th}$ space-time step.

(v) The requirements (i) - (iv) rule out the use of explicit time integration methods and the integration method in which we do not have an expression of the solution as a function of time explicitly for an increment of time.

Remarks [a.]

a. Requirements (i) - (iv) rule out methods like Runge-Kutta in which we do not have an analytical expression of solution as a function of time in a time interval as required in the $L_2$ -norm computations.

b. Wilson’s $\theta$ -method with linear acceleration and Newmark’s linear acceleration methods are suitable for second-order systems of ordinary differential equations in time.

c. In the case of first-order systems of ordinary differential equations in time, we must derive details of time integration methods based on Wilson’s $\theta$ method and Newmark method.

d. In the present work, we consider the Newmark method for both first-order and second-order systems of ordinary differential equations in time. For simplicity, we only consider linear ordinary differential equations in time. In the following, we present details of Newmark time integration method for linear first-order and second-order ordinary differential equations in time.

3.4.1. Newmark Time Integration Method for a System of Linear First Order Ordinary Differential Equations in Time

Consider the following first-order system of linear ordinary differential equations in time obtained from decoupling of space and time

$\left[C\right]\left\{\dot{\delta }\right\}+\left[K\right]\left\{\delta \right\}=\left[F\right]$ (27)

Consider a time interval $\left[t,t+\text{Δ}t\right]$ and assume that $\left\{\dot{\delta }\right\}$ is linear in the interval $\left[t,t+\text{Δ}t\right]$ (Figure 2).

Figure 2. Newmark linear method for first order ordinary differential equations.

Let $\tau$ be the time measured from $t$ . Let ${\left\{\dot{\delta }\right\}}_t$ and ${\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}$ be values of $\left\{\dot{\delta }\right\}$ at $t$ and $t+\text{Δ}t$ , respectively. Then

${\left\{\dot{\delta }\right\}}_{t+\tau }={\left\{\dot{\delta }\right\}}_t+\frac{\tau }{\text{Δ}t}\left({\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}-{\left\{\dot{\delta }\right\}}_t\right)$(28)

Integrating with respect to $\tau$ , we get:

${\left\{\delta \right\}}_{t+\tau }=\tau {\left\{\dot{\delta }\right\}}_t+\frac{{\tau }^2}{2\text{Δ}t}\left({\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}-{\left\{\dot{\delta }\right\}}_t\right)+C$(29)

At $\tau =0$ ,

${\left\{\delta \right\}}_{t+\tau }={\left\{\delta \right\}}_t$ (30)

Substituting in (29), we obtain:

$C={\left\{\delta \right\}}_t$ (31)

$\therefore {\left\{\delta \right\}}_{t+\tau }={\left\{\delta \right\}}_t+\tau {\left\{\dot{\delta }\right\}}_t+\frac{{\tau }^2}{2\text{Δ}t}\left({\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}-{\left\{\dot{\delta }\right\}}_t\right)$(32)

Letting $\tau =\text{Δ}t$ , we obtain:

${\left\{\delta \right\}}_{t+\text{Δ}t}={\left\{\delta \right\}}_t+\text{Δ}t{\left\{\dot{\delta }\right\}}_t+\frac{\text{Δ}t}{2}\left({\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}-{\left\{\dot{\delta }\right\}}_t\right)$(33)

From (28),

${\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}={\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}$

From (33), we solve for ${\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}$:

${\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}=\frac{2}{\text{Δ}t}\left({\left\{\delta \right\}}_{t+\text{Δ}t}-{\left\{\delta \right\}}_t\right)-{\left\{\dot{\delta }\right\}}_t$(34)

Consider (27) at $t+\text{Δ}t$ :

$\left[C\right]{\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}+\left[K\right]{\left\{\delta \right\}}_{t+\text{Δ}t}={\left\{F\right\}}_{t+\text{Δ}t}$(35)

Substitute ${\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}$ from (34) into (35):

$\left[C\right]\left(\frac{2}{\text{Δ}t}\left({\left\{\delta \right\}}_{t+\text{Δ}t}-{\left\{\delta \right\}}_t\right)-{\left\{\dot{\delta }\right\}}_t\right)+\left[K\right]{\left\{\delta \right\}}_{t+\text{Δ}t}={\left\{F\right\}}_{t+\text{Δ}t}$(36)

or

$\left(\left[K\right]+\frac{2}{\text{Δ}t}\left[C\right]\right){\left\{\delta \right\}}_{t+\text{Δ}t}={\left\{F\right\}}_{t+\text{Δ}t}+\frac{2}{\text{Δ}t}\left[C\right]{\left\{\delta \right\}}_t-\left[C\right]{\left\{\dot{\delta }\right\}}_t$(37)

or

$\left[B\right]{\left\{\delta \right\}}_{t+\text{Δ}t}={\left\{Q\right\}}_t+{\left\{F\right\}}_{t+\text{Δ}t}$ (38)

where

${\left\{Q\right\}}_t=\frac{2}{\text{Δ}t}\left[C\right]{\left\{\delta \right\}}_t-\left[C\right]{\left\{\dot{\delta }\right\}}_t$(39)

$\left[B\right]$ is called the time approximation operator, and ${\left\{Q\right\}}_t$ is called the load operator.

Knowing ${\left\{\delta \right\}}_t$ , ${\left\{\dot{\delta }\right\}}_t$ , and ${\left\{F\right\}}_{t+\text{Δ}t}$ , we can calculate ${\left\{\delta \right\}}_{t+\text{Δ}t}$ using (37), and then ${\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}$ using (34). Thus, the solution ${\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}$ and ${\left\{\delta \right\}}_{t+\text{Δ}t}$ is now known.

Initial conditions

1) Since (27) is a first-order system of ordinary differential equations, we require only one initial condition. If we begin with $t=0$ then ${\left\{\delta \right\}}_0$ is the initial condition (known). Using (27) at $t=0$ , we can write:

$\left[C\right]{\left\{\dot{\delta }\right\}}_0+\left[K\right]{\left\{\delta \right\}}_0={\left\{F\right\}}_0$ (40)

From (27), we can solve for ${\left\{\dot{\delta }\right\}}_0$ :

${\left\{\dot{\delta }\right\}}_0={\left[C\right]}^{-1}\left(-\left[K\right]{\left\{\delta \right\}}_0+{\left\{F\right\}}_0\right)$ (41)

Thus, now the complete solution is known at $t=0$ .

2) Using the initial condition in (1) at $t=0$ , we can use (34) at $t=0$ to calculate ${\left\{\delta \right\}}_{t+\text{Δ}t}$ and then use (33) to calculate ${\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}$

3) This procedure of time marching is continued until final time $t=\tau$ is reached.

Remarks

1) First, we note that for each time interval, we have an analytical expression for $\left\{\delta \left(t\right)\right\}$ and $\left\{\dot{\delta }\left(t\right)\right\}$ needed for residual computations.

2) Secondly, we note that at a time $t$ between the two time intervals, $\left\{\delta \right\}$ and $\left\{\dot{\delta }\right\}$ are continuous, implying that in this time integration scheme, $\delta \left(t\right)$ is of class $C^1$ as required by (14) to ensure that the time integrals in $L_2$ -norm are Riemann in time.

3.4.2. Newmark Time Integration Method for a System of Second Order Ordinary Differential Equations in Time

Consider a system of second-order linear ordinary differential equations in time:

$\left[M\right]\left\{\ddot{\delta }\right\}+\left[C\right]\left\{\dot{\delta }\right\}+\left[K\right]\left\{\delta \right\}=\left\{F\right\}$ (42)

These ordinary differential equations in time result when decoupling space and time in balance of linear momenta in linear elasticity and when constructing integral forms over spatial discretization using Galerkin method with weak form (GM/WF). $\left[M\right]$ , $\left[C\right]$ , and $\left[K\right]$ are mass, damping, and stiffness matrices, and $\left\{\delta \right\}$ , $\left\{\dot{\delta }\right\}$ , and $\left\{\ddot{\delta }\right\}$ are nodal degrees of freedom, their first and second time derivatives, when the degrees of freedom $\left\{\delta \right\}$ are displacements (C0 Lagrange family finite elements), then $\left\{\dot{\delta }\right\}$ and $\left\{\ddot{\delta }\right\}$ are velocities and accelerations of the grid points of the spatial discretization.

We consider Newmark’s linear acceleration method for integrating ordinary differential Equations (42) in time. In this method, we assume that acceleration $\left\{\ddot{\delta }\right\}$ is linear in the time interval $\left[t,t+\text{Δ}t\right]$ (Figure 3).

Figure 3. Acceleration $\left\{\ddot{\delta }\right\}$ versus time t (Newmark linear acceleration method).

Complete details of this method can be found in [1]. Important steps and the calculation procedure are described in the following. If we measure time $\tau$ from $t$ , then

${\left\{\ddot{\delta }\right\}}_{t+\tau }={\left\{\ddot{\delta }\right\}}_t+\frac{\tau }{\text{Δ}t}\left({\left\{\ddot{\delta }\right\}}_{t+\text{Δ}t}-{\left\{\ddot{\delta }\right\}}_t\right)$(43)

Integrating with respect to $\tau$ and using $\left\{\dot{\delta }\right\}$ at $\tau =0$ equal as ${\left\{\dot{\delta }\right\}}_t$ , to evaluate constant of integration, we obtain

${\left\{\dot{\delta }\right\}}_{t+\tau }=\tau {\left\{\ddot{\delta }\right\}}_t+{\left\{\dot{\delta }\right\}}_t+\frac{{\tau }^2}{2\text{Δ}t}\left({\left\{\ddot{\delta }\right\}}_{t+\text{Δ}t}-{\left\{\ddot{\delta }\right\}}_t\right)$(44)

Integrating with respect to $\tau$ and using ${\left\{\delta \right\}}_{t+\tau }$ at $\tau =0$ , as ${\left\{\delta \right\}}_t$ to evaluate the constant of integration, we obtain

${\left\{\delta \right\}}_{t+\tau }={\left\{\delta \right\}}_t+\frac{{\tau }^2}{2}{\left\{\ddot{\delta }\right\}}_t+\frac{{\tau }^3}{6\text{Δ}t}\left({\left\{\ddot{\delta }\right\}}_{t+\text{Δ}t}-{\left\{\ddot{\delta }\right\}}_t\right)$(45)

Using (43), (44), and (45), we obtain expressions for ${\left\{\delta \right\}}_{t+\text{Δ}t}$ , ${\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}$ and ${\left\{\ddot{\delta }\right\}}_{t+\text{Δ}t}$ by letting $\tau =\text{Δ}t$ . From these, we solve for ${\left\{\ddot{\delta }\right\}}_{t+\text{Δ}t}$ and ${\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}$ in terms of ${\left\{\delta \right\}}_t$ , ${\left\{\dot{\delta }\right\}}_t$ , ${\left\{\ddot{\delta }\right\}}_t$ , ${\left\{\delta \right\}}_{t+\text{Δ}t}$ and ${\left\{F\right\}}_{t+\text{Δ}t}$ and substitute these into (43) at $t+\text{Δ}t$ .

$\left[M\right]{\left\{\ddot{\delta }\right\}}_{t+\text{Δ}t}+\left[C\right]{\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}+\left[K\right]{\left\{\delta \right\}}_{t+\text{Δ}t}={\left\{F\right\}}_{t+\text{Δ}t}$(46)

We express resulting (46) as a system of linear simultaneous algebraic equations in unknown ${\left\{\delta \right\}}_{t+\text{Δ}t}$ (see [1] for details), Equation (47). Final expression used for computations are given in the following:

$\begin{array}{l}\left(\frac{6}{\text{Δ}{t}^2}\left[M\right]+\frac{3}{\text{Δ}t}\left[C\right]+\left[K\right]\right){\left\{\delta \right\}}_{t+\text{Δ}t}\\ ={\left\{F\right\}}_{t+\text{Δ}t}+\left(\frac{6}{\text{Δ}{t}^2}\left[M\right]+\frac{3}{\text{Δ}t}\left[C\right]\right){\left\{\delta \right\}}_t+\left(\frac{6}{\text{Δ}t}\left[M\right]+2\left[C\right]\right){\left\{\dot{\delta }\right\}}_t\\ +\left(2\left[M\right]+\frac{\text{Δ}t}{2}\left[C\right]\right){\left\{\ddot{\delta }\right\}}_t\end{array}$(47)

and

${\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}=\frac{3}{\text{Δ}t}\left({\left\{\delta \right\}}_{t+\text{Δ}t}-{\left\{\delta \right\}}_t\right)-2{\left\{\ddot{\delta }\right\}}_t-\frac{\text{Δ}t}{2}{\left\{\ddot{\delta }\right\}}_t$(48)

${\left\{\ddot{\delta }\right\}}_{t+\text{Δ}t}=\frac{6}{\text{Δ}{t}^2}\left({\left\{\delta \right\}}_{t+\text{Δ}t}-{\left\{\delta \right\}}_t\right)-\frac{6}{\text{Δ}t}{\left\{\dot{\delta }\right\}}_t-2{\left\{\ddot{\delta }\right\}}_t$(49)

We can write (47) as

$\left[B\right]{\left\{\delta \right\}}_{t+\text{Δ}t}={\left\{Q\right\}}_t+{\left\{F\right\}}_{t+\text{Δ}t}$

where

$\begin{array}{c}{\left\{Q\right\}}_t=\left(\frac{6}{\text{Δ}{t}^2}\left[M\right]+\frac{3}{\text{Δ}t}\left[C\right]\right){\left\{\delta \right\}}_t+\left(\frac{6}{\text{Δ}t}\left[M\right]+2\left[C\right]\right){\left\{\dot{\delta }\right\}}_t\\ +\left(2\left[M\right]+\frac{\text{Δ}t}{2}\left[C\right]\right){\left\{\ddot{\delta }\right\}}_t\end{array}$(51)

Initial Conditions:

Since (42) is a system of linear second-order ordinary differential equations in time, it requires two initial conditions. If we assume that the evolution commences at $t=0$ , then we could have ${\left\{\delta \right\}}_0$ and ${\left\{\dot{\delta }\right\}}_0$ as initial conditions. Using (42) at $t=0$ , $\left[M\right]{\left\{\ddot{\delta }\right\}}_0+\left[C\right]{\left\{\dot{\delta }\right\}}_0+\left[K\right]{\left\{\delta \right\}}_0={\left\{F\right\}}_0$ (52)

From (52), we can determine ${\left\{\ddot{\delta }\right\}}_0$ .

${\left\{\ddot{\delta }\right\}}_0={\left[M\right]}^{-1}\left({\left\{f\right\}}_0-\left[C\right]{\left\{\dot{\delta }\right\}}_0-\left[K\right]\left\{{\delta }_0\right\}\right)$

Thus, ${\left\{\delta \right\}}_0$ , ${\left\{\dot{\delta }\right\}}_0$ , and ${\left\{\ddot{\delta }\right\}}_0$ are all known at time $t_0=0$ .

Remarks

1) If the evolution commences at $t=0$ , then ${\left\{\delta \right\}}_0$ , ${\left\{\dot{\delta }\right\}}_0$ , and ${\left\{\ddot{\delta }\right\}}_0$ are all known at $t=0$ .

2) We use $t=0$ in (45) to calculate ${\left\{\delta \right\}}_{\text{Δ}t}$ followed by the calculation of ${\left\{\dot{\delta }\right\}}_{t+\text{Δ}t}$ and ${\left\{\ddot{\delta }\right\}}_{t+\text{Δ}t}$ using (48) and (49). Thus, the solution at $t=\text{Δ}t$ is now known.

3) Using the solution at $t=\text{Δ}t$ , the solution at $t=2\text{Δ}t$ i.e. ${\left\{\delta \right\}}_{2\text{Δ}t}$ , ${\left\{\dot{\delta }\right\}}_{2\text{Δ}t}$, and ${\left\{\ddot{\delta }\right\}}_{2\text{Δ}t}$ can be calculated using (47) - (49).

4) This time marching process is continued until the final time $t=\tau$ is reached.

3.5. A Posteriori Error Computations in Space-Time Decoupled Methods for a First-Order System of Linear Ordinary Differential Equations in Time

Consider the following system of linear first order ordinary differential equations in time:

$\left[C\right]\left\{\dot{\delta }\right\}+\left[H\right]\left\{\delta \right\}=\left\{F\right\}$ (53)

Obtained by decoupling space and time in the initial value problem using a finite element discretization ${\left({\overline{\Omega }}_{xt}\right)}^T={\displaystyle \underset{e}{\cup }{\overline{\Omega }}_x^e}$ in space of the spatial domain ${\overline{\Omega }}_x\subset {ℝ}^1$ or ${ℝ}^2$ or ${ℝ}^3$ ; ${\overline{\Omega }}_x^e$ could be a line, an area or a volume finite element in space,with p-version hierarchical with higher orders global differentiability in space.

For illustration, consider the following time dependent diffusion equation

$c\frac{\partial \varphi }{\partial t}-k\frac{{\partial }^2\varphi }{\partial x^2}-f=0,\forall x,t\in {\Omega }_x\times {\Omega }_t=\left(0,1\right)\times \left(0,\tau \right)$ (54)

Consider discretization of ${\overline{\Omega }}_x^T$ of the spatial domain ${\overline{\Omega }}_x$

${\overline{\Omega }}_x^T={\displaystyle \underset{e}{\cup }{\overline{\Omega }}_x^e}$ (55)

${\overline{\Omega }}_x^e$ is a three-node p-version hierarchical finite element in the spatial domain $x$ with higher-order global differentiability local approximation. Let ${\varphi }_h^e\left(x,t\right)$ be the approximation of $\varphi$ over ${\overline{\Omega }}_x^e$ , then in space-time decoupled finite element method, we write:

${\varphi }_h^e\left(x,t\right)=\left[N\left(x\right)\right]\left\{{\delta }^e\left(t\right)\right\}$ (56)

in which ${\delta }^e\left(t\right)$ is the computed solution for time integration method for an increment of time.

We substitute (56) in (54) to obtain residual function $E_{xt}^e$ . ${\underset{˜}{\overline{\Omega }}}_{xt}^e={\overline{\Omega }}_x^e\times \left[\left(n-1\right)t\le t\le nt\right]$ , for the $n^{th}$ time interval

$E_{xt}^e\left(x,t\right)=C\left[N\left(x\right)\right]\left\{\dot{\delta }\left(t\right)\right\}-k\left[\frac{{\partial }^{2}N\left(x\right)}{\partial x^2}\right]\left\{\delta \left(t\right)\right\}-f\left(t\right)$ (57)

and the residual functional $I_{xt}^e$ is given by

$I_{xt}^e={\left(E_{xt}^e\left(x,t\right),E_{xt}^e\left(x,t\right)\right)}_{{\underset{˜}{\overline{\Omega }}}_{xt}^e}$ (58)

and for the the whole spatial discretization

$I_{xt}={\displaystyle \sum }_e{I}_{xt}^e$ (59)

We note that $N\left(x\right)$ are approximation functions for ${\overline{\Omega }}_x^e$ mapped in natural coordinate space $\xi$ . Thus,

$E_{\xi t}^e\left(\xi ,t\right)=C\left(\left[N\left(\xi \right)\right]\right)\left\{\dot{\delta }\left(t\right)\right\}+K\left[\frac{{\partial }^{2}N\left(\xi \right)}{\partial x^2}\right]\left\{\delta \left(t\right)\right\}-f\left(t\right)$ (60)

$I_{\xi t}={\displaystyle \sum }_e{I}_{\xi t}^e={\displaystyle \sum }_e{\left(E_{\xi t}^e\left(\xi ,t\right),E_{\xi t}^e\left(\xi ,t\right)\right)}_{{\underset{˜}{\overline{\Omega }}}_{\xi t}^e}$(61)

We map time $t$ for an increment of time into natural coordinate $\nu$ in a two-unit length. Consider the $j^{th}$ time increment $\left[t^j,t^j+\text{Δ}t\right]=\left[t_1^j,t_2^j\right]$ mapped into a two-unit length $\left[{\nu }_1,{\nu }_2\right]=\left[-1,1\right]$ . We can write the following:

$t\left(\nu \right)=\left(\frac{1-\nu }{2}\right)t_1^j+\left(\frac{1+\nu }{2}\right)t_2^j$ (62)

and

$\tau \left(\nu \right)=\left(\frac{1+\nu }{2}\right)\text{Δ}t\text{and}\text{d}\tau =\frac{\text{Δ}t}{2}\text{d}\nu$ (63)

time $\tau$ is measured from $t=t_1^j$ .

Recall ${\left\{\delta \right\}}_{t+\tau }$ and ${\left\{\dot{\delta }\right\}}_{t+\tau }$ (from Equations (29) and (28)), we can express these as functions of $\nu$ :

${\left\{\delta \right\}}_{t_1+\tau \left(\nu \right)}=\tau \left(\nu \right){\left\{\dot{\delta }\right\}}_{t_1}^j+\frac{{\left(\tau \left(\nu \right)\right)}^2}{2\text{Δ}t}\left({\left\{\dot{\delta }\right\}}_{t_1^j+\text{Δ}t}-{\left\{\dot{\delta }\right\}}_{t_1^j}\right)$(64)

and

${\left\{\dot{\delta }\right\}}_{t_1^j+\tau \left(\nu \right)}={\left\{\dot{\delta }\right\}}_{t_1^j}+\frac{\tau \left(\nu \right)}{\text{Δ}t}\left({\left\{\dot{\delta }\right\}}_{t_1^j+\text{Δ}t}-{\left\{\dot{\delta }\right\}}_{t_1^j}\right)$(65)

where

$t_2^j=t_1^j+\text{Δ}t$ (66)

Using (62) in (61) , we obtain $E_{\xi \nu }^e\left(\xi ,\nu \right)$ and:

$I_{\xi \nu }={\displaystyle \sum }_e{I}_{\xi \nu }^e={\displaystyle \sum }_e{\left(E_{\left(\xi \nu \right)}^e\left(\xi ,\nu \right),E_{\xi \nu }^e\left(\xi ,\nu \right)\right)}_{{\underset{˜}{\overline{\Omega }}}_{\xi \nu }^e}$(67)

$I_{\xi \nu }^e$ is calculated using Gauss quadrature.

Proximity of $I_{\xi \nu }^e$ to zero is a measure of the accuracy of the solution in ${\underset{˜}{\overline{\Omega }}}_{\xi \nu }^e$ i.e. element e in space for a time increment and proximity of $I_{\xi \nu }$ to zero is a measure of accuracy of the solution over the whole spatial discretization for an increment of time.

3.6. A Posteriori Error Computations in Space-Time Decoupled Methods for Second-Order Systems of Linear Ordinary Differential Equations in Time

In linear elasticity, the balance of linear momenta is a system of second order partial differential equations in space and time in displacement when the constitutive theory for stresses is substituted in them. We discussed space-time decoupled method for obtaining solutions of such initial value problems.

Consider the initial value problem in linear elasticity, axial deformation of rod made of viscoelastic material. The rod has elasticity and dissipation mechanisms. Balance of linear momenta and constitutive theory are given by (neglecting equilibrium stress):

${\rho }_0\frac{{\partial }^{2}u}{\partial t^2}-{\rho }_0{F}_1^b-\frac{\partial {\sigma }_{11}}{\partial x}=0\forall x\in {\Omega }_x\subset {ℝ}^1$ (68)

${\sigma }_{11}=E\frac{\partial u}{\partial x}+c\frac{\partial }{\partial t}\left(\frac{\partial u}{\partial x}\right)$ (69)

By substituting (69) in (68), we obtain a single partial differential equation in $u$ :

${\rho }_0\frac{{\partial }^{2}u}{\partial t^2}-{\rho }_0{F}_1^b-E\frac{{\partial }^{2}u}{\partial x^2}-c\frac{{\partial }^{2}v}{\partial x^2}=0$ (70)

or $Au-f=0$ in which $v$ is velocity.

Consider space-time decoupled formulation of (70) with GM/WF in space followed by time integration of the resulting second order ordinary differential equations in time.

Let ${\overline{\Omega }}_x^T={\displaystyle \underset{e}{\cup }{\overline{\Omega }}_x^e}$ be discretization of ${\overline{\Omega }}_x$ in which ${\overline{\Omega }}_x^e$ is a finite element in space (1D in this case).

Let $u_h^e\left(x,t\right)={\displaystyle \sum }_{i=1}^{n}N\left(x\right){\delta }_i^e\left(t\right)$ be local approximation of $u$ over ${\overline{\Omega }}_x^e$ and $u_h\left(x,t\right)$ be approximation of $u$ over ${\overline{\Omega }}_x^T$ , then

$u_h\left(x,t\right)={\displaystyle \underset{e}{\cup }{u}_h^e\left(x,t\right)}$ (71)

Using fundamental lemma of the calculus of variations (using test function $w$ ),

${\left(A{u}_h-f,w\right)}_{{\overline{\Omega }}_x^T}=0;w=\delta u_h\text{in}\text{GM/WF}$(72)

${\left(A{u}_h-f,w\right)}_{{\overline{\Omega }}_x^T}={\displaystyle \sum }_e{\left(A{u}_h^e-f,w\right)}_{{\overline{\Omega }}_x^e}=0$ (73)

and $\left\{\delta \right\}={\displaystyle \underset{e}{\cup }\left\{{\delta }_e\right\}}$ are the degrees of freedom for ${\overline{\Omega }}_x^T$

Consider

$\begin{array}{l}\left(A{u}_h^e-f,w\right)={\left({\rho }_0\frac{{\partial }^2{u}_h^e}{\partial t^2}-{\rho }_0{F}_1^b-E\frac{{\partial }^2{u}_h^e}{\partial x^2}-c\frac{{\partial }^2{v}_h^e}{\partial x^2},w\right)}_{{\overline{\Omega }}_x^e}\\ =\underset{{\overline{\Omega }}^e}{{\displaystyle \int }}{\rho }_0\frac{{\partial }^2{u}_h^e}{\partial t^2}w\text{d}\Omega -\underset{{\overline{\Omega }}^e}{{\displaystyle \int }}{\rho }_0{F}_1^{b}w\text{d}\Omega -\underset{{\overline{\Omega }}^e}{{\displaystyle \int }}E\frac{{\partial }^2{u}_h^e}{\partial x^2}w\text{d}\Omega -\underset{{\overline{\Omega }}^e}{{\displaystyle \int }}c\frac{{\partial }^2{v}_h^e}{\partial x^2}w\text{d}\Omega \end{array}$(74)

We perform integration by part once in the last two term

$\begin{array}{l}=\underset{{\overline{\Omega }}_x^e}{{\displaystyle \int }}{\rho }_0\frac{{\partial }^2{u}_h^e}{\partial t^2}w\text{d}\Omega -\underset{{\overline{\Omega }}_x^e}{{\displaystyle \int }}{\rho }_0{F}_1^{b}w\text{d}\Omega +\underset{{\overline{\Omega }}^e}{{\displaystyle \int }}E\frac{\partial w}{\partial x}\frac{\partial u}{\partial x}w\text{d}\Omega \\ +\underset{{\overline{\Omega }}_x^e}{{\displaystyle \int }}c\frac{\partial w}{\partial x}\frac{\partial v_h^e}{\partial x}w\text{d}\Omega -{\left[w\left(E\frac{\partial u_h^e}{\partial x}\right)\right]}_{x_e}^{x_{e+1}}-{\left[wc\frac{\partial v_h^e}{\partial x}\right]}_{x_e}^{x_{e+1}}\end{array}$(75)

Substituting (68), $w=\delta u_h^e=N_j;j=1,2,\cdots ,n$ and $v_h^e={\displaystyle \sum }_{i=1}^n{N}_i\left(x\right){\dot{\delta }}_i\left(t\right)$ in (75). Boundary terms lead to secondary variables (standard procedure) $\left\{P^e\right\}$ for an element $e$ .

$\begin{array}{l}{\left(A{u}_h^e-f,w\right)}_{{\overline{\Omega }}_x^e}=\underset{{\overline{\Omega }}^e}{{\displaystyle \int }}{\rho }_0\left({\displaystyle \sum }_{i=1}^n{N}_i{\ddot{\delta }}_i^e\right)N_j\text{d}\Omega -\underset{{\overline{\Omega }}_x^e}{{\displaystyle \int }}{\rho }_0{F}_1^b{N}_j\text{d}\Omega \\ -\underset{{\overline{\Omega }}_x^e}{{\displaystyle \int }}E\frac{\partial N_j}{\partial x}{\displaystyle \sum }_{i=1}^n\frac{\partial N_i}{\partial x}{\delta }_i^e\text{d}\Omega +\underset{{\overline{\Omega }}_x^e}{{\displaystyle \int }}c\frac{\partial N_j}{\partial x}{\displaystyle \sum }_{i=1}^n\frac{\partial N_i}{\partial x}{\dot{\delta }}_i^e\text{d}\Omega \\ -{\left[N_j\left(E\frac{\partial u_h^e}{\partial x}\right)\right]}_{x_e}^{x_{e+1}}-{\left[N_{j}c\frac{\partial v_h^e}{\partial x}\right]}_{x_e}^{x_{e+1}},j=1,2,\cdots ,n\end{array}$(76a)

in which

$M_{ij}^e={\rho }_0\underset{{\overline{\Omega }}_x^e}{{\displaystyle \int }}{N}_j{N}_i\text{d}\Omega ;C_{ij}^e=\underset{{\overline{\Omega }}_x^e}{{\displaystyle \int }}\frac{\partial N_j}{\partial x}\frac{\partial N_i}{\partial x}\text{d}\Omega ;K_{ij}^e=\underset{{\overline{\Omega }}_x^e}{{\displaystyle \int }}E\frac{\partial N_j}{\partial x}\frac{\partial N_i}{\partial x}\text{d}\Omega$(76b)

We can write (76a) in the matrix and vector form

${\left(A{u}_h^e-f,w\right)}_{{\overline{\Omega }}_x^e}=\left[M^e\right]\left\{{\ddot{\delta }}^e\right\}+\left[C^e\right]\left\{{\dot{\delta }}^e\right\}+\left[K^e\right]\left\{{\delta }^e\right\}-\left\{F^e\right\}-\left\{P^e\right\}$ (77)

Using (77) in (73) i.e. assembly of element expressions we obtain:

$\begin{array}{l}{\left(A{u}_h-f,v\right)}_{{\overline{\Omega }}_x^T}={\displaystyle \sum }_e{\left(A{u}_h^e-f,v\right)}_{{\overline{\Omega }}_x^e}\\ ={\displaystyle \sum }_e\left[M^e\right]\left\{{\ddot{\delta }}^e\right\}+{\displaystyle \sum }_e\left[C^e\right]\left\{{\dot{\delta }}^e\right\}+{\displaystyle \sum }_e\left[K^e\right]\left\{{\delta }^e\right\}-{\displaystyle \sum }_e\left\{F^e\right\}-{\displaystyle \sum }_e\left\{P^e\right\}=0\end{array}$ (78)

or

$\left[M\right]\left\{\ddot{\delta }\right\}+\left[C\right]\left\{\dot{\delta }\right\}+\left[K\right]\left\{\delta \right\}=\left\{F\right\}+\left\{P\right\}$ (79)

in which

$\begin{array}{l}\left[M\right]={\displaystyle \sum }_e\left[M^e\right],\left[C\right]={\displaystyle \sum }_e\left[C^e\right],\left[K\right]={\displaystyle \sum }_e\left[K^e\right],\left\{F\right\}={\displaystyle \sum }_e\left\{F^e\right\},\\ \left\{P\right\}={\displaystyle \sum }_e\left\{P^e\right\},\left\{\delta \right\}={\displaystyle \underset{e}{\cup }\left\{{\delta }^e\right\}},\left\{\dot{\delta }\right\}={\displaystyle \underset{e}{\cup }\left\{{\dot{\delta }}^e\right\}}\text{and}\left\{\ddot{\delta }\right\}={\displaystyle \underset{e}{\cup }\left\{{\ddot{\delta }}^e\right\}}\end{array}$(80)

Equations (79) are a system of second-order ordinary differential equations in time.

Let ${\overline{\Omega }}_x^e$ be a 3-node p-version hierarchical element in space with higher order global differentiability

$u_h^e\left(x,t\right)={\displaystyle \sum }_{i=1}^n{N}_i\left(x\right){\delta }_i^e\left(t\right)$ (81)

We substitute (81) in (68) to obtain residual function $E_{xt}^e$ $\forall$ $x,t\in {\underset{˜}{\overline{\Omega }}}_{xt}^e={\overline{\Omega }}_x^e\times \left[\left(n-1\right)t,nt\right]$ , for the $n^{th}$ space-time strip:

$E_{xt}^e={\rho }_0\left[N\left(x\right)\right]\left\{{\ddot{\delta }}^e\right\}-c\left[\frac{{\partial }^{2}N}{\partial x^2}\right]\left\{{\dot{\delta }}^e\right\}-E\left[\frac{{\partial }^{2}N}{\partial x^2}\right]\left\{{\delta }^e\right\}-{\rho }_0{F}_1^b$ (82)

and the residual functional $I^e$ for ${\underset{˜}{\overline{\Omega }}}_{xt}^e$ is given by

$I_{xt}^e={\left(E_{xt}^e,E_{xt}^e\right)}_{{\underset{˜}{\overline{\Omega }}}_{xt}^e}$ (83)

and for the whole spatial discretization and time $t\in \left[\left(n-1\right)t,nt\right]$ ,

$I_{xt}={\displaystyle \sum }_e{I}_{xt}^e$ (84)

We note that $N_i\left(x\right)$ are approximation functions for ${\overline{\Omega }}_x^e$ mapped with natural coordinate space $\xi \in \left[-1,1\right]$ . Thus,

$E_{\xi t}^e={\rho }_0\left[N\left(\xi \right)\right]\left\{{\ddot{\delta }}^e\left(t\right)\right\}-c\left[\frac{{\partial }^{2}N\left(\xi \right)}{\partial x^2}\right]\left\{{\dot{\delta }}^e\left(t\right)\right\}-E\left[\frac{{\partial }^{2}N\left(\xi \right)}{\partial x^2}\right]\left\{{\delta }^e\left(t\right)\right\}-{\rho }_0{F}^b$ (85)

and

$I_{\xi t}={\displaystyle \sum }_e{I}_{\xi t}^e={\displaystyle \sum }_e{\left(E_{\xi t}^e\left(\xi ,t\right),E_{\xi t}^e\left(\xi ,t\right)\right)}_{{\underset{˜}{\overline{\Omega }}}_{\xi t}^e}$(86)

We map $t$ with natural coordinate space $\nu$ in a two-unit length. Consider the $j^{th}$ time increment $\left[t^j,t^j+\text{Δ}t\right]=\left[t_1^j,t_2^j\right]$ mapped into a two-unit length $\left[{\nu }_1,{\nu }_2\right]=\left[-1,1\right]$ , we can write the following:

$t\left(\nu \right)=\left(\frac{1-\nu }{2}\right)t_1^j+\left(\frac{1+\nu }{2}\right)t_2^j$ (87)

and

$\tau \left(\nu \right)=\left(\frac{1+\nu }{2}\right)\text{Δ}t\text{and}\text{d}\tau =\frac{\text{Δ}t}{2}\text{d}\nu$ (88)

Recall ${\left\{\delta \right\}}_{t+\tau }$ , ${\left\{\dot{\delta }\right\}}_{t+\tau }$ and ${\left\{\ddot{\delta }\right\}}_{t+\tau }$ , (Equation (43), (44) and (45)) we can write these as functions of $\nu$ using (88) for $\tau \left(\nu \right)$ :

${\left\{\delta \right\}}_{t_1^j+\tau \left(\nu \right)}={\left\{\delta \right\}}_{t_1^j}+\frac{{\tau }^2\left(\nu \right)}{2}{\left\{\dot{\delta }\right\}}_{t_1}^j+\frac{{\left(\tau \left(\nu \right)\right)}^3}{6\text{Δ}t}\left({\left\{\ddot{\delta }\right\}}_{t_1^j+\text{Δ}t}-{\left\{\ddot{\delta }\right\}}_{t_1^j}\right)$(89)

${\left\{\dot{\delta }\right\}}_{t_1^j+\tau \left(\nu \right)}=\tau {\left\{\ddot{\delta }\right\}}_{t_1}^j+\frac{{\left(\tau \left(\nu \right)\right)}^2}{2\text{Δ}t}\left({\left\{\ddot{\delta }\right\}}_{t_1^j+\text{Δ}t}-{\left\{\ddot{\delta }\right\}}_{t_1^j}\right)$(90)

${\left\{\ddot{\delta }\right\}}_{t_1^j+\tau \left(\nu \right)}={\left\{\ddot{\delta }\right\}}_{t_1^j}+\frac{\tau \left(\nu \right)}{\text{Δ}t}\left({\left\{\ddot{\delta }\right\}}_{t_1^j+\text{Δ}t}-{\left\{\ddot{\delta }\right\}}_{t_1^j}\right)$(91)

in which

$t_2^j=t_1^j+\text{Δ}t$ (92)

Using (86) - (88) in (82), we can obtain $E_{\xi \nu }^e\left(\xi ,\nu \right)$ and

$I_{\xi \nu }={\displaystyle \sum }_e{I}_{\xi \nu }^e={\displaystyle \sum }_e{\left(E_{\left(\xi \nu \right)}^e\left(\xi ,\nu \right),E_{\xi \nu }^e\left(\xi ,\nu \right)\right)}_{{\underset{˜}{\overline{\Omega }}}_{\xi \nu }^e}$(93)

$I_{\xi \nu }^e$ is calculated using Gauss quadrature.

Proximity of $I_{\xi \nu }^e$ to zero is a measure of the accuracy of the solution in ${\underset{˜}{\overline{\Omega }}}_{\xi \gamma }^e$ i.e. for an element in space and for a time increment how close $I_{\xi \nu }$ is to zero is a measure of the accuracy of the solution or the error in the solution for the spatial discretization for an increment of time.

4. Model Problem Studies

In this section, we consider two model problems: Transient Convection Diffusion Equation and Transient Pure Advection. Evolutions for both problems are obtained using space-time coupled methods as well as space-time decoupled methods. In the space-time coupled method, converged solutions are obtained for the first space-time step, and then time is marched to obtain the evolution up to desired values of time. Solutions of class $C^1$ and $C^2$ are considered for both in space and time. In the case of space-time decoupled method, integral form in space is constructed using Galerkin method with weak form. Element equations in space are constructed and assembled to yield a system of ordinary differential equations in time that are integrated using Newmark’s linear method. Influence of p-levels in space and time as well as the order of the approximation space are investigated in the case of space-time coupled finite element methods. In space-time decoupled methods, the influence of p-levels and the order of the approximation space in space are investigated in conjunction with varying integration time steps in Newmark’s linear method.

4.1. Transient Convection Diffusion Equation

Consider a 1D transient convection diffusion equation given by:

$\frac{\partial \varphi }{\partial t}+\frac{\partial \varphi }{\partial x}-\frac{1}{Pe}\frac{{\partial }^2\varphi }{\partial x^2}=f\left(x,t\right)\forall x,t\in {\Omega }_{xt}={\Omega }_x\times {\Omega }_t=\left(0,1\right)\times \left(0,\tau \right)$ (94)

where Pe is Peclet number.

with the following boundary conditions and initial conditions:

$\{\begin{array}{l}\text{BCs}:\varphi \left(0,t\right)=0,{\frac{\partial \varphi }{\partial x}|}_{x=0}=0\forall t\in {\overline{\Omega }}_t=\left[0,\tau \right],\\ \text{ICs}:\varphi \left(x,0\right)=g\left(x\right)\forall x\in {\overline{\Omega }}_x=\left[0,1\right].\end{array}$(95)

$g\left(x\right)$ defines the initial condition on $\varphi \left(x,0\right)$ at $t=0$ . We can write (94) as:

$A\varphi -f=0\forall x,t\in {\Omega }_{x,t};A=\frac{\partial }{\partial t}+\frac{\partial }{\partial x}-\frac{1}{Pe}\frac{{\partial }^2}{\partial x^2}$ (96)

4.1.1. Space-Time Coupled Finite Element Method Using Space-Time Residual Functional

Let be a discretization of ${\overline{\Omega }}_{xt}^T=\left[0,1\right]\times \left[0,\Delta t\right]$ ${\overline{\Omega }}_{xt}={\Omega }_{xt}\cup \Gamma$ , $\text{Γ}$ is boundary of ${\Omega }_{xt}$ and ${\overline{\Omega }}_{xt}^e$ is a space-time element of the discretization ${\overline{\Omega }}_{xt}^T$ of the first space-time strip. Let ${\varphi }_h^e\left(x,t\right)$ be the local approximation of $\varphi \left(x,t\right)$ over ${\overline{\Omega }}_{xt}^e$ and let ${\varphi }_h\left(x,t\right)$ be global approximation of $\varphi \left(x,t\right)$ over ${\overline{\Omega }}_{xt}^T$ , such that:

${\varphi }_h\left(x,t\right)={\displaystyle \underset{e}{\cup }{\varphi }_h^e}\left(x,t\right);{\varphi }_h^e\left(x,t\right)={\displaystyle \sum }_{i=1}^n{N}_i\left(x,t\right){\delta }_i^e;\left\{\delta \right\}={\displaystyle \underset{e}{\cup }\left\{{\delta }^e\right\}}$ (97)

$\left\{{\delta }^e\right\}$ are degrees of freedom for ${\overline{\Omega }}_{xt}^e$ and $\left\{\delta \right\}$ are degrees of freedom for ${\overline{\Omega }}_{xt}^T$

The residual function $E^e$ for ${\overline{\Omega }}_{xt}^e$ is given by:

$E^e=A{\varphi }_h^e-f=\frac{\partial {\varphi }_h^e}{\partial t}+\frac{\partial {\varphi }_h^e}{\partial x}-\frac{1}{Pe}\frac{{\partial }^2{\varphi }_h^e}{\partial x^2}\forall \left(x,t\right)\in {\overline{\Omega }}_{xt}^e$(98)

and space-time functional $I$ for discretization ${\overline{\Omega }}_{xt}^T$ is defined by:

$I^e={\displaystyle \sum }_e{I}^e={\displaystyle \sum {\left(E^e,E^e\right)}_{{\overline{\Omega }}_x^e}}\forall \left(x,t\right)\in {\overline{\Omega }}_{xt}^T$(99)

in which $I^e$ is the space-time residual functional over ${\overline{\Omega }}_{xt}^e$ , space-time domain of space-time element $e$ .

$\{\begin{array}{l}{I}^e=I^e\left\{{\delta }^e\right\},\\ I=I\left\{\delta \right\}.\end{array}$ (100)

If $I(\cdot )$ is differentiable in its arguments, then $\delta I$ is unique and $\delta I=0$ is a necessary condition for extremum of $I$ :

$\delta I={\displaystyle \sum }_e\delta I^e={\displaystyle \sum }_{e}2\left(E^e,\delta E^e\right)={\displaystyle \sum }_e\left\{g^e\right\}=\left\{g\right\}=0$ (101)

or

${\delta }^{2}I={\displaystyle \sum }_{e}2\left(\delta E^e,\delta E^e\right)>0\forall \delta E^e$ (102)

then (102) is unique extremum principle, and the space-time integral form (100) is space-time variationally consistent.

$\delta I^e=2\left(E^e,\delta E^e\right)=\left[K^e\right]\left\{{\delta }^e\right\}-\left[F^e\right]$ (103)

in which

$K_{ij}^e={\left(\frac{\partial N_j}{\partial t}+\frac{\partial N_j}{\partial x}-\frac{1}{Pe}\frac{{\partial }^2{N}_j}{\partial x^2},\frac{\partial N_i}{\partial t}+\frac{\partial N_i}{\partial x}-\frac{1}{Pe}\frac{{\partial }^2{N}_i}{\partial x^2}\right)}_{{\overline{\Omega }}_{xt}^e};i,j=1,2,\cdots ,$(104)

and

$F_j^e={\left(f\left(x,t\right),\frac{\partial N_j}{\partial t}+\frac{\partial N_j}{\partial x}-\frac{1}{Pe}\frac{{\partial }^2{N}_j}{\partial x^2}\right)}_{{\overline{\Omega }}_{xt}^e};j=1,2,\cdots ,$(105)

$\therefore \delta I={\displaystyle \sum }_e\delta I^e=\left[K\right]\left\{\delta \right\}-\left\{F\right\}=0$ (106)

in which

$\left[K\right]={\displaystyle \sum }_e\left[K^e\right];\left\{F\right\}={\displaystyle \sum }_e\left\{F^e\right\}$ (107)

Equation (106) is the algebraic system for the first (or any) space-time strip. We calculate $\left\{\delta \right\}$ from (106) after imposing boundary conditions and initial conditions. In the numerical studies, choice of ${\varphi }_h^e\left(x,t\right)$ of class $C^1$ in space and time would yield Lebesgue integral in space but Riemann in time, whereas ${\varphi }_h^e\left(x,t\right)$ of class $C^2$ in space and time would ensure that all space-time integrals over ${\overline{\Omega }}_{xt}^T$ are Riemann. In the numerical computations we consider both $C^1$ as well as $C^2$ approximations in space and time. ${\Vert E\Vert }_{L_2}$ for ${\overline{\Omega }}_{xt}^T$ is given by ${\Vert E\Vert }_{L_2}=\sqrt{I}=\sqrt{{\displaystyle \sum }_e{\left(E^e,E^e\right)}_{{\overline{\Omega }}_{xt}^e}}$ in which $E^e$ is given by (98).

4.1.2. Space-Time Decoupled Finite Element Method

Let ${\overline{\Omega }}_x^T={\displaystyle \underset{e}{\cup }{\overline{\Omega }}_x^e}$ be a discretization of ${\overline{\Omega }}_x=\left[0,1\right]$ in which ${\overline{\Omega }}_x^e$ is a three-node p-version hierarchical 1D element in space with higher-order global differentiability.

Let

${\varphi }_h^e\left(x,t\right)={\displaystyle \sum }_{i=1}^n{N}_i\left(x\right){\delta }_i^e\left(t\right)$ (108)

be the local approximation of $\varphi \left(x,t\right)$ over ${\overline{\Omega }}_x^e$ , in which $N_i\left(x\right)$ are the p-version hierarchical local approximations with higher-order global differentiability in space $x$ and ${\delta }_i^e\left(t\right)$ are nodal degrees of freedom that are functions of time. Equation (108) causes decoupling of space and time. Using (108), we can write:

${\varphi }_h\left(x,t\right)={\displaystyle \underset{e}{\cup }{\varphi }_h^e\left(x,t\right)}\text{and}\left\{\delta \left(t\right)\right\}={\displaystyle \underset{e}{\cup }\left\{{\delta }^e\left(t\right)\right\}}$ (109)

in which $\left\{{\delta }^e\left(t\right)\right\}$ are the degrees of freedom for an element ${\overline{\Omega }}_x^e$ and $\left\{\delta \left(t\right)\right\}$ are the degrees of freedom for the discretization ${\overline{\Omega }}_x^T$ of the spatial domain ${\overline{\Omega }}_x$ . Consider the integral form of (94) over ${\overline{\Omega }}_x^T$ using Galerkin methods with weak form:

${\left(A{\varphi }_h-f,v\left(x\right)\right)}_{{\overline{\Omega }}_x^T}={\displaystyle \sum }_e{\left(A{\varphi }_h^e-f,v\left(x\right)\right)}_{{\overline{\Omega }}_x^e}=0$ (110)

Consider ${\left(A{\varphi }_h^e-f,v\right)}_{{\overline{\Omega }}_x^e}$ in which $v\left(x\right)=\delta {\varphi }_h^e$ (Galerkin method)

$\begin{array}{c}{\left(A{\varphi }_h^e-f,v\right)}_{{\overline{\Omega }}_x^e}={\left(\frac{\partial {\varphi }_h^e}{\partial t}+\frac{\partial {\varphi }_h^e}{\partial x}-\frac{1}{Pe}\frac{{\partial }^2{\varphi }_h^e}{\partial x^2}-f,v\left(x\right)\right)}_{{\overline{\Omega }}_x^e}\\ =\underset{{\overline{\Omega }}_x^e}{{\displaystyle \int }}\left(\frac{\partial {\varphi }_h^e}{\partial t}v+\frac{\partial {\varphi }_h^e}{\partial x}v-\frac{1}{Pe}\frac{{\partial }^2{\varphi }_h^e}{\partial x^2}v-fv\right)\text{d}x\end{array}$(111)

We perform integration by parts in the third term to transfer one order of differentiation from ${\varphi }_h^e$ to $v$

$\left(A{\varphi }_h^e-f,v\right)=\underset{{\overline{\Omega }}_x^e}{{\displaystyle \int }}\left(\frac{\partial {\varphi }_h^e}{\partial t}v+\frac{\partial {\varphi }_h^e}{\partial x}v+\frac{1}{Pe}\frac{\partial v}{\partial x}\frac{\partial {\varphi }_h^e}{\partial x}-fv\right)\text{d}x-{\left[v\left(x\right)\left(\frac{1}{Pe}\frac{\partial {\varphi }_h^e}{\partial x}\right)\right]}_{x_e}^{x_{e+1}}$ (112)

$x_e$ and $x_{e+1}$ are the coordinates of the end points of the element.

Using (108) and

$v\left(x\right)=\delta {\varphi }_h^e\left(x,t\right)=N_j\left(x\right);j=1,2,\cdots ,n.$ (113)

in (112), we obtain:

$\begin{array}{l}{\left(A{\varphi }_h^e-f,v\right)}_{{\overline{\Omega }}_x^e}=\underset{{\overline{\Omega }}_x^e}{{\displaystyle \int }}({\displaystyle \sum }_{i=1}^n{N}_i\left(x\right){\dot{\delta }}_i^e\left(t\right)N_j\left(x\right)+{\displaystyle \sum }_{i=1}^n\frac{\partial N_i}{\partial x}{\delta }_i^e\left(t\right)N_j\left(x\right)\\ +\frac{1}{Pe}\frac{\partial N_j}{\partial x}\frac{\partial N_i}{\partial x}{\delta }^e\left(t\right)-f{N}_j\left(x\right))\text{d}x\\ -{\left[N_j\frac{1}{Pe}\frac{\partial {\varphi }_h^e}{\partial x}\right]}_{x_e}^{x_{e+1}};j=1,2,\cdots ,n.\end{array}$(114)

Defining:

$P_1^e=-\frac{1}{Pe}{\left(\frac{\partial {\varphi }_h^e}{\partial x}\right)}_{x_e}\text{and}{P}_2^e=\frac{1}{Pe}{\left(\frac{\partial {\varphi }_h^e}{\partial x}\right)}_{x_{e+1}}$ (115)

We can write (114) as:

${\left(A{\varphi }_h^e-f,v\right)}_{{\overline{\Omega }}_x^e}=\left[C^e\right]\left\{{\dot{\delta }}^e\right\}+\left[K^e\right]\left\{{\delta }^e\right\}-\left\{F^e\right\}-\left\{P^e\right\}$ (116)

in which:

$\begin{array}{l}{C}_{ij}^e={\left(N_j,N_i\right)}_{{\overline{\Omega }}_x^e},\\ \left[K_{ij}^e\right]={\left(\frac{\partial N_j}{\partial x},N_i\left(x\right)\right)}_{{\overline{\Omega }}_x^e}+\frac{1}{Pe}{\left(\frac{\partial N_j}{\partial x},\frac{\partial N_i}{\partial x}\right)}_{{\overline{\Omega }}_x^e};i,j=1,2,\cdots ,n\\ \left\{F_i^e\right\}={\left(f\left(x,t\right),N_i\left(x\right)\right)}_{{\overline{\Omega }}_x^e}\\ \left\{p^e\right\}=\left[\begin{array}{l}{P}_1^e\hfill \\ 0\hfill \\ 0\hfill \\ P_2^e\hfill \end{array}\right]\end{array}$(117)

Substituting (116) into (110), we obtain:

${\left(A{\varphi }_h^e-f,v\right)}_{{\overline{\Omega }}_x^T}=\left[C\right]\left\{\dot{\delta }\right\}+\left[K\right]\left\{\delta \right\}-\left\{F\right\}-\left\{P\right\}=0$ (118)

in which:

$\left[C\right]={\displaystyle \sum }_e\left[C^e\right],\left[K\right]={\displaystyle \sum }_e\left[K^e\right],\left\{F\right\}={\displaystyle \sum }_e\left\{F^e\right\},\left\{p\right\}={\displaystyle \sum }_e\left\{p^e\right\}$ (119)

Equation (118) are a system of first-order ordinary differential equations that are integrated using Newmark’s linear method.

$L_2$ -norm of $E$ i.e. ${\Vert E\Vert }_{L_2}$ is given by first substituting ${\varphi }_h^e$ from (107) in initial value problem and obtaining $E^e$ $\forall$ $x,t\in {\overline{\Omega }}_{xt}^e$ and then calculating $L_2$ -norm of $E^e$ followed by sum over the element in the spatial discretization

$E^e={\displaystyle \sum }_{i=1}^n{N}_i{\dot{\delta }}_i^e-{\displaystyle \sum }_{i=1}^n\frac{\partial N_i}{\partial x}{\delta }_i^e-\frac{1}{Pe}{\displaystyle \sum }_{i=1}^n\frac{{\partial }^2{N}_i}{\partial x^2}{\delta }_i^e$(120)

and

${\Vert E\Vert }_{L_2}=\sqrt{I}=\sqrt{{\displaystyle \sum }_e{\left(E^e,E^e\right)}_{{\overline{\Omega }}_{xt}^e}}$ (121)

${\overline{\Omega }}_{xt}^e$ is mapped into a two-unit square in a natural coordinate space $\xi ,\nu$ space.

4.1.3. Computation of Solutions and a Posteriori Errors

In both space-time coupled and space-time decoupled methods, we consider the following $g\left(x\right)$ i.e. $\varphi \left(x,0\right)$ defining initial conditions (Figure 4):

Figure 4. Initial condition.

$\{\begin{array}{l}\text{At}x=0.1,t=0.0,\varphi \left(0.1,0\right)=0,\frac{\partial \varphi }{\partial x}=0,\\ \text{At}x=0.2,t=0.0,\frac{\partial \varphi }{\partial x}=0,\varphi =1\\ \text{At}x=0.3,t=0.0,\varphi \left(0.3,0\right)=0,\frac{\partial \varphi }{\partial x}=0,\\ \varphi \left(x,0\right)=0,\forall x\le 0.1\text{and}t=0.0\\ \varphi \left(x,0\right)=0,\forall x\ge 0.3\text{and}t=0.0\end{array}$ (122)

We consider a ten-element uniform discretization in space-time coupled as well as space-time decoupled methods. The values of $\text{Δ}t$ used for various studies are either shown on the graphs or are described with the results. We choose $Pe=30$ . For such low $Pe$ , the physical diffusion is significant, thus substantial base elongation and amplitude decay of the applied pulse is expected during evolution.

We consider solutions of class $C^1$ and $C^2$ in space and time for space-time coupled finite element method, and local approximations in class $C^1$ and $C^2$ in space for space-time decoupled finite element method with time integration of ordinary differential equations using Newmark’s linear method. In both, space-time coupled finite element method and space-time decoupled finite element method, solutions are calculated for $\text{Δ}t=0.05$ with p-levels 5, 7, 9 and 11 in space and time i.e. in space-time decoupled methods we also choose $\text{Δ}t=0.05$ with p-levels 5, 7, 9 and 11 for spatial finite element discretization. Solutions are computed for first time step and then time-marched in both space-time coupled and space-time decoupled methods, keeping $\text{Δ}t=0.05$ constant.

Figure 5 shows graphs of ${\Vert E\Vert }_{L_2}$ or $\sqrt{I}$ vs degrees of freedom using log-log scale for solutions of class $C^1$ and $C^2$ for the first space-time strip $\left(0\le t\le 0.05\right)$ obtained by using the solution from the space-time coupled finite element method. The data points on the graphs correspond to p-levels of 5, 7, 9, 11, both in space and time. Progressively increasing p-levels results in progressively reduced ${\Vert E\Vert }_{L_2}$ , hence progressively improved accuracy for both solutions of class $C^1$ and $C^2$ . Solutions of $C^2$ yield almost an order of magnitude reduction in ${\Vert E\Vert }_{L_2}$ corresponding to each p-level. Slight improvement in convergence rate of ${\Vert E\Vert }_{L_2}$ is observed for the solutions of class $C^2$ . This behavior of ${\Vert E\Vert }_{L_2}$ is what is expected according to the theory of space-time coupled finite element processes based on space-time residual functional. Our objective is to see if the solution of similar accuracy are possible in space-time decoupled methods for same discretization, p-levels and $\text{Δ}t$ .

Figure 5. $L_2$ -norm of $E$ versus degrees of freedom for the solution of class $C^1$ and $C^2$ ; space-time coupled method.

Figure 6 shows plots of ${\Vert E\Vert }_{L_2}$ versus degrees of freedom using log-log scale for space-time decoupled method for the first time increment obtained using the solutions of class $C^1$ and $C^2$ in space, for p-levels of $5,7,9,11$ and $\text{Δ}t=0.05$ in Newmark’s linear time integration method. Lower values of ${\Vert E\Vert }_{L_2}$ are observed in class $C^2$ in space compared to class $C^1$ . For the same computational resources (discretization, p-levels, solution class, and $\text{Δ}t$ ), the ${\Vert E\Vert }_{L_2}$ values are of orders of magnitude higher in space-time decoupled methods compared to space-time coupled method (Figure 5 and Figure 6). Another disturbing aspect of the graphs in Figure 6 for space-time decoupled methods is that with increasing p-levels, ${\Vert E\Vert }_{L_2}$ increases indicating progressively reducing aaccuracy of the solution. This holds true for both solutions of classes $C^1$ and $C^2$ in space. We provide two explanations for this behavior in the following.

Figure 6. $L_2$ -norm of $E$ versus degrees of freedom for solution of class $C^1$ and $C^2$ ; space-time decoupled method.

The first possible explanation is due to time approximation operator [B] in the time integration process. In space-time decoupled method we are time marching the solution using Equation (38):

$\left[B\right]{\left\{\delta \right\}}_{t+\text{Δ}t}={\left\{Q\right\}}_t+{\left\{F\right\}}_{t+\text{Δ}t}$ (123)

The matrix $\left[B\right]$ has eigen-spectrum $\left({\omega }_i,{\left\{\varphi \right\}}_i\right)$ , where $i=1,2,\cdots ,n$ for some p-level say $p$ for a fixed order of space. For a fixed discretization, a change in the p-level say $p_1>p$ results in changes in the eigen-spectrum of $\left[B\right]$ . Hence, we now have $\left({\omega }_i,{\left\{\varphi \right\}}_i\right)$ , where $i=1,2,\cdots ,n$ instead of $\left({\omega }_i,{\varphi }_i\right),i=1,2,\cdots ,n$ . In case of $\left({\omega }_i,{\left\{\varphi \right\}}_i\right)$ for p-level of $p$ , $\text{Δ}t=0.05$ produces results with some fixed accuracy i.e. a value of ${\Vert E\Vert }_{L_2}$ . With $p_1>p$ , increased ${\Vert E\Vert }_{L_2}$ , hence increased error in the solution implies that for this p-level of $p_1$ , $\text{Δ}t=0.05$ is not adequate to integrate Equations (123). This is possible only if $\left[B\right]$ at p-level of $p_1$ is a stiffer system with higher frequency contents than at p-level of $p$ contributing to the solution. If this is indeed true than some $\text{Δ}t$ lower than 0.05 can be used to integrate the higher frequency contents accuracy, thus achieve accuracy similar to p-level of p. We present numerical study to demonstrate this.

Figure 7 shows graphs of ${\Vert E\Vert }_{L_2}$ vs degrees of freedom using log-log scale for solutions of class $C^1$ in space and time using space-time coupled methods for p-levels of 5, 7, 9, 11 with $\text{Δ}t=0.05$ . This figure also shows graphs ${\Vert E\Vert }_{L_2}$ vs degrees of freedom for space-time decoupled methods with local approximation of $C^1$ in space, p-levels of 5, 7, 9, 11, and $\text{Δ}t=0.05,0.01,0.005,0.0005$ and 0.00005. We note that for fixed $\text{Δ}t$ (regardless of it’s value) increase in ${\Vert E\Vert }_{L_2}$ with increasing p-level persists. But a decrease in $\text{Δ}t$ with increase p-level has dramatic effect in reducing ${\Vert E\Vert }_{L_2}$ as shown in Figure 7 (dashed line graphs). Even with $\text{Δ}t=0.0005$ , i.e., 100 time steps in space-time decoupled methods in the time interval $\left[0,0.05\right]$ compared to only one time step for space-time coupled method, the space-time decoupled method struggles in approaching the results obtained for ${\Vert E\Vert }_{L_2}$ versus degrees of freedom for space-time coupled finite element method. A further reduction in $\text{Δ}t$ to 0.00005 (1000 time steps in $t=0$ to 0.005), ${\Vert E\Vert }_{L_2}$ versus dofs improves and is much closer to ${\Vert E\Vert }_{L_2}$ versus dofs for space-time coupled methods, but not better.

Figure 7. $L_2$ -norm of $E$ versus degrees of freedom for solution of class $C^1$ ; space-time coupled and decoupled methods.

The second explanation for increasing ${\Vert E\Vert }_{L_2}$ with increasing p-level is as follows. We note that with progressively reducing $\text{Δ}t$ , increase in ${\Vert E\Vert }_{L_2}$ with increasing p-level reduces in magnitude. At $\text{Δ}t=0.00005$ , ${\Vert E\Vert }_{L_2}$ versus degrees of freedom for the first time shows progressively reducing ${\Vert E\Vert }_{L_2}$ with increasing p-levels, which is what is expected. Another observation is that all ${\Vert E\Vert }_{L_2}$ values at all p-levels for $\text{Δ}t>0.005$ are greater than $\mathcal{O}\left({10}^{-1}\right)$ . Only for $\text{Δ}t=0.0005$ , but more so far $\text{Δ}t=0.00005$ , ${\Vert E\Vert }_{L_2}$ values are within the 3 × 10⁻¹ to 5 × 10⁻² range for which progressively reducing ${\Vert E\Vert }_{L_2}$ with progressively increasing p-levels is observed, which is the expected behavior, i.e., increasing p-levels must result in better accuracy. Comparisons of ${\Vert E\Vert }_{L_2}$ versus dofs for $\text{Δ}t>0.005$ and $\text{Δ}t=0.00005$ suggest that perhaps ${\Vert E\Vert }_{L_2}$ of the order $\mathcal{O}\left({10}^{-1}\right)$ are not low enough for the computed solutions to be reliable, thus ${\Vert E\Vert }_{L_2}$ computations are not with sufficient accuracy, hence the reason for increasing ${\Vert E\Vert }_{L_2}$ with increasing p-level for $\text{Δ}t>0.005$ . It is further disturbing to note that what is achievable in space-time coupled methods in one increment of time with $\text{Δ}t=0.05$ , cannot be achieved even with 1000 time steps between $t=0$ to $t=0.05$ . This study undoubtedly establishes that ${\Vert E\Vert }_{L_2}$ values of the order $\mathcal{O}\left({10}^{-1}\right)$ correspond to solutions that are not reliable, hence cannot be used to judge solution error or useful in adaptive processes. Much lower values ($\mathcal{O}\left({10}^{-3}\right)$ or lower) of ${\Vert E\Vert }_{L_2}$ are needed to be reliable measures.

Figure 8 shows graphs of ${\Vert E\Vert }_{L_2}$ versus degrees of freedom for the solution of class $C^2$ in space and time obtained from the space-time coupled method for p-levels 5, 7, 9, and 11, with $\text{Δ}t=0.05$ . This figure also shows graphs of ${\Vert E\Vert }_{L_2}$ versus degrees of freedom for space-time decoupled method with p-levels 5, 7, 9, and 11 and $\text{Δ}t=0.05,0.01,0.005,0.0005$ , and 0.00005. Here also, we note that for fixed $\text{Δ}t$ (regardless of its value), increases in ${\Vert E\Vert }_{L_2}$ with increasing p-levels persists indicating inaccurale ${\Vert E\Vert }_{L_2}$ values. But, a decrease in $\text{Δ}t$ with in increase in p-level has a dramatic effect in reducing ${\Vert E\Vert }_{L_2}$ , as shown in Figure 8 (dashed line graphs). Even with $\text{Δ}t=0.0005$ i.e. 100 time steps in the space-time decoupled method in the time interval $\left[0,0.05\right]$ compared to one time step in space-time coupled finite element method, the space-time decoupled method struggles in approaching the results obtained for ${\Vert E\Vert }_{L_2}$ versus degrees of freedom for space-time coupled finite element methods. Further reduction in $\text{Δ}t$ to 0.00005 improves ${\Vert E\Vert }_{L_2}$ (more improvement at higher p-levels) but is unable to match ${\Vert E\Vert }_{L_2}$ versus dofs for space-time coupled method. Here also, we observe that for $\text{Δ}t>0.005$ , ${\Vert E\Vert }_{L_2}$ versus dofs shows increasing values for increasing p-levels and for this range of $\text{Δ}t$ , ${\Vert E\Vert }_{L_2}$ values are of the order of $\mathcal{O}\left({10}^{-1}\right)$ , same as for solutions of class $C^1$ (shown in Figure 7) indicating that these do not have enough accuracy. However, for $\text{Δ}t=0.0005$ and 0.00005, we clearly observe decreasing ${\Vert E\Vert }_{L_2}$ with increasing p-levels, expected behavior. ${\Vert E\Vert }_{L_2}$ are $\mathcal{O}\left({10}^{-2}\right)$ , close enough to be somewhat reliable. This study for solution of class $C^2$ in space also confirms that ${\Vert E\Vert }_{L_2}$ must be $\mathcal{O}\left({10}^{-2}\right)$ or lower for it to be meaningful in judging accuracy of the solution and for it to be of any benefit in adaptivity. Results presented in Figure 7 and Figure 8 demonstrate that space-time decoupled finite element processes cannot be competitive with space-time coupled finite element processes.

Figure 8. $L_2$ -norm of $E$ versus degrees of freedom for solution of class $C^2$ ; space-time coupled and decoupled methods.

Figure 9 shows plots of solution $\varphi$ versus $x$ for $t=0.05$ and $t=0.25$ , for solutions of class $C^1$ in space and time, $\text{Δ}t=0.05$ and p-level of 11 in space and time for space-time coupled finite element method. Solutions obtained for space-time decoupled method with local approximation of class $C^1$ in space with p-level 11 for varying values of $\text{Δ}t$ ($\text{Δ}t=0.05,0.01,0.005,0.0005$ ), are also shown in Figure 9 for $t=\text{Δ}t$ and $t=5\text{Δ}t$ . We note from the solution plotted in Figure 9 at $t=\text{Δ}t$ that a good agreement of the space-time decoupled solution with space-time coupled finite element solution is only observed for $\text{Δ}t=0.0005$ .

Figure 10 shows graphs of the solution $\varphi$ versus the distance $x$ for $t=0.05$ and $t=0.25$ for solutions of class $C^2$ in space and time, $\text{Δ}t=0.05$ , and p-level of 11 in space and time for space-time coupled finite element method. Solutions obtained from space-time decoupled method with local approximation of class $C^2$ in space and p-level of 11 for $\text{Δ}t=0.05,0.01,0.005,0.0005$ are also shown in Figure 10 for $t=\text{Δ}t$ and at $t=5\text{Δ}t$ . We note from the solutions shown in Figure 10 at $t=\text{Δ}t$ that a good agreement with space-time coupled method solutions is only observed for $\text{Δ}t=0.0005$ . These plots of the solution confirm that the ${\Vert E\Vert }_{L_2}$ -norm is meaningful for $\text{Δ}t=0.0005$ . In some instances, it is entirely possible that ${\Vert E\Vert }_{L_2}$ has acceptable value but the solution may have some spurious behavior in partial or the entire domain. We illustrate this point in the next model problem.

Figure 9. Solution $\varphi$ versus distance $x$ for the solution of class $C^1$ ; space-time coupled and decoupled methods.

Figure 10. Solution $\varphi$ versus distance $x$ for the solution of class $C^2$ ; space-time coupled and decoupled methods.

4.2. Transient Advection Equation

We consider a 1D transient advection equation given by:

$\frac{\partial \varphi }{\partial t}+\frac{\partial \varphi }{\partial x}=g\left(x,t\right),\forall x,t\in {\Omega }_{xt}={\Omega }_x\times {\Omega }_t=\left(0,1\right)\times \left(0,\tau \right)$ (124)

with the following boundary conditions and initial conditions:

$\{\begin{array}{l}\text{BCs}:\varphi \left(0,t\right)=0,\forall t\in \left[0,\tau \right]\\ \text{ICs}:\varphi \left(x,0\right)=g\left(x\right),\forall x\in \left[0,1\right]\end{array}$ (125)

$g\left(x\right)$ defines initial conditions on $\varphi \left(x,t\right)$ at $t=0$

We can write (124):

$A\varphi -f=0\forall x,t\in {\Omega }_{xt}$ (126)

in which:

$A=\frac{\partial }{\partial t}+\frac{\partial }{\partial x}$ (127)

4.3. Space-Time Coupled Finite Element Method Based on Space-Time Residual Function

Let ${\overline{\Omega }}_{xt}={\displaystyle \underset{e}{\cup }{\overline{\Omega }}_{xt}^e}$ be the discretization of ${\overline{\Omega }}_{xt}={\Omega }_{xt}\cup \Gamma =\left[0,1\right]\times \left(0,\tau \right)$ , $\text{Γ}$ is the boundary of ${\Omega }_{xt}$ and ${\overline{\Omega }}_{xt}^e$ is a space-time element of the discretization ${\overline{\Omega }}_{xt}^T$ of the first space-time strip. Let ${\varphi }_h^e\left(x,t\right)$ be local approximation of $\varphi \left(x,t\right)$ over ${\overline{\Omega }}_{xt}^e$ and let ${\varphi }_h\left(x,t\right)$ be global approximation of $\varphi \left(x,t\right)$ over ${\overline{\Omega }}_{xt}^T$ such that:

$\{\begin{array}{l}{\varphi }_h\left(x,t\right)={\displaystyle \underset{e}{\cup }{\varphi }_h^e\left(x,t\right)}\\ {\varphi }_h^e\left(x,t\right)={\displaystyle \sum }_{i=1}^n{N}_i\left(x,t\right){\delta }_i^e,\left\{\delta \right\}={\displaystyle \underset{e}{\cup }\left\{{\delta }_i^e\right\}}\end{array}$ (128)

Consider discretization ${\overline{\Omega }}_{xt}^T={\displaystyle \underset{e}{\cup }{\overline{\Omega }}_{xt}^e}$ exactly similar to Section 4.1.1., then space-time residual functional $E^e$ over ${\overline{\Omega }}_{xt}^e$ is given by:

$\frac{\partial {\varphi }_h^e}{\partial t}+\frac{\partial {\varphi }_h^e}{\partial x}-f=E^e\left(x,t\right),\forall x,t\in {\overline{\Omega }}_{xt}^e.$(129)

Following the procedure described in Section 4.1.1 for transient convection diffusion we can obtain:

$\left[K\right]\left\{\delta \right\}=\left\{F\right\}$ (130)

in which

$\left[K\right]={\displaystyle \sum }_e\left[K^e\right],\left\{F\right\}={\displaystyle \sum }_e\left\{F^e\right\},\left\{\delta \right\}={\displaystyle \underset{e}{\cup }\left\{{\delta }^e\right\}}$ (131)

$K_{ij}^e={\left(\frac{\partial N_j}{\partial t}+\frac{\partial N_j}{\partial x},\frac{\partial N_i}{\partial t}+\frac{\partial N_i}{\partial x}\right)}_{{\overline{\Omega }}_{xt}^e};i,j=1,2,\cdots ,n$(132)

$F_i^e={\left(f,\frac{\partial N_i}{\partial t}+\frac{\partial N_j}{\partial x}\right)}_{{\overline{\Omega }}_{xt}^e}$ (133)

${\Vert E\Vert }_{L_2}=\sqrt{{\displaystyle \sum }_e{\left(E^e,E^e\right)}_{{\overline{\Omega }}_{xt}^e}}$ (134)

Equations (130) are a system of linear algebraic equations for the first (or any) space-time strip. We calculate $\left\{\delta \right\}$ using (122) after imposing boundary conditions and initial conditions. Since the initial value problem, in this case is first order space and also first order in time, local approximation of class $C^1$ in space and time will maintain space-time integrals over ${\overline{\Omega }}_{xt}^T$ in Riemann sense.

4.3.1. Space-Time Decoupled Finite Element Method

Let ${\overline{\Omega }}_x^T={\displaystyle \underset{e}{\cup }{\overline{\Omega }}_x^e}$ be the discretization of spatial domain ${\overline{\Omega }}_x$ in which ${\overline{\Omega }}_x^e$ is three node p-version hierarchical 1D element in space with higher-order global differentiability.

Let

${\varphi }_h^e\left(x,t\right)={\displaystyle \sum }_{i=1}^n{N}_i\left(x\right){\delta }_i^e\left(t\right)\forall x\in {\Omega }_x^e,\forall t\in \left(0,t\right)$ (135)

and let ${\varphi }_h\left(x,t\right)={\displaystyle \underset{e}{\cup }{\varphi }_h^e\left(x,t\right)}$ be the approximation of $\varphi \left(x,t\right)$ over the space-time domain ${\Omega }_{xt}^T={\Omega }_x^T\times \left(0,t\right)$ and $t=\text{Δ}t$ .

Consider the integral form of (125) based on the fundamental lemma over ${\overline{\Omega }}_x^T$ using Galerkin method with weak form:

${\left(A{\varphi }_h-f,v\right)}_{{\overline{\Omega }}_x^T}={\displaystyle \sum }_e{\left(A{\varphi }_h^e-f,v\left(x\right)\right)}_{{\overline{\Omega }}_x^e}$ (136)

Consider ${\left(A{\varphi }_h^e-f,v\left(x\right)\right)}_{{\overline{\Omega }}_x^e}$ , in which $v\left(x\right)=\delta {\varphi }_h^e=N_j\left(x\right);j=1,2,\cdots ,n$ :

$\left(A{\varphi }_h^e-f,v\right)={\left(\frac{\partial {\varphi }_h^e}{\partial t}+\frac{\partial {\varphi }_h^e}{\partial x}-f,v\left(x\right)\right)}_{{\overline{\Omega }}_x^e}$ (137)

Using (135):

$\begin{array}{l}\frac{\partial {\varphi }_h^e}{\partial t}={\displaystyle \sum }_{i=1}^n{N}_i\left(x\right){\dot{\delta }}_i^e\left(t\right)\\ \frac{\partial {\varphi }_h^e}{\partial x}={\displaystyle \sum }_{i=1}^n\frac{\partial N_i\left(x\right)}{\partial x}{\delta }_i^e\left(t\right)\\ \text{and}v\left(x\right)=N_j\left(x\right);j=1,2,\cdots ,n\end{array}$(138)

Substituting (138) into (137), we can write:

${\left(A{\varphi }_h^e-f,v\right)}_{{\overline{\Omega }}_x^e}=\left[C^e\right]\left\{{\dot{\delta }}^e\right\}+\left[K^e\right]\left\{{\delta }^e\right\}-\left\{F^e\right\}=0$ (19)

And using (139) into (136), we obtain:

${\left(A{\varphi }_h-f,v\right)}_{{\overline{\Omega }}_x^T}={\displaystyle \sum }_e\left[C^e\right]\left\{{\dot{\delta }}^e\right\}+{\displaystyle \sum }_e\left[K^e\right]\left\{{\delta }^e\right\}-{\displaystyle \sum }_e\left\{F^e\right\}=0$ (140)

$=\left[C\right]\left\{\dot{\delta }\right\}+\left[K\right]\left\{\delta \right\}-\left\{F\right\}=0$ (141)

in which:

$\begin{array}{l}\left[C\right]={\displaystyle \sum }_e\left[C^e\right],\left[K\right]={\displaystyle \sum }_e\left[K^e\right],\left\{F\right\}={\displaystyle \sum }_e\left\{F^e\right\}\\ \left\{\delta \right\}={\displaystyle \underset{e}{\cup }\left\{{\delta }^e\right\}}\\ C_{ij}^e={\left(N_j,N_i\right)}_{{\overline{\Omega }}_x^e}\\ K_{ij}^e={\left(\frac{\partial N_j\left(x\right)}{\partial x},\frac{\partial N_i\left(x\right)}{\partial x}\right)}_{{\overline{\Omega }}_x^e}\\ F_i^e={\left(f,N_i\left(x\right)\right)}_{{\overline{\Omega }}_x^e};\text{for}i,j=1,2,\cdots ,n.\end{array}$(142)

$L_2$ -norm of $E$ i.e., ${\Vert E\Vert }_{L_2}$ , is given by first substituting (142) in (131) to obtain $E^e$ , the residual function for ${\overline{\Omega }}_x^e$ :

$E^e={\displaystyle \sum }_{i=1}^n{N}_i\left(x\right){\dot{\delta }}^e\left(t\right)+{\displaystyle \sum }_e\frac{\partial N_i\left(x\right)}{\partial x}{\delta }_i^e\left(t\right)$(143)

Now,

${\Vert E\Vert }_{L_2}=\sqrt{{\displaystyle \sum }_e{\left(E^e,E^e\right)}_{{\overline{\Omega }}_{xt}}}.$ (144)

${\overline{\Omega }}_{xt}^e$ is a two-square unit in the natural coordinate space $\xi ,\nu$ .

4.3.2. Computation of Solutions and a Posteriori Errors

We consider the initial conditions shown in Figure 4, same mesh used as in Section 4.1.3 for transient convection-diffusion equations. In this initial value problem, there is no physical diffusion, hence we expect the initial pulse to remain the same but advance in the x-direction as time elapses. Amplitude decay and base elongation in this model problem can only happen if the numerical computational process has numerical dispersion (artificial inherent diffusion).

We consider solutions of class $C^1$ and $C^2$ in space and time for space-time coupled finite element methods and local approximations of class $C^1$ and $C^2$ in space for space-time decoupled finite element methods with time integration of resulting ordinary differential equations in time using Newmark’s linear method. In both space-time coupled and space-time decoupled methods, solutions are calculated for $\text{Δ}t=0.05$ with p-levels of 5, 7, 9, and 11 in space and time i.e. in space-time decoupled method with p-levels of 5, 7, 9, and 11 for spatial finite element discretization. Solutions are computed for the first time step and then time-marched in both space-time coupled and space-time decoupled methods, keeping $\text{Δ}t=0.05$ constant.

Figure 11 shows graphs of ${\Vert E\Vert }_{L_2}$ or $\sqrt{I}$ versus degrees of freedom using log-log scale for solutions of class $C^1$ and $C^2$ for the first space-time strip $\left(\text{Δ}t=0.05\right)$ obtained by using the solutions from the space-time coupled finite element method. The data points on the graphs correspond to p-levels of 5, 7, 9, and 11, both in space and time. Progressively increasing p-levels results in progressively reduced ${\Vert E\Vert }_{L_2}$ , hence progressively improved accuracy of both solutions of class $C^1$ and $C^2$ . Solutions of class $C^2$ result in reduction in ${\Vert E\Vert }_{L_2}$ at p-level 9 and 11. Slight improvement in the convergence rate of ${\Vert E\Vert }_{L_2}$ is also observed for solutions of class $C^2$ beyond p-level 7. This behavior of ${\Vert E\Vert }_{L_2}$ versus degrees of freedom is generally expected.

Figure 11. $L_2$ -norm of $E$ versus degrees of freedom for solution of class $C^1$ and $C^2$ ; space-time coupled method.

Figure 12 show plots ${\Vert E\Vert }_{L_2}$ versus degrees of freedom using log-log scale for the first time increment obtained from the solutions of class $C^1$ and $C^2$ in space for p-levels of 5, 7, 9, and 11, at $\text{Δ}t=0.05$ in Newmark’s linear time integration method. Slightly higher values of ${\Vert E\Vert }_{L_2}$ are observed for solutions of $C^2$ in space compared to solutions of $C^1$ in space. This is an indication that for the model problem, ${\Vert E\Vert }_{L_2}$ values of $\mathcal{O}\left({10}^0\right)$ are not sufficiently accurate to draw conlcusions for class $C^1$ and $C^2$ solutions. Note in Figure 11 that ${\Vert E\Vert }_{L_2}$ for space-time coupled methods is between $\mathcal{O}\left({10}^{-2}\right)-\mathcal{O}\left({10}^{-3}\right)$ , the range we need to achieve for space-time decoupled method to have meaningful values of ${\Vert E\Vert }_{L_2}$ and accurate solution. Thus, any conclusions based on graphs in Figure 12 are not meaningful as the ${\Vert E\Vert }_{L_2}$ values are not reliable.

Figure 13 shows plots of ${\Vert E\Vert }_{L_2}$ versus degrees of freedom for solutions of class $C^1$ in space and time for the space-time coupled method with p-levels of 5, 7, 9, and 11 using $\text{Δ}t=0.05$ . This figure also shows graphs of ${\Vert E\Vert }_{L_2}$ versus degrees of freedom for the space-time decoupled method with local approximation of class $C^1$ in space, p-levels of 5, 7, 9, and 11, and $\text{Δ}t=0.05,0.01,0.005,0.0005$ and 0.00005. We note that for fixed $\text{Δ}t$ (regardless of its value), increase in ${\Vert E\Vert }_{L_2}$ with increasing p-level persists. But, a decrease in $\text{Δ}t$ with an increase in $p$ -level has dramatic effects on reducing ${\Vert E\Vert }_{L_2}$ , as shown in Figure 13 (graphs with dashed lines). Even with $\text{Δ}t=0.0005$ i.e., 100 time steps and $\text{Δ}t=0.00005$ i.e. 1000 time steps in the space-time decoupled methods in the time interval $\left[0,0.05\right]$ , the increase in ${\Vert E\Vert }_{L_2}$ with increasing p-levels persists. ${\Vert E\Vert }_{L_2}$ values for all time steps and p-levels of the order of $\mathcal{O}\left({10}^{-1}\right)$ confirm the lack of accuracy of the computed solution, resulting in higher values of ${\Vert E\Vert }_{L_2}$ in the range that is not reliable.

Figure 12. $L_2$ -norm of $E$ versus degrees of freedom for solution of class $C^1$ and $C^2$ ; space-time decoupled methods.

Figure 13. $L_2$ -norm of $E$ versus degrees of freedom for solution of class $C^1$ ; space-time coupled and decoupled method.

Figure 14 shows plots of ${\Vert E\Vert }_{L_2}$ versus degrees of freedom for solutions of class $C^2$ in space and time obtained from the space-time coupled method for p-levels of 5, 7, 9, and 11 with $\text{Δ}t=0.05$ . This figure also shows graphs of ${\Vert E\Vert }_{L_2}$ versus degrees of freedom for the space-time decoupled method for p-levels of 5, 7, 9, and 11 with $\text{Δ}t=0.05,0.01,0.005,0.0005$ and 0.00005. First thing we observe that for $\text{Δ}t=0.05$ and 0.01, there is virtually no change in ${\Vert E\Vert }_{L_2}$ with increasing p-level, i.e., we neither observe an increase in ${\Vert E\Vert }_{L_2}$ with increasing p-levels nor decreases for these values of $\text{Δ}t$ . For $\text{Δ}t\le 0.005$ , ${\Vert E\Vert }_{L_2}$ progressively reduces with increasing p-level, expected behavior, and further reduction in $\text{Δ}t$ results in even reduced ${\Vert E\Vert }_{L_2}$ for the same p-level, also expected correct behavior of ${\Vert E\Vert }_{L_2}$ versus dofs. ${\Vert E\Vert }_{L_2}$ values in the range $\mathcal{O}\left({10}^{-1}\right)-\mathcal{O}\left({10}^{-2}\right)$ for $\text{Δ}t\le 0.005$ confirm that we should expect solutions with reasonable accuracy. We confirm this in Figure 15.

Figure 14. $L_2$ -norm of $E$ versus degrees of freedom for solution of class $C^2$ ; space-time coupled and decoupled methods.

Figure 15 shows plots of solution $\varphi$ of class $C^1$ versus distance $x$ from space-time coupled as well as space-time decoupled methods at $t=0.05$ and $t=0.25$ , i.e., for the first and the fifth time steps. Even for the first time step, the space-time decoupled method solution shows oscillations, some amplitude decay, and minor base elongation. At $t=0.25$ , we observe significant oscillations, amplitude decay, and base elongation in the space-time decoupled solution, while in the solution from the space-time coupled methods, the initial pulse propagates without oscillations, base elongation, and amplitude decay.

Figure 16 shows plots of $\varphi$ versus dofs also for solutions of class $C^1$ for space-time coupled and space-time decoupled methods at $t=0.05$ and $t=0.25$ , computed using $\text{Δ}t=0.05$ for space-time coupled methods and $\text{Δ}t=0.0005$ for space-time decoupled methods. Also, Figure 16 shows $\varphi$ versus $x$ at $t=0.05$ and $t=0.25$ obtained for space-time coupled methods ($C^1$ in space and time), using p-level 11 and $\text{Δ}t=0.05$ (same solution as in Figure 15) and a comparison with the solutions obtained from space-time decoupled methods at $t=0.05$ and $t=0.25$ obtained using solution of class $C^1$ in space and integration interval $\text{Δ}t=0.0005$ . Dramatic improvement is obvious, still some minor oscillations of extremely small amplitude are observed but the solution at $t=0.05$ and $t=0.25$ (at the end of 500-time steps) are much improved. Extremely poor performance of space-time decoupled methods is rather obvious compared to space-time coupled methods for identical computational resources i.e. same discretization, same p-level in space and $\text{Δ}t$ (Figure 15).

Figure 15. Solution $\varphi$ versus distance $x$ for solution of class $C^1$ for space-time coupled and decoupled methods.

Figure 16. Solution $\varphi$ versus distance $x$ for solution of class $C^1$ in space-time coupled and decoupled methods.

Figure 17 and Figure 18 show parallel studies to those presented in Figure 15 and Figure 16 but using solutions of class $C^2$ in space and time for space-time coupled methods and solution of class $C^2$ in space for space-time decoupled methods. In Figure 17, where $\text{Δ}t=0.05$ for space-time coupled and space-time decoupled methods, we do not observe any visible improvement in the solutions from space-time decoupled method compared to solutions in Figure 15 for class $C^1$ space. However, in Figure 18, with solutions of class $C^2$ , $\text{Δ}t=0.0005$ , the space-time decoupled solutions are virtually oscillation and much improved (compared to those in Figure 15 for class $C^1$ with $\text{Δ}t=0.0005$ ) illustrating the influence of higher order space.

Figure 17. Solution $\varphi$ versus distance $x$ for solution of class $C^2$ in space-time coupled and decoupled methods.

Figure 18. Solution $\varphi$ versus distance $x$ for solution of class $C^2$ for space-time coupled and decoupled methods.

Figure 19 shows plots of $\varphi$ versus $x$ for the first six-time steps with solutions of class $C^1$ in space, p = 11 and $\text{Δ}t=0.05$ obtained using space-time decoupled method. We have already seen that this value of $\text{Δ}t$ is too large for accurate evoution of $\varphi$ , but here we demonstrate a different aspect of the study. Figure 19 shows progressively more oscillations, progressively decaying amplitudes, and progressively increasing base with increasing time. Intuitively, we expect progressively increasing ${\Vert E\Vert }_{L_2}$ for each subsequent time step as the solution is deteriorating.

Figure 19. Solution $\varphi$ versus distance $x$ : Solution of class $C^1$ , $p=11$ , $\text{Δ}t=0.05$ ; space-time decoupled method: six time-steps.

Figure 20 shows plots of ${\Vert E\Vert }_{L_2}$ versus dofs for solutions of class $C^1$ with $\text{Δ}t=0.05$ for the first, sixth, and tenth time steps. We observe no significant deterioration in ${\Vert E\Vert }_{L_2}$ versus dofs as time elapses. Thus, based on the graphs in Figure 20, solution for all ten time steps are of same quality and accuracy. This is obviously not true, as shown in Figure 19.

Figure 20. $L_2$ -norm of $E$ versus DOFs for solution of class $C^1$ ; space-time decoupled methods.

Figure 19 clearly shows deteriorating solutions as time elapses but Figure 20 suggests otherwise. The important conclusion from this study is that may be ${\Vert E\Vert }_{L_2}$ alone is not an absolute measure of the accuracy of the solution when its computed value is outside of the threshold of computed error. High magnitude ${\Vert E\Vert }_{L_2}$ in Figure 20 for all p-levels indicates that the solution is not good. What is not clear is that, based on Figure 20, we could conclude that all solutions for each of the six-time steps are of equal accuracy, which is not true. A sure way to avoid this trap is to make sure that ${\Vert E\Vert }_{L_2}$ is sufficienty close to zero (in this study $<\mathcal{O} \left({10}^{-2}\right)$ is sufficient), then these issues will not arise. Unfortunately, this level of accuracy is quite difficult to achieve in space-time decoupled methods. Surana et al. have shown that ${\Vert E\Vert }_{L_2}<\mathcal{O}\left({10}^{-6}\right)$ generally indicates solutions of good accuracy in dimensionless mathematical models in space-time coupled method based on residual functional. Since this level of accuracy is difficult to achieve in space-time decoupled methods, what magnitude of ${\Vert E\Vert }_{L_2}$ is a measure of good accuracy of solutions requires convergence studies possibly for each initial value problem of interest.

5. Summary and Conclusions

An important aspect of the a posteriori computations presented in this paper is that they are based on space-time residual functional obtained by substituting the element local approximations in the space-time differential equation(s) describing the initial value problem for both space-time coupled and space-time decoupled finite element processes. The necessitates that the local approximations in both space-time coupled and space-time decoupled methods must be in higher order minimally conforming spaces. Choice of the orders of approximation space in space and time for space-time coupled methods must ensure that the integrals over the space-time discretization are Riemann. If the theoretical solution is smooth, the integrals in Lebesgue sense over the discretization may suffice. In case of space-time decoupled finite element method, the order of the approximation space in space must be the same as in coupled methods. Since space-time decoupled methods, ordinary differential equations in time are integrated using time integration methods, the selection of time integration method, must ensure: (i) desired global differentiability of approximation in time over the entire discretized time domain, (ii) that the explicit expression for the approximation of the solution and its derivatives as a function of time is obtainable. In the absence of either these two aspects, the local approximations are not admissible in the mathematical model. We make some remarks in the following and draw some conclusions.

1. We note that accurate and meaningful a posteriori computations for initial value problems proposed here are not possible in space-time coupled or space-time decoupled processes without the use of “$k$ ”, higher-order approximation spaces in space and time i.e. $k$ -version of the finite element method.

2. In space-time coupled methods, the local approximation functions are functions of space and time and dofs are constants to be determined. Thus, ${\Vert E\Vert }_{L_2}$ for space-time element, hence for the entire space-time strip is straightforward to complete with minimally conforming or appropriate approximation spaces.

3. In space-time decoupled methods, the approximation functions are function of space coordinates and dofs are function of time, but both are in minimally conforming spaces. Hence, using this approximation, ${\Vert E\Vert }_{L_2}$ for an element in space and for an increment of time can be calculated in a similar manner as in space-time coupled methods, followed by ${\Vert E\Vert }_{L_2}$ for the spatial discretization for an increment of time.

4. When ${\varphi }_h\left(x,t\right)$ , the finite element solution for the spce-time discretization in either methods, approaches the theoretical solution $\varphi$ , then ${\Vert E\Vert }_{L_2}=\sqrt{I}\to 0$ , thus proximity of ${\Vert E\Vert }_{L_2}$ to zero is a measure of accuracy or error in the computed solution in both space-time coupled and space-time decoupled methods.

5. ${\Vert E^e\Vert }_{L_2}$ is space-time finite element or for elements in space for an increment of time can be used as a means for designing adaptive refinements. Elements with the ${\Vert E^e\Vert }_{L_2}$ above a threshold can be candidates for $h,p$ refinements.

6. The a posteriori framework proposed here is general framework and is applicable to space-time coupled as well as space-time decoupled finite element computations regardless of the nature of the space-time differential operators but requires local approximations in $h,p$ and $k$ scalar product spaces, i.e. $k$ -version of finite element method is an essential requirement.

7. In the finite difference and the finite volume methods, the a posteriori computations presented here, based on space-time residual functional, are not possible as in these methods, the computed solution has discrete values at the nodes or at the grid points, hence lacks the description similar to the local approximation that is required for the computation of space-time residual functional.

8. Even in space-time coupled and space-time decoupled finite element processes, if the local approximations are of class $C^0$ in space and time, such computations are not possible if the space-time integrals over the discretization domain are to be in Riemann sense, necessary for assured accuracy.

9. A significant aspect of this work is that it provides a quantitative measure of the error in the computed solution through a posteriori computations without the knowledge of the theoretical solution, which is necessary in real applications to quantitatively determine errors in the computed solution.

10. Model problems studied, presented using transient convection-diffusion equations and transient advection equations, reveal additional insights into the a posteriori computations presented in this paper.

(a) In space-time coupled methods, an increase in p-level for a fixed $k$ results in reduced ${\Vert E\Vert }_{L_2}$ . Also, for fixed p-level, an increase in $k$ results in reduced ${\Vert E\Vert }_{L_2}$ . Increasing p-level for fixed $k$ results in an increasing rate of convergence. However, for fixed $h$ and $p$ , an increase in $k$ results in better accuracy i.e. lower ${\Vert E\Vert }_{L_2}$ but convergence rate remains unaffected. These findings hold regardless of the nature of differential operators in the initial value problem, regardless of complexity of the mathematical model of the initial value problem, and regardless of the application.

(b) In the case of space-time decoupled methods, in both model problem studies, we find that for a larger $\text{Δ}t$ , increasing p-level in space results in increasing ${\Vert E\Vert }_{L_2}$ , which is counter to space-time coupled methods. There are two possible explanations for this behavior: first is due to the eigen spectrum of time approximation operator $\left[B\right]$ . An increase in p-level changes the eigen-spectrum of $\left[B\right]$ , in which high-frequencies begin to contribute to the time response. Integrating these correctly would require lower $\text{Δ}t$ , but if we use same $\text{Δ}t$ as for lower p-level than accuracy of ${\Vert E\Vert }_{L_2}$ will deteriorate. In other words, in space-time decoupled methods an increase in p-level in space must be followed by a decrease in the time integration interval $\text{Δ}t$ to maintain accuracy of ${\Vert E\Vert }_{L_2}$ . The second possible explanation is that, in the range of $\text{Δ}t$ in which ${\Vert E\Vert }_{L_2}$ increases with increasing p-level, ${\Vert E\Vert }_{L_2}$ is of the order $\mathcal{O}\left({10}^{-1}\right)$ , which is not a good measure of computed zero, hence the solutions corresponding to this range of ${\Vert E\Vert }_{L_2}$ are not reliable. When $\text{Δ}t=0.00005$ , ${\Vert E\Vert }_{L_2}$ versus dofs is monotonically decreasing with increasing p-levels, and computed ${\Vert E\Vert }_{L_2}$ is of the order of $\mathcal{O}\left({10}^{-2}\right)$ or lower. This explanation has clearly evidence in the graphs. We keep in mind that ${\Vert E\Vert }_{L_2}$ must be closer to zero for the computed solutions that satisfy partial differential equations, in the mathematical model, hence ${\Vert E\Vert }_{L_2}$ of the order of $\mathcal{O}\left({10}^{-1}\right)$ undoubtedly corresponds to incorrect solutions as partial differential equations are not satisfied accurately. This study demonstrates that to obtain monotonically decreasing ${\Vert E\Vert }_{L_2}$ with increasing p-level, essential for correct convergence behavior of ${\Vert E\Vert }_{L_2}$ , $\text{Δ}t$ of the order of $\mathcal{O}\left({10}^{-4}\right)$ or $\mathcal{O}\left({10}^{-5}\right)$ is necessary, clearly demonstrating serious weaknesses of space-time decoupled methods.

(c) An increase in the order of space for spatial discretization has a beneficial effect on ${\Vert E\Vert }_{L_2}$ in both methods.

(d) Studies presented here show that ${\Vert E\Vert }_{L_2}$ values must be below a certain threshold for them to be a good measures of accuracy. In the model problems, ${\Vert E\Vert }_{L_2}\le \mathcal{O}\left({10}^{-2}\right)$ is found to be adequate when the dimensionless form of the mathematical model is used in computations with this value of ${\Vert E\Vert }_{L_2}$ computed solutions have sufficient accuracy in the model problems shown here.

(e) Our studies show that accuracy ${\Vert E\Vert }_{L_2}$ comparable to space-time comparable

to space-time coupled methods is not achievable in space-time decoupled methods even after reducing $\text{Δ}t$ by a factor of 1000. This indicates that the damage done by decoupling space and time is quite difficult to recover in space-time decoupled methods, even after extreme refinements in the integration time step.

(f) This is hardly a surprise that space-time decoupled methods are inferior to space-time coupled methods, however necessity of using them for 3D initial value problem is almost unavoidable.

(g) The work presented in this paper provides a framework and a systematic procedure for computing a posteriori errors in space-time coupled and decoupled finite element methods for initial value problems without the knowledge of theoretical solutions. This framework is applicable to all initial value problems provided we follow the guidelines presented in this paper, the most important of them being use of the k-version of finite element method.

Acknowledgements

The authors are grateful for the facilities provided by the Computational Mechanics Laboratory of the Department of Mechanical Engineering. The first author is grateful to his Deane E. Ackers funds for providing support to the second author. The financial assistance provided to the second author by the Department of Mechanical Engineering is also acknowledged.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Surana, K.S. and Reddy, J.N. (2018) The Finite Element Method for Initial Value Problems. CRC/Taylor and Francis.
[2]	Shilov, G.E. (1977) Mathematical Analysis: A Special Course. Elsevier.
[3]	Surana, K.S., Ahmadi, A.R. and Reddy, J.N. (2004) The k-Version of Finite Element Method for Nonlinear Operators in BVP. International Journal of Computational Engineering Science, 5, 133-207. https://doi.org/10.1142/s1465876304002307
[4]	Surana, K.S., Ahmadi, A.R. and Reddy, J.N. (2003) The k-Version of Finite Element Method for Non-Self-Adjoint Operators in BVP. International Journal of Computational Engineering Science, 4, 737-812. https://doi.org/10.1142/s1465876303002179
[5]	Surana, K.S., Ahmadi, A.R. and Reddy, J.N. (2002) The k-Version of Finite Element Method for Self-Adjoint Operators in BVP. International Journal of Computational Engineering Science, 3, 155-218. https://doi.org/10.1142/s1465876302000605
[6]	Surana, K.S., Carranza, C.H. and Mathi, S.S.C. (2021) K-Version of Finite Element Method for BVPs and IVPs. Mathematics, 9, Article 1333. https://doi.org/10.3390/math9121333
[7]	Maduri, R.K., Surana, K.S., Romkes, A. and Reddy, J.N. (2009) Higher Order Global Differentiability Local Approximations for 2-D Distorted Triangular Elements. International Journal for Computational Methods in Engineering Science and Mechanics, 10, 20-26. https://doi.org/10.1080/15502280802572270
[8]	Ahmadi, A., Surana, K.S., Maduri, R.K., Romkes, A. and Reddy, J.N. (2009) Higher Order Global Differentiability Local Approximations for 2-D Distorted Quadrilateral Elements. International Journal for Computational Methods in Engineering Science and Mechanics, 10, 1-19. https://doi.org/10.1080/15502280802572262
[9]	Surana, K.S. and Reddy, J.N. (2016) The Finite Element Method for Boundary Value Problems. CRC/Taylor and Francis.
[10]	Ainsworth, M. and Oden, J.T. (2011) A Posteriori Error Estimation in Finite Element Analysis. Wiley.
[11]	Surana, K.S., Joy, A.D. and Reddy, J.N. (2016) Error Estimations, Error Computations, and Convergence Rates in FEM for BVPs. Applied Mathematics, 7, 1359-1407. https://doi.org/10.4236/am.2016.712120
[12]	Surana, K.S., Euler, L., Reddy, J.N. and Romkes, A. (2011) Methods of Approximation in hpk Framework for ODEs in Time Resulting from Decoupling of Space and Time in IVPs. American Journal of Computational Mathematics, 1, 83-103. https://doi.org/10.4236/ajcm.2011.12009