1 Introduction↩︎

In the last years, there has been a growing interest in simultaneous space-time boundary element methods (BEM) for the heat equation [1]–[14]. In contrast to the differential-operator-based variational formulation on the space-time cylinder, the variational formulation corresponding to space-time BEM is coercive [15], [16] so that the discretized version always has a unique solution regardless of the chosen trial space, which is even quasi-optimal in the natural energy norm. Moreover, it is naturally applicable on unbounded domains and only requires a mesh of the lateral boundary of the space-time cylinder (as well as a mesh of the spatial domain in case of nonhomogeneous initial data) resulting in a dimension reduction. The potential disadvantage that discretizations lead to dense matrices due to the nonlocality of the boundary integral operators has been tackled, e.g., in [2]–[4], [10], [14], [17]–[19] via wavelets, the fast multipole method, and \(\mathcal{H}\)-matrices.

Two often mentioned advantages of simultaneous space-time methods are their potential for massive parallelization as well as their potential for fully adaptive refinement to resolve singularities local in both space and time. While the first advantage has been investigated in, e.g., [7], [9], the latter requires suitable a posteriori computable error estimators, which are so far only available for the heat equation with Dirichlet boundary conditions [20], [21]. In the latter work, the Faermann error estimator [22], [23] for the Laplace problem is generalized.

In the present manuscript, we extend the weighted-residual estimators [24]–[26] for the Laplace problem to the heat equation with Neumann boundary conditions: Let \(\Omega\subset\mathbb{R}^{d}\), \(d=2,3\), be a Lipschitz domain with connected boundary \(\Gamma:=\partial\Omega\) and \(T>0\) a given end time point with corresponding time interval \(I:=(0,T)\). While extensions to piecewise smooth curved boundaries are possible, for the ease of presentation, we assume that \(\Omega\) is polygonal for \(d=2\) and polyhedral for \(d=3\). We abbreviate the space-time cylinder \(Q:=I\times\Omega\) with lateral boundary \(\Sigma:=I\times\Gamma\) and corresponding outer normal vector \({\boldsymbol{n}}\in\mathbb{R}^{d}\). With the heat kernel \[\begin{align} G(t,{\boldsymbol{x}}) := \begin{cases} \frac{1}{(4\pi t)^{d/2}} \, e^{-\frac{|{\boldsymbol{x}}|^2}{4t}} \quad &\text{for }(t,{\boldsymbol{x}})\in (0,\infty)\times\mathbb{R}^{d},\\ 0 \quad &\text{else}, \end{cases} \end{align}\] and a given function \(f:\Sigma\to\mathbb{R}\), we consider the boundary integral equation \[\begin{align} \label{eq:hypersing32operator} (\mathscr{W} u)(t,{\boldsymbol{x}}) := -\partial_{{\boldsymbol{n}}({\boldsymbol{x}})} \int_\Sigma \partial_{{\boldsymbol{n}}({\boldsymbol{y}})}G(t-s,{\boldsymbol{x}}-{\boldsymbol{y}}) u(s,{\boldsymbol{y}}) \,{\rm d}{\boldsymbol{y}}\,{\rm d}s = f(t,{\boldsymbol{x}}) \quad\text{for }(t,{\boldsymbol{x}})\in \Sigma. \end{align}\tag{1}\] Here, \(\mathscr{W}\) is the hyper-singular operator. Note that the integral kernel is given as \[\label{eq:partialnyG} \partial_{{\boldsymbol{n}}({\boldsymbol{y}})}G(t-s,{\boldsymbol{x}}-{\boldsymbol{y}})=({\boldsymbol{x}}-{\boldsymbol{y}})\cdot {\boldsymbol{n}}({\boldsymbol{y}})\frac{1}{2(t-s)}G(t-s,{\boldsymbol{x}}-{\boldsymbol{y}}).\tag{2}\] For given initial condition \(U_0:\Omega\to\mathbb{R}\) and Neumann datum \(\phi:\Sigma\to\mathbb{R}\), such equations arise from the heat equation \[\begin{align} \label{eq:interior} \begin{array}{rcll} \partial_t U - \Delta_{\boldsymbol{x}}U & = & 0 & \text{ on } Q,\\ \partial_{\boldsymbol{n}}U & = & \phi & \text{ on }\Sigma,\\ U(0,\cdot) & = & U_0 & \text{ on }\Omega; \end{array} \end{align}\tag{3}\] see Section 2.3.

Let \(\mathcal{P}_h\) be a mesh of the space-time boundary \(\Sigma\) consisting of prismatic elements \(J\times K\) with intervals \(J\subseteq \overline{I}\) and simplices \(K\subseteq \Gamma\), and \(\widetilde{S}^{p_t,p_{\boldsymbol{x}}}(\mathcal{P}_h)\) the associated space of continuous piecewise polynomials of some fixed degree in time and space that vanish at \(t=0\). Further, let \(u_h\in\widetilde{S}^{p_t,p_{\boldsymbol{x}}}(\mathcal{P}_h)\) be the Galerkin approximation of \(u\) with respect to the coercive bilinear form induced by \(\mathscr{W}\). As \(\mathscr{W}\) is an isomorphism from the anisotropic Sobolev space \(H^{1/2,1/4}(\Sigma):= L^2 (I; H^{1/2}(\Gamma)) \cap H^{1/4}(I; L^2(\Gamma))\) to its dual space \(H^{-1/2,-1/4}(\Sigma):=H^{1/2,1/4}(\Sigma)'\), the discretization error \(\|u-u_h\|_{H^{1/2,1/4}(\Sigma)}\) is equivalent to the norm of the residual \(\|f-\mathscr{W}u_h\|_{H^{-1/2,-1/4}(\Sigma)}\). Under certain regularity assumptions on the mesh and the local parabolic scaling \(|J|\eqsim\mathop{\mathrm{diam}}(K)^2\) for all \(J\times K\in\mathcal{P}_h\), we show that \[\label{eq:intreliability} \|u-u_h\|_{H^{1/2,1/4}(\Sigma)}^2 \eqsim \|f-\mathscr{W}u_h\|_{H^{-1/2,-1/4}(\Sigma)}^2 \lesssim \sum_{J\times K\in\mathcal{P}_h} \mathop{\mathrm{diam}}(K) \|f-\mathscr{W}u_h\|_{L^2(J\times K)}^2.\tag{4}\] Assuming \(f\in L^2(\Sigma)\), the right-hand side only makes sense if \(\mathscr{W}u_h\in L^2(\Sigma)\), which follows from the regularity \(\mathscr{W}:\widetilde{H}^{1,1/2}(\Sigma) \to L^2(\Sigma)\) in conjunction with the inclusion \(\widetilde{S}^{p_t,p_{\boldsymbol{x}}}(\mathcal{P}_h)\subset \widetilde{H}^{1,1/2}(\Sigma):=\big\{v \in L^2(I;H^1(\Gamma)) \cap H^{1/2}(I;L^2(\Gamma))\,:\,v(0,\cdot) = 0\big\}\). Such a regularity result is valid for piecewise smooth (including in particular polygonal and polyhedral) domains \(\Omega\); see [16].

Building on the derived a posteriori computable error estimator 4 , we further propose an adaptive algorithm and numerically investigate it. In all considered experiments (for \(d=2\)), reliability of the estimator 4 is empirically confirmed and also efficiency, i.e., the converse inequality, is observed. Compared to uniform refinement, the adaptive algorithm yields significantly faster convergence, achieving near-optimal rates with respect to the number of degrees of freedom, even in the presence of strong singularities. We provide full details on the numerical implementation of the proposed algorithm, including the computation of the Galerkin system and the error estimator. In particular, for arbitrary polynomial degree \(p_t\) in time, we give explicit formulas for the time integrals in the boundary integral operators applied to and tested with piecewise polynomials of that degree, which are so far only available for \(p_t = 0\); see, e.g., [9], [16], [27].

Outline↩︎

The remainder of this work is organized as follows: Section 2 summarizes the general principles of the space-time boundary element method for the heat equation, including prismatic boundary meshes and piecewise polynomial ansatz spaces along with a basis for the lowest-order case \(p_t = p_x = 1\). Section 3 establishes Poincaré-type inequalities for anisotropic Sobolev spaces (Section 3.1) and constructs a Clément-type interpolation operator (Section 3.2). The latter is used to derive our main result, reliability of a weighted-residual estimator (Theorem 3.6). Finally, Section 4 introduces an adaptive algorithm which is based on the derived error estimator. This algorithm is subsequently applied for \(d=2\) to several concrete examples with typical singularities in space and time. The stable implementation is discussed in Appendix A.

2 Preliminaries↩︎

2.1 General notation↩︎

Throughout and without any ambiguity, \(|\cdot|\) denotes the absolute value of scalars, the Euclidean norm of vectors in \(\mathbb{R}^n\), or the the measure of a set in \(\mathbb{R}^n\), e.g., the length of an interval or the area of a surface in \(\mathbb{R}^3\). We write \(A\lesssim B\) to abbreviate \(A\le CB\) with some generic constant \(C>0\), which is clear from the context. Moreover, \(A\eqsim B\) abbreviates \(A\lesssim B\lesssim A\).

2.2 Anisotropic Sobolev spaces↩︎

For a measurable \(d\)-dimensional \(\omega \subseteq \Omega\) or \((d-1)\)-dimensional \(\omega \subseteq \Gamma\), and \(\mu\in(0,1]\), we first recall the Sobolev space \[\begin{align} H^{\mu}(\omega):=\big\{v\in L^2(\omega)\,:\,\|v\|_{H^{\mu}(\omega)} < \infty\big\} \end{align}\] associated with the Sobolev–Slobodeckij norm \[\begin{align} \|v\|_{H^{\mu}(\omega)}^2 := \|v\|_{L^2(\omega)}^2 + |v|_{H^{\mu}(\omega)}^2, \quad |v|_{H^{\mu}(\omega)}^2 :=\begin{cases} \int_\omega\int_\omega \frac{|v({\boldsymbol{x}})-v({\boldsymbol{y}})|^2}{|{\boldsymbol{x}}-{\boldsymbol{y}}|^{{\rm dim}(\omega)+2\mu}}\,{\rm d}{\boldsymbol{y}}\,{\rm d}{\boldsymbol{x}}&\text{ if }\mu\in(0,1), \\ \|\nabla_\omega v\|_{L^2(\omega)}^2 \quad&\text{ if }\mu=1, \end{cases} \end{align}\] where \({\rm dim}(\omega)\) denotes the dimension of \(\omega\), i.e., \(d\) or \(d-1\), and \(\nabla_\omega\) denotes the (weak) gradient on \(\omega\), i.e., the standard gradient or the surface gradient.

Moreover, we define for any subinterval \(J\subseteq\overline{I}\), \(\nu\in(0,1]\), and any Banach space \(X\), \[\begin{align} H^\nu(J;X) := \big\{v\in L^2(J;X)\,:\,\|v\|_{H^{\nu}(J;X)} < \infty\big\} \end{align}\] associated with the norm \[\begin{align} \|v\|_{H^{\nu}(J;X)}^2 := \|v\|_{L^2(J;X)}^2 + |v|_{H^{\nu}(J;X)}^2, \quad |v|_{H^{\nu}(J;X)}^2 :=\begin{cases} \int_J\int_J \frac{\|v(t)-v(s)\|_{X}^2}{|t-s|^{1+2\nu}}\,{\rm d}s\,{\rm d}t &\text{ if }\nu\in(0,1), \\ \|\partial_t v\|_{L^2(J;X)}^2 &\text{ if }\nu=1, \end{cases} \end{align}\] where \(\partial_t\) denotes the (weak) time derivative. If \(X=\mathbb{R}\), we simply write \(H^\nu(J)\), \(\|v\|_{H^\nu(J)}\), and \(|v|_{H^\nu(J)}\). We recall the anisotropic Sobolev space \[\begin{align} H^{\mu,\nu}(J\times\omega) := L^2(J;H^{\mu}(\omega)) \cap H^{\nu}(J;L^2(\omega)) \end{align}\] with corresponding norm \[\begin{align} \|v\|_{H^{\mu,\nu}(J\times\omega)}^2 := \|v\|_{L^2(J;H^{\mu}(\omega))}^2 + \|v\|_{H^{\nu}(J;L^2(\omega))}^2 \quad\text{for all } v\in H^{\mu,\nu}(J\times\omega). \end{align}\] We will sometimes use the abbreviation \[\begin{align} |v|_{L^2(J;H^{\mu}(\omega))}^2 := \int_J |v(t,\cdot)|_{H^\mu(\omega)}^2 \,{\rm d}t \quad\text{for all } v\in L^2(J;H^\mu(\omega)). \end{align}\]

For \(\omega\in\{\Omega,\Gamma\}\), we denote by \(H^{-\mu,-\nu}(I\times\omega)\) the dual space of \(H^{\mu,\nu}(I\times\omega)\) with duality pairing \(\langle\cdot\,,\,\cdot\rangle_{I\times\omega}\). We interpret \(L^2(I\times\omega)\) as subspace of \(H^{-\mu,-\nu}(I\times\omega)\) via \[\begin{align} \langle v\,,\,\psi\rangle_{I\times\omega} := \int_I\int_\omega v(t,{\boldsymbol{x}}) \psi(t,{\boldsymbol{x}}) \,{\rm d}{\boldsymbol{x}}\,{\rm d}t \quad\text{for all }v\in H^{\mu,\nu}(I\times\omega) \text{ and }\psi\in L^2(I\times\omega). \end{align}\] We further recall the space \[\begin{align} \widetilde{H}^{\mu,\nu}(I\times\omega) := \big\{v\in H^{\mu,\nu}(I\times\omega)\,:\,v(0,\cdot)=0\big\}, \end{align}\] which coincides with \(H^{\mu,\nu}(I\times\omega)\) for \(\nu<1/2\), and \[\begin{align} H^{1,1/2}{(Q;\partial_t-\Delta_{\boldsymbol{x}})}:=\big\{v\in H^{1,1/2}(Q)\,:\,\partial_tv-\Delta_{\boldsymbol{x}}v\in L^2(Q)\big\}, \end{align}\] which is equipped with the canonical graph norm. For more details on the defined anisotropic Sobolev spaces, we refer, e.g., to [16].

2.3 Boundary integral equations↩︎

For given initial condition \(U_0:\Omega\to\mathbb{R}\) and Neumann datum \(\phi:\Sigma\to\mathbb{R}\), we consider the heat equation 3 . It is well-known that for \(U_0\in L^2(\Omega)\) and \(\phi\in H^{-1/2,-1/4}(\Sigma)\), the heat equation 3 admits a unique solution \(u\in H^{1,1/2}(Q)\). With the lateral trace \(u:=U|_\Sigma \in H^{1/2,1/4}(\Sigma)\), \(U\) satisfies the representation formula \[\begin{align} \label{eq:representation} U = \widetilde{\mathscr{M}}_0 U_0 + \widetilde{\mathscr{V}}\phi - \widetilde{\mathscr{K}}u, \end{align}\tag{5}\] where \[\begin{align} \tag{6} (\widetilde{\mathscr{M}}_0 U_0)(t,{\boldsymbol{x}})&:=\int_\Omega G(t,{\boldsymbol{x}}-{\boldsymbol{y}}) U_0({\boldsymbol{y}}) \,{\rm d}{\boldsymbol{y}}\quad\text{for all }(t,{\boldsymbol{x}})\in Q \intertext{denotes the initial potential,} \tag{7} (\widetilde{\mathscr{V}}\phi)(t,{\boldsymbol{x}})&:=\int_\Sigma G(t-s,{\boldsymbol{x}}-{\boldsymbol{y}}) \phi(s, {\boldsymbol{y}}) \,{\rm d}{\boldsymbol{y}}\,{\rm d}s \quad\text{for all }(t,{\boldsymbol{x}})\in Q \intertext{denotes the single-layer potential, and } \tag{8} (\widetilde{\mathscr{K}}u)(t,{\boldsymbol{x}})&:=\int_\Sigma \partial_{{\boldsymbol{n}}({\boldsymbol{y}})}G(t-s,{\boldsymbol{x}}-{\boldsymbol{y}}) u(s, {\boldsymbol{y}}) \,{\rm d}{\boldsymbol{y}}\,{\rm d}s \quad\text{for all }(t,{\boldsymbol{x}})\in Q \end{align}\] denotes the double-layer potential. These linear operators satisfy the mapping properties \(\widetilde{\mathscr{M}}_0:L^2(\Omega)\to H^{1,1/2}(Q;\partial_t-\Delta_{\boldsymbol{x}})\), \(\widetilde{\mathscr{V}}:H^{-1/2,-1/4}(\Sigma)\to H^{1,1/2}(Q;\partial_t-\Delta_{\boldsymbol{x}})\), and \(\widetilde{\mathscr{K}}:H^{1/2,1/4}(\Sigma)\to H^{1,1/2}(Q;\partial_t-\Delta_{\boldsymbol{x}})\). We denote the normal derivative \(\partial_{\boldsymbol{n}}(\cdot):H^{1,1/2}(Q;\partial_t-\Delta_{\boldsymbol{x}})\to H^{-1/2,-1/4}(\Sigma)\) of these potentials as follows \[\begin{align} \partial_{\boldsymbol{n}}(\widetilde{\mathscr{M}}_0 U_0) =: \mathscr{M}_1 U_0, \quad \partial_{\boldsymbol{n}}(\widetilde{\mathscr{V}}\phi) - \phi/2 =: {\mathscr N} \phi, \quad -\partial_{\boldsymbol{n}}(\widetilde{\mathscr{K}}u) =: \mathscr{W}u . \end{align}\] Applying the normal derivative to 5 thus results in \[\begin{align} \label{eq:direct} \mathscr{W}u = (1/2-{\mathscr N}) \phi - \mathscr{M}_1 U_0 , \end{align}\tag{9}\] i.e., 1 with \(f:= (1/2-{\mathscr N}) \phi - \mathscr{M}_1 U_0\). For \(U_0\in L^2(\Omega)\) and sufficiently smooth \(\phi\), one has the identities \[(\mathscr{M}_1 U_0)(t,{\boldsymbol{x}})=\int_\Omega \partial_{{\boldsymbol{n}}({\boldsymbol{x}})}G(t,{\boldsymbol{x}}-{\boldsymbol{y}}) U_0({\boldsymbol{y}}) \,{\rm d}{\boldsymbol{y}}\quad \text{for all } (t,{\boldsymbol{x}}) \in \Sigma\] and \[({\mathscr N} \phi)(t,{\boldsymbol{x}})= \int_\Sigma \partial_{{\boldsymbol{n}}({\boldsymbol{x}})}G(t-s,{\boldsymbol{x}}-{\boldsymbol{y}})\phi(s,{\boldsymbol{y}})\,{\rm d}{\boldsymbol{y}}\,{\rm d}s \quad \text{for all } (t,{\boldsymbol{x}}) \in \Sigma.\] As the hyper-singular operator \(\mathscr{W}\) is coercive, i.e., \[\begin{align} \label{eq:coercivity} \langle\mathscr{W}v\,,\,v\rangle_\Sigma \ge c_{\rm coe} \|v\|_{H^{1/2,1/4}(\Sigma)}^2 \quad \text{for all }v\in H^{1/2,1/4}(\Sigma) \end{align}\tag{10}\] with some constant \(c_{\rm coe}>0\), 9 is uniquely solvable and the solution \(u\) is just the missing lateral trace \(U|_\Sigma\) to compute \(U\) via the representation formula 5 .

Alternatively, one can make the ansatz \(U=\widetilde{\mathscr{M}}_0 U_0 - \widetilde{\mathscr{K}}u\). Indeed, both \(\widetilde{\mathscr{M}}_0 U_0\) and \(\widetilde{\mathscr{K}}u\) satisfy the heat equation, where \(\widetilde{\mathscr{M}}_0U_0\) restricted to \(\{0\}\times\Omega\) coincides with \(U_0\) and \(\widetilde{\mathscr{K}}u\) vanishes there. To satisfy the Dirichlet boundary conditions, one has to solve \[\begin{align} \label{eq:indirect} \mathscr{W}u = \phi - \mathscr{M}_1 U_0, \end{align}\tag{11}\] i.e., 1 with \(f:=\phi - \mathscr{M}_1 U_0\). While 9 is called direct method, as it directly provides the physically relevant quantity \(u=U|_\Sigma\), 11 is called indirect method.

For more details and proofs, we refer to the seminal works [15], [16], [28], which consider \(U_0=0\), and to [6], [29] for the general case.

We finally mention the additional regularities \(\mathscr{W}:\widetilde{H}^{1,1/2}(\Sigma)\to L^2(\Sigma)\) and \(\mathscr{N}:L^2(\Sigma)\to L^2(\Sigma)\), stated in [16] for piecewise smooth domains \(\Omega\); see also [30] for the additional regularity of the spatial trace operator. The mapping property \(\mathscr{W}^{-1}: L^2(\Sigma) \to \widetilde{H}^{1,1/2}(\Sigma)\) is even satisfied for general Lipschitz domains \(\Omega\). Moreover, replacing \(\Omega\) by \(\mathbb{R}^d\) and considering \((t,{\boldsymbol{x}})\in I\times \mathbb{R}^d\) in 6 , we extend \(\widetilde{\mathscr{M}}_0\) to an operator on \(H^1(\mathbb{R}^d)\). Parabolic regularity gives that \(\widetilde{\mathscr{M}}_0:H^1(\mathbb{R}^d)\to H^{2,1}(I\times \mathbb{R}^d)\) with \(H^{2,1}(I\times \mathbb{R}^d)\) defined analogously as in Section 2.2. (Indeed, on \(\mathbb{R}^d\), this is a simple consequence of the characterization of Sobolev norms in terms of the Fourier transform, e.g.,  [31].) Extending functions \(U_0\in H_0^1(\Omega)\) by zero to \(H^1(\mathbb{R}^d)\) and restricting \(\widetilde{\mathscr{M}}_0 U_0\) to \(Q\), we further see that \(\widetilde{\mathscr{M}}_0:H_0^1(\Omega)\to H^{2,1}(Q)\) with \(H^{2,1}(Q)\) defined as in Section 2.2. Since \(\partial_{\boldsymbol{n}}(\cdot): H^{2,1}(Q) \to L^2(\Sigma)\), we conclude that \(\mathscr{M}_1:H_0^1(\Omega) \to L^2(\Sigma)\).

2.4 Boundary meshes↩︎

Throughout this work, we consider prismatic meshes \(\mathcal{P}_h\) of \(\Sigma\). This means that \(\mathcal{P}_h\) is a set of prisms of the form \(P = J\times K\) with closed time intervals \(J \subset \overline{I}\) and closed simplices \(K \subset \Gamma\) that forms a partition of \(\Sigma\) in the sense that \(\overline{\Sigma} = \bigcup_{P\in \mathcal{P}_h} P\) and \(P \cap P'\) has \(d\)-dimensional measure zero for all \(P, P' \in \mathcal{P}_h\) with \(P \neq P'\). In addition, we assume that \(\mathcal{P}_h\) satisfies the following three properties.

• If the intersection \(P \cap P'\) of any \(P,P'\in\mathcal{P}_h\) has positive \((d-1)\)-dimensional measure, then \(P \cap P'\) is a hyperface of \(P\) or \(P'\).

• If \(J'\times K' = P'\in \mathcal{P}_h\) and \(J''\times K'' = P'' \in \mathcal{P}_h\) are temporal neighbors of \(J\times K=P\in \mathcal{P}_h\) on the same side, i.e., if \(P \cap P'\) and \(P \cap P''\) are \((d-1)\)-dimensional subsets of the same hyperface \(\{t\}\times K\) of \(P\), then \(J' = J''\).

• For almost every \(t\in I\), the spatial triangulation \(\big\{K\,:\,J\times K\in \mathcal{P}_h, t\in J\big\}\) is conforming (which is trivially satisfied if \(d=2\)).

We further suppose the existence of \(\emph{uniform}\) constants \(C_\text{shape}, C_\text{lqu}^t, C_\text{lqu}^{\boldsymbol{x}}\geq 1\) such that:

• the spatial parts in \(\mathcal{P}_h\) are shape-regular, i.e., \[\label{Cshape} C_\text{shape}^{-1} |K| \le \mathop{\mathrm{diam}}(K)^{d-1} \le C_\text{shape} |K| \quad \text{for all } J \times K \in \mathcal{P}_h,\tag{12}\] (which is trivially satisfied if \(d=2\));

• \(\mathcal{P}_h\) is locally quasi-uniform, i.e., \[\begin{align} \label{Cquasiunif} |J| \le C_\text{lqu}^t |J'| \quad \text{and} \quad \mathop{\mathrm{diam}}(K) \le C_\text{lqu}^{\boldsymbol{x}}\mathop{\mathrm{diam}}(K')\quad \end{align}\tag{13}\] for all \(J \times K, J' \times K' \in \mathcal{P}_h\) with \((J \times K) \cap (J' \times K') \neq \emptyset\).

In particular, the reliability constant of 4 will depend on these constants; see Theorem 4 below. A mesh refinement strategy that, starting from a tensor mesh, guarantees all these properties is given in Section 4.1.

2.5 Boundary element method↩︎

Given a finite-dimensional subspace \(X_h \subset H^{1/2,1/4}(\Sigma)\), let \(u_h\in X_h\) denote the Galerkin discretization of the solution \(u\) of the boundary integral equation 1 , i.e., \[\begin{align} \label{eq:Galerkin} \langle\mathscr{W}u_h\,,\,v_h\rangle_\Sigma = \langle f\,,\,v_h\rangle_\Sigma \quad \text{for all }v_h\in X_h, \end{align}\tag{14}\] which is equivalent to the Galerkin orthogonality \[\begin{align} \label{eq:orthogonality} \langle\mathscr{W}(u-u_h)\,,\,v_h\rangle_\Sigma = 0 \quad \text{for all }v_h\in X_h. \end{align}\tag{15}\] Note that coercivity 10 guarantees unique solvability of the latter equations, and the Céa lemma applies, i.e., \[\begin{align} \label{eq:cea} \|u-u_h\|_{H^{1/2,1/4}(\Sigma)} \le \frac{\|\mathscr{W}\|}{c_{\rm coe}}\,\min_{v_h\in X_h} \|u-v_h\|_{H^{1/2,1/4}(\Sigma)}, \end{align}\tag{16}\] where \(\|\mathscr{W}\|\) denotes the operator norm of \(\mathscr{W}:H^{1/2,1/4}(\Sigma) \to H^{-1/2,-1/4}(\Sigma)\).

Given a prismatic mesh as in Section 2.4, a natural choice of \(X_h\) would be the space of all \(\mathcal{P}_h\)-piecewise polynomials of some fixed degree \(p_t \in \mathbb{N}_0\) in time and \(p_x \in \mathbb{N}\) in space that are continuous in space. However, in order to ensure \(\mathscr{W}u_h \in L^2(\Sigma)\) for the envisaged a posteriori error estimate 4 , we require \(X_h \subset \widetilde{H}^{1,1/2}(\Sigma)\); see Section 2.3. With the space of all continuous \(\mathcal{P}_h\)-piecewise polynomials \(S^{p_t,p_x}(\mathcal{P}_h)\) of fixed degree \(p_t\in \mathbb{N}\) in time and \(p_x \in \mathbb{N}\) in space, we hence consider \[X_h := \widetilde{S}^{p_t,p_{\boldsymbol{x}}}(\mathcal{P}_h) := \big\{v_h \in S^{p_t,p_{\boldsymbol{x}}}(\mathcal{P}_h)\,:\,v_h(0,\cdot) = 0\big\}.\] Note that for \(f \in L^2(\Sigma)\), e.g., as in 9 or 11 with \(\phi\in L^2(\Sigma)\) and \(U_0 \in H^1_0(\Omega)\), we have the additional regularity \(u \in \widetilde{H}^{1,1/2}(\Sigma)\) (see Section 2.3) so that it makes sense to approximate \(u\) by functions vanishing at \(t = 0\).

If \(\mathcal{P}_h = \big\{J\times K\,:\,J\in \mathcal{P}_{h_t}, K\in\mathcal{P}_{h_{\boldsymbol{x}}}\big\}\) is a full tensor mesh corresponding to a mesh \(\mathcal{P}_{h_{\boldsymbol{x}}}\) of \(\Gamma\) with uniform mesh size \(h_{\boldsymbol{x}}\eqsim \mathop{\mathrm{diam}}(K)\) for all \(K\in\mathcal{P}_{h_{\boldsymbol{x}}}\) and a mesh \(\mathcal{P}_{h_t}\) of \(\overline{I}\) with uniform step size \(h_t\eqsim h_{\boldsymbol{x}}^{\sigma}\) for some \(\sigma>0\), then, [16]¹ suggests the error decay rate \[\begin{align} \label{eq:rates} \min_{v_h\in\widetilde{S}^{p_t,p_{\boldsymbol{x}}}(\mathcal{P}_h)} \|u-v_h\|_{H^{1/2,1/4}(\Sigma)} \lesssim N_h^{-\frac{\min\{p_{\boldsymbol{x}}+1/2,(p_t+3/4)\sigma\}}{d-1+\sigma}} \quad\text{for all smooth }u \text{ with }u(0,\cdot)=0; \end{align}\tag{17}\] see also [5] for Dirichlet boundary value problems. Here, \(N_h\eqsim h_{\boldsymbol{x}}^{-(d-1)} h_t^{-1} = h_{\boldsymbol{x}}^{-d+1-\sigma}\) denotes the number of degrees of freedom in \(X_h\). The optimal grading parameter is thus given by \(\sigma = (p_{\boldsymbol{x}}+\tfrac12)/(p_t+\tfrac34)\) with resulting rate \({\mathcal{O}}\big(N_h^{-\frac{p_{\boldsymbol{x}}+1/2}{d-1+\sigma}}\big)\). Similarly, we expect the decay rate \({\mathcal{O}}\big(N_h^-{\frac{\min\{p_{\boldsymbol{x}}+1, (p_t+1)\sigma\}}{d-1+\sigma}}\big)\) for the minimal \(L^2(\Sigma)\)-error.

2.6 Partition of unity↩︎

We introduce a basis of the spaces \(S^{1,1}(\mathcal{P}_h)\) and \(\widetilde{S}^{1,1}(\mathcal{P}_h)\). The sets of free nodes for \(S^{1,1}(\mathcal{P}_h)\) and \(\widetilde{S}^{1,1}(\mathcal{P}_h)\) are defined as \[\mathcal{N}_h := \{{\boldsymbol{z}}\in \overline{\Sigma} : {\boldsymbol{z}}\text{ is a vertex of every } P\in\mathcal{P}_h \text{ with } {\boldsymbol{z}}\in P\} \quad \text{and}\quad \widetilde{\mathcal{N}_h}:=\mathcal{N}_h \setminus (\{0\}\times\Gamma),\] respectively. Any other vertex that is not in \(\mathcal{N}_h\) is called hanging node, and the corresponding set is denoted by \(\mathcal{N}_h^\perp\). For each free node \({\boldsymbol{z}}\in\mathcal{N}_h\), the nodal basis function \(\varphi_{h,{\boldsymbol{z}}} \in S^{1,1}(\mathcal{P}_h)\) is characterized by \[\label{eq:kronecker} \varphi_{h,{\boldsymbol{z}}}({\boldsymbol{z}}^\prime)=\delta_{{\boldsymbol{z}},{\boldsymbol{z}}^\prime} \quad\text{ for all } {\boldsymbol{z}}^\prime\in \mathcal{N}_h.\tag{18}\] Proposition 2 states that \(\varphi_{h,{\boldsymbol{z}}}\) is indeed well defined by 18 and provides important properties of both \(\varphi_{h,{\boldsymbol{z}}}\) as well as its corresponding support \(\omega_h({\boldsymbol{z}}):=\mathop{\mathrm{supp}}(\varphi_{h,{\boldsymbol{z}}})\). The proof is given in Appendix 5.

3 A posteriori error estimation↩︎

The proof of the a posteriori error estimate 4 requires some preparations. In Section 3.1, we prove Poincaré-type inequalities. These are subsequently used to derive an approximation property of the interpolation operator introduced in Section 3.2. With this interpolation operator, we conclude the proof of 4 in Section 3.3.

3.1 Poincaré-type inequalities↩︎

The following lemma is implicitly stated in the proof of [16], while the proposition itself states the result for the full \(H^{1/2,1/4}(\Sigma)\)-norm. We give a detailed proof for the sake of completeness.

In order to prove a suitable approximation property for the interpolation operator constructed in Section 3.2 below, we require a Poincaré-type inequality on more general local subsets of \(\Sigma\). We follow the arguments of [33], where such a generalization is provided for the space \(H^{1,1/2}(Q)\) instead of \(H^{1/2,1/4}(\Sigma)\). For a measurable set \(\omega \subset \Sigma\) and \(v \in H^{1/2,1/4}(\Sigma)\), we define the localized seminorms \[\begin{align} |v|_{L^2H^{1/2}(\omega)}^2 &:=\int_I\int_\Gamma\int_\Gamma \chi_{\omega}(s,{\boldsymbol{x}})\,\chi_{\omega}(s,{\boldsymbol{y}})\, \frac{|v(s,{\boldsymbol{x}})-v(s,{\boldsymbol{y}})|^2}{|{\boldsymbol{x}}-{\boldsymbol{y}}|^d} \,{\rm d}{\boldsymbol{y}}\,{\rm d}{\boldsymbol{x}}\,{\rm d}s,\\[4pt] |v|_{H^{1/4}L^2(\omega)}^2 &:=\int_I\int_I\int_\Gamma \chi_{\omega}(t,{\boldsymbol{x}})\,\chi_{\omega}(s,{\boldsymbol{x}})\, \frac{|v(t,{\boldsymbol{x}})-v(s,{\boldsymbol{x}})|^2}{|t-s|^{3/2}}\,{\rm d}{\boldsymbol{x}}\,{\rm d}s\,{\rm d}t, \end{align}\] where \(\chi_{\omega}\colon\Sigma\to\{0,1\}\) denotes the characteristic function on \(\omega\). For all \(P\in \mathcal{P}\), we consider the following local subsets \[\label{eq:patchomegaP} \omega_h(P) := \bigcup_{{\boldsymbol{z}}\in \mathcal{N}_h: P\subseteq \mathop{\mathrm{supp}}(\varphi_{h,{\boldsymbol{z}}})} \omega_h({\boldsymbol{z}}).\tag{19}\] From the Poincaré-type inequality on \(\omega_h(P)\), we then also derive a local estimate for \(\| v \|_{L^2(\omega_h(P))}\) if \(P\) is close to \(\{0\}\times \Gamma\).

3.2 Interpolation operator↩︎

The crucial ingredient in the proof of the a posteriori error estimate 4 is a suitable interpolation operator \(\mathcal{I}_h: H^{1/2,1/4}(\Sigma) \to \widetilde{S}^{1,1}(\mathcal{P}_h)\). We introduce a Clément–type interpolation operator \(\mathcal{I}_h\colon L^2(\Sigma)\to \widetilde{S}^{1,1}(\mathcal{P}_h)\) by \[\label{Interpolationop} \mathcal{I}_h v:= \sum_{{\boldsymbol{z}}\in\widetilde{\mathcal{N}}_h} \left( \frac{1}{|\omega_h({\boldsymbol{z}})|}\int_{\omega_h({\boldsymbol{z}})}v(t,{\boldsymbol{x}})\,{\rm d}{\boldsymbol{x}}\,{\rm d}t \right)\varphi_{h,{\boldsymbol{z}}}.\tag{25}\]

For the crucial local approximation property of Lemma 3 below, we first need to show the following elementary properties of \(\mathcal{I}_h\).

For the following approximation property of \(\mathcal{I}_h\), we additionally require local parabolic scaling of \(\mathcal{P}_h\) in the sense that there exists a uniform constant \(C_\text{par} \ge 1\) such that \[\label{eq:localparabolicscaling} C_\text{par}^{-1} \mathop{\mathrm{diam}}(K)^2 \le |J| \le C_\text{par} \mathop{\mathrm{diam}}(K)^2 \quad \text{for all } J \times K \in \mathcal{P}_h.\tag{26}\] To ease notation, we introduce the local spatial mesh-size function \(h_{\boldsymbol{x}}\in L^\infty(\Sigma)\), which is defined \(\mathcal{P}_h\)-piecewise via \[h_{\boldsymbol{x}}|_{J \times K} := \mathop{\mathrm{diam}}(K) \quad \text{for all } J \times K \in \mathcal{P}_h.\]

3.3 A posteriori error estimate↩︎

With the preparations of the previous sections, we are finally in the position to prove our main result.

4 Numerical experiments↩︎

In this section, we employ the error estimator \[\label{eq:estimator} \eta_h^2:=\sum_{J\times K\in \mathcal{P}_h} \eta_h(J\times K)^2 \quad \text{with} \quad \eta_h(J\times K):=\text{diam}(K)^{1/2}\|f-\mathscr{W}u_h\|_{L^2(J\times K)}\tag{29}\] within an adaptive algorithm and numerically investigate the resulting convergence rates. We restrict ourselves to the case \(d=2\), with \(\Gamma=\partial \Omega\) being the boundary of a polygonal domain in \(\mathbb{R}^2\). Details on the numerical computation of all involved (singular) integrals are given in Appendix 6 for \(d\in\{2,3\}\).

4.1 Adaptive algorithm↩︎

The following adaptive algorithm is applied.

To ensure reliability of the estimator, it is crucial that the refinement step (iv) ensures parabolic scaling 26 and uniformly preserves the mesh assumptions of Section 2.4, which also requires refining non-marked elements. We start with an arbitrary tensor mesh \(\mathcal{P}_0\) and assign level 0 to every element in it. With each parabolic refinement (Figure 2), the level of the new elements is increased by one. We choose \(\mathcal{P}_{\ell+1}\) in the refinement step (iv) as the minimal refinement of \(\mathcal{P}_\ell\) such that all elements in \(\mathcal{M}_\ell\) are refined and such that the level difference between two elements sharing an edge is at most one. It follows immediately that all conditions from Section 2.4 are fulfilled. The constants \(C_\text{lqu}^t\) and \(C_\text{lqu}^{\boldsymbol{x}}\) from 13 can be chosen as \(C_\text{lqu}^t=16\) and \(C_\text{lqu}^{\boldsymbol{x}}=4\).

4.2 Reference for exact error↩︎

Since the exact error \(\|u-u_h\|_{H^{1/2,1/4}(\Sigma)}\) is not available in practice, either because its computation is too expensive or because the exact solution is unknown, we use the following \((h-h/2)\)-estimator \(\zeta_h\) as a reference quantity instead. For a prismatic mesh \(\mathcal{P}_h\), let \(\widehat{\mathcal{P}}_h\) denote the uniform isotropic refinement of \(\mathcal{P}_h\), i.e., \(\widehat \mathcal{P}_h\) results from \(\mathcal{P}_h\) by bisecting all elements once in space and once in time. Let \(\widehat{u}_h\) be the corresponding Galerkin approximation of the exact solution \(u\). The \((h-h/2)\)-estimator is defined by \[\begin{align} \label{eq:hhalbe} \|u_h-\widehat u_h\|_{\mathscr{W}} := \Bigl(\langle \mathscr{W}(u_h-\widehat u_h),\, u_h-\widehat u_h\rangle_\Sigma\Bigr)^{1/2}. \end{align}\tag{30}\] Under the saturation assumption \[\label{eq:qsat} \|u-\widehat u_h\|_{H^{1/2,1/4}(\Sigma)} \le q_{\rm sat}\, \|u-u_h\|_{H^{1/2,1/4}(\Sigma)}, \qquad 0<q_{\rm sat}<1,\tag{31}\] the triangle inequality and boundedness and coercivity of \(\mathscr{W}\) show that this estimator is equivalent to the error \(\|u-u_h\|_{H^{1/2,1/4}(\Sigma)}\). Note that the saturation assumption is indeed satisfied under the realistic (asymptotic) assumption that \(\|u-u_h\|_{H^{1/2,1/4}(\Sigma)} = \mathcal{O}((\# \mathcal{P}_h)^{-s})\) for some arbitrary rate \(s > 0\).

If the exact solution \(u\) is known, we additionally compute the \(L^2\)-error \(\|u-u_h\|_{L^2(\Sigma)}\).

4.3 Experiments↩︎

In the following numerical experiments, we consider the heat equation 3 with end time \(T=1\) on different two-dimensional domains \(\Omega\) and given initial condition \(U_0 \in H_0^1(\Omega)\) and Neumann datum \(\phi \in L^2(\Sigma)\). We approximate the corresponding Dirichlet datum \(u\), being the exact solution of the boundary integral equation 9 , by the adaptive Algorithm 1 with \(\theta = 0.5\). Note that the estimator is indeed well-defined by the regularity of the data \(U_0\) and \(\phi\) and the approximations \(u_\ell\); see Section 2.4. For comparison, we also consider \(\theta = 1\), which results in uniform refinement.

In each example, we choose the initial mesh \(\mathcal{P}_0\) consisting of rectangular elements of size 0.5 in both space and time. The employed refinement strategy guarantees parabolically scaled meshes, i.e., \(\sigma = 2\), which is required to guarantee reliability of the estimator; see Theorem 4. In case of smooth solutions \(u\), Section 2.5 yields that for \(p_t=1\), the optimal choice of \(p_{\boldsymbol{x}}\) is 3; see also Table 1. We fix the temporal polynomial degree as \(p_t:=1\) and consider the lowest-order choice \(p_{\boldsymbol{x}}:=1\) and the optimal choice \(p_{\boldsymbol{x}}:=3\) in space.

4.3.1 Smooth problem on the square↩︎

We prescribe the exact solution \[\begin{align} \label{problem:cos} U(t,{\boldsymbol{x}}):=1-e^{-4\pi^2 t}(\cos(2\pi x_1)+\cos(2\pi x_2))+e^{-8\pi^2 t}\cos(2\pi x_1)\cos(2\pi x_2) \end{align}\tag{32}\] on the square \(\Omega := (0,1)^2\) and choose the initial condition \(U_0\) and the Neumann datum \(\phi\) accordingly. It is noteworthy that the Dirichlet trace converges fast to \(1\) as \(t\) increases. In the numerical results, displayed in Figure 3, we observe optimal rates for both adaptive and uniform refinement. However, owing to the mentioned behavior of the solution \(U\) (and its Dirichlet trace \(u\)), adaptive refinement achieves much better results. In each case, we see a close resemblance between the estimator \(\eta_\ell =\|h_{{\boldsymbol{x}}}^{1/2}(f-\mathscr{W}u_\ell)\|_{L^2(\Sigma)}\) and the (\(h-h/2\))-estimator \(\|u_\ell - \widehat u_\ell\|_{\mathscr{W}}\), which numerically indicates both reliability and efficiency of \(\eta_\ell\). Note that in this plot and in all following plots, \(\|u_\ell-\widehat u_\ell\|_{\mathscr{W}}\) is not always computed over the same range as \(\eta_\ell\), owing to its high computational cost.

4.3.2 Almost singular problem on the square↩︎

We next prescribe the exact solution as shifted heat kernel \[\label{problem:heatkernel} U(t,{\boldsymbol{x}}) := G(t,{\boldsymbol{x}}-{\boldsymbol{z}})\tag{33}\] with a fixed singularity \({\boldsymbol{z}}\) lying outside of \(\overline{\Omega}\) with \(\Omega = (0,1)^2\). The initial condition \(U_0\) and the Neumann datum \(\phi\) are chosen accordingly. Note that \(U_0 = 0\). While the Dirichlet datum \(u\) is smooth, depending on the choice of \({\boldsymbol{z}}\), significant differences between adaptive and uniform refinement can be observed. We examine the problem with the singularity \({\boldsymbol{z}}= (0.5,-0.25)\). In the numerical results, displayed in Figure 4, we observe optimal rates for adaptive refinement and uniform refinement with \(p_{\boldsymbol{x}}= 1\), while uniform refinement with \(p_{\boldsymbol{x}}=3\) still seems to be in a pre-asymptotic regime. Again, the displayed curves indicate reliability and efficiency of the estimator \(\eta_\ell\). The strong adaptive refinement towards the singularity \({\boldsymbol{z}}\) can be seen in Figure 5. The smallest elements are of the size \(h_t = 2^{-11}\) and \(h_{\boldsymbol{x}}= 2^{-6}\).

4.3.3 Singular problem on the square↩︎

We consider the heat equation 3 with given data \[\begin{align} \label{problem:singular} U_0({\boldsymbol{x}})=0\quad\text{and}\quad \phi(t,{\boldsymbol{x}})=1 \end{align}\tag{34}\] on the square \(\Omega = (0,1)^2\). Due to the incompatibility between \(\phi\) and \(U_0\) in the sense that \(\phi(0,{\boldsymbol{x}}) \neq \partial_{\boldsymbol{n}}U_0({\boldsymbol{x}})\), we expect the unknown solution \(U\) (and its Dirichlet trace \(u\)) to be singular at \(t = 0\). In the numerical results, displayed in Figure 6, we observe convergence rates of order \(\mathcal{O}(N_\ell^{-1/2})\) for uniform refinement, which is suboptimal for \(p_{\boldsymbol{x}}\neq 1\). In contrast, for \(p_{\boldsymbol{x}}= 3\), adaptive refinement yields nearly optimal rates, namely approximately \(\mathcal{O}(N_\ell^{-1})\), in contrast to \(\mathcal{O}(N_\ell^{-7/6})\) expected for smooth solutions.

4.3.4 Singular problem on the L-shape↩︎

We consider the heat equation 3 with given data 34 on the L-shaped domain \(\Omega = (-1,1)^2\setminus [0,1]^2\). Due to the incompatibility between \(\phi\) and \(U_0\) in the sense that \(\phi(0,{\boldsymbol{x}}) \neq \partial_{\boldsymbol{n}}U_0({\boldsymbol{x}})\), we expect the unknown solution \(U\) (and its Dirichlet trace \(u\)) to be singular at \(t = 0\). In the numerical results, displayed in Figure 7, we observe convergence rates of order \(\mathcal{O}(N_\ell^{-1/2})\) for uniform refinement, which is suboptimal for \(p_{\boldsymbol{x}}\neq 1\). In contrast, for \(p_{\boldsymbol{x}}= 3\), adaptive refinement yields nearly optimal rates, namely approximately \(\mathcal{O}(N_\ell^{-1})\), in contrast to \(\mathcal{O}(N_\ell^{-7/6})\) expected for smooth solutions.

5 Proof of Proposition 2↩︎

In this section, we prove Proposition 2. We start with the following lemma.

6 Numerical computation↩︎

This section describes the implementation of the Galerkin system 14 with \(f\) given as in 9 and of the estimator 29 for \(d\in\{2,3\}\). Throughout, let \(\mathcal{P}_h\) be a prismatic mesh of \(\Sigma\).

6.1 Right-hand side↩︎

We discuss the discretization of the two parts \((1/2-\mathscr{N})\phi_h\) and \(\mathscr{M}_1 U_0\) of the right-hand side \(f\) given as in 9 .

6.1.1 Neumann datum↩︎

For our computations, we replace the Neumann datum \(\phi\) by its \(L^2(\Sigma)\)-orthogonal projection \(\phi_h := \Pi_h \phi\) onto \(\mathcal{P}_h\)-piecewise polynomials of degree \(p_t\) in time and \(p_{\boldsymbol{x}}\) in space. In particular, we replace the right-hand side \(f\) of 9 by \[f_h := (1/2-\mathscr{N})\phi_h - \mathscr{M}_1 U_0.\] To ensure that the best possible convergence rate \({\mathcal{O}}\big(N_h^{-\frac{\min\{p_{\boldsymbol{x}}+1/2,2p_t+3/2\}}{d+1}}\big)\) for parabolically scaled meshes \(\mathcal{P}_h\) (see 17 ) is not deteriorated, we must control the resulting error \(\| u_h - \widetilde{u}_h\|_{H^{1/2,1/4}(\Sigma)}\), where \(\widetilde{u}_h\) denotes the Galerkin solution of 9 with \(f\) replaced by \(f_h\). By stability of the Galerkin projection and the operator \(\mathscr{N}\), we see that

\[\|u_h - \widetilde{u}_h\|_{H^{1/2,1/4}(\Sigma)} \lesssim \|(1/2-\mathscr{N})(\phi-\phi_h)\|_{H^{-1/2,-1/4}(\Sigma)} \lesssim \|\phi-\phi_h\|_{H^{-1/2,-1/4}(\Sigma)}.\] Using the characterization of the dual norm and the \(L^2(\Sigma)\)-orthogonality of \(\Pi_h\), we get that \[\begin{align} \|\phi-\phi_h\|_{H^{-1/2,-1/4}(\Sigma)} &= \sup_{v\in H^{1/2,1/4}(\Sigma)\setminus\{0\}} \frac{\bigl|\langle (1-\Pi_h)\phi, v\rangle_{\Sigma}\bigr|}{\|v\|_{H^{1/2,1/4}(\Sigma)}}\\ &= \sup_{v\in H^{1/2,1/4}(\Sigma)\setminus\{0\}} \frac{\bigl|\langle h_{\boldsymbol{x}}^{1/2}(1-\Pi_h)\phi,\,h_{\boldsymbol{x}}^{-1/2}(1-\Pi_h)v\rangle_{\Sigma}\bigr|}{\|v\|_{H^{1/2,1/4}(\Sigma)}}\\ &\leq \| h_{\boldsymbol{x}}^{1/2}(1-\Pi_h)\phi\|_{L^2(\Sigma)} \sup_{v\in H^{1/2,1/4}(\Sigma)\setminus\{0\}} \frac{\|h_{\boldsymbol{x}}^{-1/2}(1-\Pi_h)v\|_{L^2(\Sigma)}}{\|v\|_{H^{1/2,1/4}(\Sigma)}}. \end{align}\] Since \(\Pi_h\) is a local projection, we may replace it by the interpolation operator \(\mathcal{I}_h\) and we conclude by Proposition 3 that the second factor is indeed finite. For uniform parabolically scaled meshes \(\mathcal{P}_h\), standard approximation results show for the first factor that \[\| h_{\boldsymbol{x}}^{1/2}(1-\Pi_h)\phi\|_{L^2(\Sigma)} = {\mathcal{O}}\big(\|h_{\boldsymbol{x}}\|_{L^\infty(\Sigma)}^{\min \{ p_{\boldsymbol{x}}+3/2,2p_t+5/2 \} }\big)= {\mathcal{O}}\big(N_h^{-\frac{\min \{ p_{\boldsymbol{x}}+3/2,2p_t+5/2 \} }{d+1}}\big).\] We conclude that the approximation error introduced by replacing \(\phi\) by \(\phi_h\) is of higher order.

We next discuss the computaton of \(\mathscr{N} \phi_h\) as well as \(\langle \mathscr{N} \phi_h,v_h\rangle\). For \((t,{\boldsymbol{x}}) \in \Sigma\), we can write

\[\begin{align} \label{eq:appB1} (\mathscr{N} \phi_h)(t,{\boldsymbol{x}}) \overset{\eqref{eq:partialnyG}}{=} \sum_{J\times K \in \mathcal{P}_h} \int_K \int_J \frac{({\boldsymbol{y}}-{\boldsymbol{x}})\cdot{\boldsymbol{n}}({\boldsymbol{x}})}{2(t-s)}G(t-s,{\boldsymbol{x}}-{\boldsymbol{y}})\phi_h(s,{\boldsymbol{y}})\,{\rm d}s \,{\rm d}{\boldsymbol{y}}. \end{align}\tag{36}\] These evaluations of 36 are required for the numerical quadrature of the error estimator 29 .

For piecewise polynomial test functions \(v_h\), we can write

\[\begin{align} \label{eq:appB2} &\langle \mathscr{N}\phi_h,v_h\rangle =\\ \notag & \sum_{J_1\times K_1 \in \mathcal{P}_h} \sum_{J_2\times K_2 \in \mathcal{P}_h} \int_{K_1} \int_{K_2} \int_{J_1} \int_{J_2} \frac{({\boldsymbol{y}}-{\boldsymbol{x}})\cdot{\boldsymbol{n}}({\boldsymbol{x}})}{2(t-s)}G(t-s,{\boldsymbol{x}}-{\boldsymbol{y}})\phi_h(s,{\boldsymbol{y}})v_h(t,{\boldsymbol{x}})\,{\rm d}s \,{\rm d}t\,{\rm d}{\boldsymbol{y}}\,{\rm d}{\boldsymbol{x}}. \end{align}\tag{37}\] These duality products of 37 are required for the numerical quadrature of the Galerkin system 14 .

In Section 6.3, we will further simplify the expressions 36 –37 by analytically computing the time integrals, and in Section 6.4, we will address the numerical computation of the space integrals.

6.1.2 Initial condition↩︎

To compute the normal derivative of the initial potential applied to the initial condition \(U_0\), we introduce a mesh \(\mathcal{T}_h\) of \(\Omega\). For \((t,{\boldsymbol{x}}) \in \Sigma\), we can write \[\begin{align} \label{eq:initialcond1} (\mathscr{M}_1 U_0)(t,{\boldsymbol{x}})&=\sum_{K \in \mathcal{T}_h} \int_K ({\boldsymbol{y}}-{\boldsymbol{x}})\cdot{\boldsymbol{n}}({\boldsymbol{x}})\frac{1}{2t} G(t,{\boldsymbol{x}}-{\boldsymbol{y}}) U_0({\boldsymbol{y}}) \,{\rm d}{\boldsymbol{y}}. \end{align}\tag{38}\] Note that the integrands are smooth for \(t>0\) so that we can use standard quadrature. These evaluations of 38 are required for the numerical quadrature of the error estimator 29 .

For piecewise polynomial test functions \(v_h\), we can write

\[\begin{align} \label{eq:initialcond2} \langle \mathscr{M}_1 U_0,v_h \rangle &= \sum_{J_1 \times K_1 \in \mathcal{P}_h} \sum_{K \in \mathcal{T}_h} \int_{K_1} \int_K \int_{J_1}({\boldsymbol{y}}-{\boldsymbol{x}})\cdot{\boldsymbol{n}}({\boldsymbol{x}})\frac{1}{2t} G(t,{\boldsymbol{x}}-{\boldsymbol{y}}) U_0({\boldsymbol{y}}) v_h(t,{\boldsymbol{x}}) \,{\rm d}t\,{\rm d}{\boldsymbol{y}}\,{\rm d}{\boldsymbol{x}}. \end{align}\tag{39}\] These duality products of 39 are required for the numerical quadrature of the Galerkin system 14 .

In Section 6.3, we will further simplify the expression 38 by analytically computing the time integral, and in Section 6.4, we will address the numerical computation of the space integrals.

6.2 Left-hand side↩︎

Similarly as in the stationary case [35], for sufficiently smooth functions \(u\), \(v\) on \(\Sigma\), the duality product \(\langle \mathscr{W}u,v\rangle_\Sigma\) for the hypersingular integral operator \(\mathscr{W}= -\partial_{\boldsymbol{n}}\widetilde{\mathscr{K}}\colon H^{1/2,1/4}(\Sigma)\to H^{-1/2,-1/4}(\Sigma)\) can be rewritten in terms of the single-layer operator \(\mathscr{V}:=(\cdot)|_\Sigma \circ \widetilde{\mathscr{V}}\). For sufficiently smooth \(\phi\), \(\mathscr{V}\phi\) is given as

\[(\mathscr{V}\phi)(t,{\boldsymbol{x}}):=\int_\Sigma G(t-s,{\boldsymbol{x}}-{\boldsymbol{y}})\phi (s,{\boldsymbol{y}})\,{\rm d}{\boldsymbol{y}}\,{\rm d}s \quad \text{for all } (t,{\boldsymbol{x}}) \in \Sigma.\] The following integration-by-parts formula is found in [16] for \(d=2\) and in [3] for \(d=3\); see also [11] for a detailed proof.

Using another integration by parts [35], we obtain the following explicit formula for the hypersingular operator \[\label{eq:hypersingularopvar} \mathscr{W}u=- D_\Gamma'(\mathscr{V} (D_\Gamma u))+{\boldsymbol{n}}\cdot \partial_t \mathscr{V}(u \,{\boldsymbol{n}})\tag{40}\] with \[D_\Gamma' (\cdot):=\begin{cases}D_\Gamma (\cdot) & \quad \text{if } d=2,\\ {\boldsymbol{n}}\cdot D_\Gamma (\cdot) & \quad \text{if } d=3.\end{cases}\] For a continuous piecewise polynomial \(u_h\) with \(u_h(0,\cdot) = 0\) and \((t,{\boldsymbol{x}}) \in \Sigma\), we can write the first summand in 40 without the differential operator \(D_\Gamma'\) in front as

\[\label{eq:hypersingular1} (\mathscr{V}(D_\Gamma u_h))(t,{\boldsymbol{x}}) = \sum_{J\times K \in \mathcal{P}_h} \int_K \int_J G(t-s,{\boldsymbol{x}}-{\boldsymbol{y}})D_\Gamma u_h (s,{\boldsymbol{y}})\,{\rm d}s\,{\rm d}{\boldsymbol{y}}.\tag{41}\] For the implementation, we replace the application of the differential operator \(D_\Gamma'\) to 41 by the application of \(D_\Gamma'\) to the Lagrange interpolation of 41 . We can write the second summand in 40 with integration by parts in time applied to \(\partial_t G(t-s,{\boldsymbol{x}}-{\boldsymbol{y}}) = -\partial_s G(t-s,{\boldsymbol{x}}-{\boldsymbol{y}})\) as \[\begin{align} \label{eq:hypersingular2} \begin{aligned} {\boldsymbol{n}}({\boldsymbol{x}}) \cdot(\partial_t \mathscr{V}( u_h {\boldsymbol{n}}))(t,{\boldsymbol{x}}) &= \int_\Gamma \int_0^T \partial_t G(t-s,{\boldsymbol{x}}-{\boldsymbol{y}})u_h (s,{\boldsymbol{y}}) {\boldsymbol{n}}({\boldsymbol{x}})\cdot {\boldsymbol{n}}({\boldsymbol{y}}) \,{\rm d}s\,{\rm d}{\boldsymbol{y}}\\ &= \sum_{J\times K \in \mathcal{P}_h} \int_K \int_J G(t-s,{\boldsymbol{x}}-{\boldsymbol{y}})\partial_s u_h (s,{\boldsymbol{y}}) {\boldsymbol{n}}({\boldsymbol{x}})\cdot {\boldsymbol{n}}({\boldsymbol{y}}) \,{\rm d}s\,{\rm d}{\boldsymbol{y}}.\end{aligned} \end{align}\tag{42}\] The boundary terms vanish since \(u_h(0,\cdot)=0\) and \(G(t-T,{\boldsymbol{x}}-{\boldsymbol{y}})=0\). The evaluation of 40 is required for the computation of the error estimator 29 .

For continuous piecewise polynomial test functions \(v_h\), 41 –42 imply for the two summands in ?? that

\[\label{eq:hypersingular3}\begin{align} &\langle\mathscr{V}(D_\Gamma u_h),D_\Gamma v_h\rangle\\ = &\sum_{J_1\times K_1 \in \mathcal{P}_h} \sum_{J_2\times K_2 \in \mathcal{P}_h} \int_{K_1} \int_{K_2} \int_{J_1} \int_{J_2} G(t-s,{\boldsymbol{x}}-{\boldsymbol{y}})D_\Gamma u_h(s,{\boldsymbol{y}})D_\Gamma v_h (t,{\boldsymbol{x}}) \,{\rm d}s\,{\rm d}t\,{\rm d}{\boldsymbol{y}}\,{\rm d}{\boldsymbol{x}}, \end{align}\tag{43}\] and

\[\begin{align} \label{eq:hypersingular4} &\langle\partial_t \mathscr{V}(u_h {\boldsymbol{n}}), v_h {\boldsymbol{n}}\rangle\\\notag = &\sum_{J_1\times K_1 \in \mathcal{P}_h} \sum_{J_2\times K_2 \in \mathcal{P}_h} \int_{K_1} \int_{K_2} \int_{J_1} \int_{J_2} G(t-s,{\boldsymbol{x}}-{\boldsymbol{y}})\partial_s u_h (s,{\boldsymbol{y}}) v_h(t,{\boldsymbol{x}}){\boldsymbol{n}}({\boldsymbol{x}})\cdot {\boldsymbol{n}}({\boldsymbol{y}})\,{\rm d}s\,{\rm d}t\,{\rm d}{\boldsymbol{y}}\,{\rm d}{\boldsymbol{x}}. \end{align}\tag{44}\] The duality products of ?? are required for the computation of the Galerkin system 14 .

In Section 6.3, we will further simplify the expressions 41 –44 by analytically computing the time integrals, and in Section 6.4, we will address the numerical computation of the space integrals.

6.3 Analytical computation of time integrals.↩︎

It is well-known that integration in time for the considered boundary integral operators can be performed exactly in case that the ansatz and test functions are piecewise constant in time; see, e.g., [9], [16], [21], [27]. We generalize these explicit formulas to arbitrary polynomial degree in time.

6.3.1 Single time integrals↩︎

As we have seen in the previous sections, the time integrals of \(\langle \mathscr{M}_1 U_0,v_h\rangle_\Sigma\) in the Galerkin system 14 and of the residual in the estimator 29 can be reduced to \[\begin{align} \label{formelintegral0} &\int_a^b \frac{1}{(t-s)^{1-m}}\,G(t-s,{\boldsymbol{z}})\,s^j\,{\rm d}s \end{align}\tag{45}\] with \(m\in\mathbb{Z}\) (more precisely only \(m\in\{0,1\}\)), a time interval \([a,b]\subseteq \overline{I}\), a time point \(t\in I\), \({\boldsymbol{z}}\in \mathbb{R}^d\setminus\{0\}\), and \(j\in\mathbb{N}_0\). For \(\alpha\in\mathbb{R}\) and \(r>0\), we introduce the incomplete gamma function \[\begin{align} \label{ableitungGamma} \Gamma(\alpha,r):=\int_{r}^\infty s^{\alpha-1}e^{-s}~\,{\rm d}s\quad\text{ with derivative }\quad\frac{\,{\rm d}}{\,{\rm d}r}\Gamma(\alpha,r)=-r^{\alpha-1}e^{-r}. \end{align}\tag{46}\] We abbreviate \(\delta:=\frac{d}{2}\) and \(\zeta:=\frac{|{\boldsymbol{z}}|^2}{4}\). For \(t>b\), integration by substitution \((s\mapsto t-\frac{\zeta}{s})\) and the binomial formula show that \[\begin{align} \label{eq:innerintegral} &\int_{a}^{b} \frac{1}{(t-s)^{1-m}}G(t-s,{\boldsymbol{z}})s^j \,{\rm d}s=\frac{1}{(4\pi)^\delta}\int_{a}^{b} \frac{1}{(t-s)^{1-m+\delta}}e^{\frac{-\zeta}{t-s}}s^j \,{\rm d}s\\ \notag & \quad=\frac{1}{(4\pi)^\delta}\int^{\frac{\zeta}{t-b}}_{\frac{\zeta}{t-a}} \Bigl( \frac{s}{\zeta}\Bigr)^{1-m+\delta} e^{-s}\frac{\zeta}{s^2}\Bigl(t-\frac{\zeta}{s}\Bigr)^j~\,{\rm d}s\\ \notag & \quad =\frac{1}{(4\pi)^\delta}\sum_{k=0}^j \binom{j}{k} t^{j-k} (-1)^k\zeta^{m-\delta+k} \int^{\frac{\zeta}{t-b}}_{\frac{\zeta}{t-a}} s^{\delta-k-m-1}e^{-s}~\,{\rm d}s\\ \notag & \quad =\frac{1}{(4\pi)^\delta}\sum_{k=0}^j\binom{j}{k} t^{j-k} (-1)^k \zeta^{m-\delta+k}\biggl[ \Gamma\Big(\delta-k-m,\frac{\zeta}{t-a}\Big)-\Gamma\Big(\delta-k-m,\frac{\zeta}{t-b}\Big)\biggr]. \end{align}\tag{47}\] For general \(t\geq 0\), the piecewise definition of \(G\) yields that \[\label{formelintegral} \int_{a}^{b} \frac{1}{(t-s)^{1-m}}G(t-s,{\boldsymbol{z}})s^j \,{\rm d}s=\mathfrak{g}^{j,m}_{t,a}({\boldsymbol{z}})-\mathfrak{g}^{j,m}_{t,b}({\boldsymbol{z}})\tag{48}\] with \[\begin{align} \mathfrak{g}_{t,\sigma}^{j,m}({\boldsymbol{z}}):=\begin{cases}\frac{1}{(4\pi)^\delta}\sum_{k=0}^j \binom{j}{k}(-1)^k \zeta^{m-\delta+k}t^{j-k} \Gamma(\delta-k-m,\frac{\zeta}{t-\sigma})&\text{if } t-\sigma>0, \\ 0&\text{else,}\end{cases} \end{align}\] \(\sigma\in \{a,b\}\).

6.3.2 Double time integrals↩︎

As we have seen in the previous sections, the time integrals for the duality products in the Galerkin system 14 can be reduced to \[\int_{a_1}^{b_1}\int_{a_2}^{b_2}\frac{1}{(t-s)^{1-m}}G(t-s,{\boldsymbol{z}})s^jt^i\,{\rm d}s\,{\rm d}t=\int_{a_1}^{b_1}(\mathfrak{g}^{j,m}_{t,a_2}({\boldsymbol{z}})-\mathfrak{g}^{j,m}_{t,b_2}({\boldsymbol{z}}))t^i\,{\rm d}t\] with \(m\in\mathbb{Z}\) (more precisely only \(m\in\{0,1\}\)), time intervals \([a_1,b_1],[a_2,b_2]\subseteq \overline{I}\), \({\boldsymbol{z}}\in \mathbb{R}^d\setminus\{0\}\), and \(i,j\in\mathbb{N}_0\). Recalling 48 , it remains to calculate the integral of \(t^{i+j-k}\Gamma(\delta-k-m,\frac{\zeta}{t-\sigma})\) with respect to \(t\). We abbreviate \(\beta := \delta -k -m\). Then, the substitution \(t \mapsto t+\sigma\) gives that \[\int_{a_1}^{b_1} t^{i+j-k}\Gamma\Big(\beta,\frac{\zeta}{t-\sigma}\Big)\,{\rm d}t =\int_{a_1-\sigma}^{b_1-\sigma} \sum_{\ell=0}^{i+j-k}\binom{i+j-k}{\ell}\sigma^{i+j-k-\ell}t^\ell\Gamma\Big(\beta,\frac{\zeta}{t}\Big)\,{\rm d}t\] For \(t \in \{a_1,b_1\}\), integration by parts (using 46 to get \(\frac{\,{\rm d}}{\,{\rm d}t}\Gamma(\beta,\frac{\zeta}{t})=\zeta^\beta t^{-1-\beta} e^{-\zeta/t}\)) and integration by substitution \((t\mapsto \frac{\zeta}{t})\) give that \[\begin{align} \int_{a_1-\sigma}^{b_1-\sigma} t^\ell \Gamma \Big(\beta,\frac{\zeta}{t}\Big)\,{\rm d}t&=\biggl[\frac{t^{\ell+1}}{\ell+1}\Gamma\Big(\beta,\frac{\zeta}{t}\Big) \biggr]_{a_1-\sigma}^{b_1-\sigma}-\int_{a_1-\sigma}^{b_1-\sigma} \frac{t^{\ell+1}}{\ell+1}\zeta^\beta t^{-1-\beta}e^{-\zeta/t}\,{\rm d}t\\ &=\biggl[\frac{t^{\ell+1}}{\ell+1}\Gamma\Big(\beta,\frac{\zeta}{t}\Big) \biggr]_{a_1-\sigma}^{b_1-\sigma}-\frac{1}{\ell+1}\int_{\zeta/(a_1-\sigma)}^{\zeta/(b_1-\sigma)} \Big(\frac{\zeta}{t}\Big)^{\ell-\beta}\zeta^\beta \frac{-\zeta}{t^2}e^{-t}\,{\rm d}t\\ &=\frac{1}{\ell+1}\Bigl[t^{\ell+1}\Gamma\Big(\beta,\frac{\zeta}{t}\Big)-\zeta^{\ell+1}\Gamma\Big(\beta-\ell-1,\frac{\zeta}{t}\Big)\Bigr]_{a_1-\sigma}^{b_1-\sigma}=:\mathfrak{H}^{\ell,\beta}_{b_1-\sigma}({\boldsymbol{z}})-\mathfrak{H}^{\ell,\beta}_{a_1-\sigma}({\boldsymbol{z}}). \end{align}\]

The piecewise definition of \(\mathfrak{g}_{t,\sigma}^{j,m}\) yields that \[\label{formeldoppelintegral} \int_{a_1}^{b_1}\int_{a_2}^{b_2} \frac{1}{(t-s)^{1-m}}G(t-s,{\boldsymbol{z}})s^jt^i \,{\rm d}s\,{\rm d}t=\mathfrak{G}^{i,j,m}_{a_1,a_2}({\boldsymbol{z}})-\mathfrak{G}^{i,j,m}_{b_1,a_2}({\boldsymbol{z}})-\mathfrak{G}^{i,j,m}_{a_1,b_2}({\boldsymbol{z}})+\mathfrak{G}^{i,j,m}_{b_1,b_2}({\boldsymbol{z}}),\tag{49}\] with \[\begin{align} \mathfrak{G}^{i,j,m}_{\tau,\sigma}({\boldsymbol{z}}):=\begin{cases} \frac{1}{(4\pi)^\delta}\sum_{k=0}^j \binom{j}{k}(-1)^k\zeta^{m-\delta+k} & \\ \quad\quad \sum_{\ell=0}^{i+j-k}\binom{i+j-k}{\ell}\sigma^{i+j-k-\ell}\mathfrak{H}^{\ell,\delta-k-m}_{\tau-\sigma}({\boldsymbol{z}})&\text{if } \tau-\sigma>0,\\ 0&\text{else}, \end{cases} \end{align}\] \(\tau \in \{a_1,b_1\}\), \(\sigma\in\{a_2,b_2\}\).

6.4 Approximate computation of space integrals↩︎

It remains to discuss the numerical quadrature of the singular integrals

\[\label{Einfachintegral} \int_{K_2} [({\boldsymbol{y}}- {\boldsymbol{x}})\cdot{\boldsymbol{n}}({\boldsymbol{x}})]^{1-m} \mathfrak{g}_{t,\sigma}^{j,m}({\boldsymbol{x}}-{\boldsymbol{y}}) u_{h_{\boldsymbol{x}}}({\boldsymbol{y}}) \,{\rm d}{\boldsymbol{y}}\tag{50}\] and \[\label{Doppelintegral} \int_{K_1} \int_{K_2} [({\boldsymbol{y}}- {\boldsymbol{x}})\cdot {\boldsymbol{n}}({\boldsymbol{x}})]^{1-m} \mathfrak{G}_{\tau,\sigma}^{i,j,m}({\boldsymbol{x}}-{\boldsymbol{y}}) u_{h_{\boldsymbol{x}}}({\boldsymbol{y}}) v_{h_{\boldsymbol{x}}}({\boldsymbol{x}}) \,{\rm d}{\boldsymbol{y}}\,{\rm d}{\boldsymbol{x}}\tag{51}\] for \((t,{\boldsymbol{x}})\in \Sigma\), \(J_1 \times K_1, J_2 \times K_2 \in \mathcal{P}_h\), \(\tau,\sigma\in \overline{I}\) with \(\tau>\sigma\) , \(i, j \in \mathbb{N}_0\), \(m \in \{0,1\}\), and piecewise polynomials \(u_{h_{\boldsymbol{x}}}\) and \(v_{h_{\boldsymbol{x}}}\) on \(\Gamma\). For the term \(\langle \mathscr{M}_1 U_0,v_h\rangle_\Sigma\), we additionally need to consider \[\label{Doppelintegral2} \int_{K_1}\int_K [({\boldsymbol{y}}- {\boldsymbol{x}})\cdot {\boldsymbol{n}}({\boldsymbol{x}})]^{1-m} \mathfrak{g}_{t,\sigma}^{i,m}({\boldsymbol{x}}-{\boldsymbol{y}}) U_0({\boldsymbol{y}}) v_{h_{\boldsymbol{x}}}({\boldsymbol{x}}) \,{\rm d}{\boldsymbol{y}}\,{\rm d}{\boldsymbol{x}}\tag{52}\] for \(K \subset \Omega\).

Since the boundary \(\Gamma\) is piecewise smooth, it holds that \(({\boldsymbol{y}}-{\boldsymbol{x}})\cdot {\boldsymbol{n}}({\boldsymbol{x}})={\mathcal{O}}(|{\boldsymbol{x}}-{\boldsymbol{y}}|^2)\) for \({\boldsymbol{x}},{\boldsymbol{y}}\) on the same smooth part of \(\Gamma\) and \(|{\boldsymbol{x}}-{\boldsymbol{y}}|\to 0\); see, e.g., [30]. (In fact, for piecewise flat \(\Gamma\), we even have that \(({\boldsymbol{y}}-{\boldsymbol{x}})\cdot {\boldsymbol{n}}({\boldsymbol{x}}) = 0\).) It remains to determine the strength of the singularity of \(\mathfrak{g}^{j,m}_{t,\sigma}(\mathbf{z})\) and \(\mathfrak{G}^{\,i,j,m}_{\tau,\sigma}(\mathbf{z})\) at \({\boldsymbol{z}}=0\). For \(\mathfrak{g}^{j,m}_{t,\sigma}(\mathbf{z})\) from 48 , the only term that requires a detailed study is \[\zeta^{\,m-\delta+k}\, \Gamma\Bigl(\delta-k-m,\frac{\zeta}{t-\sigma}\Bigr)\] for \(k \in \{0,\ldots,j\}\). Here, we use again the abbreviations \(\delta=d/2\) and \(\zeta=|\mathbf{z}|^2/4\). Similarly, for \(\mathfrak{G}^{\,i,j,m}_{\tau,\sigma}(\mathbf{z})\) from 49 , it is sufficient to study the term \[\zeta^{m-\delta+k}\Bigl(\zeta^{\ell+1}\Gamma(\delta-k-m-\ell-1,\tfrac{\zeta}{\tau-\sigma}) - (\tau-\sigma)^{\ell+1}\Gamma(\delta-k-m,\tfrac{\zeta}{\tau-\sigma})\Bigr)\] for \(k \in \{0,\ldots,j\}\) and \(\ell \in \{0,\ldots,i+j-k\}\). Abbreviating \(\xi:=\zeta/(\tau - \sigma)\), all these terms are of the form \[\label{eq:reduced-term} \xi^{n-\delta}\,\Gamma(\delta-n,\xi)\tag{53}\] for some \(n \in \mathbb{N}_0\) (either \(n = k+m\) or \(n = k+m+\ell+1\)). The following lemma characterizes the singularity of this term. Here, \(\Gamma(\cdot)\) denotes the (complete) gamma function.

We conclude that the dominant singularity in our case occurs for \(n = m\). If \(m=1\) (which occurs for the operator \(\mathscr{W}\), see 40 –44 ), the terms \(\mathfrak{g}_{t,\sigma}^{j,m}({\boldsymbol{x}}-{\boldsymbol{y}})\) and \(\mathfrak{G}_{\tau,\sigma}^{i,j,m}({\boldsymbol{x}}-{\boldsymbol{y}})\) have the singularity \(\ln|{\boldsymbol{x}}-{\boldsymbol{y}}|\) for \(d=2\) and \(|{\boldsymbol{x}}-{\boldsymbol{y}}|^{-1}\) for \(d=3\). If \(m=0\) (which occurs for the operator \(\mathscr{N}\) and \(\mathscr{M}_1\), see 36 , 37 , and 39 ), and \({\boldsymbol{x}},{\boldsymbol{y}}\) lie on the same smooth part of \(\Gamma\), the terms \(({\boldsymbol{y}}-{\boldsymbol{x}})\cdot {\boldsymbol{n}}({\boldsymbol{x}}) \mathfrak{g}_{t,\sigma}^{j,m}({\boldsymbol{x}}-{\boldsymbol{y}})\) and \(({\boldsymbol{y}}-{\boldsymbol{x}})\cdot{\boldsymbol{n}}({\boldsymbol{x}}) \mathfrak{G}_{\tau,\sigma}^{i,j,m}({\boldsymbol{x}}-{\boldsymbol{y}})\) are regular for \(d=2\) and have the singularity \(|{\boldsymbol{x}}-{\boldsymbol{y}}|^{-1}\) for \(d=3\), since \(({\boldsymbol{y}}-{\boldsymbol{x}})\cdot {\boldsymbol{n}}({\boldsymbol{x}})={\mathcal{O}}(|{\boldsymbol{x}}-{\boldsymbol{y}}|^2)\) for \({\boldsymbol{x}}\to {\boldsymbol{y}}\).

The integrals 50 –52 can now be computed as in [21] (in case of a triangular \(K \subset \Omega\)), using Duffy transformations as in [36] together with the quadrature rules from [37] for logarithmic singularities on intervals.

References↩︎

---

1. To be precise, the given reference only applies for \(X_h = S^{p_t,p_{\boldsymbol{x}}}(\mathcal{P}_h)\), but the proof directly extends to \(X_h = \widetilde{S}^{p_t,p_{\boldsymbol{x}}}(\mathcal{P}_h)\).↩︎