Skip to content
BY 4.0 license Open Access Published by De Gruyter November 11, 2020

An Exact Realization of a Modified Hilbert Transformation for Space-Time Methods for Parabolic Evolution Equations

  • Marco Zank ORCID logo EMAIL logo

Abstract

We present different possibilities of realizing a modified Hilbert type transformation as it is used for Galerkin–Bubnov discretizations of space-time variational formulations for parabolic evolution equations in anisotropic Sobolev spaces of spatial order 1 and temporal order 1 2 . First, we investigate the series expansion of the definition of the modified Hilbert transformation, where the truncation parameter has to be adapted to the mesh size. Second, we introduce a new series expansion based on the Legendre chi function to calculate the corresponding matrices for piecewise polynomial functions. With this new procedure, the matrix entries for a space-time finite element method for parabolic evolution equations are computable to machine precision independently of the mesh size. Numerical results conclude this work.

MSC 2010: 65M12; 65M15; 65M60

1 Introduction

For the discretization of parabolic evolution equations, the classical approaches are time stepping schemes together with finite element methods in space. An alternative is to discretize the parabolic problem without separating the temporal and spatial variables, i.e. space-time methods. In general, the main advantages of space-time methods are space-time adaptivity, space-time parallelization and the treatment of moving boundaries. However, space-time approximation methods depend strongly on the space-time variational formulations on the continuous level. On the one hand, there are space-time discretizations of parabolic evolution equations based on the variational formulations in Bochner–Sobolev spaces, see, e.g., [1, 5, 8, 10, 11, 16, 17, 18, 21, 24]. On the other hand, discretizations of variational formulations in anisotropic Sobolev spaces of spatial order 1 and temporal order 1 2 became quite attractive recently, see, e.g., [4, 12, 19, 23, 25]. In this work, the approach in these anisotropic Sobolev spaces is applied. This type of space-time variational formulations allows the complete analysis of inhomogeneous Dirichlet or Neumann conditions and is used for the analysis of the resulting boundary integral operators, see [3]. Hence, discretizations of variational formulations in these anisotropic Sobolev spaces can be used for the interior problems of FEM-BEM couplings for transmission problems.

The model problem for a parabolic evolution equation is the homogeneous Dirichlet problem of the heat equation

(1.1) { t u ( x , t ) - Δ x u ( x , t ) = f ( x , t ) for  ( x , t ) Q := Ω × ( 0 , T ) , u ( x , t ) = 0 for  ( x , t ) Σ := Ω × ( 0 , T ) , u ( x , 0 ) = 0 for  x Ω ,

where Ω d , d = 1 , 2 , 3 , is a bounded Lipschitz domain with boundary Ω , T > 0 is a given terminal time and f is a given right-hand side. Next, we consider the space-time variational formulation of (1.1) to find u H 0 ;  0 , 1 , 1 / 2 ( Q ) such that

(1.2) a ( u , v ) = f , v Q

for all v H 0 ; , 0 1 , 1 / 2 ( Q ) , where f [ H 0 ; , 0 1 , 1 / 2 ( Q ) ] is a given right-hand side. Here, the bilinear form

a ( , ) : H 0 ;  0 , 1 , 1 / 2 ( Q ) × H 0 ; , 0 1 , 1 / 2 ( Q )

is defined by

a ( u , v ) := t u , v Q + x u , x v L 2 ( Q ) , u H 0 ;  0 , 1 , 1 / 2 ( Q ) , v H 0 ; , 0 1 , 1 / 2 ( Q ) ,

and , Q denotes the duality pairing as extension of the inner product in L 2 ( Q ) . Furthermore, the anisotropic Sobolev spaces

H 0 ;  0 , 1 , 1 / 2 ( Q ) := H 0 , 1 / 2 ( 0 , T ; L 2 ( Ω ) ) L 2 ( 0 , T ; H 0 1 ( Ω ) ) ,
H 0 ; , 0 1 , 1 / 2 ( Q ) := H , 0 1 / 2 ( 0 , T ; L 2 ( Ω ) ) L 2 ( 0 , T ; H 0 1 ( Ω ) )

are endowed with the Hilbertian norms

v H 0 ;  0 , 1 , 1 / 2 ( Q ) := v H 0 , 1 / 2 ( 0 , T ; L 2 ( Ω ) ) 2 + x v L 2 ( Q ) 2 ,
w H 0 ; , 0 1 , 1 / 2 ( Q ) := w H , 0 1 / 2 ( 0 , T ; L 2 ( Ω ) ) 2 + x w L 2 ( Q ) 2

with the usual Bochner–Sobolev norms

v H 0 , 1 / 2 ( 0 , T ; L 2 ( Ω ) ) := v H 1 / 2 ( 0 , T ; L 2 ( Ω ) ) 2 + 0 T v ( , t ) L 2 ( Ω ) 2 t d t ,
w H , 0 1 / 2 ( 0 , T ; L 2 ( Ω ) ) := w H 1 / 2 ( 0 , T ; L 2 ( Ω ) ) 2 + 0 T w ( , t ) L 2 ( Ω ) 2 T - t d t

and the usual Bochner–Sobolev spaces

H 0 , 1 / 2 ( 0 , T ; L 2 ( Ω ) ) := { v H 1 / 2 ( 0 , T ; L 2 ( Ω ) ) : v H 0 , 1 / 2 ( 0 , T ; L 2 ( Ω ) ) < } ,
H , 0 1 / 2 ( 0 , T ; L 2 ( Ω ) ) := { w H 1 / 2 ( 0 , T ; L 2 ( Ω ) ) : w H , 0 1 / 2 ( 0 , T ; L 2 ( Ω ) ) < } ,

see [15, 14, 23, 25] for more details. Moreover, the dual space [ H 0 ; , 0 1 , 1 / 2 ( Q ) ] is characterized as completion of L 2 ( Q ) with respect to the Hilbertian norm

f [ H 0 ; , 0 1 , 1 / 2 ( Q ) ] := sup 0 w H 0 ; , 0 1 , 1 / 2 ( Q ) | f , w Q | w H 0 ; , 0 1 , 1 / 2 ( Q ) .

With these Sobolev spaces, the bilinear form a ( , ) is bounded, i.e. there exists a constant C > 0 such that

| a ( u , v ) | C u H 0 ;  0 , 1 , 1 / 2 ( Q ) v H 0 ; , 0 1 , 1 / 2 ( Q ) for all  u H 0 ;  0 , 1 , 1 / 2 ( Q )  and all  v H 0 ; , 0 1 , 1 / 2 ( Q ) ,

which is proven by a density argument, see [3, Lemma 2.6, p. 505] and [25, Remark 3.3.1, p. 65]. In [3], the following existence and uniqueness theorem is proven by a transposition and interpolation argument as in [15, 14], see also [9].

Theorem 1.1.

Let the right-hand side f [ H 0 ; , 0 1 , 1 / 2 ( Q ) ] be given. Then the variational formulation (1.2) has a unique solution u H 0 ;  0 , 1 , 1 / 2 ( Q ) satisfying

u H 0 ;  0 , 1 , 1 / 2 ( Q ) C f [ H 0 ; , 0 1 , 1 / 2 ( Q ) ]

with a constant C > 0 . Furthermore, the solution operator

: [ H 0 ; , 0 1 , 1 / 2 ( Q ) ] H 0 ;  0 , 1 , 1 / 2 ( Q ) , f := u ,

is an isomorphism.

For a discretization scheme, let the bounded Lipschitz domain Ω d be an interval Ω = ( 0 , L ) for d = 1 , or polygonal for d = 2 , or polyhedral for d = 3 . For a tensor-product ansatz, we consider admissible decompositions

Q ¯ = Ω ¯ × [ 0 , T ] = i = 1 N x ω i ¯ × = 1 N t [ t - 1 , t ]

with N := N x N t space-time elements, where the time intervals ( t - 1 , t ) with mesh sizes h t , = t - t - 1 are defined via the decomposition

(1.3) 0 = t 0 < t 1 < t 2 < < t N t - 1 < t N t = T

of the time interval ( 0 , T ) . The maximal and the minimal time mesh sizes are denoted by

h t := h t , max := max h t , and h t , min := min h t , .

For the spatial domain Ω, we consider a shape-regular sequence ( 𝒯 ν ) ν of admissible decompositions 𝒯 ν := { ω i d : i = 1 , , N x } of Ω into finite elements ω i d with mesh sizes h x , i and the maximal mesh size h x := max i h x , i . The spatial elements ω i are intervals for d = 1 , triangles or quadrilaterals for d = 2 , and tetrahedra or hexahedra for d = 3 . Next, we introduce the finite element space

(1.4) Q h 1 ( Q ) := V h x , 0 1 ( Ω ) S h t 1 ( 0 , T )

of piecewise multilinear, continuous functions, i.e.

V h x , 0 1 ( Ω ) = span { ψ j 1 } j = 1 M x H 0 1 ( Ω ) , S h t 1 ( 0 , T ) = span { φ 1 } = 0 N t H 1 ( 0 , T ) .

In fact, V h x , 0 1 ( Ω ) is either the space S h x 1 ( Ω ) H 0 1 ( Ω ) of piecewise linear, continuous functions on intervals ( d = 1 ), triangles ( d = 2 ), and tetrahedra ( d = 3 ), or V h x , 0 1 ( Ω ) is the space Q h x 1 ( Ω ) H 0 1 ( Ω ) of piecewise linear/bilinear/trilinear, continuous functions on intervals ( d = 1 ), quadrilaterals ( d = 2 ), and hexahedra ( d = 3 ). Analogously, for a fixed polynomial degree p , we consider the space of piecewise polynomial, continuous functions

(1.5) Q h p ( Q ) := V h x , 0 p ( Ω ) S h t p ( 0 , T ) .

By using the finite element space (1.4), it turns out that a discretization of (1.2) with the conforming ansatz space Q h 1 ( Q ) H 0 ;  0 , 1 , 1 / 2 ( Q ) and the conforming test space Q h 1 ( Q ) H 0 ; , 0 1 , 1 / 2 ( Q ) is not stable, see [25, Section 3.3]. A possible way out is the modified Hilbert transformation T defined by

( T u ) ( x , t ) := i = 1 k = 0 u i , k cos ( ( π 2 + k π ) t T ) ϕ i ( x ) , ( x , t ) Q ,

where the given function u L 2 ( Q ) is represented by

u ( x , t ) = i = 1 k = 0 u i , k sin ( ( π 2 + k π ) t T ) ϕ i ( x ) , ( x , t ) Q ,

with the eigenfunctions ϕ i H 0 1 ( Ω ) and eigenvalues μ i , satisfying

- Δ ϕ i = μ i ϕ i in  Ω , ϕ i = 0 on  Ω , ϕ i L 2 ( Ω ) = 1 , i .

This approach was introduced in [23] and [25, Section 3.4], see also [4, 6, 7, 12] for analogous considerations for an infinite time interval ( 0 , ) with the classical Hilbert transformation, which is related to . The map

T : H 0 ;  0 , 1 , 1 / 2 ( Q ) H 0 ; , 0 1 , 1 / 2 ( Q )

is norm preserving, bijective and fulfills the coercivity property

(1.6) t v , T v Q v H 0 , 1 / 2 ( 0 , T ; L 2 ( Ω ) ) 2 for all  v H 0 , 1 / 2 ( 0 , T ; L 2 ( Ω ) ) .

Moreover, the relation

(1.7) v , T v L 2 ( Q ) 0 for all  v L 2 ( Q )

holds true. With the modified Hilbert transformation T , the variational formulation (1.2) is equivalent to find u H 0 ;  0 , 1 , 1 / 2 ( Q ) such that

(1.8) a ( u , T v ) = f , T v Q for all  v H 0 ;  0 , 1 , 1 / 2 ( Q ) .

Hence, unique solvability of the variational formulation (1.8) follows from the unique solvability of (1.2). Thus, Theorem 1.1 and the properties of T give the stability estimate

c u H 0 ;  0 , 1 , 1 / 2 ( Q ) sup 0 v H 0 ;  0 , 1 , 1 / 2 ( Q ) | a ( u , T v ) | v H 0 ;  0 , 1 , 1 / 2 ( Q ) for all  u H 0 ;  0 , 1 , 1 / 2 ( Q )

with a constant c > 0 . When using some conforming space-time finite element space 𝒱 h H 0 ;  0 , 1 , 1 / 2 ( Q ) , the Galerkin variational formulation of (1.8) is to find u h 𝒱 h such that

(1.9) a ( u h , T v h ) = f , T v h Q for all  v h 𝒱 h .

Note that ansatz and test spaces are equal. With the coercivity property (1.6) and property (1.7), there exists a constant c > 0 such that

(1.10) a ( v h , T v h ) c v h H 0 , 1 / 2 ( 0 , T ; L 2 ( Ω ) ) 2 for all  v h 𝒱 h ,

which leads to the following theorem, where its proof is contained in [25].

Theorem 1.2.

Let V h H 0 ;  0 , 1 , 1 / 2 ( Q ) be a conforming space-time finite element space and let f [ H 0 ; , 0 1 , 1 / 2 ( Q ) ] be a given right-hand side. Then, a unique solution u h V h of the Galerkin–Bubnov variational formulation (1.9) exists. If, in addition, the right-hand side fulfills f [ H , 0 1 / 2 ( 0 , T ; L 2 ( Ω ) ) ] [ H 0 ; , 0 1 , 1 / 2 ( Q ) ] , then the stability estimate

u h H 0 , 1 / 2 ( 0 , T ; L 2 ( Ω ) ) c f [ H , 0 1 / 2 ( 0 , T ; L 2 ( Ω ) ) ]

is true with a constant c > 0 .

Theorem 1.2 states that, under the assumption f [ H , 0 1 / 2 ( 0 , T ; L 2 ( Ω ) ) ] , any conforming space-time finite element space 𝒱 h H 0 ;  0 , 1 , 1 / 2 ( Q ) leads to an unconditionally stable method, i.e. no CFL condition is required. For the choice of the tensor-product space-time finite element space

𝒱 h = Q h p ( Q ) H 0 ;  0 , 1 , 1 / 2 ( Q )

with (1.5), the Galerkin–Bubnov variational formulation (1.9) to find u h Q h p ( Q ) H 0 ;  0 , 1 , 1 / 2 ( Q ) such that

(1.11) a ( u h , T v h ) = f , T v h Q for all  v h Q h p ( Q ) H 0 ;  0 , 1 , 1 / 2 ( Q )

fulfills the space-time error estimates

u - u h H 0 , 1 / 2 ( 0 , T ; L 2 ( Ω ) ) c h p + 1 / 2 ,
(1.12) u - u h L 2 ( Q ) c h p + 1 ,
(1.13) | u - u h | H 1 ( Q ) c h p

with h = max { h t , h x } and with a constant c > 0 , see [23, 25]. Here, we have to assume that the solution u of (1.2) is sufficiently smooth and that Ω is sufficiently regular, e.g., convex, such that the extended H 0 1 ( Ω ) projection Q h x p : L 2 ( 0 , T ; H 0 1 ( Ω ) ) V h x , 0 p ( Ω ) L 2 ( 0 , T ) , defined for a function w L 2 ( 0 , T ; H 0 1 ( Ω ) ) by

x Q h x p w , x w h x L 2 ( Q ) = x w , x w h x L 2 ( Q ) for all  w h x V h x , 0 p ( Ω ) L 2 ( 0 , T ) ,

fulfills the standard error estimate

v - Q h x p v L 2 ( Q ) c h x p + 1 v L 2 ( 0 , T ; H p + 1 ( Ω ) ) for all  v L 2 ( 0 , T ; H 0 1 ( Ω ) H p + 1 ( Ω ) )

with a constant c > 0 . For the H 1 ( Q ) error estimate (1.13), the sequence ( 𝒯 ν ) ν of decompositions of Ω is additionally assumed to be globally quasi-uniform, see [23, 25] for details.

In the remainder of this work, we consider p = 1 , i.e. the tensor-product space of piecewise multilinear, continuous functions 𝒱 h = Q h 1 ( Q ) H 0 ;  0 , 1 , 1 / 2 ( Q ) , where analogous results hold true for an arbitrary polynomial degree p > 1 . Furthermore, for f L 2 ( Q ) , we approximate the right-hand side f by

(1.14) f Q h 0 f = j = 1 N x = 1 N t f j , ψ j 0 φ 0 S h x 0 ( Ω ) S h t 0 ( 0 , T )

with coefficients f j , , where Q h 0 : L 2 ( Q ) S h x 0 ( Ω ) S h t 0 ( 0 , T ) is the L 2 ( Q ) projection on the piecewise constant functions S h x 0 ( Ω ) S h t 0 ( 0 , T ) with

S h x 0 ( Ω ) = span { ψ j 0 } j = 1 N x and S h t 0 ( 0 , T ) = span { φ 0 } = 1 N t .

So, we consider the perturbed variational formulation to find u ~ h Q h 1 ( Q ) H 0 ;  0 , 1 , 1 / 2 ( Q ) such that

(1.15) a ( u ~ h , T v h ) = Q h 0 f , T v h L 2 ( Q ) for all  v h Q h 1 ( Q ) H 0 ;  0 , 1 , 1 / 2 ( Q ) .

The discrete variational formulation (1.15) is equivalent to the global linear system

(1.16) K h u ¯ ~ = F ¯ ~ T

with the system matrix

K h = A h t T M h x + M h t T A h x N t M x × N t M x ,

where M h x M x × M x and A h x M x × M x denote spatial mass and stiffness matrices given by

M h x [ i , j ] = ψ j 1 , ψ i 1 L 2 ( Ω ) , A h x [ i , j ] = x ψ j 1 , x ψ i 1 L 2 ( Ω ) , i , j = 1 , , M x ,

and M h t T N t × N t and A h t T N t × N t are defined by

(1.17) M h t T [ , k ] := φ k 1 , T φ 1 L 2 ( 0 , T ) , A h t T [ , k ] := t φ k 1 , T φ 1 L 2 ( 0 , T )

for , k = 1 , , N t . Here, the modified Hilbert transformation T : L 2 ( 0 , T ) L 2 ( 0 , T ) , acting on solely time-dependent functions, is given as

(1.18) ( T w ) ( t ) := η = 0 w η cos ( ( π 2 + η π ) t T ) , t ( 0 , T ) ,

where the given function w L 2 ( 0 , T ) is represented by

(1.19) w ( t ) = η = 0 w η sin ( ( π 2 + η π ) t T ) , t ( 0 , T ) ,

with the coefficients

w η = 2 T 0 T w ( s ) sin ( ( π 2 + η π ) s T ) d s .

We use the same notation T for solely time-dependent functions and functions, which depend on ( x , t ) , since for a function u L 2 ( Q ) with u ( x , t ) = v ( x ) w ( t ) , v L 2 ( Ω ) , w L 2 ( 0 , T ) , the equality

T u ( x , t ) = v ( x ) T w ( t ) , ( x , t ) Q ,

holds true. Next, we state some properties of T , acting on solely time-dependent functions. For the Sobolev spaces

H 0 , 1 / 2 ( 0 , T ) = { v H 1 / 2 ( 0 , T ) : v H 0 , 1 / 2 ( 0 , T ) < } ,
H , 0 1 / 2 ( 0 , T ) = { w H 1 / 2 ( 0 , T ) : w H , 0 1 / 2 ( 0 , T ) < }

with the Hilbertian norms

v H 0 , 1 / 2 ( 0 , T ) := v H 1 / 2 ( 0 , T ) 2 + 0 T | v ( t ) | 2 t d t ,
w H , 0 1 / 2 ( 0 , T ) := w H 1 / 2 ( 0 , T ) 2 + 0 T | w ( t ) | 2 T - t d t ,

the modified Hilbert transformation T , as given in (1.18), is an isomorphism

T : H 0 , 1 / 2 ( 0 , T ) H , 0 1 / 2 ( 0 , T )

or an isomorphism

T : H 0 , 1 ( 0 , T ) H , 0 1 ( 0 , T ) ,

where the latter spaces are defined accordingly. As shown in [23, 25], for all u , v H 0 , 1 / 2 ( 0 , T ) with expansion coefficients u η , v η as in (1.19),

t u , T v ( 0 , T ) = 1 2 η = 0 ( π 2 + η π ) u η v η = : u , v H 0 , 1 / 2 ( 0 , T ) , F

defines an inner product in H 0 , 1 / 2 ( 0 , T ) , where , ( 0 , T ) denotes the duality pairing in [ H , 0 1 / 2 ( 0 , T ) ] and H , 0 1 / 2 ( 0 , T ) as extension of the inner product in L 2 ( 0 , T ) . With this notation, the coercivity property

(1.20) t v , T v ( 0 , T ) = v H 0 , 1 / 2 ( 0 , T ) , F 2 v H 0 , 1 / 2 ( 0 , T ) 2 for all  v H 0 , 1 / 2 ( 0 , T )

is proven in [23, 25], where the norm H 0 , 1 / 2 ( 0 , T ) , F is induced by the inner product , H 0 , 1 / 2 ( 0 , T ) , F . Additionally, the relation

(1.21) v , T v L 2 ( 0 , T ) 0 for all  v L 2 ( 0 , T )

holds true.

With the coercivity property (1.20), the matrix A h t T in (1.17) is symmetric and positive definite, whereas the matrix M h t T is nonsymmetric and positive semi-definite, see (1.21). Furthermore, the vector of the right-hand side in (1.16) is given by

F ¯ ~ T := ( f ¯ ~ 1 , , f ¯ ~ N t ) N t M x

with the vectors f ¯ ~ k M x , k = 1 , , N t , where, with the help of (1.14),

f ¯ ~ k [ i ] := Q h 0 f , ψ i 1 T φ k 1 L 2 ( Q ) = j = 1 N x = 1 N t f j , ψ j 0 , ψ i 1 L 2 ( Ω ) φ 0 , T φ k 1 L 2 ( 0 , T ) , i = 1 , , M x .

To assemble the vector of the right-hand side in (1.16), the relation

f ¯ ~ k [ i ] := F ~ [ i , k ] , i = 1 , , M x , k = 1 , , N t ,

holds true with

F ~ := M h x 1 , 0 F ( C h t T ) M x × N t ,

where

M h x 1 , 0 [ i , j ] := ψ j 0 , ψ i 1 L 2 ( Ω ) , i = 1 , , M x , j = 1 , , N x ,
F [ j , ] := f j , , j = 1 , , N x , = 1 , , N t ,
(1.22) C h t T [ k , ] := φ 0 , T φ k 1 L 2 ( 0 , T ) , k = 1 , , N t , = 1 , , N t .

Hence, to assemble the global linear system (1.16), a realization of the modified Hilbert transformation T acting on piecewise linear, continuous functions, is needed, i.e. the assembling of the matrices M h t T , A h t T and C h t T . One possibility is introduced in [22], where the modified Hilbert transformation T is represented as an integral operator and the matrices M h t T , A h t T , C h t T are assembled by the use of numerical integration, i.e. quadrature errors perturb the entries of the matrices M h t T , A h t T , C h t T . Another possibility of a realization of the modified Hilbert transformation T is to truncate the series expansion (1.18) of the definition of T , where the convergence is slow and depends on the time mesh size h t , min , see numerical examples in Section 3. In this work, we propose a new series representation of the modified Hilbert transformation T for piecewise polynomial functions, which converges very fast independently of the time mesh size h t , min , i.e. only a few terms in this new series are needed to calculate the matrix entries of M h t T , A h t T and C h t T to machine precision.

The rest of the paper is organized as follows: In Section 2 we show different possibilities to realize the modified Hilbert transformation T . First, we use the series expansion of the definition of T for assembling M h t T , A h t T and C h t T . Second, a new series representation of the entries of the matrices M h t T , A h t T and C h t T is derived and analyzed. In Section 3 numerical examples show the quality of the new assembling method of M h t T , A h t T and C h t T . In addition, a numerical example for the heat equation in a two-dimensional spatial domain is shown. In Section 4 we give some conclusions.

2 Realizations of the Modified Hilbert Transformation T

In this section, we consider realizations of the modified Hilbert transformation T to compute the matrices A h t T , M h t T , C h t T in (1.17) and (1.22). Hence, only the temporal part of the global linear system (1.16) is investigated. In this section, the space S h t 1 ( 0 , T ) H 0 , 1 / 2 ( 0 , T ) = span { φ k 1 } k = 1 N t of the piecewise linear basis functions

(2.1) φ k 1 ( t ) = { 1 h t , k ( t - t k - 1 ) , t ( t k - 1 , t k ] , 1 h t , k + 1 ( t k + 1 - t ) , t ( t k , t k + 1 ) , 0 , otherwise ,

for k = 1 , , N t - 1 and

(2.2) φ N t 1 ( t ) = { 1 h t , N t ( t - t N t - 1 ) , t ( t N t - 1 , t N t ] , 0 , otherwise ,

is considered as a special case, where a generalization to a polynomial degree p > 1 or high-order splines is straightforward.

2.1 Series Representation of the Definition

In this subsection, we use the definition (1.18) of T to calculate approximations of the matrices A h t T , M h t T , C h t T in (1.17) and (1.22). For this purpose, we represent the piecewise linear basis functions (2.1), (2.2) by

(2.3) φ k 1 ( t ) = η = 0 φ k , η sin ( ( π 2 + η π ) t T ) , k = 1 , , N t ,

with the expansion coefficients

φ k , η = 2 T 0 T φ k 1 ( t ) sin ( ( π 2 + η π ) t T ) d t .

For η 0 , the expansion coefficients are given by

φ k , η = 2 T 0 T φ k 1 ( t ) sin ( ( π 2 + η π ) t T ) d t
= 8 T ( - ( h t , k + h t , k + 1 ) sin ( π ( 2 η + 1 ) t k 2 T ) + h t , k + 1 sin ( π ( 2 η + 1 ) t k - 1 2 T ) + h t , k sin ( π ( 2 η + 1 ) t k + 1 2 T ) ) - ( 2 π η + π ) 2 h t , k h t , k + 1

for k = 1 , , N t - 1 and by

φ N t , η = 2 T 0 T φ N t 1 ( t ) sin ( ( π 2 + η π ) t T ) d t = - 8 T ( sin ( π ( 2 η + 1 ) t N t - 1 2 T ) + ( - 1 ) η + 1 ) ( 2 π η + π ) 2 h t , N t .

Analogously, for the expansion

(2.4) φ k 1 ( t ) = η = 0 φ ^ k , η cos ( ( π 2 + η π ) t T ) , k = 1 , , N t ,

the expansion coefficients

φ ^ k , η = 2 T 0 T φ k 1 ( t ) cos ( ( π 2 + η π ) t T ) d t

are

φ ^ k , η = 8 T ( - ( h t , k + h t , k + 1 ) cos ( π ( 2 η + 1 ) t k 2 T ) + h t , k + 1 cos ( π ( 2 η + 1 ) t k - 1 2 T ) + h t , k cos ( π ( 2 η + 1 ) t k + 1 2 T ) ) - ( 2 π η + π ) 2 h t , k h t , k + 1

for k = 1 , , N t - 1 , and

φ ^ N t , η = - 8 T cos ( π ( 2 η + 1 ) t N t - 1 2 T ) + 4 h t , N t ( 2 π η + π ) ( - 1 ) η ( 2 π η + π ) 2 h t , N t .

Using the expansions (2.3), (2.4), the matrix entries (1.17) are given by

(2.5) A h t T [ k , ] = η = 0 a , η a k , η and M h t T [ k , ] = η = 0 b , η a k , η

for k , = 1 , , N t , where

a k , η := 2 π η + π 2 φ k , η and b k , η := T 2 π η + π φ ^ k , η .

For the matrix C h t T in (1.22), analogous results hold true, which we omit here since A h t T and C h t T are related, see (2.15) and (2.17).

2.1.1 Approximation of the Matrices A h t T and M h t T

With the help of the series representation (2.5), approximations of the matrices A h t T and M h t T are given by

(2.6) A ~ h t T [ k , ] := η = 0 ρ a , η a k , η and M ~ h t T [ k , ] := η = 0 ρ b , η a k , η

for k , = 1 , , N t with the truncation parameter ρ 0 . Note that the approximate matrix A ~ h t T is also symmetric. For each matrix entry (2.6) of A ~ h t T , the errors are estimated by

| A h t T [ k , ] - A ~ h t T [ k , ] | = | η = ρ + 1 a , η a k , η | η = ρ + 1 8 T ( h t , k + h t , k + 1 ) ( 2 π η + π ) 3 / 2 h t , k h t , k + 1 8 T ( h t , + h t , + 1 ) ( 2 π η + π ) 3 / 2 h t , h t , + 1
= ( h t , k + h t , k + 1 ) ( h t , + h t , + 1 ) h t , k h t , k + 1 h t , h t , + 1 η = ρ + 1 64 T 2 ( 2 π η + π ) 3
(2.7) ( h t , k + h t , k + 1 ) ( h t , + h t , + 1 ) h t , k h t , k + 1 h t , h t , + 1 16 T 2 π 3 ( 2 ρ + 1 ) 2

for k , = 1 , , N t - 1 , and

(2.8) | A h t T [ N t , N t ] - A ~ h t T [ N t , N t ] | 1 h t , N t 2 16 T 2 π 3 ( 2 ρ + 1 ) 2 ,
(2.9) | A h t T [ N t , ] - A ~ h t T [ N t , ] | = | A h t T [ , N t ] - A ~ h t T [ , N t ] | h t , + h t , + 1 h t , N t h t , h t , + 1 16 T 2 π 3 ( 2 ρ + 1 ) 2

for = 1 , , N t - 1 . Analogously, for M ~ h t T , the errors are estimated by

(2.10) | M h t T [ k , ] - M ~ h t T [ k , ] | ( h t , k + h t , k + 1 ) ( h t , + h t , + 1 ) h t , k h t , k + 1 h t , h t , + 1 64 T 3 3 π 4 ( 2 ρ + 1 ) 3

for k , = 1 , , N t - 1 , and

(2.11) | M h t T [ N t , N t ] - M ~ h t T [ N t , N t ] | 8 T 2 1 h t , N t 2 3 h t , N t ( 2 π ρ + π ) + 4 T 3 π 4 ( 2 ρ + 1 ) 3 ,
(2.12) | M h t T [ , N t ] - M ~ h t T [ , N t ] | 8 T 2 h t , + h t , + 1 h t , N t h t , h t , + 1 3 h t , N t ( 2 π ρ + π ) + 4 T 3 π 4 ( 2 ρ + 1 ) 3 ,
(2.13) | M h t T [ N t , ] - M ~ h t T [ N t , ] | h t , + h t , + 1 h t , h t , + 1 h t , N t 64 T 3 3 π 4 ( 2 ρ + 1 ) 3

for = 1 , , N t - 1 .

Since the stability and error analysis of the Galerkin–Bubnov scheme (1.15) is based on the coercivity property (1.10), which is related to the positive definiteness of A h t T , we have to prove the positive definiteness of the approximate matrix A ~ h t T , which is done in the following. First, the estimate of the approximation error

| ( ( A h t T - A ~ h t T ) u ¯ , v ¯ ) | k = 1 N t = 1 N t | A h t T [ k , ] - A ~ h t T [ k , ] | | u | | v k | 192 T 3 h t , min 4 1 π 3 ( 2 ρ + 1 ) 2 u h t L 2 ( 0 , T ) v h t L 2 ( 0 , T )

holds true for all N t u ¯ , v ¯ u h t , v h t S h t 1 ( 0 , T ) H 0 , 1 / 2 ( 0 , T ) , where the error estimates (2.7), (2.8), (2.9) and the Cauchy–Schwarz inequality are used. With this last estimate, the coercivity property (1.20) and the Poincaré inequality [25, Lemma 3.4.5], it follows that

( A ~ h t T u ¯ , u ¯ ) = ( A h t T u ¯ , u ¯ ) + ( ( A h t T - A ~ h t T ) u ¯ , u ¯ )
= t u h t , T u h t L 2 ( 0 , T ) + ( ( A h t T - A ~ h t T ) u ¯ , u ¯ )
u h t H 0 , 1 / 2 ( 0 , T ) , F 2 - 192 T 3 h t , min 4 1 π 3 ( 2 ρ + 1 ) 2 u h t L 2 ( 0 , T ) 2
( 1 - 384 T 4 h t , min 4 1 π 4 ( 2 ρ + 1 ) 2 ) u h t H 0 , 1 / 2 ( 0 , T ) , F 2

for all N t u ¯ u h t S h t 1 ( 0 , T ) H 0 , 1 / 2 ( 0 , T ) , i.e. we can prove that the approximate matrix A ~ h t T is positive definite for

(2.14) ρ > 4 6 T 2 π 2 h t , min 2 - 1 2 .

Numerical examples show that estimate (2.14) is not sharp. It seems that the truncation parameter ρ depends linearly on the time mesh size h t , min , see Table 4. However, in the next subsection, we propose a new possibility to avoid the approximate matrices (2.6), which allows to calculate the matrix entries of the matrices A h t T , M h t T , C h t T in (1.17) and (1.22) to machine precision independently of the mesh size h t , min .

2.2 New Series Representation via the Legendre Chi Function

In this subsection, we introduce a new possibility to calculate the matrices A h t T , M h t T , C h t T in (1.17) and (1.22). For this purpose, the piecewise linear basis functions (2.1), (2.2) are represented as

φ k 1 ( t ) = 1 h t , k ( α 1 ( t ) - t k - 1 ) φ k 0 ( t ) + 1 h t , k + 1 ( t k + 1 - α 1 ( t ) ) φ k + 1 0 ( t )

for k = 1 , , N t - 1 and

φ N t 1 ( t ) = 1 h t , N t ( α 1 ( t ) - t N t - 1 ) φ N t 0 ( t )

with the piecewise constant functions φ 0 , = 1 , , N t , where α 1 ( t ) := t . Analogously, the derivative of the functions φ k 1 is

t φ k 1 ( t ) = 1 h t , k φ k 0 ( t ) - 1 h t , k + 1 φ k + 1 0 ( t )

for k = 1 , , N t - 1 and

t φ N t 1 ( t ) = 1 h t , N t φ N t 0 ( t ) .

With these representations, the matrices A h t T , M h t T , C h t T in (1.17) and (1.22) are given by

(2.15) A h t T = [ Z h t , 1 A h t 1 , 0 + Z h t , 0 A h t 0 , 0 ] Z h t , 1 ,
(2.16) M h t T = [ Z h t , 1 A h t 1 , 1 + Z h t , 0 A h t 0 , 1 ] Z h t , 1 + [ Z h t , 1 A h t 1 , 0 + Z h t , 0 A h t 0 , 0 ] Z h t , 0 ,
(2.17) C h t T = Z h t , 1 A h t 1 , 0 + Z h t , 0 A h t 0 , 0

with the assembling matrices

(2.18) Z h t , 1 := ( 1 h t , 1 - 1 h t , 2 1 h t , 2 - 1 h t , 3 1 h t , N t - 1 - 1 h t , N t 1 h t , N t ) N t × N t ,
(2.19) Z h t , 0 := ( - t 0 h t , 1 t 2 h t , 2 - t 1 h t , 2 t 3 h t , 3 - t N t - 2 h t , N t - 1 t N t h t , N t - t N t - 1 h t , N t ) N t × N t .

Here, the auxiliary matrix A h t r , q N t × N t is defined by

(2.20) A h t r , q [ k , ] := α q φ 0 , T ( α r φ k 0 ) L 2 ( 0 , T ) = t - 1 t t q T ( α r φ k 0 ) ( t ) d t

for k = 1 , , N t and = 1 , , N t , where α r ( t ) := t r and α q ( t ) := t q are monomials of degrees r 0 and q 0 . The following theorem states a new representation of the entries (2.20) with the help of the Legendre chi function χ ν : { z : | z | 1 } of order ν , ν 2 , given by the series

(2.21) χ ν ( z ) = η = 0 z 2 η + 1 ( 2 η + 1 ) ν , z  with  | z | 1 ,

see [2, 13] for more details. In this work, z and z are the real and imaginary part of a complex number z , and ι denotes the imaginary unit.

Theorem 2.1.

Let the polynomial degrees r N 0 and q N 0 be fixed. Then the matrix entries (2.20) of the auxiliary matrix A h t r , q R N t × N t are given by

A h t r , q [ k , ] = n = 0 r m = 0 q ( 2 T π ) n + m + 2 - 1 T ι ( m + n ) mod 2 ι n + m { - 1 , 1 } r ! ( r - n ) ! q ! ( q - m ) !
( ι ( m + n ) mod 2 { t k r - n t q - m [ ( - 1 ) n χ n + m + 2 ( e ι π ( t - t k ) 2 T ) + χ n + m + 2 ( e ι π ( t + t k ) 2 T ) ]
          - t k - 1 r - n t q - m [ ( - 1 ) n χ n + m + 2 ( e ι π ( t - t k - 1 ) 2 T ) + χ n + m + 2 ( e ι π ( t + t k - 1 ) 2 T ) ]
          - t k r - n t - 1 q - m [ ( - 1 ) n χ n + m + 2 ( e ι π ( t - 1 - t k ) 2 T ) + χ n + m + 2 ( e ι π ( t - 1 + t k ) 2 T ) ]
          + t k - 1 r - n t - 1 q - m [ ( - 1 ) n χ n + m + 2 ( e ι π ( t - 1 - t k - 1 ) 2 T ) + χ n + m + 2 ( e ι π ( t - 1 + t k - 1 ) 2 T ) ] } )

for k = 1 , , N t and = 1 , , N t .

Proof.

Let k , { 1 , , N t } be fixed. With the help of the representations (1.18), (1.19), it follows that

A h t r , q [ k , ] = 2 T η = 0 t k - 1 t k sin ( π T ( 1 2 + η ) s ) s r d s t - 1 t cos ( π T ( 1 2 + η ) t ) t q d t
= 2 4 T ι η = 0 t k - 1 t k ( e Θ ( η ) ι s - e - Θ ( η ) ι s ) s r d s t - 1 t ( e Θ ( η ) ι t + e - Θ ( η ) ι t ) t q d t

with Θ ( η ) := π T ( 1 2 + η ) . Using the integration rule

t r e c t d t = e c t ( n = 0 r ( - 1 ) n c - n - 1 r ! ( r - n ) ! t r - n )

with a constant 0 c gives

A h t r , q [ k , ] = - 1 2 T ι η = 0 n = 0 r m = 0 q ( - 1 ι ) n + m Θ ( η ) - n - m - 2 r ! ( r - n ) ! q ! ( q - m ) !
{ t k r - n e Θ ( η ) ι t k - t k - 1 r - n e Θ ( η ) ι t k - 1 - ( - 1 ) - n - 1 t k r - n e - Θ ( η ) ι t k + ( - 1 ) - n - 1 t k - 1 r - n e - Θ ( η ) ι t k - 1 }
{ t q - m e Θ ( η ) ι t - t - 1 q - m e Θ ( η ) ι t - 1 + ( - 1 ) - m - 1 t q - m e - Θ ( η ) ι t - ( - 1 ) - m - 1 t - 1 q - m e - Θ ( η ) ι t - 1 }

and so, we have

A h t r , q [ k , ] = - 1 2 T ι η = 0 n = 0 r m = 0 q ( - 1 ι ) n + m Θ ( η ) - n - m - 2 r ! ( r - n ) ! q ! ( q - m ) !
{ t k r - n [ e Θ ( η ) ι t k + ( - 1 ) - n e - Θ ( η ) ι t k ] - t k - 1 r - n [ e Θ ( η ) ι t k - 1 + ( - 1 ) - n e - Θ ( η ) ι t k - 1 ] }
{ t q - m [ e Θ ( η ) ι t - ( - 1 ) - m e - Θ ( η ) ι t - t - 1 q - m [ e Θ ( η ) ι t - 1 - ( - 1 ) - m e - Θ ( η ) ι t - 1 ] } .

Splitting n 0 and m 0 into odd and even indices yields four cases:

  1. n 0 is even and m 0 is even:

    - 1 2 T ι η = 0 1 ι n + m Θ ( η ) - n - m - 2 r ! ( r - n ) ! q ! ( q - m ) ! { t k r - n [ e Θ ( η ) ι t k + e - Θ ( η ) ι t k ] - t k - 1 r - n [ e Θ ( η ) ι t k - 1 + e - Θ ( η ) ι t k - 1 ] }
    { t q - m [ e Θ ( η ) ι t - e - Θ ( η ) ι t ] - t - 1 q - m [ e Θ ( η ) ι t - 1 - e - Θ ( η ) ι t - 1 ] }
    = - 2 T η = 0 1 ι n + m Θ ( η ) - n - m - 2 r ! ( r - n ) ! q ! ( q - m ) ! { t k r - n cos ( Θ ( η ) t k ) - t k - 1 r - n cos ( Θ ( η ) t k - 1 ) }
    { t q - m sin ( Θ ( η ) t ) - t - 1 q - m sin ( Θ ( η ) t - 1 ) }
    = - 1 T η = 0 1 ι n + m Θ ( η ) - n - m - 2 r ! ( r - n ) ! q ! ( q - m ) !
    { t k r - n t q - m [ sin ( Θ ( η ) ( t - t k ) ) + sin ( Θ ( η ) ( t + t k ) ) ]
    - t k - 1 r - n t q - m [ sin ( Θ ( η ) ( t - t k - 1 ) ) + sin ( Θ ( η ) ( t + t k - 1 ) ) ]
    - t k r - n t - 1 q - m [ sin ( Θ ( η ) ( t - 1 - t k ) ) + sin ( Θ ( η ) ( t - 1 + t k ) ) ]
    + t k - 1 r - n t - 1 q - m [ sin ( Θ ( η ) ( t - 1 - t k - 1 ) ) + sin ( Θ ( η ) ( t - 1 + t k - 1 ) ) ] }
    = ( 2 T π ) n + m + 2 - 1 T 1 ι n + m { - 1 , 1 } r ! ( r - n ) ! q ! ( q - m ) !
    { t k r - n t q - m [ χ n + m + 2 ( e ι π ( t - t k ) 2 T ) + χ n + m + 2 ( e ι π ( t + t k ) 2 T ) ]
    - t k - 1 r - n t q - m [ χ n + m + 2 ( e ι π ( t - t k - 1 ) 2 T ) + χ n + m + 2 ( e ι π ( t + t k - 1 ) 2 T ) ]
    - t k r - n t - 1 q - m [ χ n + m + 2 ( e ι π ( t - 1 - t k ) 2 T ) + χ n + m + 2 ( e ι π ( t - 1 + t k ) 2 T ) ]
    + t k - 1 r - n t - 1 q - m [ χ n + m + 2 ( e ι π ( t - 1 - t k - 1 ) 2 T ) + χ n + m + 2 ( e ι π ( t - 1 + t k - 1 ) 2 T ) ] } .

  2. n 0 is even and m 0 is odd:

    ( 2 T π ) n + m + 2 - 1 T 1 ι n + m - 1 { - 1 , 1 } r ! ( r - n ) ! q ! ( q - m ) !
    { t k r - n t q - m [ χ n + m + 2 ( e ι π ( t - t k ) 2 T ) + χ n + m + 2 ( e ι π ( t + t k ) 2 T ) ]
    - t k - 1 r - n t q - m [ χ n + m + 2 ( e ι π ( t - t k - 1 ) 2 T ) + χ n + m + 2 ( e ι π ( t + t k - 1 ) 2 T ) ]
    - t k r - n t - 1 q - m [ χ n + m + 2 ( e ι π ( t - 1 - t k ) 2 T ) + χ n + m + 2 ( e ι π ( t - 1 + t k ) 2 T ) ]
    + t k - 1 r - n t - 1 q - m [ χ n + m + 2 ( e ι π ( t - 1 - t k - 1 ) 2 T ) + χ n + m + 2 ( e ι π ( t - 1 + t k - 1 ) 2 T ) ] } .

  3. n 0 is odd and m 0 is even:

    ( 2 T π ) n + m + 2 1 T 1 ι n + m - 1 { - 1 , 1 } r ! ( r - n ) ! q ! ( q - m ) !
    { t k r - n t q - m [ χ n + m + 2 ( e ι π ( t - t k ) 2 T ) - χ n + m + 2 ( e ι π ( t + t k ) 2 T ) ]
    - t k - 1 r - n t q - m [ χ n + m + 2 ( e ι π ( t - t k - 1 ) 2 T ) - χ n + m + 2 ( e ι π ( t + t k - 1 ) 2 T ) ]
    - t k r - n t - 1 q - m [ χ n + m + 2 ( e ι π ( t - 1 - t k ) 2 T ) - χ n + m + 2 ( e ι π ( t - 1 + t k ) 2 T ) ]
    + t k - 1 r - n t - 1 q - m [ χ n + m + 2 ( e ι π ( t - 1 - t k - 1 ) 2 T ) - χ n + m + 2 ( e ι π ( t - 1 + t k - 1 ) 2 T ) ] } .

  4. n 0 is odd and m 0 is odd:

    ( 2 T π ) n + m + 2 1 T 1 ι n + m { - 1 , 1 } r ! ( r - n ) ! q ! ( q - m ) !
    { t k r - n t q - m [ χ n + m + 2 ( e ι π ( t - t k ) 2 T ) - χ n + m + 2 ( e ι π ( t + t k ) 2 T ) ]
    - t k - 1 r - n t q - m [ χ n + m + 2 ( e ι π ( t - t k - 1 ) 2 T ) - χ n + m + 2 ( e ι π ( t + t k - 1 ) 2 T ) ]
    - t k r - n t - 1 q - m [ χ n + m + 2 ( e ι π ( t - 1 - t k ) 2 T ) - χ n + m + 2 ( e ι π ( t - 1 + t k ) 2 T ) ]
    + t k - 1 r - n t - 1 q - m [ χ n + m + 2 ( e ι π ( t - 1 - t k - 1 ) 2 T ) - χ n + m + 2 ( e ι π ( t - 1 + t k - 1 ) 2 T ) ] } .

Finally, the assertion follows by using the modulo operator ( ) mod ( ) . ∎

For the matrices A h t T , C h t T in (2.15), (2.17), only the auxiliary matrices A h t 0 , 0 and A h t 1 , 0 are needed, which are given in the next corollaries. For the matrix M h t T in (2.16), analogous corollaries hold true for the auxiliary matrices A h t 0 , 1 and A h t 1 , 1 , which are omitted here.

Corollary 2.1.

Let the polynomial degrees be r = q = 0 . Then the matrix entries (2.20) of the auxiliary matrix A h t 0 , 0 R N t × N t are given by

A h t 0 , 0 [ k , ] = ( 2 T π ) 2 - 1 T ( χ 2 ( e ι π ( t - t k ) 2 T ) + χ 2 ( e ι π ( t + t k ) 2 T ) - χ 2 ( e ι π ( t - t k - 1 ) 2 T ) - χ 2 ( e ι π ( t + t k - 1 ) 2 T ) - χ 2 ( e ι π ( t - 1 - t k ) 2 T )
- χ 2 ( e ι π ( t - 1 + t k ) 2 T ) + χ 2 ( e ι π ( t - 1 - t k - 1 ) 2 T ) + χ 2 ( e ι π ( t - 1 + t k - 1 ) 2 T ) )

for k = 1 , , N t and = 1 , , N t .

Corollary 2.2.

Let the polynomial degrees be r = 1 and q = 0 . Then the matrix entries (2.20) of the auxiliary matrix A h t 1 , 0 R N t × N t are given by

A h t 1 , 0 [ k , ] = ( 2 T π ) 2 - 1 T ( t k - 1 [ χ 2 ( e ι π ( t - 1 - t k - 1 ) 2 T ) + χ 2 ( e ι π ( t - 1 + t k - 1 ) 2 T ) - χ 2 ( e ι π ( t - t k - 1 ) 2 T ) - χ 2 ( e ι π ( t + t k - 1 ) 2 T ) ]
+ t k [ χ 2 ( e ι π ( t - t k ) 2 T ) + χ 2 ( e ι π ( t + t k ) 2 T ) - χ 2 ( e ι π ( t - 1 - t k ) 2 T ) - χ 2 ( e ι π ( t - 1 + t k ) 2 T ) ] )
+ ( 2 T π ) 3 - 1 T ( - χ 3 ( e ι π ( t - t k ) 2 T ) + χ 3 ( e ι π ( t + t k ) 2 T ) + χ 3 ( e ι π ( t - t k - 1 ) 2 T ) - χ 3 ( e ι π ( t + t k - 1 ) 2 T )
+ χ 3 ( e ι π ( t - 1 - t k ) 2 T ) - χ 3 ( e ι π ( t - 1 + t k ) 2 T ) - χ 3 ( e ι π ( t - 1 - t k - 1 ) 2 T ) + χ 3 ( e ι π ( t - 1 + t k - 1 ) 2 T ) )

for k = 1 , , N t and = 1 , , N t .

Remark 2.1.

To calculate the entries for several auxiliary matrices, e.g., A h t 0 , 0 , A h t 1 , 0 , A h t 0 , 1 and A h t 1 , 1 as needed in (2.15), (2.16), (2.17), it is possible to compute the quantities

(2.22) χ ν ( e ι β ) or χ ν ( e ι β ) , ν 0 , ν 2 , β [ - π , π ] ,

only once and to reuse them. Another possibility is to derive representations of the entries of the matrices A h t T , M h t T , C h t T in (1.17) and (1.22) as linear combinations of quantities (2.22), i.e. without using the assembling matrices (2.18), (2.19). For example, the entries of the matrix A h t T admit a representation of the form

A h t T [ , k ] = j [ c , k , j χ 2 ( e ι β , k , j ) + c ^ , k , j χ 3 ( e ι β ^ , k , j ) ]

for , k = 1 , , N t with coefficients c , k , j , c ^ , k , j and values β , k , j , β ^ , k , j , which depend on the nodes of the time mesh (1.3).

Remark 2.2.

In this section, we consider only piecewise linear basis functions. Since Theorem 2.1 holds true for arbitrary polynomial degrees r , q 0 , a realization of the discretization (1.11) for any p or high-order splines is straightforward.

2.2.1 Evaluation of the Legendre Chi Function

In this subsection, the advantage of the representation of the matrix entries (2.20) via the Legendre chi function χ ν in (2.21) is stated, where the order ν 0 , ν 2 , is fixed. Here, we apply the approach in [2] of using a special expansion for χ ν to numerically evaluate

χ ν ( e ι β ) for even  ν

and

χ ν ( e ι β ) for odd  ν

for β [ - π , π ] . Since the symmetry properties

χ ν ( e ι β ) = χ ν ( e ι ( π - β ) ) and χ ν ( e ι β ) = - χ ν ( e ι ( π - β ) ) , β [ - π , π ] ,

hold true, see [2, 13], it is sufficient to calculate χ ν ( e ι β ) and χ ν ( e ι β ) for β [ - π 2 , π 2 ] only. Using the expansion given in [2, (2.5), p. 159], we get

(2.23) χ ν ( e ι β ) for  ν  is even, χ ν ( e ι β ) for  ν  is odd } = j = 0 ν - 2 ( 1 - 2 - ν + j ) ζ ( ν - j ) β j j ! ( ι j + ( ν mod 2 ) ) { - 1 , 0 , 1 } + ι ν ι ν mod 2 { - 1 , 1 } { η = 0 ( 1 - 2 - 2 η - 1 ) ζ ( 2 η + 2 ) ( 2 η + 2 ) ( 2 η + 3 ) ( 2 η + ν + 1 ) β ν + 1 π 2 ( β π ) 2 η + β ν - 1 2 ( ν - 1 ) ! [ Γ ( 1 2 ) Γ ( 1 2 ) - Γ ( ν ) Γ ( ν ) + ln ( 2 | β | ) ] }

for β [ - π 2 , π 2 ] with the Gamma function Γ ( ) and the Riemann zeta function ζ ( ) , where for β = 0 , the term β ν - 1 ln ( 2 | β | ) disappears. In the worst case, i.e. β = ± π 2 , the expansion (2.23) converges like a geometric series of ratio β = ± 1 4 , see [2, (2.7), p. 160]. Hence, truncating the expansion (2.23) at ρ χ 0 , i.e.

(2.24) χ ν ( e ι β ) for  ν  is even, χ ν ( e ι β ) for  ν  is odd } j = 0 ν - 2 ( 1 - 2 - ν + j ) ζ ( ν - j ) β j j ! ( ι j + ( ν mod 2 ) ) { - 1 , 0 , 1 } + ι ν ι ν mod 2 { - 1 , 1 } { η = 0 ρ χ ( 1 - 2 - 2 η - 1 ) ζ ( 2 η + 2 ) ( 2 η + 2 ) ( 2 η + 3 ) ( 2 η + ν + 1 ) β ν + 1 π 2 ( β π ) 2 η + β ν - 1 2 ( ν - 1 ) ! [ Γ ( 1 2 ) Γ ( 1 2 ) - Γ ( ν ) Γ ( ν ) + ln ( 2 | β | ) ] } ,

gives an approximation of (2.23) for any β [ - π 2 , π 2 ] , where the accuracy is controlled by the truncation parameter ρ χ 0 . Note that the quantities Γ ( ν ) Γ ( ν ) , ζ ( ν - j ) and ζ ( 2 η + 2 ) in (2.24) can be precomputed and reused.

3 Numerical Examples

In this section, we give numerical examples for the assembling of the matrices A h t T and M h t T in (1.17), where the procedures of Section 2 are applied. In addition, a numerical example for the heat equation in a spatially two-dimensional domain is given.

3.1 Evaluation of the Legendre Chi Function

In this subsection, we investigate the errors of evaluating the Legendre chi function χ ν as proposed in Section 2.2.1. In Table 1 the truncation errors for the approximation (2.24) of the expansion (2.23) are given for the worst case β = π 2 for the orders ν = 2 , 3 , 4 , 5 , 6 , where we observe that the truncation errors converge very fast to 0. In addition, the truncation error decreases when the order ν increases. Hence, using the truncated series (2.24) for calculating the matrix entries (2.20) of the auxiliary matrix A h t r , q , the truncation parameter ρ χ 0 can be chosen rather small and independently of the time mesh (1.3).

Table 1

Truncation errors for the approximation (2.24) of the expansion (2.23) for β = π 2 .

ρ χ ν = 2 ν = 3 ν = 4 ν = 5 ν = 6
0 5.33e - 03 1.35e - 03 2.95e - 04 5.70e - 05 9.83e - 06
2 1.02e - 04 1.55e - 05 2.16e - 06 2.77e - 07 3.30e - 08
4 3.02e - 06 3.30e - 07 3.38e - 08 3.26e - 09 2.96e - 10
6 1.10e - 07 9.38e - 09 7.60e - 10 5.86e - 11 4.32e - 12
8 4.51e - 09 3.15e - 10 2.11e - 11 1.36e - 12 8.42e - 14
10 1.99e - 10 1.18e - 11 6.73e - 13 3.72e - 14 1.99e - 15
12 9.22e - 12 4.75e - 13 2.37e - 14 1.15e - 15 5.39e - 17
14 4.45e - 13 2.03e - 14 8.96e - 16 3.86e - 17 1.62e - 18
16 2.21e - 14 9.02e - 16 3.59e - 17 1.39e - 18 5.27e - 20
18 1.13e - 15 4.16e - 17 1.50e - 18 5.30e - 20 1.83e - 21
20 5.85e - 17 1.97e - 18 6.52e - 20 2.11e - 21 6.70e - 23

3.2 Approximations A ~ h t T A h t T and M ~ h t T M h t T

In this subsection, we investigate numerical examples, regarding the quality of the approximations A ~ h t T and M ~ h t T in (2.6) of the matrices A h t T , M h t T in (1.17) via truncating the series expansion of the definition (1.18) of the modified Hilbert transformation T . For this purpose, the temporal domain ( 0 , 1 2 ) = ( 0 , T ) is decomposed into nonuniform elements with the nodes

(3.1) t 0 = 0.0 , t 1 = 1 32 , t 2 = 1 16 , t 3 = 1 8 , t 4 = 1 2 = T ,

or we use a graded mesh with nodes

(3.2) t = ( N t ) q T , = 0 , , N t ,

with a grading parameter q 1 . First, in Table 2 we show the relative error of the matrix entries A h t T , M h t T when approximating them by (2.6) for increasing truncation parameter ρ and the fixed time mesh (3.1).

Table 2

Relative errors for the approximate matrices A ~ h t T , M ~ h t T in (2.6) with increasing truncation parameter ρ andfixed time mesh (3.1) with N t = 4 .

ρ max k , | A h t T [ k , ] - A ~ h t T [ k , ] | max k , | A h t T [ k , ] | eoc max k , | M h t T [ k , ] - M ~ h t T [ k , ] | max k , | M h t T [ k , ] | eoc
16 7.7e - 02 4.0e - 03
32 1.6e - 02 2.3 1.7e - 04 4.6
64 4.4e - 03 1.9 4.2e - 05 2.0
128 1.1e - 03 2.0 1.1e - 05 2.0
256 2.8e - 04 2.0 2.7e - 06 2.0
512 7.0e - 05 2.0 6.7e - 07 2.0
1024 1.8e - 05 2.0 1.7e - 07 2.0
Table 3

Errors for the approximate matrix A ~ h t T in (2.6) with mesh refinements and a fixed truncation parameter ρ = 10 5 . In the second and third column, uniform refinements are applied to the mesh (3.1). In the fourth and fifth column, we use the graded mesh (3.2) with grading parameter q = 1.5 .

N t max k , | A h t T [ k , ] - A ~ h t T [ k , ] | eoc max k , | A h t T [ k , ] - A ~ h t T [ k , ] | eoc
4 2.1e - 09 2.8e - 10
8 9.9e - 09 - 2.3 2.2e - 09 - 3.0
16 4.0e - 08 - 2.0 1.8e - 08 - 3.0
32 1.6e - 07 - 2.0 1.4e - 07 - 3.0
64 6.4e - 07 - 2.0 1.1e - 06 - 3.0
128 2.5e - 06 - 2.0 9.1e - 06 - 3.0
256 1.0e - 05 - 2.0 7.4e - 05 - 3.0

Second, in Table 3 we give the errors for the approximate matrix A ~ h t T in (2.6) when the number of elements N t doubles, i.e. refining the mesh, and the truncation parameter ρ = 10 5 is fixed. For the second and third column in Table 3, we refine the time mesh (3.1) uniformly. For the fourth and fifth column in Table 3, we use the graded mesh (3.2) with the grading parameter q = 1.5 . In all cases of Table 2 and Table 3, we use the procedure of Section 2.2, i.e. (2.24) with ρ χ = 20 , to calculate the reference values of the matrix entries of A h t T and M h t T . The behavior of the errors in Table 2 and Table 3 is as expected from the error estimates (2.7), (2.8), (2.9) and (2.10), (2.11), (2.12), (2.13). Note that we have, e.g., h t , 1 , h t , 2 N t - q for the graded mesh (3.2), which explains the results of the fourth and fifth column in Table 3.

Last, in Table 4 the smallest eigenvalue λ min of the approximate matrix A ~ h t T in (2.6) is given with the truncation parameter ρ for the time mesh (3.1) with uniform refinements. In Table 4 we see that the choice of the truncation parameter ρ, such that the approximate matrix A ~ h t T is positive definite, depends linearly on the mesh size h t , min , i.e. estimate (2.14) is not sharp. However, an accurate approximation of the matrix A h t T is needed, which is possible independently of the mesh size h t , min by using the procedure proposed in Section 2.2.

Table 4

Smallest eigenvalue λ min of the approximate matrix A ~ h t T in (2.6) with truncation parameter ρ for the time mesh (3.1) with uniform refinements.

ρ N t = 4 N t = 8 N t = 16 N t = 32 N t = 64
1 - 3.1e - 17 - 2.2e - 17 - 4.9e - 17 - 2.6e - 17 - 6.2e - 18
2 - 1.5e - 17 - 8.9e - 17 - 6.0e - 17 - 4.5e - 17 - 1.6e - 17
4 3.0e - 03 - 1.4e - 17 - 5.4e - 17 - 6.8e - 17 - 7.4e - 17
8 3.6e - 01 5.5e - 11 - 2.4e - 17 - 1.0e - 16 - 1.4e - 16
16 4.8e - 01 3.5e - 04 - 3.0e - 17 - 6.2e - 17 - 2.2e - 16
32 4.8e - 01 3.1e - 01 4.5e - 12 - 1.3e - 16 - 1.8e - 16
64 4.8e - 01 3.1e - 01 1.7e - 01 - 3.9e - 17 - 3.4e - 16
128 4.8e - 01 3.1e - 01 1.7e - 01 8.4e - 02 - 3.6e - 16
256 4.8e - 01 3.1e - 01 1.7e - 01 8.4e - 02 4.2e - 02
512 4.8e - 01 3.1e - 01 1.7e - 01 8.4e - 02 4.2e - 02
1024 4.8e - 01 3.1e - 01 1.7e - 01 8.4e - 02 4.2e - 02

3.3 Heat Equation

In this subsection, a numerical example for the Galerkin–Bubnov finite element method (1.15) is given. For this purpose, the two-dimensional spatial L-shaped domain

(3.3) Ω := ( - 1 , 1 ) 2 ( [ 0 , 1 ] × [ - 1 , 0 ] ) 2

and the terminal time T = 1 2 are considered for the solution

u ( x 1 , x 2 , t ) = 5 2 π t e - ( x 1 - 1 4 ) 2 - ( x 2 + 1 4 ) 2 4 t sin ( π x 1 x 2 ) , ( x 1 , x 2 , t ) Q = Ω × ( 0 , T ) .

The inhomogeneous Dirichlet boundary condition is treated via homogenization as for the elliptic case, see [20, p. 246]. In greater detail, as there exists an extension u g H ;  0 , 1 , 1 / 2 ( Q ) := H 0 , 1 / 2 ( 0 , T ; L 2 ( Ω ) ) L 2 ( 0 , T ; H 1 ( Ω ) ) with u g | Σ = g for any Dirichlet data g H 1 / 2 , 1 / 4 ( Σ ) , see [3], we introduce the splitting u 0 = u - u g H 0 ;  0 , 1 , 1 / 2 ( Q ) . Hence, we consider the variational formulation to find u 0 H 0 ;  0 , 1 , 1 / 2 ( Q ) such that

a ( u 0 , T v ) = f , T v Q - a ( u g , T v ) for all  v H 0 ;  0 , 1 , 1 / 2 ( Q ) .

Assuming that u g C ( Q ¯ ) , we can replace u g with its nodal interpolant I h u g . Thus, with the help of the Galerkin–Bubnov finite element method (1.15), we approximate u 0 u ~ 0 , h Q h 1 ( Q ) H 0 ;  0 , 1 , 1 / 2 ( Q ) . This leads to an approximation u u ~ h := u ~ 0 , h + I h u g . Note that the extension u g is not needed explicitly for the implementation, i.e. it is sufficient to know the Dirichlet data g on the lateral boundary Σ to calculate u ~ h , see the splitting in [20, (11.8), p. 246].

The spatial domain Ω is decomposed into uniform triangles with uniform mesh size h x as given in Figure 1 for level 0. The temporal domain ( 0 , 1 2 ) = ( 0 , T ) is decomposed into nonuniform elements with the nodes (3.1). The assembling of the global linear system (1.16) is done via the representations (2.15), (2.16), (2.17) of A h t T , M h t T , C h t T , where the auxiliary matrices A h t 0 , 0 , A h t 1 , 0 , A h t 0 , 1 , A h t 1 , 1 with entries (2.20) are calculated by the truncated series (2.24) for the truncation parameter ρ χ = 10 . The integrals to compute the projection Q h 0 f in (1.14) are calculated by using high-order quadrature rules. The global linear system (1.16) is solved by a direct solver. The numerical results for the smooth solution u, when a uniform refinement strategy is applied as in Figure 1, are given in Table 5, where unconditional stability is observed and the convergence rates in L 2 ( Q ) and | | H 1 ( Q ) are as expected from the error estimates (1.12) and (1.13).

Figure 1 
                  Uniform refinement strategy: Starting mesh, the meshes after one and two uniform refinement steps.
Figure 1 
                  Uniform refinement strategy: Starting mesh, the meshes after one and two uniform refinement steps.
Figure 1 
                  Uniform refinement strategy: Starting mesh, the meshes after one and two uniform refinement steps.
Figure 1

Uniform refinement strategy: Starting mesh, the meshes after one and two uniform refinement steps.

Table 5

Numerical results of the Galerkin–Bubnov finite element discretization (1.15) for the L-shape (3.3) and T = 1 2 for the function u for a uniform refinement strategy with truncation parameter ρ χ = 10 .

dof h x h t , max h t , min u - u ~ h L 2 ( Q ) eoc | u - u ~ h | H 1 ( Q ) eoc
20 0.3536 0.3750 0.0313 3.33e - 01 4.31e + 00
264 0.1768 0.1875 0.0156 1.09e - 01 1.3 2.70e + 00 0.5
2576 0.0884 0.0938 0.0078 3.14e - 02 1.6 1.44e + 00 0.8
22560 0.0442 0.0469 0.0039 8.31e - 03 1.8 6.98e - 01 1.0
188480 0.0221 0.0234 0.0020 2.13e - 03 1.9 3.45e - 01 1.0
1540224 0.0111 0.0117 0.0010 5.38e - 04 2.0 1.71e - 01 1.0

4 Conclusion

First, a truncated series, coming from the definition of the modified Hilbert transformation T , was used to approximate the matrix entries of the matrices A h t T and M h t T , where the choice of the truncation parameter depends on the mesh size h t , min . Second, a new possibility of realizing the modified Hilbert transformation T based on the Legendre chi function χ ν was introduced. The main advantage of this new procedure is that the matrix entries of the matrices A h t T and M h t T , which are needed for a Galerkin–Bubnov discretization of parabolic equations, can be calculated to machine precision independently of the mesh size h t , min . Moreover, since the main theorem of this new series representation is formulated for any polynomial degree, a generalization to high-order splines or piecewise polynomial functions of an arbitrary degree is straightforward.

References

[1] R. Andreev, Stability of sparse space-time finite element discretizations of linear parabolic evolution equations, IMA J. Numer. Anal. 33 (2013), no. 1, 242–260. 10.1093/imanum/drs014Search in Google Scholar

[2] J. Boersma and J. P. Dempsey, On the numerical evaluation of Legendre’s chi-function, Math. Comp. 59 (1992), no. 199, 157–163. 10.1090/S0025-5718-1992-1134715-0Search in Google Scholar

[3] M. Costabel, Boundary integral operators for the heat equation, Integral Equations Operator Theory 13 (1990), no. 4, 498–552. 10.1007/BF01210400Search in Google Scholar

[4] D. Devaud, Petrov–Galerkin space-time hp-approximation of parabolic equations in H 1 / 2 , IMA J. Numer. Anal. 40 (2020), no. 4, 2717–2745. 10.1093/imanum/drz036Search in Google Scholar

[5] D. Devaud and C. Schwab, Space-time hp-approximation of parabolic equations, Calcolo 55 (2018), no. 3, Paper No. 35. 10.1007/s10092-018-0275-2Search in Google Scholar

[6] M. Fontes, Parabolic Equations with Low Regularity, ProQuest LLC, Ann Arbor, 1996; Thesis (Takn.dr)–Lunds Universitet (Sweden). Search in Google Scholar

[7] M. Fontes, Initial-boundary value problems for parabolic equations, Ann. Acad. Sci. Fenn. Math. 34 (2009), no. 2, 583–605. Search in Google Scholar

[8] T. Führer and M. Karkulik, Space-time least-squares finite elements for parabolic equations, preprint (2019), https://arxiv.org/abs/1911.01942. 10.1016/j.camwa.2021.03.004Search in Google Scholar

[9] O. A. Ladyženskaja, V. A. Solonnikov and N. N. Ural’ceva, Linear and Quasilinear Equations of Parabolic Type, Transl. Math. Monogr. 23, American Mathematical Society, Providence, 1968. 10.1090/mmono/023Search in Google Scholar

[10] U. Langer, S. E. Moore and M. Neumüller, Space-time isogeometric analysis of parabolic evolution problems, Comput. Methods Appl. Mech. Engrg. 306 (2016), 342–363. 10.1016/j.cma.2016.03.042Search in Google Scholar

[11] S. Larsson and M. Molteni, Numerical solution of parabolic problems based on a weak space-time formulation, Comput. Methods Appl. Math. 17 (2017), no. 1, 65–84. 10.1515/cmam-2016-0027Search in Google Scholar

[12] S. Larsson and C. Schwab, Compressive space-time Galerkin discretizations of parabolic partial differential equations, preprint (2015), https://arxiv.org/abs/1501.04514. Search in Google Scholar

[13] L. Lewin, Polylogarithms and Associated Functions, North-Holland, New York, 1981. Search in Google Scholar

[14] J.-L. Lions and E. Magenes, Problèmes aux limites non homogènes et applications. Vol. 1, Trav. Rech. Math. 17, Dunod, Paris, 1968. Search in Google Scholar

[15] J.-L. Lions and E. Magenes, Problèmes aux limites non homogènes et applications. Vol. 2, Trav. Rech. Math. 18, Dunod, Paris, 1968. Search in Google Scholar

[16] C. Mollet, Stability of Petrov–Galerkin discretizations: Application to the space-time weak formulation for parabolic evolution problems, Comput. Methods Appl. Math. 14 (2014), no. 2, 231–255. 10.1515/cmam-2014-0001Search in Google Scholar

[17] M. Neumüller, Space-Time Methods: Fast Solvers and Applications, Monogr. Ser. TU Graz Comput. Eng. Sci. 20, TU Graz, Graz, 2013. Search in Google Scholar

[18] C. Schwab and R. Stevenson, Space-time adaptive wavelet methods for parabolic evolution problems, Math. Comp. 78 (2009), no. 267, 1293–1318. 10.1090/S0025-5718-08-02205-9Search in Google Scholar

[19] C. Schwab and R. Stevenson, Fractional space-time variational formulations of Navier–Stokes equations, SIAM J. Math. Anal. 49 (2017), no. 4, 2442–2467. 10.1137/15M1051725Search in Google Scholar

[20] O. Steinbach, Numerical Approximation Methods for Elliptic Boundary Value Problems. Finite and Boundary Elements, Springer, New York, 2008. 10.1007/978-0-387-68805-3Search in Google Scholar

[21] O. Steinbach, Space-time finite element methods for parabolic problems, Comput. Methods Appl. Math. 15 (2015), no. 4, 551–566. 10.1515/cmam-2015-0026Search in Google Scholar

[22] O. Steinbach and M. Zank, A note on the efficient evaluation of a modified Hilbert transformation, J. Numer. Math. (2020), 10.1515/jnma-2019-0099. 10.1515/jnma-2019-0099Search in Google Scholar

[23] O. Steinbach and M. Zank, Coercive space-time finite element methods for initial boundary value problems, Electron. Trans. Numer. Anal. 52 (2020), 154–194. 10.1553/etna_vol52s154Search in Google Scholar

[24] R. Stevenson and J. Westerdiep, Stability of Galerkin discretizations of a mixed space-time variational formulation of parabolic evolution equations, IMA J. Numer. Anal. (2020), 10.1093/imanum/drz069. 10.1093/imanum/drz069Search in Google Scholar

[25] M. Zank, Inf-Sup Stable Space-Time Methods for Time-Dependent Partial Differential Equations, Monogr. Ser. TU Graz Comput. Eng. Sci. 36, TU Graz, Graz, 2020. Search in Google Scholar

Received: 2020-03-01
Revised: 2020-09-23
Accepted: 2020-10-23
Published Online: 2020-11-11
Published in Print: 2021-04-01

© 2020 Walter de Gruyter GmbH, Berlin/Boston

This work is licensed under the Creative Commons Attribution 4.0 International License.

Downloaded on 28.3.2024 from https://www.degruyter.com/document/doi/10.1515/cmam-2020-0026/html
Scroll to top button