Integral Representations and Quadrature Schemes for the Modified Hilbert Transformation

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.

03/07/2021

Numerical results for an unconditionally stable space-time finite element method for the wave equation

In this work, we introduce a new space-time variational formulation of t...
09/24/2019

A sharp error estimate of piecewise polynomial collocation for nonlocal problems with weakly singular kernels

As is well known, using piecewise linear polynomial collocation (PLC) an...
07/18/2022

Modified BDF2 schemes for subdiffusion models with a singular source term

The aim of this paper is to study the time stepping scheme for approxima...
03/02/2022

A modified convolution quadrature combined with the method of fundamental solutions and Galerkin BEM for acoustic scattering

We describe a numerical method for the solution of acoustic exterior sca...
04/10/2020

Numerical methods for stochastic Volterra integral equations with weakly singular kernels

In this paper, we first establish the existence, uniqueness and Hölder c...
08/20/2022

A Modified Trapezoidal Rule for a Class of Weakly Singular Integrals in n Dimensions

In this paper we propose and analyze a general arbitrarily high-order mo...
07/28/2017

Two Hilbert schemes in computer vision

We study the multiview moduli problems that arise in computer vision. We...

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 , 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 , are an -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 , 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 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 afore-mentioned -FEM. In addition, for discretizations of hyperbolic problems as in [9, 10], the weakly singular integral representation of the modified Hilbert transformation , 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 for high polynomial degrees, as needed in an -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 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 , and

, we consider the ordinary differential equation (ODE)

 ∂tu(t)+μu(t)=f(t) for t∈(0,T),u(0)=0, (1)

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

 ∂ttu(t)+μu(t)=f(t) for t∈(0,T),u(0)=∂tu(0)=0, (2)

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), (2), we introduce classical Sobolev spaces, where details and further references are contained in [18]. The classical Sobolev space is endowed with the norm . Its closed subspaces

 H10,(0,T)={v∈H1(0,T):v(0)=0} and H1,0(0,T)={v∈H1(0,T):v(T)=0}

are endowed with the Hilbertian norm where is the usual norm and denotes the (weak) derivative. Additionally, we have the Sobolev spaces

 H1/20,(0,T) :={U|(0,T):U∈H1/2(−∞,T) with U(t)=0,t<0}, H1/2,0(0,T) :={U|(0,T):U∈H1/2(0,∞) with U(t)=0,t>T}

with the Hilbertian norms

 ∥v∥H1/2,0(0,T):=(∥v∥2H1/2(0,T)+∫T0|v(t)|2T−tdt)1/2,

where is a norm, e.g., the Slobodetskii norm, in the usual Sobolev space for open intervals . Further, denotes the duality pairing in and , or in and , as extension of the inner product in Last, we define the modified Hilbert transformation as

 (HTv)(t)=∞∑k=0vkcos((π2+kπ)tT),t∈(0,T), (3)

where

 vk=2T∫T0v(t)sin((π2+kπ)tT)dt

are the Fourier coefficients of the series representation

 v(t)=∞∑k=0vksin((π2+kπ)tT)dt,t∈(0,T),

when is given. The modified Hilbert transformation is a bijective mapping for with . Additional properties of are given in [15, 16, 19, 14].

With this notation, a related variational formulation to the parabolic-type problem (1) is to find such that

 ∀v∈H1/20,(0,T):⟨∂tu,HTv⟩(0,T)+μ⟨u,HTv⟩L2(0,T)=⟨f,HTv⟩(0,T) (4)

for a given right-hand side . For the hyperbolic-type problem (2), a related variational formulation is to find such that

 ∀v∈H10,(0,T):⟨HT∂tu,∂tv⟩L2(0,T)+μ⟨u,HTv⟩L2(0,T)=⟨f,HTv⟩(0,T) (5)

for a given right-hand side . The variational formulations (4), (5) are uniquely solvable, where their derivation and analysis is given in [15, 18, 9].

Next, we state conforming discretizations of the variational formulations (4), (5). For this purpose, we define a partition of for given , i.e.,

 0=t0

with elements . The mesh sizes are for , and the maximal mesh size is . On , with the distribution of polynomial degrees, we introduce the space of piecewise polynomial, continuous functions on intervals

 Sp(TN):={v∈C[0,T]:∀ℓ∈{1,…,N}:v|τℓ∈Ppℓ(τℓ)}=span{φi}Mi=1, (7)

where is the space of continuous functions and is the space of polynomials on a subset , , of global degree at most . Here,

is the number of degrees of freedom and the functions

are basis functions of , satisfying

 ∀j∈{2,…,M}:φj(0)=0, (8)

i.e., Further, we define the subspace

 Sp0,(TN):=Sp(TN)∩H10,(0,T)=span{φj+1}M−1j=1.

Thus, a conforming finite element method of the variational formulation (4) is to find such that

 ∀vh∈Sp0,(TN):⟨∂tuh,HTvh⟩L2(0,T)+μ⟨uh,HTvh⟩L2(0,T)=⟨f,HTvh⟩(0,T). (9)

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

 (˜AHT+μ˜MHT)–u=F––HT (10)

with the matrices

 ˜MHT[i,j]:= ⟨φj+1,HTφi+1⟩L2(0,T), (11) ˜AHT[i,j]:= ⟨∂tφj+1,HTφi+1⟩L2(0,T) (12)

for and the corresponding right-hand side .

Analogously, the conforming finite element method of the variational formulation (5) to find such that

 ∀vh∈Sp0,(TN):⟨HT∂tuh,∂tvh⟩L2(0,T)+μ⟨uh,HTvh⟩L2(0,T)=⟨f,HTvh⟩(0,T)

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

 ((˜BHT)⊤+μ˜MHT)u––=F––HT

with the matrix , given in (11), and the matrix

 ˜BHT[i,j]:=⟨∂tφj+1,HT∂tφi+1⟩L2(0,T) (13)

for and the right-hand side as for the linear system (10). Note that for , this right-hand side can be realized, using the matrix

 MHT[i,j]:=⟨φj,HTφi⟩L2(0,T) (14)

for , see [19, 21] for all the details. Additionally, we set

 AHT[i,j]:= ⟨∂tφj,HTφi⟩L2(0,T), (15) BHT[i,j]:= ⟨∂tφj,HT∂tφi⟩L2(0,T) (16)

for . Hence, the matrices , , , given in (11), (12), (13), are submatrices of the matrices , , , given in (14), (15), (16), respectively. Thus, we ask for an accurate realization of these matrices , , , 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 , which are the starting point of the calculation of the matrices , , . In Section 3, Gauß 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 , , to high accuracy. In Section 5, numerical examples show the quality of the new assembling method of , , . 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 for functions in and prove a new weakly singular integral representation of for functions, which are only piecewise in . These weakly singular integral representations are used for the calculation of the matrices , , , given in (14), (15), (16), respectively. We start to recall the integral representations for functions in or , which is proven in [16]:

Lemma 1 (Lemma 2.1 in [16]).

For the operator , as defined in (3), allows the integral representation

 (HTv)(t)=v.p.∫T0v(s)K(s,t)ds,t∈(0,T), (17)

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

 K(s,t):=12T⎡⎢⎣1sinπ(s+t)2T+1sinπ(s−t)2T⎤⎥⎦.

For the operator , as defined in (3), allows the integral representation

 (HTv)(t)=−2πv(0)lntanπt4T+∫T0∂tv(s)K(s,t)ds,t∈(0,T), (18)

as a weakly singular integral, where

 K(s,t):=−1πln[tanπ(s+t)4Ttanπ|t−s|4T]. (19)

The integral representation (18) in connection with specialized numerical integration gives the possibility to calculate matrices, which are based on applied to a function in , e.g., the matrices in (14) and in (15). However, the integral representation (18) is not applicable for functions, which have jumps, since these functions do not belong to . As the integral representation (17) is a Cauchy principal value integral and therefore, is complicated to use in practical implementations, we prove a new integral representation of for functions, which are piecewise in but possibly having jumps. This integral representation is used for the calculation of the matrix in (16).

Lemma 2.

Let with and a function be given. Then, for the function

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

the operator , as defined in (3), allows the integral representation

 (HTvf)(t)=−f(b)K(b,t)+f(a)K(a,t)+∫ba∂tf(s)K(s,t)ds (20)

for as a weakly singular integral, where the kernel function is given in (19).

If, in addition, , i.e., , then the representation

 (HTvf)(a)=−f(b)K(b,a)+∫ba∂tf(s)K(s,a)ds

holds true as a weakly singular integral. Analogously, if , i.e., , then we have the representation

 (HTvf)(b)=f(a)K(a,b)+∫ba∂tf(s)K(s,b)ds

as a weakly singular integral.

Proof.

Let be arbitrary but fixed. We split the proof accordingly to the relation between , and .

First, consider the case . The integral representation (17), integration by parts, the relation and the Hölder continuity of with exponent , see [11, Chapitre 2, Théorème 3.8], yield

 (HTvf)(t) =limε↘0(∫t−εaf(s)K(s,t)ds+∫bt+εf(s)K(s,t)ds) =limε↘0(−f(t−ε)K(t−ε,t)+f(a)K(a,t)+∫t−εa∂tf(s)K(s,t)ds =f(a)K(a,t)−f(b)K(b,t)+∫ba∂tf(s)K(s,t)ds,

where the last integral exists as a weakly singular integral.

Second, we examine the case when , or when . The integral representation (17) and integration by parts give

 (HTvf)(t)=limε↘0∫baf(s)K(s,t)ds=f(a)K(a,t)−f(b)K(b,t)+∫ba∂tf(s)K(s,t)ds.

Third, the case is investigated when . As in the cases before, the integral representation (17), the Hölder continuity of and integration by parts lead to

 (HTvf)(a) =limε↘0∫ba+εf(s)K(s,a)ds =limε↘0(f(a+ε)K(a+ε,a)−f(b)K(b,a)+∫ba∂tf(s)K(s,a)ds) =−f(b)K(b,a)+∫ba∂tf(s)K(s,a)ds,

where the last integral exists as a weakly singular integral.

Last, the case , when , is proven analogously. ∎

Remark 1.

Choose and in Lemma 2. Then, the representation (20) of Lemma 2 coincides with the representation (18) of Lemma 1 due to and for

Remark 2.

Consider the situation of Lemma 2. For with or , the function admits a singularity for or , respectively.

In this section, Gauß quadratures are stated, which are used in the following sections. For this purpose, we introduce Gauß quadratures on of order , see [3]. First, the classical Gauß–Legendre quadrature with weights and nodes , , fulfills

 ∀g∈P2K−1[0,1]:∫10g(t)dt=K∑ν=1ων,Kg(ξν,K), (21)

see [3, (1.6)]. This Gauß quadrature is generalized to the square domain via using the Gauß–Legendre quadrature (21

) for each coordinate direction. More precisely, the tensor Gauß quadrature

 K1∑ν1=1K2∑ν2=1ων1,K1ων2,K2G(ξν1,K1,ξν2,K2), (22)

approximates the integral

 ∫10∫10G(ξ1,ξ2)dξ1dξ2

for a function , where and are the Gauß integration weights and Gauß integration nodes of order with respect to the coordinate direction .

For singular integrals with a logarithmic term, we introduce a nonclassical Gauß–Jacobi quadrature. In greater detail, for , the weights and nodes , , of the Gauß–Jacobi quadrature for a logarithmic term satisfy

 ∀g∈P2K−1[0,1]:−∫10g(t)ln(t)dt=K∑ν=1^ων,Kg(^ξν,K), (23)

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

 ∫10∫10g(s,t)ln|s−t|dsdt=−K∑μ=1K∑ν=1^ωμ^ξμων(g((1−ξν)^ξμ,^ξμ)+g(1−(1−ξν)^ξμ,1−^ξμ))−K∑μ=1K∑ν=1ωμξμ^ων(g((1−^ξν)ξμ,ξμ)+g(1−(1−^ξν)ξμ,1−ξμ)) (24)

holds true for all , where

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

are defined in (21), (23), see [3, (2.7)], [4].

4 Quadrature Schemes for the Modified Hilbert Transformation

In this section, we describe the assembling of the matrices in (14), in (15) and in (16). The crucial point is the realization of the modified Hilbert transformation , where different possibilities exist, see [16, 19]

. In particular, for a uniform degree vector

with a fixed, low polynomial degree , e.g., or , the matrices , and in (14), (15), (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 -FEM, the degree vector is not uniform or for high polynomial degrees, it is convenient to apply numerical quadrature rules to numerically approximate the matrix entries of , and in (14), (15), (16), respectively. This is the main topic of the paper and is described in great detail in the following.

From the integral representation of in (18), we have

 MHT[i,j]=⟨φj,HTφi⟩L2(0,T)=−2πφi(0)∫T0φj(t)lntanπt4Tdt+∫T0φj(t)∫T0K(s,t)∂tφi(s)dsdt (25)

and

 (26)

for with the basis functions in (7) satisfying (8).

For the matrix in (16), the integral representation of in (20) yields

 HT∂tφi(t)=N∑k=1[−(∂tφi)k−K(tk,t)+(∂tφi)k−1+K(tk−1,t)+∫tktk−1∂tt(φi|τk)(s)K(s,t)ds]

for , where we use the notation

 vk−:=limε↘0v(tk−ε) and vk−1+:=limε↘0v(tk−1+ε),k=1,…,N,

for a sufficiently smooth function Thus, the entries of the matrix in (16) admit the representation

 BHT[i,j]= ⟨∂tφj,HT∂tφi⟩L2(0,T) = N∑k=1[−(∂tφi)k−∫T0∂tφj(t)K(tk,t)dt+(∂tφi)k−1+∫T0∂tφj(t)K(tk−1,t)dt (27)

for with the basis functions in (7) satisfying (8).

The matrix entries , , in (25), (26), (27), respectively, are computed element-wise for the partition of into intervals , . For this purpose, we assume that we use standard Lagrange finite elements as in [2, Section 6.3], or for the -FEM, we use basis functions , , as described in [13, Subsection 3.1.6]. Thus, we assume that the basis functions , , allow for a local representation on an element , by shape functions , , defined on the reference interval In greater detail, for a basis function , , there exists an element , , with a related polynomial degree such that the local representation

 ∀t∈¯¯¯τℓ:φα(m,ℓ)(t)=ψpℓm(t−tℓ−1hℓ) (28)

holds true, where is the global index related to the local index for the interval and is a shape function. Possible choices of such shape functions are the classical Lagrange polynomials with equidistant nodes or Gauß–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 , with indices and related local polynomial degrees , . We define the local matrix by

 MHTk,ℓ[n,m]= −2πφα(n,k)(0)∫tℓtℓ−1φα(m,ℓ)(t)lntanπt4Tdt +∫tℓtℓ−1φα(m,ℓ)(t)∫tktk−1K(s,t)∂tφα(n,k)(s)dsdt = −δk,1ψpkn(0)2hℓπ∫10ψpℓm(η)lntanπ(tℓ−1+ηhℓ)4Tdη +hℓ∫10ψpℓm(η)∫10K(tk−1+ξhk,tℓ−1+ηhℓ)∂tψpkn(ξ)dξdη (29)

and, analogously, the local matrix by

 AHTk,ℓ[n,m]=−δk,1ψpkn(0)2π∫10∂tψpℓm(η)lntanπ(tℓ−1+ηhℓ)4Tdη+∫10∂tψpℓm(η)∫10K(tk−1+ξhk,tℓ−1+ηhℓ)∂tψpkn(ξ)dξdη (30)

for and , where we used the property (8), the local representation (28) and is the usual Kronecker delta.

In the same way, the local matrix is defined by

 BHTk,ℓ[n,m]=−1πhk∂tψpkn(0)J0k,ℓ[m]+1πhk∂tψpkn(1)J1k,ℓ[m]+1hk∫10∂tψpℓm(η)