Direct Finite Elements in Time (DFET)

• DFET is a finite element approach that treats time as an additional dimension, enabling unified spatial-temporal discretization for matrix exponentials and nonstationary PDEs.
• It employs a fictitious-time formulation with Chebyshev polynomial-based basis functions to ensure elementwise continuity and high accuracy.
• The method supports various implementation strategies—such as all-at-once and slab solves—balancing rigorous discretization with efficient, parallelizable computation.

Searching arXiv for relevant papers on Direct Finite Elements in Time and space-time finite elements. Tool call unavailable in this environment. Proceeding with the arXiv papers explicitly provided: (0811.2612) and (Thiele, 2024). Direct Finite Elements in Time (DFET) denotes a numerical viewpoint in which temporal evolution is approximated with finite elements rather than handled solely through a derived time-stepping formula. In the matrix-exponential setting, DFET evaluates $e^A$ by introducing a fictitious time variable, defining $V(t)=e^{At}$, solving the initial value problem $\dot V(t)=A\,V(t)$ with $V(0)=I$, and recovering the desired exponential from the terminal value $V(1)=e^A$ (0811.2612). In a broader space-time finite element formulation, DFET is framed as treating time exactly like space in a variational discretization, with finite element basis functions in both space and time, particularly through tensor-product space-time slabs for nonstationary partial differential equations (Thiele, 2024).

1. Conceptual scope and historical framing

The literature represented by these two papers uses DFET at two closely related levels. In the narrower 2008 formulation, the method is a finite-element-in-time procedure for a first-order matrix ODE whose terminal value is the matrix exponential. In the 2024 formulation, DFET is the broader idea of building a space-time variational formulation and approximating the solution with finite element basis functions in both space and time. This suggests that the matrix-exponential construction can be viewed as a specialized instance of a wider space-time finite element perspective.

For nonstationary PDEs, three algorithmic viewpoints are explicitly distinguished. One is reformulation as a time-stepping scheme, in which one discretizes in space and then marches forward in time with methods such as backward Euler, Crank–Nicolson, or Runge–Kutta. A second is a native $(d+1)$-dimensional space-time discretization on a full space-time mesh in $\mathbb{R}^{d+1}$. A third, and the one emphasized for implementation, is the tensor-product space-time slab approach, in which time is discretized with one-dimensional finite elements and space with standard $d$-dimensional finite elements, and the two are combined as a tensor product (Thiele, 2024).

A recurring distinction from standard time-stepping is philosophical as well as algebraic. In the DFET/space-time FE viewpoint, temporal continuity or discontinuity is encoded by trace and jump terms in a Galerkin weak form rather than by manually deriving a one-step update or a Butcher tableau. A common misconception is therefore that DFET is simply another notation for a time-marching scheme. The cited formulation instead presents time-stepping, all-at-once space-time solves, and slab-wise solves as organizational consequences of a space-time finite element discretization, not as its defining principle.

2. Fictitious-time formulation for the matrix exponential

The 2008 construction begins from the classic computational linear algebra problem of evaluating $e^A$. Rather than approximating the exponential directly, it defines

$V(t)=e^{At}, \qquad t\in[0,1],$

so that

$V(t)=e^{At}$0

More generally, the same construction applies if one wants $V(t)=e^{At}$1 for arbitrary $V(t)=e^{At}$2, with the endpoint moved accordingly (0811.2612).

Because $V(t)=e^{At}$3, it satisfies the first-order linear matrix ODE

$V(t)=e^{At}$4

or, in component form,

$V(t)=e^{At}$5

The matrix exponential is thus recast as a fictitious-time evolution problem in which the identity matrix is propagated from $V(t)=e^{At}$6 to $V(t)=e^{At}$7.

The time interval $V(t)=e^{At}$8 is partitioned into finite elements $V(t)=e^{At}$9. On each element, the global time $\dot V(t)=A\,V(t)$0 is mapped to a local coordinate $\dot V(t)=A\,V(t)$1 by

$\dot V(t)=A\,V(t)$2

In local coordinates, the ODE takes the form

$\dot V(t)=A\,V(t)$3

equivalently written with the factor $\dot V(t)=A\,V(t)$4 moved to the other side as in the paper’s formulation (0811.2612).

The elementwise solution is decomposed as

$\dot V(t)=A\,V(t)$5

with the auxiliary function constrained by

$\dot V(t)=A\,V(t)$6

This enforces inter-element continuity: the endpoint value from element $\dot V(t)=A\,V(t)$7 becomes the left boundary value for element $\dot V(t)=A\,V(t)$8. The finite-element representation is therefore not merely local polynomial approximation; it is a propagation mechanism across temporal elements.

3. Basis construction, weak formulation, and elementwise propagation

On each time element, the auxiliary function is expanded as

$\dot V(t)=A\,V(t)$9

with basis functions built from Chebyshev polynomials of the first kind,

$V(0)=I$0

where $V(0)=I$1 denotes the Chebyshev polynomial. Since

$V(0)=I$2

the basis satisfies the required elementwise initial condition automatically, and the solution representation becomes

$V(0)=I$3

The weak formulation is obtained by inserting the expansion into the ODE, multiplying by $V(0)=I$4, and integrating over $V(0)=I$5, where the Chebyshev weight is

$V(0)=I$6

The required integrals are defined by

$V(0)=I$7

$V(0)=I$8

and

$V(0)=I$9

After substitution and rearrangement, the coefficient equations become

$V(1)=e^A$0

which the paper writes compactly as

$V(1)=e^A$1

For each time element, this is a simultaneous linear system of size $V(1)=e^A$2, where $V(1)=e^A$3 is the matrix dimension and $V(1)=e^A$4 is the number of time basis functions (0811.2612).

The numerical procedure proceeds element by element. One starts from $V(1)=e^A$5, solves the first element’s linear system for $V(1)=e^A$6, reconstructs

$V(1)=e^A$7

evaluates the endpoint $V(1)=e^A$8, and uses that value as the left boundary condition for the next element: $V(1)=e^A$9 The paper uses LU decomposition for the linear system and notes that the coefficient matrix $(d+1)$0 is independent of the step index in the implementation, so the decomposition need only be performed once; subsequent propagation is largely matrix-vector work. That structure is presented as favorable for parallel/high-performance computation.

4. Space-time variational formulation for PDEs

In the broader space-time FE framing, DFET starts from an evolution problem on a spatial domain $(d+1)$1, a time interval $(d+1)$2, and the space-time cylinder

$(d+1)$3

The functional setting is cast in terms of Bochner spaces over an evolution triple $(d+1)$4: $(d+1)$5 and, for vector-valued functions,

$(d+1)$6

A key embedding theorem quoted in the paper is

$(d+1)$7

This guarantees meaningful traces in time and underpins jump terms and weak imposition of initial conditions in dG-in-time formulations. The applicability is noted not only for $(d+1)$8, but also for $(d+1)$9 and $\mathbb{R}^{d+1}$0 (Thiele, 2024).

For the heat equation,

$\mathbb{R}^{d+1}$1

the continuous weak form is written as

$\mathbb{R}^{d+1}$2

for all $\mathbb{R}^{d+1}$3, together with $\mathbb{R}^{d+1}$4, where

$\mathbb{R}^{d+1}$5

The time interval is partitioned as

$\mathbb{R}^{d+1}$6

and the dG-in-time solution space is the broken Bochner space

$\mathbb{R}^{d+1}$7

Traces are then well defined: $\mathbb{R}^{d+1}$8 with jump

$\mathbb{R}^{d+1}$9

The dG heat-equation formulation includes both the volumetric terms and the jump terms, as well as weak enforcement of the initial condition (Thiele, 2024).

The tensor-product structure is central: $d$0 so that a space-time basis function can be written as

$d$1

with

$d$2

Time is discretized by piecewise polynomials of degree $d$3,

$d$4

and space by standard conforming Lagrangian elements

$d$5

The fully discrete tensor-product space is

$d$6

The shorthand $d$7 is used for conforming spatial degree $d$8 and discontinuous temporal degree $d$9, with the implementation focusing on dG in time (Thiele, 2024).

5. Slabs, algebraic organization, and implementation infrastructure

The algebraic organization of DFET in the space-time FE setting admits several modes. The finite element ansatz is written as

$e^A$0

or equivalently reindexed with a single global space-time index. From this viewpoint, three solve strategies are identified: an all-at-once space-time solve on the full slab, a sequential slab solve that yields a sequence of systems $e^A$1, and a compromise based on slabs consisting of several temporal elements times one spatial triangulation. The all-at-once approach is the most direct DFET interpretation but is memory-expensive; the slab compromise is presented as a way to balance memory cost and solver flexibility, and as potentially favorable for space-time multigrid or preconditioners (Thiele, 2024).

The same paper distinguishes fixed and dynamic mesh variants. In the fixed-mesh case, one forms

$e^A$2

For time-dependent spatial meshes, one may use a different triangulation $e^A$3 on each temporal interval $e^A$4, leading effectively to

$e^A$5

Conceptual support for dynamic meshes is described, but only fixed triangulations are implemented in the current version. The paper also notes that mesh-to-mesh projection between time slabs can introduce oscillations, especially for mixed methods such as Taylor–Hood Navier–Stokes where projection may destroy discrete divergence-free properties.

A distinctive implementation detail is the choice of temporal support points. The cited library supports Gauss-Lobatto, Gauss-Legendre, Gauss-Radau Left, and Gauss-Radau Right points. These choices affect temporal sparsity patterns and interface extraction. Gauss-Lobatto places support points on the boundaries, making trace values coincide with degree-of-freedom values; Gauss-Radau places one support point on a boundary and one inside; Gauss-Legendre places both inside and produces a fuller off-diagonal block structure.

The implementation layer ideal.II is designed as a C++ extension of deal.II to simplify assembly and solution of nonstationary PDEs with the tensor-product approach. Its design philosophy is to preserve the “feel” of deal.II, avoid hidden framework magic, keep application code close to stationary deal.II code, support arbitrary spatial finite elements via the FiniteElement interface, preserve extensibility and dimension-independent programming, and reuse deal.II and third-party infrastructure whenever possible. To that end it introduces space-time wrappers around Triangulation, DoFHandler, FEValues, FEFaceValues, FEJumpValues, quadrature rules, and time-iteration utilities. The DG_FiniteElement class combines a spatial FiniteElement, internally constructed temporal basis functions, a temporal polynomial degree $e^A$6, and a temporal support type, thereby creating the full space-time basis on the reference cell $e^A$7. The FEValues family provides methods such as shape_value(i,q), shape_dt(i,q), shape_space_grad(i,q), and temporal trace access through shape_value_plus() and shape_value_minus(). For global linear algebra, the ordering

$e^A$8

groups spatial degrees of freedom by temporal layer. Support for adjoint or backward-in-time solves is accommodated by storing slab collections in a std::list, and a TimeIteratorCollection keeps several iterators synchronized (Thiele, 2024).

6. Accuracy, reported behavior, advantages, and limitations

For matrix exponentials, the reported numerical results are a central part of the 2008 presentation. For a “pathological” matrix $e^A$9, with the exact exponential known, the method achieved about 13 decimal digits of accuracy using only 8 time steps and 8 basis functions. For the test matrix

$V(t)=e^{At}, \qquad t\in[0,1],$0

whose exponential is

$V(t)=e^{At}, \qquad t\in[0,1],$1

the reported accuracy was about 14 digits with the same discretization. For random real and complex test matrices $V(t)=e^{At}, \qquad t\in[0,1],$2 and $V(t)=e^{At}, \qquad t\in[0,1],$3, the results showed saturation at roughly 12–13 significant digits with moderate discretization. The tables indicate that, for the examples studied, about 5–8 time steps and basis functions were often sufficient to reach near machine-precision accuracy for many entries. The paper also notes that complex matrices may require somewhat more work, and that the basis size can matter more strongly than the number of time steps in some cases (0811.2612).

The matrix-exponential paper presents several advantages of the DFET strategy. Accuracy is systematically controllable through the number of time elements and basis functions. The method is described as robust for both diagonalizable and non-diagonalizable matrices. It has a parallel structure because each element leads to a structured linear solve, and it does not rely on eigen-decomposition, which can be problematic for defective or nearly defective matrices. More generally, the approach is framed as exploiting the mature numerical machinery available for initial-value problems in order to address the lack of a universally satisfactory matrix exponential algorithm.

In the PDE setting, the 2024 paper validates the tensor-product dG-in-time implementation with a manufactured heat-equation solution and reports expected $V(t)=e^{At}, \qquad t\in[0,1],$4-convergence under spatial refinement, expected $V(t)=e^{At}, \qquad t\in[0,1],$5-convergence under temporal refinement depending on $V(t)=e^{At}, \qquad t\in[0,1],$6, and combined-refinement behavior limited by the lower-order component. Different temporal support types produce slightly different errors, with Gauss-Radau Left often best among the tested choices. For the 2D-3 flow-around-a-cylinder Navier–Stokes benchmark with time-dependent inflow, the resulting drag and lift curves are reported as qualitatively similar to FEATFLOW, with dG(0) showing stronger damping and higher-order temporal discretizations giving less smearing (Thiele, 2024).

The broader benefits and challenges are correspondingly balanced. Benefits include a natural formulation for time-dependent PDEs as a space-time Galerkin problem, unified handling of spatial and temporal discretization, compatibility with all-at-once, slab, or hybrid solve strategies, support for adjoint-based methods and goal-oriented adaptivity, and the possibility of space-time preconditioners and multigrid. The approach also removes the need to derive time integrators by hand. Challenges include the size and memory cost of global space-time systems, the need for mesh-to-mesh transfers on dynamic meshes, the infrastructure required for traces, jumps, and support-point-dependent sparsity patterns, and the interaction of spatial and temporal approximation errors when choosing $V(t)=e^{At}, \qquad t\in[0,1],$7 and $V(t)=e^{At}, \qquad t\in[0,1],$8. A further misconception is that a space-time FE method must be an all-at-once global solve; the cited work explicitly treats sequential slab solves as a legitimate algebraic organization emerging from the same DFET weak formulation rather than as a departure from it (Thiele, 2024).