Marco Zank

# 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:

• 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 ) ] } .

• n 0 is even and m 0 is odd:

( 2 T π ) n + m + 2 - 1 T 1 ι n + m - 1 { - 1 , 1 } <