Integral Representations and Quadrature Schemes for the Modified Hilbert Transformation

Marco Zank

doi:10.1515/cmam-2022-0150

Artikel Open Access

Integral Representations and Quadrature Schemes for the Modified Hilbert Transformation

Marco Zank

Veröffentlicht/Copyright: 6. Oktober 2022

Veröffentlicht von

Veröffentlichen auch Sie bei De Gruyter Brill

Manuskript einreichen Informationen für Autor*innen Erkunden Sie dieses Fachgebiet

Aus der Zeitschrift Computational Methods in Applied Mathematics Band 23 Heft 2

Abstract

We present quadrature schemes to calculate matrices where the so-called modified Hilbert transformation is involved. These matrices occur as temporal parts of Galerkin finite element discretizations of parabolic or hyperbolic problems when the modified Hilbert transformation is used for the variational setting. This work provides the calculation of these matrices to machine precision for arbitrary polynomial degrees and non-uniform meshes. The proposed quadrature schemes are based on weakly singular integral representations of the modified Hilbert transformation. First, these weakly singular integral representations of the modified Hilbert transformation are proven. Second, using these integral representations, we derive quadrature schemes, which treat the occurring singularities appropriately. Thus, exponential convergence with respect to the number of quadrature nodes for the proposed quadrature schemes is achieved. Numerical results, where this exponential convergence is observed, conclude this work.

Keywords: Modified Hilbert Transformation; Weakly Singular Integral; Quadratures; Space-Time; Finite Element Method; Parabolic Equations; Hyperbolic Equations

MSC 2010: 65M60; 65R10

1 Introduction

For the discretization of parabolic or hyperbolic evolution equations, the classical approaches are time stepping schemes together with finite element methods in space, see, e.g., [17] for parabolic and [1] for hyperbolic problems. An alternative is to discretize the time-dependent problem without separating the temporal and spatial variables, i.e., space-time methods, see, e.g., [6]. 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. Different variational approaches for parabolic and hyperbolic problems are contained in [5, 8]. In recent years, novel variational settings, which use a so-called modified Hilbert transformation H T , have been introduced, see [15, 18] for parabolic problems, and see [9, 10] for hyperbolic problems. Conforming space-time discretizations by piecewise polynomials of these variational formulations lead to huge linear systems. In [7, 21, 20], fast space-time solvers are developed, where some of them allow for an easy parallelization in time. Further extensions for this variational setting, which includes the modified Hilbert transformation H T , are an h ⁢ p -FEM in temporal direction and a classical FEM with graded meshes in spatial direction, see [12]. All the mentioned discretization schemes ask for an accurate realization of the modified Hilbert transformation H T , which acts only on the temporal part, to calculate the corresponding matrices. There are different possibilities of computing the entries of these temporal matrices. In [16], a weakly singular integral representation of the modified Hilbert transformation H T for sufficiently smooth functions and quadrature schemes are derived, which can be used for discretizations of parabolic problems as in [15, 18]. However, in [16], singular integrals are proposed to be calculated analytically, which is not convenient for high-order approaches, as the aforementioned h ⁢ p -FEM. In addition, for discretizations of hyperbolic problems as in [9, 10], the weakly singular integral representation of the modified Hilbert transformation H T , given in [16], is not applicable. One way out is the approach in [19], which is well-suited for low polynomial degrees. Thus, accurate realizations of the modified Hilbert transformation H T for high polynomial degrees, as needed in an h ⁢ p -FEM, seem not to be available. In this work, we provide such realizations. For this purpose, we prove a new weakly singular integral representation of the modified Hilbert transformation H T for sufficiently smooth functions with jumps and derive quadrature schemes, which allow for the calculation of the involved matrices to high float point accuracy.

In greater detail, for given T > 0 , μ ≥ 0 and 𝑓, we consider the ordinary differential equation (ODE)

(1.1) ∂ t ⁡ u ⁢ ( t ) + μ ⁢ u ⁢ ( t ) = f ⁢ ( t ) for ⁢ t ∈ ( 0 , T ) , u ⁢ ( 0 ) = 0 ,

which is the temporal part of the heat equation, and the ordinary differential equation

(1.2) ∂ t ⁢ t ⁡ u ⁢ ( t ) + μ ⁢ u ⁢ ( t ) = f ⁢ ( t ) for ⁢ t ∈ ( 0 , T ) , u ⁢ ( 0 ) = ∂ t ⁡ u ⁢ ( 0 ) = 0 ,

which is the temporal part of the wave equation. Note that 𝜇 plays the role of an eigenvalue of the spatial operator, see [5, 15, 18].

To state related variational formulations to (1.1), (1.2), we introduce classical Sobolev spaces, where details and further references are contained in [18]. The classical Sobolev space H 1 ⁢ ( 0 , T ) is endowed with the norm ∥ ⋅ ∥ H 1 ⁢ ( 0 , T ) . Its closed subspaces

H 0 , 1 ⁢ ( 0 , T ) = { v ∈ H 1 ⁢ ( 0 , T ) : v ⁢ ( 0 ) = 0 } and H , 0 1 ⁢ ( 0 , T ) = { v ∈ H 1 ⁢ ( 0 , T ) : v ⁢ ( T ) = 0 }

are endowed with the Hilbertian norm | ⋅ | H 1 ⁢ ( 0 , T ) := ∥ ∂ t ⋅ ∥ L 2 ⁢ ( 0 , T ) , where ∥ ⋅ ∥ L 2 ⁢ ( 0 , T ) is the usual L 2 ⁢ ( 0 , T ) norm and ∂ t denotes the (weak) derivative. Additionally, we have the Sobolev spaces

H 0 , 1 / 2 ⁢ ( 0 , T ) := { U | ( 0 , T ) : U ∈ H 1 / 2 ⁢ ( - ∞ , T ) ⁢ with ⁢ U ⁢ ( t ) = 0 , t < 0 } , H , 0 1 / 2 ⁢ ( 0 , T ) := { U | ( 0 , T ) : U ∈ H 1 / 2 ⁢ ( 0 , ∞ ) ⁢ with ⁢ U ⁢ ( t ) = 0 , t > 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 ) 1 / 2 , ∥ v ∥ H , 0 1 / 2 ⁢ ( 0 , T ) := ( ∥ v ∥ H 1 / 2 ⁢ ( 0 , T ) 2 + ∫ 0 T | v ⁢ ( t ) | 2 T - t ⁢ d t ) 1 / 2 ,

where ∥ ⋅ ∥ H 1 / 2 ⁢ ( J ) is a norm, e.g., the Slobodetskii norm, in the usual Sobolev space H 1 / 2 ⁢ ( J ) for open intervals J ⊂ R . Further, ⟨ ⋅ , ⋅ ⟩ ( 0 , T ) denotes the duality pairing in [ H , 0 1 / 2 ⁢ ( 0 , T ) ] ′ and H , 0 1 / 2 ⁢ ( 0 , T ) , or in [ H , 0 1 ⁢ ( 0 , T ) ] ′ and H , 0 1 ⁢ ( 0 , T ) , as extension of the inner product ⟨ ⋅ , ⋅ ⟩ L 2 ⁢ ( 0 , T ) in L 2 ⁢ ( 0 , T ) . Last, we define the modified Hilbert transformation H T : L 2 ⁢ ( 0 , T ) → L 2 ⁢ ( 0 , T ) as

(1.3) ( H T ⁢ v ) ⁢ ( t ) = ∑ k = 0 ∞ v k ⁢ cos ⁡ ( ( π 2 + k ⁢ π ) ⁢ t T ) , t ∈ ( 0 , T ) ,

where

v k = 2 T ⁢ ∫ 0 T v ⁢ ( t ) ⁢ sin ⁡ ( ( π 2 + k ⁢ π ) ⁢ t T ) ⁢ d t

are the Fourier coefficients of the series representation

v ⁢ ( t ) = ∑ k = 0 ∞ v k ⁢ sin ⁡ ( ( π 2 + k ⁢ π ) ⁢ t T ) , t ∈ ( 0 , T ) ,

when v ∈ L 2 ⁢ ( 0 , T ) is given. The modified Hilbert transformation H T is a bijective mapping

H T : H 0 , ν ⁢ ( 0 , T ) → H , 0 ν ⁢ ( 0 , T ) for ⁢ ν ∈ { 0 , 1 / 2 , 1 } with ⁢ H 0 , 0 ⁢ ( 0 , T ) = H , 0 0 ⁢ ( 0 , T ) = L 2 ⁢ ( 0 , T ) .

Additional properties of H T are given in [15, 16, 19, 14, 18].

With this notation, a related variational formulation to the parabolic-type problem (1.1) is to find u ∈ H 0 , 1 / 2 ⁢ ( 0 , T ) such that

(1.4) ⟨ ∂ t ⁡ u , H T ⁢ v ⟩ ( 0 , T ) + μ ⁢ ⟨ u , H T ⁢ v ⟩ L 2 ⁢ ( 0 , T ) = ⟨ f , H T ⁢ v ⟩ ( 0 , T ) for all ⁢ v ∈ H 0 , 1 / 2 ⁢ ( 0 , T ) ,

for a given right-hand side f ∈ [ H , 0 1 / 2 ⁢ ( 0 , T ) ] ′ . For the hyperbolic-type problem (1.2), a related variational formulation is to find u ∈ H 0 , 1 ⁢ ( 0 , T ) such that

(1.5) ⟨ H T ⁢ ∂ t ⁡ u , ∂ t ⁡ v ⟩ L 2 ⁢ ( 0 , T ) + μ ⁢ ⟨ u , H T ⁢ v ⟩ L 2 ⁢ ( 0 , T ) = ⟨ f , H T ⁢ v ⟩ ( 0 , T ) for all ⁢ v ∈ H 0 , 1 ⁢ ( 0 , T ) ,

for a given right-hand side f ∈ [ H , 0 1 ⁢ ( 0 , T ) ] ′ . The variational formulations (1.4), (1.5) are uniquely solvable, where their derivation and analysis are given in [15, 18, 9].

Next, we state conforming discretizations of the variational formulations (1.4), (1.5). For this purpose, we define a partition T N = { τ ℓ } ℓ = 1 N of ( 0 , T ) for given N ∈ N , i.e.,

(1.6) 0 = t 0 < t 1 < t 2 < ⋯ < t N - 1 < t N = T

with elements τ ℓ := ( t ℓ - 1 , t ℓ ) ⊂ R . The mesh sizes are h ℓ := t ℓ - t ℓ - 1 for ℓ = 1 , … , N , and the maximal mesh size is h := max ℓ ⁡ h ℓ . On T N , with the distribution p = ( p 1 , … , p N ) ∈ N N of polynomial degrees, we introduce the space of piecewise polynomial, continuous functions on intervals

(1.7) S p ( T N ) := { v ∈ C [ 0 , T ] : v | τ ℓ ∈ P p ℓ ( τ ℓ ) for all ℓ ∈ { 1 , … , N } } = span { φ i } i = 1 M ,

where C ⁢ [ 0 , T ] is the space of continuous functions and P p ⁢ ( A ) is the space of polynomials on a subset A ⊂ R d , d ∈ N , of global degree at most p ∈ N 0 . Here, M = ∑ ℓ = 1 N ( p ℓ + 1 ) - ( N - 1 ) = 1 + ∑ ℓ = 1 N p ℓ is the number of degrees of freedom and the functions φ i : [ 0 , T ] → R are basis functions of S p ⁢ ( T N ) , satisfying

(1.8) φ j ⁢ ( 0 ) = 0 for all ⁢ j ∈ { 2 , … , M } ,

i.e., φ 1 ⁢ ( 0 ) ≠ 0 . Further, we define the subspace

S 0 , p ( T N ) := S p ( T N ) ∩ H 0 , 1 ( 0 , T ) = span { φ j + 1 } j = 1 M - 1 .

Thus, a conforming finite element method of the variational formulation (1.4) is to find u h ∈ S 0 , p ⁢ ( T N ) such that

(1.9) ⟨ ∂ t ⁡ u h , H T ⁢ v h ⟩ L 2 ⁢ ( 0 , T ) + μ ⁢ ⟨ u h , H T ⁢ v h ⟩ L 2 ⁢ ( 0 , T ) = ⟨ f , H T ⁢ v h ⟩ ( 0 , T ) for all ⁢ v h ∈ S 0 , p ⁢ ( T N ) .

The discrete variational formulation (1.9) is uniquely solvable, see [15, 18], and is equivalent to the linear system

(1.10) ( A ~ H T + μ ⁢ M ~ H T ) ⁢ u ¯ = F ¯ H T

with the matrices

(1.11) M ~ H T ⁢ [ i , j ] := ⟨ φ j + 1 , H T ⁢ φ i + 1 ⟩ L 2 ⁢ ( 0 , T ) ,

(1.12) A ~ H T ⁢ [ i , j ] := ⟨ ∂ t ⁡ φ j + 1 , H T ⁢ φ i + 1 ⟩ L 2 ⁢ ( 0 , T )

for i , j = 1 , … , M - 1 , and the corresponding right-hand side F ¯ H T .

Analogously, the conforming FEM of the variational formulation (1.5) to find u h ∈ S 0 , p ⁢ ( T N ) such that

⟨ H T ⁢ ∂ t ⁡ u h , ∂ t ⁡ v h ⟩ L 2 ⁢ ( 0 , T ) + μ ⁢ ⟨ u h , H T ⁢ v h ⟩ L 2 ⁢ ( 0 , T ) = ⟨ f , H T ⁢ v h ⟩ ( 0 , T ) for all ⁢ v h ∈ S 0 , p ⁢ ( T N )

is uniquely solvable, see [9, 10], and is equivalent to the linear system

( ( B ~ H T ) ⊤ + μ ⁢ M ~ H T ) ⁢ u ¯ = F ¯ H T

with the matrix M ~ H T , given in (1.11), and the matrix

(1.13) B ~ H T ⁢ [ i , j ] := ⟨ ∂ t ⁡ φ j + 1 , H T ⁢ ∂ t ⁡ φ i + 1 ⟩ L 2 ⁢ ( 0 , T )

for i , j = 1 , … , M - 1 , and the right-hand side F ¯ H T as for the linear system (1.10). Note that, for f ∈ L 2 ⁢ ( 0 , T ) , this right-hand side F ¯ H T can be realized, using the matrix

(1.14) M H T ⁢ [ i , j ] := ⟨ φ j , H T ⁢ φ i ⟩ L 2 ⁢ ( 0 , T )

for i , j = 1 , … , M , see [19, 21] for all the details. Additionally, we set

(1.15) A H T ⁢ [ i , j ] := ⟨ ∂ t ⁡ φ j , H T ⁢ φ i ⟩ L 2 ⁢ ( 0 , T ) ,

(1.16) B H T ⁢ [ i , j ] := ⟨ ∂ t ⁡ φ j , H T ⁢ ∂ t ⁡ φ i ⟩ L 2 ⁢ ( 0 , T )

for i , j = 1 , … , M . Hence, the matrices M ~ H T , A ~ H T , B ~ H T , given in (1.11), (1.12), (1.13), are submatrices of the matrices M H T , A H T , B H T , given in (1.14), (1.15), (1.16), respectively. Thus, we ask for an accurate realization of these matrices M H T , A H T , B H T , which is the main topic of this manuscript.

The rest of the paper is organized as follows. In Section 2, we recall well-known and prove new weakly singular integral representations of H T , which are the starting point of the calculation of the matrices M H T , A H T , B H T . In Section 3, Gauss quadratures are given, which are used in the following sections. Section 4 is the main part of this paper, where we state all quadrature schemes to calculate the matrices M H T , A H T , B H T to high accuracy. In Section 5, numerical examples show the quality of the new assembling method of M H T , A H T , B H T . In Section 6, we give some conclusions.

2 Integral Representations of the Modified Hilbert Transformation

In this section, we recall a weakly singular integral representation of H T for functions in H 1 ⁢ ( 0 , T ) and prove a new weakly singular integral representation of H T for functions, which are only piecewise in H 1 . These weakly singular integral representations are used for the calculation of the matrices M H T , A H T , B H T , given in (1.14), (1.15), (1.16), respectively. We start to recall the integral representations for functions in L 2 ⁢ ( 0 , T ) or H 1 ⁢ ( 0 , T ) , which are proven in [16].

Lemma 1

Lemma 1 ([16, Lemma 2.1])

For v ∈ L 2 ⁢ ( 0 , T ) , the operator H T , as defined in (1.3), allows the integral representation

(2.1) ( H T ⁢ v ) ⁢ ( t ) = v.p. ⁢ ∫ 0 T v ⁢ ( s ) ⁢ K ⁢ ( s , t ) ⁢ d s , t ∈ ( 0 , T ) ,

as a Cauchy principal value integral, where the kernel function is given as

K ⁢ ( s , t ) := 1 2 ⁢ T ⁢ [ 1 sin ⁡ π ⁢ ( s + t ) 2 ⁢ T + 1 sin ⁡ π ⁢ ( s - t ) 2 ⁢ T ] .

For v ∈ H 1 ⁢ ( 0 , T ) , the operator H T , as defined in (1.3), allows the integral representation

(2.2) ( H T ⁢ v ) ⁢ ( t ) = - 2 π ⁢ v ⁢ ( 0 ) ⁢ ln ⁡ tan ⁡ π ⁢ t 4 ⁢ T + ∫ 0 T ∂ t ⁡ v ⁢ ( s ) ⁢ K ⁢ ( s , t ) ⁢ d s , t ∈ ( 0 , T ) ,

as a weakly singular integral, where

(2.3) K ⁢ ( s , t ) := - 1 π ⁢ ln ⁡ [ tan ⁡ π ⁢ ( s + t ) 4 ⁢ T ⁢ tan ⁡ π ⁢ | t - s | 4 ⁢ T ] .

The integral representation (2.2) in connection with specialized numerical integration gives the possibility to calculate matrices, which are based on H T applied to a function in H 1 ⁢ ( 0 , T ) , e.g., the matrices M H T in (1.14) and A H T in (1.15). However, the integral representation (2.2) is not applicable for functions which have jumps since these functions do not belong to H 1 ⁢ ( 0 , T ) . As the integral representation (2.1) is a Cauchy principal value integral and therefore is complicated to use in practical implementations, we prove a new integral representation of H T for functions which are piecewise in H 1 but possibly having jumps. This integral representation is used for the calculation of the matrix B H T in (1.16).

Lemma 2

Let a , b ∈ [ 0 , T ] with a < b , and let a function f ∈ H 1 ⁢ ( a , b ) ⊂ C ⁢ [ a , b ] be given. Then, for the function v f ∈ L 2 ⁢ ( 0 , T ) ,

v f ⁢ ( s ) = { f ⁢ ( s ) , s ∈ [ a , b ] , 0 otherwise ,

the operator H T , as defined in (1.3), allows the integral representation

(2.4) ( H T ⁢ v f ) ⁢ ( t ) = - f ⁢ ( b ) ⁢ K ⁢ ( b , t ) + f ⁢ ( a ) ⁢ K ⁢ ( a , t ) + ∫ a b ∂ t ⁡ f ⁢ ( s ) ⁢ K ⁢ ( s , t ) ⁢ d s

for t ∈ ( 0 , T ) ∖ { a , b } as a weakly singular integral, where the kernel function K ⁢ ( s , t ) is given in (2.3).

If, in addition, f ⁢ ( a ) = 0 , i.e., f ∈ H 0 , 1 ⁢ ( a , b ) , then the representation

( H T ⁢ v f ) ⁢ ( a ) = - f ⁢ ( b ) ⁢ K ⁢ ( b , a ) + ∫ a b ∂ t ⁡ f ⁢ ( s ) ⁢ K ⁢ ( s , a ) ⁢ d s

holds true as a weakly singular integral. Analogously, if f ⁢ ( b ) = 0 , i.e., f ∈ H , 0 1 ⁢ ( a , b ) , then we have the representation

( H T ⁢ v f ) ⁢ ( b ) = f ⁢ ( a ) ⁢ K ⁢ ( a , b ) + ∫ a b ∂ t ⁡ f ⁢ ( s ) ⁢ K ⁢ ( s , b ) ⁢ d s

as a weakly singular integral.

Proof

Let t ∈ ( 0 , T ) be arbitrary but fixed. We split the proof accordingly to the relation between 𝑡, 𝑎 and 𝑏.

First, consider the case t ∈ ( a , b ) . The integral representation (2.1), integration by parts, the relation - ∂ s ⁡ K ⁢ ( s , t ) = K ⁢ ( s , t ) and the Hölder continuity of 𝑓 with exponent 1 2 , see [11, Chapitre 2, Théorème 3.8], yield

( H T ⁢ v f ) ⁢ ( t ) = lim ε ↘ 0 ⁡ ( ∫ a t - ε f ⁢ ( s ) ⁢ K ⁢ ( s , t ) ⁢ d ⁢ s + ∫ t + ε b f ⁢ ( s ) ⁢ K ⁢ ( s , t ) ⁢ d ⁢ s ) = lim ε ↘ 0 ( - f ( t - ε ) K ( t - ε , t ) + f ( a ) K ( a , t ) + ∫ a t - ε ∂ t f ( s ) K ( s , t ) d s - f ( b ) K ( b , t ) + f ( t + ε ) K ( t + ε , t ) + ∫ t + ε b ∂ t f ( s ) K ( s , t ) d s ) = f ⁢ ( a ) ⁢ K ⁢ ( a , t ) - f ⁢ ( b ) ⁢ K ⁢ ( b , t ) + ∫ a b ∂ t ⁡ f ⁢ ( s ) ⁢ K ⁢ ( s , t ) ⁢ d s ,

where the last integral exists as a weakly singular integral.

Second, we examine the case t ∈ ( 0 , a ) when a > 0 , or t ∈ ( b , T ) when b < T . The integral representation (2.1) and integration by parts give

( H T ⁢ v f ) ⁢ ( t ) = lim ε ↘ 0 ⁡ ∫ a b f ⁢ ( s ) ⁢ K ⁢ ( s , t ) ⁢ d s = f ⁢ ( a ) ⁢ K ⁢ ( a , t ) - f ⁢ ( b ) ⁢ K ⁢ ( b , t ) + ∫ a b ∂ t ⁡ f ⁢ ( s ) ⁢ K ⁢ ( s , t ) ⁢ d s .

Third, the case t = a is investigated when f ⁢ ( a ) = 0 . As in the cases before, the integral representation (2.1), the Hölder continuity of 𝑓 and integration by parts lead to

( H T ⁢ v f ) ⁢ ( a ) = lim ε ↘ 0 ⁡ ∫ a + ε b f ⁢ ( s ) ⁢ K ⁢ ( s , a ) ⁢ d s = lim ε ↘ 0 ⁡ ( f ⁢ ( a + ε ) ⁢ K ⁢ ( a + ε , a ) - f ⁢ ( b ) ⁢ K ⁢ ( b , a ) + ∫ a b ∂ t ⁡ f ⁢ ( s ) ⁢ K ⁢ ( s , a ) ⁢ d ⁢ s )

= - f ⁢ ( b ) ⁢ K ⁢ ( b , a ) + ∫ a b ∂ t ⁡ f ⁢ ( s ) ⁢ K ⁢ ( s , a ) ⁢ d s ,

where the last integral exists as a weakly singular integral.

Last, the case t = b , when f ⁢ ( b ) = 0 , is proven analogously. ∎

Remark 1

Choose a = 0 and b = T in Lemma 2. Then representation (2.4) of Lemma 2 coincides with representation (2.2) of Lemma 1 due to K ⁢ ( 0 , t ) = - 2 π ⁢ ln ⁡ tan ⁡ π ⁢ t 4 ⁢ T and K ⁢ ( T , t ) = 0 for t ∈ ( 0 , T ) .

Remark 2

Consider the situation of Lemma 2. For f ∈ H 1 ⁢ ( a , b ) with f ⁢ ( a ) ≠ 0 or f ⁢ ( b ) ≠ 0 , the function H T ⁢ v f ∈ L 2 ⁢ ( 0 , T ) admits a singularity for t = a or t = b , respectively.

3 Gauss Quadratures

In this section, Gauss quadratures are stated, which are used in the following sections. For this purpose, we introduce Gauss quadratures on [ 0 , 1 ] of order K ∈ N , see [3]. First, the classical Gauss–Legendre quadrature with weights ω ν , K ∈ R and nodes ξ ν , K ∈ [ 0 , 1 ] , ν = 1 , … , K , fulfills

(3.1) ∫ 0 1 g ⁢ ( t ) ⁢ d t = ∑ ν = 1 K ω ν , K ⁢ g ⁢ ( ξ ν , K ) for all ⁢ g ∈ P 2 ⁢ K - 1 ⁢ [ 0 , 1 ] ,

see [3, equation (1.6)]. This Gauss quadrature is generalized to the square domain [ 0 , 1 ] × [ 0 , 1 ] via using the Gauss–Legendre quadrature (3.1) for each coordinate direction. More precisely, the tensor Gauss quadrature

(3.2) ∑ ν 1 = 1 K 1 ∑ ν 2 = 1 K 2 ω ν 1 , K 1 ⁢ ω ν 2 , K 2 ⁢ G ⁢ ( ξ ν 1 , K 1 , ξ ν 2 , K 2 )

approximates the integral

∫ 0 1 ∫ 0 1 G ⁢ ( ξ 1 , ξ 2 ) ⁢ d ξ 1 ⁢ d ξ 2

for a function G : [ 0 , 1 ] × [ 0 , 1 ] → R , where ω ν i , K i ∈ R and ξ ν i , K i ∈ [ 0 , 1 ] are the Gauss integration weights and Gauss integration nodes of order K i ∈ N with respect to the coordinate direction i ∈ { 1 , 2 } .

For singular integrals with a logarithmic term, we introduce a nonclassical Gauss–Jacobi quadrature. In greater detail, for K ∈ N , the weights ω ^ ν , K ∈ R and nodes ξ ^ ν , K ∈ [ 0 , 1 ] , ν = 1 , … , K , of the Gauss–Jacobi quadrature for a logarithmic term satisfy

(3.3) - ∫ 0 1 g ⁢ ( t ) ⁢ ln ⁡ ( t ) ⁢ d t = ∑ ν = 1 K ω ^ ν , K ⁢ g ⁢ ( ξ ^ ν , K ) for all ⁢ g ∈ P 2 ⁢ K - 1 ⁢ [ 0 , 1 ] ,

see [3, equation (1.7)]. With quadrature (3.3), the equality

(3.4) ∫ 0 1 ∫ 0 1 g ⁢ ( s , t ) ln ⁡ | s - t | ⁢ d ⁢ s ⁢ d ⁢ t = - ∑ μ = 1 K ∑ ν = 1 K ω ^ μ ⁢ ξ ^ μ ⁢ ω ν ⁢ ( g ⁢ ( ( 1 - ξ ν ) ⁢ ξ ^ μ , ξ ^ μ ) + g ⁢ ( 1 - ( 1 - ξ ν ) ⁢ ξ ^ μ , 1 - ξ ^ μ ) ) - ∑ μ = 1 K ∑ ν = 1 K ω μ ⁢ ξ μ ⁢ ω ^ ν ⁢ ( g ⁢ ( ( 1 - ξ ^ ν ) ⁢ ξ μ , ξ μ ) + g ⁢ ( 1 - ( 1 - ξ ^ ν ) ⁢ ξ μ , 1 - ξ μ ) )

holds true for all g ∈ P 2 ⁢ K - 2 ⁢ ( [ 0 , 1 ] × [ 0 , 1 ] ) , where

ω ν = ω ν , K , ω ^ ν = ω ^ ν , K , ξ ν = ξ ν , K , ξ ^ ν = ξ ^ ν , K

are defined in (3.1), (3.3), see [3, equation (2.7)], [4].

4 Quadrature Schemes for the Modified Hilbert Transformation

In this section, we describe the assembling of the matrices M H T in (1.14), A H T in (1.15) and B H T in (1.16). The crucial point is the realization of the modified Hilbert transformation H T , where different possibilities exist, see [16, 19]. In particular, for a uniform degree vector p = ( p , p , … , p ) with a fixed, low polynomial degree p ∈ N , e.g., p = 1 or p = 2 , the matrices M H T , A H T and B H T in (1.14), (1.15), (1.16), respectively, can be calculated using a series expansion based on the Legendre chi function, which converges very fast, independently of the mesh sizes, see [19, Subsection 2.2]. As for an h ⁢ p -FEM, the degree vector 𝒑 is not uniform, or for high polynomial degrees, it is convenient to apply numerical quadrature rules to approximate the matrix entries of M H T , A H T and B H T in (1.14), (1.15), (1.16), respectively. This is the main topic of the paper and is described in great detail in the following.

From the integral representation of H T in (2.2), we have

(4.1) M H T ⁢ [ i , j ] = ⟨ φ j , H T ⁢ φ i ⟩ L 2 ⁢ ( 0 , T ) = - 2 π ⁢ φ i ⁢ ( 0 ) ⁢ ∫ 0 T φ j ⁢ ( t ) ⁢ ln ⁡ tan ⁡ π ⁢ t 4 ⁢ T ⁢ d ⁢ t + ∫ 0 T φ j ⁢ ( t ) ⁢ ∫ 0 T K ⁢ ( s , t ) ⁢ ∂ t ⁡ φ i ⁢ ( s ) ⁢ d s ⁢ d t ,

(4.2) A H T ⁢ [ i , j ] = ⟨ ∂ t ⁡ φ j , H T ⁢ φ i ⟩ L 2 ⁢ ( 0 , T ) = - 2 π ⁢ φ i ⁢ ( 0 ) ⁢ ∫ 0 T ∂ t ⁡ φ j ⁢ ( t ) ⁢ ln ⁡ tan ⁡ π ⁢ t 4 ⁢ T ⁢ d ⁢ t + ∫ 0 T ∂ t ⁡ φ j ⁢ ( t ) ⁢ ∫ 0 T K ⁢ ( s , t ) ⁢ ∂ t ⁡ φ i ⁢ ( s ) ⁢ d s ⁢ d t

for i , j = 1 , … , M with the basis functions φ i , φ j in (1.7) satisfying (1.8).

For the matrix B H T in (1.16), the integral representation of H T in (2.4) yields

H T ⁢ ∂ t ⁡ φ i ⁢ ( t ) = ∑ k = 1 N [ - ( ∂ t ⁡ φ i ) - k ⁢ K ⁢ ( t k , t ) + ( ∂ t ⁡ φ i ) + k - 1 ⁢ K ⁢ ( t k - 1 , t ) + ∫ t k - 1 t k ∂ t ⁢ t ⁡ ( φ i | τ k ) ⁡ ( s ) ⁢ K ⁢ ( s , t ) ⁢ d s ]

for t ∈ ( 0 , T ) ∖ { t k : k = 1 , … , N - 1 } , i = 1 , … , M , where we use the notation

v - k := lim ε ↘ 0 ⁡ v ⁢ ( t k - ε ) and v + k - 1 := lim ε ↘ 0 ⁡ v ⁢ ( t k - 1 + ε ) , k = 1 , … , N ,

for a sufficiently smooth function v : ( 0 , T ) → R . Thus, the entries of the matrix in (1.16) admit the representation

(4.3) B H T ⁢ [ i , j ] = ⟨ ∂ t ⁡ φ j , H T ⁢ ∂ t ⁡ φ i ⟩ L 2 ⁢ ( 0 , T ) = ∑ k = 1 N [ - ( ∂ t φ i ) - k ∫ 0 T ∂ t φ j ( t ) K ( t k , t ) d t + ( ∂ t φ i ) + k - 1 ∫ 0 T ∂ t φ j ( t ) K ( t k - 1 , t ) d t + ∫ 0 T ∂ t φ j ( t ) ∫ t k - 1 t k ∂ t ⁢ t ( φ i | τ k ) ( s ) K ( s , t ) d s d t ]

for i , j = 1 , … , M with the basis functions φ i , φ j in (1.7) satisfying (1.8).

The matrix entries M H T ⁢ [ i , j ] , A H T ⁢ [ i , j ] , B H T ⁢ [ i , j ] in (4.1), (4.2), (4.3), respectively, are computed element-wise for the partition T N = { τ ℓ } ℓ = 1 N of ( 0 , T ) into intervals τ ℓ = ( t ℓ - 1 , t ℓ ) ⊂ ( 0 , T ) , ℓ = 1 , … , N . For this purpose, we assume that we use standard Lagrange finite elements as in [2, Section 6.3], or for the h ⁢ p -FEM, we use basis functions φ i , i = 1 , … , M , as described in [13, Subsection 3.1.6]. Thus, we assume that the basis functions φ i : [ 0 , T ] → R , i = 1 , … , M , allow for a local representation on an element τ ℓ , ℓ = 1 , … , N , by shape functions ψ m p ℓ : [ 0 , 1 ] → R , m = 1 , … , p ℓ + 1 , defined on the reference interval [ 0 , 1 ] . In greater detail, for a basis function φ i , i = 1 , … , M , there exists an element τ ℓ = ( t ℓ - 1 , t ℓ ) , ℓ = 1 , … , N , with a related polynomial degree p ℓ ∈ N such that the local representation

(4.4) φ α ⁢ ( m , ℓ ) ⁢ ( t ) = ψ m p ℓ ⁢ ( t - t ℓ - 1 h ℓ ) for all ⁢ t ∈ τ ¯ ℓ

holds true, where α ⁢ ( m , ℓ ) = i ∈ { 1 , 2 , … , M } is the global index related to the local index 𝑚 for the interval τ ℓ and ψ m p ℓ is a shape function. Possible choices of such shape functions are the classical Lagrange polynomials with equidistant nodes or Gauss–Lobatto nodes, see [2, Chapter 6], or Lobatto polynomials (integrated Legendre polynomials), see [13, Subsection 3.1.4].

With this notation, we fix two intervals τ ℓ = ( t ℓ - 1 , t ℓ ) , τ k = ( t k - 1 , t k ) with indices ℓ , k ∈ { 1 , … , N } and related local polynomial degrees p ℓ , p k ∈ N . We define the local matrix M k , ℓ H T ∈ R ( p k + 1 ) × ( p ℓ + 1 ) by

(4.5) M k , ℓ H T ⁢ [ n , m ] = - 2 π ⁢ φ α ⁢ ( n , k ) ⁢ ( 0 ) ⁢ ∫ t ℓ - 1 t ℓ φ α ⁢ ( m , ℓ ) ⁢ ( t ) ⁢ ln ⁡ tan ⁡ π ⁢ t 4 ⁢ T ⁢ d ⁢ t + ∫ t ℓ - 1 t ℓ φ α ⁢ ( m , ℓ ) ⁢ ( t ) ⁢ ∫ t k - 1 t k K ⁢ ( s , t ) ⁢ ∂ t ⁡ φ α ⁢ ( n , k ) ⁢ ( s ) ⁢ d s ⁢ d t = - δ k , 1 ⁢ ψ n p k ⁢ ( 0 ) ⁢ 2 ⁢ h ℓ π ⁢ ∫ 0 1 ψ m p ℓ ⁢ ( η ) ⁢ ln ⁡ tan ⁡ π ⁢ ( t ℓ - 1 + η ⁢ h ℓ ) 4 ⁢ T ⁢ d ⁢ η + h ℓ ⁢ ∫ 0 1 ψ m p ℓ ⁢ ( η ) ⁢ ∫ 0 1 K ⁢ ( t k - 1 + ξ ⁢ h k , t ℓ - 1 + η ⁢ h ℓ ) ⁢ ∂ t ⁡ ψ n p k ⁢ ( ξ ) ⁢ d ξ ⁢ d η

and, analogously, the local matrix A k , ℓ H T ∈ R ( p k + 1 ) × ( p ℓ + 1 ) by

(4.6) A k , ℓ H T ⁢ [ n , m ] = - δ k , 1 ⁢ ψ n p k ⁢ ( 0 ) ⁢ 2 π ⁢ ∫ 0 1 ∂ t ⁡ ψ m p ℓ ⁢ ( η ) ⁢ ln ⁡ tan ⁡ π ⁢ ( t ℓ - 1 + η ⁢ h ℓ ) 4 ⁢ T ⁢ d ⁢ η + ∫ 0 1 ∂ t ⁡ ψ m p ℓ ⁢ ( η ) ⁢ ∫ 0 1 K ⁢ ( t k - 1 + ξ ⁢ h k , t ℓ - 1 + η ⁢ h ℓ ) ⁢ ∂ t ⁡ ψ n p k ⁢ ( ξ ) ⁢ d ξ ⁢ d η

for n = 1 , … , p k + 1 and m = 1 , … , p ℓ + 1 , where we used property (1.8), the local representation (4.4) and δ k , 1 is the usual Kronecker delta.

In the same way, the local matrix B k , ℓ H T ∈ R ( p k + 1 ) × ( p ℓ + 1 ) is defined by

(4.7) B k , ℓ H T ⁢ [ n , m ] = - 1 π ⁢ h k ⁢ ∂ t ⁡ ψ n p k ⁢ ( 0 ) ⁢ J k , ℓ 0 ⁢ [ m ] + 1 π ⁢ h k ⁢ ∂ t ⁡ ψ n p k ⁢ ( 1 ) ⁢ J k , ℓ 1 ⁢ [ m ] + 1 h k ⁢ ∫ 0 1 ∂ t ⁡ ψ m p ℓ ⁢ ( η ) ⁢ ∫ 0 1 K ⁢ ( t k - 1 + ξ ⁢ h k , t ℓ - 1 + η ⁢ h ℓ ) ⁢ ∂ t ⁢ t ⁡ ψ n p k ⁢ ( ξ ) ⁢ d ξ ⁢ d η

for n = 1 , … , p k + 1 and m = 1 , … , p ℓ + 1 with the integrals

(4.8) J k , ℓ 0 ⁢ [ m ] = - π ⁢ ∫ 0 1 ∂ t ⁡ ψ m p ℓ ⁢ ( η ) ⁢ K ⁢ ( t k - 1 , t ℓ - 1 + η ⁢ h ℓ ) ⁢ d η ,

(4.9) J k , ℓ 1 ⁢ [ m ] = - π ⁢ ∫ 0 1 ∂ t ⁡ ψ m p ℓ ⁢ ( η ) ⁢ K ⁢ ( t k , t ℓ - 1 + η ⁢ h ℓ ) ⁢ d η .

Note that J k , ℓ 0 ⁢ [ m ] = J k - 1 , ℓ 1 ⁢ [ m ] , i.e., investigating only one term is another possibility. To keep the same cases as for the matrices M H T , A H T and due to implementation issues, we give the quadrature schemes for J k , ℓ 0 ⁢ [ m ] and J k , ℓ 1 ⁢ [ m ] .

These local matrices M k , ℓ H T , A k , ℓ H T , B k , ℓ H T are used for the calculation of the matrices M H T , A H T and B H T , respectively. Thus, the assembling of the matrices M H T , A H T and B H T is realized element-wise, i.e., by two loops via the elements.

Set M H T = 0 M , A H T = 0 M and B H T = 0 M , where 0 M ∈ R M × M is the zero matrix.
For k , ℓ = 1 , … , N ,
- compute the matrices M k , ℓ H T , A k , ℓ H T , B k , ℓ H T ∈ R ( p k + 1 ) × ( p ℓ + 1 ) ,
- and set
  M H T ⁢ [ α ⁢ ( n , k ) , α ⁢ ( m , ℓ ) ] = M H T ⁢ [ α ⁢ ( n , k ) , α ⁢ ( m , ℓ ) ] + M k , ℓ H T ⁢ [ n , m ] , A H T ⁢ [ α ⁢ ( n , k ) , α ⁢ ( m , ℓ ) ] = A H T ⁢ [ α ⁢ ( n , k ) , α ⁢ ( m , ℓ ) ] + A k , ℓ H T ⁢ [ n , m ] , B H T ⁢ [ α ⁢ ( n , k ) , α ⁢ ( m , ℓ ) ] = B H T ⁢ [ α ⁢ ( n , k ) , α ⁢ ( m , ℓ ) ] + B k , ℓ H T ⁢ [ n , m ]
  for n = 1 , … , p k + 1 , m = 1 , … , p ℓ + 1 .

It remains to calculate the local matrices M k , ℓ H T , A k , ℓ H T , B k , ℓ H T ∈ R ( p k + 1 ) × ( p ℓ + 1 ) , which is the main purpose of this work and the content of the following subsections, Subsection 4.1, 4.2, 4.3, respectively.

4.1 Local Matrix M k , ℓ H T

In this subsection, the computation of the local matrix M k , ℓ H T in (4.5) is investigated for fixed elements τ k , τ ℓ and local polynomial degrees p k , p ℓ with k , ℓ ∈ { 1 , … , N } . We apply the strategy of [16, Subsection 3.1], i.e., the integrals in (4.5) are split into regular and singular parts. As in [16, Subsection 3.1], we distinguish three different cases for the element indices 𝑘, ℓ, which correspond to the singularities of the integrands in (4.5). For this purpose, let n ∈ { 1 , … , p k + 1 } and m ∈ { 1 , … , p ℓ + 1 } be fixed. We investigate the case k = ℓ = 1 in great detail, whereas we state only the final results for the remaining cases since they can be treated in a simpler or similar way as the case k = ℓ = 1 . These final results are formulated as integrals such that these integrals can be directly calculated or approximated by quadratures (3.1), (3.2), (3.3), (3.4).

4.1.1 Case k = ℓ = 1

We have

M 1 , 1 H T ⁢ [ n , m ] = - ψ n p 1 ⁢ ( 0 ) ⁢ 2 ⁢ h 1 π ⁢ ∫ 0 1 ψ m p 1 ⁢ ( η ) ⁢ F 1 , 1 1 ⁢ ( η ⁢ h 1 ) ⁢ d η - ψ n p 1 ⁢ ( 0 ) ⁢ 2 ⁢ h 1 π ⁢ ∫ 0 1 ψ m p 1 ⁢ ( η ) ⁢ ln ⁡ η ⁢ d ⁢ η ⏟ = ⁣ : I 1 , 1 1 - ψ n p 1 ⁢ ( 0 ) ⁢ 2 ⁢ h 1 π ⁢ ln ⁡ h 1 ⁢ ∫ 0 1 ψ m p 1 ⁢ ( η ) ⁢ d ⁢ η - h 1 π ⁢ ∫ 0 1 ψ m p 1 ⁢ ( η ) ⁢ ∫ 0 1 F 1 , 1 2 ⁢ ( ξ ⁢ h 1 , η ⁢ h 1 ) ⁢ ∂ t ⁡ ψ n p 1 ⁢ ( ξ ) ⁢ d ⁢ ξ ⁢ d ⁢ η - h 1 π ⁢ ∫ 0 1 ψ m p 1 ( η ) ∫ 0 1 ln | s - t | | s = ξ h 1 , t = η h 1 ∂ t ψ n p 1 ( ξ ) d ξ d η ⏟ = ⁣ : I 1 , 1 2 - h 1 π ⁢ ∫ 0 1 ψ m p 1 ( η ) ∫ 0 1 ln ( s + t ) | s = ξ h 1 , t = η h 1 ∂ t ψ n p 1 ( ξ ) d ξ d η ⏟ = ⁣ : I 1 , 1 3

with

F 1 , 1 1 ⁢ ( t ) := ln ⁡ tan ⁡ π ⁢ t 4 ⁢ T t , F 1 , 1 2 ⁢ ( s , t ) := ln ⁡ [ tan ⁡ π ⁢ ( s + t ) 4 ⁢ T s + t ⁢ tan ⁡ π ⁢ | t - s | 4 ⁢ T | s - t | ] .

Proceeding as in [16, equation (3.10)], i.e., the square domain of integration is split into two triangles, and on each of these triangles a Duffy transformation is applied, yields

M 1 , 1 H T ⁢ [ n , m ] = - ψ n p 1 ⁢ ( 0 ) ⁢ 2 ⁢ h 1 π ⁢ I 1 , 1 1 - h 1 π ⁢ ( I 1 , 1 2 + I 1 , 1 3 ) - ψ n p 1 ⁢ ( 0 ) ⁢ 2 ⁢ h 1 π ⁢ ∫ 0 1 ψ m p 1 ⁢ ( η ) ⁢ F 1 , 1 1 ⁢ ( η ⁢ h 1 ) ⁢ d ⁢ η - ψ n p 1 ⁢ ( 0 ) ⁢ 2 ⁢ h 1 π ⁢ ln ⁡ h 1 ⁢ ∫ 0 1 ψ m p 1 ⁢ ( η ) ⁢ d ⁢ η - h 1 π ⁢ ∫ 0 1 ∫ 0 1 ψ m p 1 ⁢ ( η ) ⁢ F 1 , 1 2 ⁢ ( ( 1 - ξ ) ⁢ η ⁢ h 1 , η ⁢ h 1 ) ⁢ η ⁢ ∂ t ⁡ ψ n p 1 ⁢ ( ( 1 - ξ ) ⁢ η ) ⏟ = ⁣ : G 1 , 1 1 ⁢ ( ξ , η ) ⁢ d ξ ⁢ d η - h 1 π ⁢ ∫ 0 1 ∫ 0 1 ψ m p 1 ⁢ ( ( 1 - η ) ⁢ ξ ) ⁢ F 1 , 1 2 ⁢ ( ξ ⁢ h 1 , ( 1 - η ) ⁢ ξ ⁢ h 1 ) ⁢ ξ ⁢ ∂ t ⁡ ψ n p 1 ⁢ ( ξ ) ⏟ = ⁣ : G 1 , 1 2 ⁢ ( ξ , η ) ⁢ d η ⁢ d ξ .

Further, applying the Gauss–Legendre quadratures (3.1), (3.2) of order K ∈ N gives the approximation

M 1 , 1 H T ⁢ [ n , m ] ≈ - ψ n p 1 ⁢ ( 0 ) ⁢ 2 ⁢ h 1 π ⁢ I 1 , 1 1 - h 1 π ⁢ ( I 1 , 1 2 + I 1 , 1 3 ) - ψ n p 1 ⁢ ( 0 ) ⁢ 2 ⁢ h 1 π ⁢ ∑ ν = 1 K ω ν , K ⁢ ψ m p 1 ⁢ ( ξ ν , K ) ⁢ ( F 1 , 1 1 ⁢ ( ξ ν , K ⁢ h 1 ) + ln ⁡ h 1 ) - h 1 π ⁢ ∑ ν 1 = 1 K ∑ ν 2 = 1 K ω ν 1 , K ⁢ ω ν 2 , K ⁢ ( G 1 , 1 1 ⁢ ( ξ ν 1 , K , ξ ν 2 , K ) + G 1 , 1 2 ⁢ ( ξ ν 1 , K , ξ ν 2 , K ) ) .

Using [16, Corollary 3.1] for the two-dimensional parts and similar arguments for the one-dimensional part results in exponential convergence of the applied Gauss–Legendre quadratures with respect to the number 𝐾 of Gauss integration nodes. It remains to calculate the singular integrals I 1 , 1 1 , I 1 , 1 2 and I 1 , 1 3 . In [16, Subsection 3.1], an analytic integration is proposed, which is practical only for constant low polynomial degrees, e.g., p k = p ℓ ∈ { 1 , 2 } . For arbitrary high polynomial degrees p k , p ℓ , quadrature rules of order adapted to p k , p ℓ are easier to implement. This approach is not contained in [16] and thus is investigated in the following.

For the integral I 1 , 1 1 , the Gauss–Jacobi quadrature (3.3) for a logarithmic term yields

I 1 , 1 1 = ∫ 0 1 ψ m p 1 ⁢ ( η ) ⁢ ln ⁡ η ⁢ d ⁢ η = - ∑ ν = 1 K log ω ^ ν , K log ⁢ ψ m p 1 ⁢ ( ξ ^ ν , K log )

with p 1 + 1 2 ≤ K log ∈ N .

For the integral I 1 , 1 2 , we apply the Gauss–Legendre quadrature (3.1) and quadrature (3.4), which lead to

I 1 , 1 2 = ln ⁡ h 1 ⁢ ∫ 0 1 ψ m p 1 ⁢ ( η ) ⁢ ∫ 0 1 ∂ t ⁡ ψ n p 1 ⁢ ( ξ ) ⁢ d ξ ⁢ d η + ∫ 0 1 ∫ 0 1 ψ m p 1 ⁢ ( η ) ⁢ ∂ t ⁡ ψ n p 1 ⁢ ( ξ ) ⏟ = H 1 , 1 ⁢ ( ξ , η ) ⁢ ln ⁡ | ξ - η | ⁢ d ξ ⁢ d η = ln ⁡ h 1 ⁢ ( ψ n p 1 ⁢ ( 1 ) - ψ n p 1 ⁢ ( 0 ) ) ⁢ ∑ ν = 1 K reg ω ν , K reg ⁢ ψ m p 1 ⁢ ( ξ ν , K reg ) - ∑ μ = 1 K log ∑ ν = 1 K log ω ^ μ ⁢ ξ ^ μ ⁢ ω ν ⁢ ( H 1 , 1 ⁢ ( ( 1 - ξ ν ) ⁢ ξ ^ μ , ξ ^ μ ) + H 1 , 1 ⁢ ( 1 - ( 1 - ξ ν ) ⁢ ξ ^ μ , 1 - ξ ^ μ ) ) - ∑ μ = 1 K log ∑ ν = 1 K log ω μ ⁢ ξ μ ⁢ ω ^ ν ⁢ ( H 1 , 1 ⁢ ( ( 1 - ξ ^ ν ) ⁢ ξ μ , ξ μ ) + H 1 , 1 ⁢ ( 1 - ( 1 - ξ ^ ν ) ⁢ ξ μ , 1 - ξ μ ) )

with p 1 + 1 2 ≤ K reg ∈ N , p 1 + 1 2 ≤ K log ∈ N and the function H 1 , 1 ⁢ ( ξ , η ) := ψ m p 1 ⁢ ( η ) ⁢ ∂ t ⁡ ψ n p 1 ⁢ ( ξ ) .

For the integral I 1 , 1 3 , the square domain of integration of the first line is split into two triangles, and on each of these triangles, a Duffy transformation is applied, which yields

I 1 , 1 3 = ln ⁡ h 1 ⁢ ∫ 0 1 ψ m p 1 ⁢ ( η ) ⁢ ∫ 0 1 ∂ t ⁡ ψ n p 1 ⁢ ( ξ ) ⁢ d ξ ⁢ d η + ∫ 0 1 ∫ 0 1 ψ m p 1 ⁢ ( η ) ⁢ ln ⁡ ( ξ + η ) ⁢ ∂ t ⁡ ψ n p 1 ⁢ ( ξ ) ⁢ d ξ ⁢ d η = ln ⁡ h 1 ⁢ ∫ 0 1 ψ m p 1 ⁢ ( η ) ⁢ ∫ 0 1 ∂ t ⁡ ψ n p 1 ⁢ ( ξ ) ⁢ d ξ ⁢ d η + ∫ 0 1 ∫ 0 1 ψ m p 1 ⁢ ( η ) ⁢ η ⁢ ln ⁡ ( η ⁢ ξ + η ) ⏟ = ln ⁡ η + ln ⁡ ( 1 + ξ ) ⁢ ∂ t ⁡ ψ n p 1 ⁢ ( ξ ⁢ η ) ⁢ d ξ ⁢ d η + ∫ 0 1 ∫ 0 1 ψ m p 1 ⁢ ( ξ ⁢ η ) ⁢ ξ ⁢ ln ⁡ ( ξ + ξ ⁢ η ) ⏟ = ln ⁡ ξ + ln ⁡ ( 1 + η ) ⁢ ∂ t ⁡ ψ n p 1 ⁢ ( ξ ) ⁢ d ξ ⁢ d η .

Further, using quadratures (3.1), (3.3), (3.2) gives

I 1 , 1 3 ≈ ln ⁡ h 1 ⁢ ( ψ n p 1 ⁢ ( 1 ) - ψ n p 1 ⁢ ( 0 ) ) ⁢ ∑ ν = 1 K 1 ω ν , K 1 ⁢ ψ m p 1 ⁢ ( ξ ν , K 1 ) - ∑ μ = 1 K 2 ∑ ν = 1 K 3 ω μ , K 2 ⁢ ω ^ ν , K 3 ⁢ H 1 3 ⁢ ( ξ μ , K 2 , ξ ^ ν , K 3 ) - ∑ μ = 1 K 4 ∑ ν = 1 K 5 ω ^ μ , K 4 ⁢ ω ν , K 5 ⁢ H 2 3 ⁢ ( ξ ^ μ , K 4 , ξ ν , K 5 ) + ∑ ν 1 = 1 K ∑ ν 2 = 1 K ω ν 1 , K ⁢ ω ν 2 , K ⁢ H 3 3 ⁢ ( ξ ν 1 , K , ξ ν 2 , K )

with the functions

H 1 3 ⁢ ( ξ , η ) := ψ m p 1 ⁢ ( η ) ⁢ η ⁢ ∂ t ⁡ ψ n p 1 ⁢ ( ξ ⁢ η ) , H 2 3 ⁢ ( ξ , η ) := ψ m p 1 ⁢ ( ξ ⁢ η ) ⁢ ξ ⁢ ∂ t ⁡ ψ n p 1 ⁢ ( ξ ) , H 3 3 ⁢ ( ξ , η ) := ψ m p 1 ⁢ ( η ) ⁢ η ⁢ ln ⁡ ( 1 + ξ ) ⁢ ∂ t ⁡ ψ n p 1 ⁢ ( ξ ⁢ η ) + ψ m p 1 ⁢ ( ξ ⁢ η ) ⁢ ξ ⁢ ln ⁡ ( 1 + η ) ⁢ ∂ t ⁡ ψ n p 1 ⁢ ( ξ )

and the number of integration nodes

p 1 + 1 2 ≤ K 1 ∈ N , p 1 2 ≤ K 2 ∈ N , p 1 + 1 2 ≤ K 3 ∈ N , p 1 + 1 2 ≤ K 4 ∈ N , p 1 + 1 2 ≤ K 5 ∈ N

and K ∈ N sufficiently large. Note that all integrals of I 1 , 1 3 are calculated exactly except for the integral with the integrand H 3 3 , where arguments similar to [16, Corollary 3.1] yield exponential convergence of the applied Gauss–Legendre quadratures with respect to the number 𝐾 of Gauss integration nodes.

4.1.2 Case k = ℓ = N

For this and the remaining cases, we state only the integrals, where the quadratures of Section 3 have to be applied. We calculate

M N , N H T ⁢ [ n , m ] = - h N π ⁢ ∫ 0 1 ∫ 0 1 ψ m p N ⁢ ( η ) ⁢ G N , N 1 ⁢ ( ξ , η ) ⁢ ∂ t ⁡ ψ n p N ⁢ ( ( 1 - ξ ) ⁢ η ) ⁢ d ξ ⁢ d η - h N π ⁢ ∫ 0 1 ∫ 0 1 ψ m p N ⁢ ( ( 1 - η ) ⁢ ξ ) ⁢ G N , N 2 ⁢ ( ξ , η ) ⁢ ∂ t ⁡ ψ n p N ⁢ ( ξ ) ⁢ d η ⁢ d ξ - h N π ⁢ ∫ 0 1 ∫ 0 1 ψ m p N ⁢ ( η ) ⁢ ∂ t ⁡ ψ n p N ⁢ ( ξ ) ⁢ ln ⁡ | ξ - η | ⁢ d ξ ⁢ d η + h N π ⁢ ∫ 0 1 ∫ 0 1 ψ m p N ⁢ ( 1 - η ) ⁢ ∂ t ⁡ ψ n p N ⁢ ( 1 - η ⁢ ξ ) ⁢ [ ln ⁡ η + ln ⁡ ( 1 + ξ ) ] ⁢ η ⁢ d ξ ⁢ d η + h N π ⁢ ∫ 0 1 ∫ 0 1 ψ m p N ⁢ ( 1 - ξ ⁢ η ) ⁢ ∂ t ⁡ ψ n p N ⁢ ( 1 - ξ ) ⁢ [ ln ⁡ ξ + ln ⁡ ( 1 + η ) ] ⁢ ξ ⁢ d ξ ⁢ d η

with the functions

F N , N ⁢ ( s , t ) := ln ⁡ [ tan ⁡ π ⁢ ( s + t ) 4 ⁢ T ⁢ ( 2 ⁢ T - s - t ) ⁢ tan ⁡ π ⁢ | t - s | 4 ⁢ T | s - t | ] , G N , N 1 ⁢ ( ξ , η ) := F N , N ⁢ ( t N - 1 + ( 1 - ξ ) ⁢ η ⁢ h N , t N - 1 + η ⁢ h N ) ⁢ η , G N , N 2 ⁢ ( ξ , η ) := F N , N ⁢ ( t N - 1 + ξ ⁢ h N , t N - 1 + ( 1 - η ) ⁢ ξ ⁢ h N ) ⁢ ξ ,

which do not have any singularity in the domain of integration. For the first and the second integral, we apply the tensor Gauss quadrature (3.2), for the third integral, quadrature (3.4), and for the last two integrals, quadratures (3.1), (3.3), as shown in Subsection 4.1.1.

4.1.3 Cases Excluding k = ℓ = 1 and k = ℓ = N

We compute

M k , ℓ H T ⁢ [ n , m ] = - h ℓ π ⁢ ∫ 0 1 ∫ 0 1 ψ m p ℓ ⁢ ( η ) ⁢ G k , ℓ 1 ⁢ ( ξ , η ) ⁢ ∂ t ⁡ ψ n p k ⁢ ( ( 1 - ξ ) ⁢ η ) ⁢ d ⁢ ξ ⁢ d ⁢ η - h ℓ π ⁢ I k , ℓ 2 - h ℓ π ⁢ ∫ 0 1 ∫ 0 1 ψ m p ℓ ⁢ ( ( 1 - η ) ⁢ ξ ) ⁢ G k , ℓ 2 ⁢ ( ξ , η ) ⁢ ∂ t ⁡ ψ n p k ⁢ ( ξ ) ⁢ d η ⁢ d ξ - δ k , 1 ⁢ ψ n p k ⁢ ( 0 ) ⁢ 2 ⁢ h ℓ π ⁢ I ℓ 1

with the integrals

I ℓ 1 := ∫ 0 1 ψ m p ℓ ⁢ ( η ) ⁢ ln ⁡ tan ⁡ π ⁢ ( t ℓ - 1 + η ⁢ h ℓ ) 4 ⁢ T ⁢ d ⁢ η , I k , ℓ 2 := ∫ 0 1 ∫ 0 1 ψ m p ℓ ⁢ ( η ) ⁢ ∂ t ⁡ ψ n p k ⁢ ( ξ ) ⁢ ln ⁡ | ( t k - 1 + ξ ⁢ h k ) - ( t ℓ - 1 + η ⁢ h ℓ ) | ⁢ d ξ ⁢ d η

and with the functions

F k , ℓ ⁢ ( s , t ) := ln ⁡ [ tan ⁡ π ⁢ ( s + t ) 4 ⁢ T ⁢ tan ⁡ π ⁢ | t - s | 4 ⁢ T | s - t | ] , G k , ℓ 1 ⁢ ( ξ , η ) := F k , ℓ ⁢ ( t k - 1 + ( 1 - ξ ) ⁢ η ⁢ h k , t ℓ - 1 + η ⁢ h ℓ ) ⁢ η , G k , ℓ 2 ⁢ ( ξ , η ) := F k , ℓ ⁢ ( t k - 1 + ξ ⁢ h k , t ℓ - 1 + ( 1 - η ) ⁢ ξ ⁢ h ℓ ) ⁢ ξ ,

which do not have any singularity in the domain of integration, i.e., the tensor Gauss quadrature (3.2) is applied to these integrands. The singular parts I ℓ 1 , I k , ℓ 2 are treated differently corresponding to the indices k , ℓ .

First, we consider the part I ℓ 1 . The case ℓ = 1 is excluded since I ℓ 1 contributes only for k = 1 due to δ k , 1 = 0 for k > 1 . For ℓ > 1 , the term I ℓ 1 is not singular. Thus, we apply the Gauss–Legendre quadrature (3.1) to

I ℓ 1 = ∫ 0 1 ψ m p ℓ ⁢ ( η ) ⁢ ln ⁡ tan ⁡ π ⁢ ( t ℓ - 1 + η ⁢ h ℓ ) 4 ⁢ T ⁢ d ⁢ η

for ℓ > 1 .

Second, we consider the part I k , ℓ 2 . We distinguish four cases.

Case k = ℓ : We compute
I k , ℓ 2 = ln ⁡ h ℓ ⁢ ( ψ n p ℓ ⁢ ( 1 ) - ψ n p ℓ ⁢ ( 0 ) ) ⁢ ∫ 0 1 ψ m p ℓ ⁢ ( η ) ⁢ d η + ∫ 0 1 ∫ 0 1 ψ m p ℓ ⁢ ( η ) ⁢ ∂ t ⁡ ψ n p ℓ ⁢ ( ξ ) ⁢ ln ⁡ | ξ - η | ⁢ d ξ ⁢ d η ,
where quadratures (3.1), (3.4) are applied.
Case k = ℓ + 1 : The equality
I k , ℓ 2 = ∫ 0 1 ∫ 0 1 ψ m p ℓ ⁢ ( 1 - η ) ⁢ ∂ t ⁡ ψ n p k ⁢ ( η ⁢ ξ ) ⁢ η ⁢ ( ln ⁡ η + ln ⁡ ( ξ ⁢ h k + h ℓ ) ) ⁢ d ⁢ ξ ⁢ d ⁢ η + ∫ 0 1 ∫ 0 1 ψ m p ℓ ⁢ ( 1 - ξ + η ⁢ ξ ) ⁢ ∂ t ⁡ ψ n p k ⁢ ( ξ ) ⁢ ξ ⁢ ( ln ⁡ ξ + ln ⁡ ( ( 1 - η ) ⁢ h ℓ + h k ) ) ⁢ d ξ ⁢ d η
holds true, where we use quadratures (3.1), (3.3) and (3.2).
Case k + 1 = ℓ : We have that
I k , ℓ 2 = ∫ 0 1 ∫ 0 1 ψ m p ℓ ⁢ ( ( 1 - ξ ) ⁢ η ) ⁢ ∂ t ⁡ ψ n p k ⁢ ( 1 - η ) ⁢ η ⁢ ( ln ⁡ η + ln ⁡ ( h k + ( 1 - ξ ) ⁢ h ℓ ) ) ⁢ d ⁢ ξ ⁢ d ⁢ η + ∫ 0 1 ∫ 0 1 ψ m p ℓ ⁢ ( η ) ⁢ ∂ t ⁡ ψ n p k ⁢ ( 1 - η + ξ ⁢ η ) ⁢ η ⁢ ( ln ⁡ η + ln ⁡ ( ( 1 - ξ ) ⁢ h k + h ℓ ) ) ⁢ d ξ ⁢ d η
holds true, where we use again quadratures (3.1), (3.3) and (3.2).
Otherwise, we apply the tensor Gauss quadrature (3.2) to
I k , ℓ 2 = ∫ 0 1 ∫ 0 1 ψ m p ℓ ⁢ ( η ) ⁢ ∂ t ⁡ ψ n p k ⁢ ( ξ ) ⁢ ln ⁡ | ( t k - 1 + ξ ⁢ h k ) - ( t ℓ - 1 + η ⁢ h ℓ ) | ⁢ d ξ ⁢ d η ,
as the integral is not singular.

4.2 Local Matrix A k , ℓ H T

In this subsection, the computation of the local matrix A k , ℓ H T in (4.6) is investigated for fixed elements τ k , τ ℓ and related polynomial degrees p k , p ℓ with k , ℓ ∈ { 1 , … , N } . To calculate the local matrix A k , ℓ H T , replace in Subsection 4.1 the function ψ m p ℓ ⁢ ( ⋅ ) with the function 1 h ℓ ⁢ ∂ t ⁡ ψ m p ℓ ⁢ ( ⋅ ) in all occurring quantities.

4.3 Local Matrix B k , ℓ H T

In this subsection, the computation of the local matrix B k , ℓ H T in (4.7) is investigated for fixed elements τ k , τ ℓ and related polynomial degrees p k , p ℓ with k , ℓ ∈ { 1 , … , N } . To calculate the second line in (4.7), replace in Subsection 4.1 the function ψ m p ℓ ⁢ ( ⋅ ) with the function 1 h ℓ ⁢ ∂ t ⁡ ψ m p ℓ ⁢ ( ⋅ ) and the function ∂ t ⁡ ψ n p k ⁢ ( ⋅ ) with the function 1 h k ⁢ ∂ t ⁢ t ⁡ ψ n p k ⁢ ( ⋅ ) in all occurring quantities which are related to the last line in (4.5). Thus, it remains to compute the integrals J k , ℓ 0 ⁢ [ m ] and J k , ℓ 1 ⁢ [ m ] given in (4.8) and (4.9), respectively, where m = 1 , … , p ℓ + 1 .

As in Subsection 4.1, we distinguish three different cases for the element indices 𝑘, ℓ, which correspond to the singularities of the integrands in (4.8) and (4.9). For this purpose, let m ∈ { 1 , … , p ℓ + 1 } be fixed. We state only the final results since their derivation is given in a simpler or similar way as in Subsection 4.1. These final results are formulated as integrals such that these integrals can be directly calculated or approximated by the Gauss–Legendre quadrature (3.1) and the Gauss–Jacobi quadrature (3.3).

4.3.1 Case k = ℓ = 1

We calculate

J k , ℓ 0 ⁢ [ m ] = ∫ 0 1 ∂ t ⁡ ψ m p 1 ⁢ ( η ) ⁢ F 1 , 1 ⁢ ( 0 , η ⁢ h 1 ) ⁢ d η + 2 ⁢ ∫ 0 1 ∂ t ⁡ ψ m p 1 ⁢ ( η ) ⁢ ln ⁡ η ⁢ d ⁢ η + 2 ⁢ ln ⁡ h 1 ⁢ ( ψ m p 1 ⁢ ( 1 ) - ψ m p 1 ⁢ ( 0 ) ) ,

J k , ℓ 1 ⁢ [ m ] = ∫ 0 1 ∂ t ⁡ ψ m p 1 ⁢ ( η ) ⁢ F 1 , 1 ⁢ ( t 1 , η ⁢ h 1 ) ⁢ d η + ∫ 0 1 ∂ t ⁡ ψ m p 1 ⁢ ( η ) ⁢ ln ⁡ ( 1 + η ) ⁢ d η + ∫ 0 1 ∂ t ⁡ ψ m p 1 ⁢ ( 1 - η ) ⁢ ln ⁡ η ⁢ d ⁢ η + 2 ⁢ ln ⁡ h 1 ⁢ ( ψ m p 1 ⁢ ( 1 ) - ψ m p 1 ⁢ ( 0 ) )

with the function

F 1 , 1 ⁢ ( s , t ) := ln ⁡ [ tan ⁡ π ⁢ ( s + t ) 4 ⁢ T s + t ⁢ tan ⁡ π ⁢ | t - s | 4 ⁢ T | s - t | ] ,

which does not have any singularity in the domain of integration. For these integrals, we apply quadratures (3.1), (3.3), as shown in Subsection 4.1.1.

4.3.2 Case k = ℓ = N

We compute

J k , ℓ 0 ⁢ [ m ] = ∫ 0 1 ∂ t ⁡ ψ m p N ⁢ ( η ) ⁢ F N , N ⁢ ( t N - 1 , t N - 1 + η ⁢ h N ) ⁢ d η + ∫ 0 1 ∂ t ⁡ ψ m p N ⁢ ( η ) ⁢ ln ⁡ η ⁢ d ⁢ η - ∫ 0 1 ∂ t ⁡ ψ m p N ⁢ ( η ) ⁢ ln ⁡ ( 2 - η ) ⁢ d η

and J k , ℓ 1 ⁢ [ m ] = 0 with the function

F N , N ⁢ ( s , t ) := ln ⁡ [ tan ⁡ π ⁢ ( s + t ) 4 ⁢ T ⁢ ( 2 ⁢ T - s - t ) ⁢ tan ⁡ π ⁢ | t - s | 4 ⁢ T | s - t | ] ,

which does not have any singularity in the domain of integration. For these integrals, we apply again quadratures (3.1), (3.3), as shown in Subsection 4.1.1.

4.3.3 Cases Excluding k = ℓ = 1 and k = ℓ = N

We have

J k , ℓ 0 ⁢ [ m ] = ∫ 0 1 ∂ t ⁡ ψ m p ℓ ⁢ ( η ) ⁢ F k , ℓ ⁢ ( t k - 1 , t ℓ - 1 + η ⁢ h ℓ ) ⁢ d η + J k , ℓ 0 , sing ⁢ [ m ] , J k , ℓ 1 ⁢ [ m ] = ∫ 0 1 ∂ t ⁡ ψ m p ℓ ⁢ ( η ) ⁢ F k , ℓ ⁢ ( t k , t ℓ - 1 + η ⁢ h ℓ ) ⁢ d η + J k , ℓ 1 , sing ⁢ [ m ]

with the function

F k , ℓ ⁢ ( s , t ) := ln ⁡ [ tan ⁡ π ⁢ ( s + t ) 4 ⁢ T ⁢ tan ⁡ π ⁢ | t - s | 4 ⁢ T | s - t | ] ,

which does not have any singularity in the domain of integration, and the terms

J k , ℓ 0 , sing ⁢ [ m ] = { ln ⁡ h ℓ ⁢ ( ψ m p ℓ ⁢ ( 1 ) - ψ m p ℓ ⁢ ( 0 ) ) + ∫ 0 1 ∂ t ⁡ ψ m p ℓ ⁢ ( η ) ⁢ ln ⁡ η ⁢ d ⁢ η , k = ℓ , ln ⁡ h ℓ ⁢ ( ψ m p ℓ ⁢ ( 1 ) - ψ m p ℓ ⁢ ( 0 ) ) + ∫ 0 1 ∂ t ⁡ ψ m p ℓ ⁢ ( 1 - η ) ⁢ ln ⁡ η ⁢ d ⁢ η , k = ℓ + 1 , ∫ 0 1 ∂ t ⁡ ψ m p ℓ ⁢ ( η ) ⁢ ln ⁡ ( h k + η ⁢ h ℓ ) ⁢ d η , k + 1 = ℓ , ∫ 0 1 ∂ t ⁡ ψ m p ℓ ⁢ ( η ) ⁢ ln ⁡ | t k - 1 - ( t ℓ - 1 + η ⁢ h ℓ ) | ⁢ d η otherwise ,

and

J k , ℓ 1 , sing ⁢ [ m ] = { ln ⁡ h ℓ ⁢ ( ψ m p ℓ ⁢ ( 1 ) - ψ m p ℓ ⁢ ( 0 ) ) + ∫ 0 1 ∂ t ⁡ ψ m p ℓ ⁢ ( 1 - η ) ⁢ ln ⁡ η ⁢ d ⁢ η , k = ℓ , ∫ 0 1 ∂ t ⁡ ψ m p ℓ ⁢ ( η ) ⁢ ln ⁡ ( h k + ( 1 - η ) ⁢ h ℓ ) ⁢ d η , k = ℓ + 1 , ln ⁡ h ℓ ⁢ ( ψ m p ℓ ⁢ ( 1 ) - ψ m p ℓ ⁢ ( 0 ) ) + ∫ 0 1 ∂ t ⁡ ψ m p ℓ ⁢ ( η ) ⁢ ln ⁡ η ⁢ d ⁢ η , k + 1 = ℓ , ∫ 0 1 ∂ t ⁡ ψ m p ℓ ⁢ ( η ) ⁢ ln ⁡ | t k - ( t ℓ - 1 + η ⁢ h ℓ ) | ⁢ d η otherwise .

Again, we apply quadratures (3.1), (3.3) to approximate J k , ℓ 0 ⁢ [ m ] and J k , ℓ 1 ⁢ [ m ] .

4.4 Exponential Convergence of the Proposed Quadrature Schemes

In this subsection, we summarize the quality of the quadrature schemes proposed in Subsection 4.1, 4.2, 4.3.

Theorem 1

Let the mesh (1.6) fulfill the assumption max ℓ ⁡ h ℓ ≤ T / 2 . Further, choose all integration orders of the Gauss–Jacobi quadratures (3.3), (3.4) such that the related integrals in Subsections 4.1, 4.2, 4.3 are calculated in an exact way. Then all Gauss–Legendre quadratures (3.1), (3.2) applied in Subsections 4.1, 4.2, 4.3 converge exponentially with respect to the number of Gauss integration nodes. In other words, the entries of the matrices M H T , A H T , B H T are computable to high float point accuracy.

Proof

We apply the same arguments as in [16, Corollary 3.1]. ∎

5 Numerical Examples

In this section, we give numerical examples for the assembling of the matrices M H T , A H T , B H T , given in (1.14), (1.15), (1.16), respectively. Further, we give a numerical example for the ordinary differential equation (1.1) related to the heat equation, using an ℎ- and h ⁢ p -FEM as discretization schemes.

For this purpose, we fix the choice of the shape functions. In this section, we use the Lobatto polynomials (or integrated Legendre polynomials) as hierarchical shape functions on the reference element [ 0 , 1 ] . In greater detail, for a given local polynomial degree p ∈ N , we set

ψ 1 p ⁢ ( ξ ) = 1 - ξ , ψ 2 p ⁢ ( ξ ) = ξ , ψ m p ⁢ ( ξ ) = ∫ 0 ξ L m - 2 ⁢ ( ζ ) ⁢ d ζ for ⁢ m ≥ 3 , ξ ∈ [ 0 , 1 ] ,

where L m denotes the 𝑚-th Legendre polynomial on [ 0 , 1 ] , see [13, Subsection 3.1.4].

5.1 Numerical Integration

In this subsection, we show a numerical example for the quadrature schemes presented in Section 4 to calculate the entries of the matrices M H T , A H T , B H T . For this purpose, for T = 10 , we fix the non-uniform mesh by t 0 = 0 , t ℓ = 2 ℓ - N ⁢ T for ℓ = 1 , … , N with N = 6 and the polynomial degree vector by p = ( 2 , 2 , 2 , 2 , 2 , 2 ) , i.e., the degree vector is uniform. Thus, the number of the degrees of freedom is M = 13 . In this situation, we are able to calculate the matrices M H T , A H T , B H T as proposed in [19, Subsection 2.2], which serve as reference values for these matrix entries. We compare these reference values with the values, calculated as proposed in Section 4, when increasing the number of integration nodes K = K 1 = K 2 of the Gauss–Legendre quadratures (3.1), (3.2). Here, we choose the orders of the Gauss–Jacobi quadratures (3.3), (3.4) such that all related integrals are computed in an exact way, i.e., no additional approximation occurs. The errors are measured in ∥ ⋅ ∥ max defined by ∥ A ∥ max := max i , j ⁡ | A ⁢ [ i , j ] | for a matrix 𝐴 with entries A ⁢ [ i , j ] .

In Figure 1, we observe exponential convergence with respect to the number of Gauss integration nodes 𝐾 for the matrix entries of M H T , A H T , B H T calculated as proposed in Section 4. Note that this exponential convergence is in accordance with Theorem 1.

$Figure 1 Numerical results of the quadrature schemes of Section 4 with 𝐾 Gauss–Legendre points per coordinate direction with a non-uniform mesh for the matrices M H T M^{\mathcal{H}_{T}} , A H T A^{\mathcal{H}_{T}} , B H T B^{\mathcal{H}_{T}} for piecewise quadratic functions.$

Figure 1

Numerical results of the quadrature schemes of Section 4 with 𝐾 Gauss–Legendre points per coordinate direction with a non-uniform mesh for the matrices M H T , A H T , B H T for piecewise quadratic functions.

5.2 ODE Related to the Heat Equation

In this subsection, we give a numerical example for the ordinary differential equation (1.1), which is related to the heat equation. We use the discrete variational formulation (1.9) for different choices of the meshes T N and polynomial degrees 𝒑.

First, we apply an ℎ-FEM. For given N ∈ N , N ≥ 2 , we consider the uniform mesh defined by t ℓ = T ⁢ ℓ N , ℓ = 0 , … , N , i.e., a uniform refinement strategy is used when 𝑁 doubles. The vector of the polynomial degrees is uniform, i.e., p = ( p , … , p ) with p = 2 or p = 10 .

Second, we apply an h ⁢ p -FEM. For given N ∈ N , N ≥ 2 , the geometric mesh is defined by t 0 = 0 , t ℓ = T ⁢ σ N - ℓ for ℓ ∈ { 1 , … , N } with the grading parameter σ = 0.17 , see [13, p. 96]. The vector of the polynomial degrees p = ( p 1 , … , p N ) is chosen as p ℓ = ℓ for ℓ ∈ { 1 , … , N } .

Further, we set K = K 1 = K 2 = 20 for the number of integration nodes of the Gauss–Legendre quadratures (3.1), (3.2). In addition, we choose the orders of the Gauss–Jacobi quadratures (3.3), (3.4) such that all related integrals are computed in an exact way.

Last, we measure the error in the norm ∥ ⋅ ∥ H 0 , 1 / 2 ⁢ ( 0 , T ) , which is hardly computable. Thus, we use the approximation

[ v ] H 0 , 1 / 2 ⁢ ( 0 , T ) := ∥ v ∥ L 2 ⁢ ( 0 , T ) ⁢ ∥ ∂ t ⁡ v ∥ L 2 ⁢ ( 0 , T ) , v ∈ H 0 , 1 ⁢ ( 0 , T ) ,

which is an upper bound for C ⁢ ∥ ⋅ ∥ H 0 , 1 / 2 ⁢ ( 0 , T ) with C > 0 , due to the interpolation inequality, see [8, p. 23].

In Figure 2, we state the numerical results of the discrete variational formulation (1.9) for the exact solution u ⁢ ( t ) = t 3 / 4 , t ∈ [ 0 , T ] for T = 1 and μ = 10 . For the ℎ-FEM with p = 2 or p = 10 , we observe the reduced convergence due to the low regularity u ∈ H 5 / 4 - ε ⁢ ( 0 , T ) with ε ∈ ( 0 , 1 ) . For the h ⁢ p -FEM, exponential convergence with respect to the number of degrees of freedom 𝑀 is achieved. Note that this exponential convergence is proven in [12]. These numerical examples show that the proposed quadrature schemes in Section 4 provide the matrices M H T , A H T , B H T in high float point accuracy.

$Figure 2 Numerical results of the finite element discretization (1.9) for an ℎ-FEM with p = 1 p=1 or p = 10 p=10 , and for an h ⁢ p hp -FEM for the exact solution u ⁢ ( t ) = t 3 / 4 u(t)=t^{3/4} , t ∈ [ 0 , T ] t\in[0,T] , with T = 1 T=1 and μ = 10 \mu=10 .$

Figure 2

Numerical results of the finite element discretization (1.9) for an ℎ-FEM with p = 1 or p = 10 , and for an h ⁢ p -FEM for the exact solution u ⁢ ( t ) = t 3 / 4 , t ∈ [ 0 , T ] , with T = 1 and μ = 10 .

6 Conclusion

The matrices M H T , A H T , B H T occur for finite element discretizations of the heat or wave equation. The aim of this work was to realize these matrices for arbitrary polynomial degrees, which are used for, e.g., a temporal h ⁢ p -FEM. First, we stated weakly singular integral representations of the modified Hilbert transformation H T . These integral representations formed the basis of the quadrature schemes to calculate the matrices M H T , A H T , B H T . All occurring singular integrals were reformulated such that they can be calculated to machine precision. In the last part, exponential convergence with respect to the number of Gauss integration nodes was observed in numerical examples. Moreover, an h ⁢ p -FEM for an ordinary differential equation, which is related to the heat equation, showed the potential of the proposed quadrature schemes.

References

[1] W. Bangerth, M. Geiger and R. Rannacher, Adaptive Galerkin finite element methods for the wave equation, Comput. Methods Appl. Math. 10 (2010), no. 1, 3–48. 10.2478/cmam-2010-0001Suche in Google Scholar

[2] A. Ern and J.-L. Guermond, Finite Elements I—Approximation and interpolation, Texts Appl. Math. 72, Springer, Cham, 2021. 10.1007/978-3-030-56341-7Suche in Google Scholar

[3] W. Gautschi, Numerical integration over the square in the presence of algebraic/logarithmic singularities with an application to aerodynamics, Numer. Algorithms 61 (2012), no. 2, 275–290. 10.1007/s11075-012-9611-9Suche in Google Scholar

[4] W. Gautschi, Erratum to: “Numerical integration over the square in the presence of algebraic/logarithmic singularities with an application to aerodynamics”, Numer. Algorithms 64 (2013), no. 4, 759–759. 10.1007/s11075-013-9787-7Suche in Google Scholar

[5] O. A. Ladyzhenskaya, The Boundary Value Problems of Mathematical Physics, Appl. Math. Sci. 49, Springer, New York, 1985. 10.1007/978-1-4757-4317-3Suche in Google Scholar

[6] U. Langer and O. Steinbach, Space-Time Methods—Applications to Partial Differential Equations, Radon Ser. Comput. Appl. Math.25, De Gruyter, Berlin, 2019. 10.1515/9783110548488Suche in Google Scholar

[7] U. Langer and M. Zank, Efficient direct space-time finite element solvers for parabolic initial-boundary value problems in anisotropic Sobolev spaces, SIAM J. Sci. Comput. 43 (2021), no. 4, A2714–A2736. 10.1137/20M1358128Suche in Google Scholar

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

[9] R. Löscher, O. Steinbach and M. Zank, Numerical results for an unconditionally stable space-time finite element method for the wave equation. Vol. 26, Domain Decomposition Methods in Science and Engineering, Lect. Notes Comput. Sci. Eng. 145, Springer, Cham, (2022), 587–594. Suche in Google Scholar

[10] R. Löscher, O. Steinbach and M. Zank, An unconditionally stable space-time finite element method for the wave equation, in preparation, 2022. 10.1007/978-3-030-95025-5_68Suche in Google Scholar

[11] J. Nečas, Les méthodes directes en théorie des équations elliptiques, Masson et Cie, Paris, 1967. Suche in Google Scholar

[12] I. Perugia, C. Schwab and M. Zank, Exponential convergence of h ⁢ p -time-stepping in space-time discretizations of parabolic PDEs, preprint (2022), https://arxiv.org/abs/2203.11879; to appear in ESAIM Math. Model. Numer. Anal. 10.1051/m2an/2022081Suche in Google Scholar

[13] C. Schwab, 𝑝- and h ⁢ p -Finite Element Methods, Numer. Math. Sci. Comput., The Clarendon, New York, 1998. Suche in Google Scholar

[14] O. Steinbach and A. Missoni, A note on a modified Hilbert transform, Appl. Anal. (2022), 10.1080/00036811.2022.2030725. 10.1080/00036811.2022.2030725Suche in Google Scholar

[15] 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_vol52s154Suche in Google Scholar

[16] O. Steinbach and M. Zank, A note on the efficient evaluation of a modified Hilbert transformation, J. Numer. Math. 29 (2021), no. 1, 47–61. Suche in Google Scholar

[17] V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, 2nd ed., Springer Ser. Comput. Math. 25, Springer, Berlin, 2006. Suche in Google Scholar

[18] M. Zank, Inf-Sup Stable Space-Time Methods for Time-Dependent Partial Differential Equations, Monographic Series TU Graz 36, TU Graz, Graz, 2020. Suche in Google Scholar

[19] M. Zank, An exact realization of a modified Hilbert transformation for space-time methods for parabolic evolution equations, Comput. Methods Appl. Math. 21 (2021), no. 2, 479–496. 10.1515/cmam-2020-0026Suche in Google Scholar

[20] M. Zank, Efficient direct space-time finite element solvers for the wave equation in second-order formulation, in preparation, 2022. Suche in Google Scholar

[21] M. Zank, High-order discretisations and efficient direct space-time finite element solvers for parabolic initial-boundary value problems, to appear in Proceedings of Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2020+1. Suche in Google Scholar

Received: 2022-07-25

Accepted: 2022-09-06

Published Online: 2022-10-06

Published in Print: 2023-04-01

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

Artikel in diesem Heft

https://doi.org/10.1515/cmam-2022-0150

Schlagwörter für diesen Artikel

Modified Hilbert Transformation; Weakly Singular Integral; Quadratures; Space-Time; Finite Element Method; Parabolic Equations; Hyperbolic Equations

Creative Commons

BY 4.0