A posteriori error estimates for hypersingular integral equation on spheres with spherical splines

TL;DR

This paper develops and validates a posteriori error estimates for solving hypersingular integral equations on spheres using spherical splines, enabling adaptive mesh refinement to improve computational efficiency.

Contribution

It introduces new a posteriori residual and hierarchical error bounds for spherical splines on the sphere, facilitating adaptive algorithms for hypersingular integral equations.

Findings

Error estimates are proven and validated numerically.

Adaptive refinement reduces computational complexity.

Results demonstrate improved accuracy and efficiency.

Abstract

A posteriori residual and hierarchical upper bounds for the error estimates were proved when solving the hypersingular integral equation on the unit sphere by using the Galerkin method with spherical splines. Based on these a posteriori error estimates, adaptive mesh refining procedures are used to reduce complexity and computational cost of the discrete problems. Numerical experiments illustrate our theoretical results.

Tables

Table 1. Table 1 : Errors vs degrees of freedom for f 1 subscript 𝑓 1 f_{1}

Uniform		Residual		hierarchical
DoFs	Error	DoFs	Error	DoFs	Error
6	0.77566	6	0.77566	6	0.77566
18	0.38229	26	0.43544	14	0.68900
66	0.16686	78	0.07714	95	0.18822
258	0.09537	102	0.04493	119	0.07424
1026	0.05792	128	0.03864	141	0.04222
4098	0.03564	211	0.03495	170	0.03574

Table 2. Table 2 : Degrees of freedom and accumulating computation time for f 1 subscript 𝑓 1 f_{1}

Uniform		Residual		hierarchical
DoFs	Comp. time	DoFs	Comp. time	DoFs	Comp. time
6	1.58	6	1.58	6	2.54
18	7.09	26	11.07	14	9.60
66	30.12	78	53.41	95	125.18
258	192.91	102	91.39	119	245.08
1026	2654.11	128	144.25	141	401.22
4098	38754.89	211	259.89	170	612.70

Table 3. Table 3 : Errors vs degrees of freedom for f 2 subscript 𝑓 2 f_{2}

Uniform		Residual		hierarchical
DoFs	Error	DoFs	Error	DoFs	Error
6	0.78050	6	0.78050	6	0.78050
18	0.36153	40	0.38340	54	0.38262
66	0.15705	151	0.06762	153	0.16873
258	0.09356	199	0.04232	199	0.06693
1026	0.05826	253	0.03668	247	0.04151
4098	0.03682	448	0.03269	302	0.03606

Table 4. Table 4 : Degrees of freedom and accumulating computation time for f 2 subscript 𝑓 2 f_{2}

Uniform		Residual		hierarchical
DoFs	Comp. time	DoFs	Comp. time	DoFs	Comp. time
6	1.67	6	2.01	6	3.68
18	7.49	40	27.59	54	88.60
66	31.44	151	176.24	153	346.21
258	184.11	199	311.59	199	722.11
1026	2421.76	253	509.00	247	1242.09
4098	35351.71	448	1051.12	302	1968.70

Equations

$$-Nu+{\omega }^2{\int }_{\mathbb{S}}ud\sigma =f\text{on }\mathbb{S},$$

$$Nv(\mathbf{x}):=\frac{1}{4\pi }\frac{\partial }{\partial {\nu }_{\mathbf{x}}}{\int }_{\mathbb{S}}v(\mathbf{y})\frac{\partial }{\partial {\nu }_{\mathbf{y}}}\frac{1}{\left|\mathbf{x}-\mathbf{y}\right|}d{\sigma }_{\mathbf{y}},$$

$${\Vert u-u_{\Delta }\Vert }_{H^{1/2}(\mathbb{S})}\le C{h}_{\Delta }^{s-1/2}{\Vert u\Vert }_{H^s(\mathbb{S})},$$

$${\Vert u-u_{\Delta }\Vert }_{H^s(\mathbb{S})}\le C{\left(\sum _{\tau \in \Delta }{h}_{\tau }^{2-2s}{\Vert f+N{u}_{\Delta }-{\omega }^2⟨u_{\Delta },1⟩\Vert }_{L_2(\tau )}^2\right)}^{1/2},$$

$${\Vert u-u_{\Delta }\Vert }_{H^{1/2}(\mathbb{S})}^2\le C\sum _{\tau \in \Delta }\sum _{\genfrac{}{}{0ex}{}{{\mathbf{v}}_i\in V_{{\Delta }^{\prime }}}{{\mathbf{v}}_i\in \tau }}{\left(\frac{⟨f+N{u}_{\Delta }-{\omega }^2⟨u_{\Delta },1⟩,B_{{\mathbf{v}}_i}^{\prime }⟩}{{\Vert B_{{\mathbf{v}}_i}^{\prime }\Vert }_{H^{1/2}(\mathbb{S})}}\right)}^2,$$

$$T=\{\mathbf{v}\in {\mathbb{R}}^3:\mathbf{v}=b_1{\mathbf{v}}_1+b_2{\mathbf{v}}_2+b_3{\mathbf{v}}_3\mathrm{with}{b}_i\ge 0,i=1,2,3\}.$$

$$S_d^r(\Delta )=\{s\in C^r(\mathbb{S}):s{|}_{\tau }\in {\Pi }_d,\tau \in \Delta \}.$$

$$\{\begin{array}{ll}d\ge 3r+2& \text{if }r\ge 1\\ d\ge 1& \text{if }r=0\end{array}$$

$$\mathbf{v}=b_{1,\tau }(\mathbf{v}){\mathbf{v}}_1+b_{2,\tau }(\mathbf{v}){\mathbf{v}}_2+b_{3,\tau }(\mathbf{v}){\mathbf{v}}_3.$$

$$\left(\begin{array}{ccc}{v}_1^x& v_2^x& v_3^x\\ v_1^y& v_2^y& v_3^y\\ v_1^z& v_2^z& v_3^z\end{array}\right)\left(\begin{array}{c}{b}_{1,\tau }\\ b_{2,\tau }\\ b_{3,\tau }\end{array}\right)=\left(\begin{array}{c}{v}^x\\ v^y\\ v^z\end{array}\right).$$

$$b_{1,\tau }(\mathbf{v})=\frac{\det \text{/}(\mathbf{v},{\mathbf{v}}_2,{\mathbf{v}}_3)}{\det \text{/}({\mathbf{v}}_1,{\mathbf{v}}_2,{\mathbf{v}}_3)},b_{2,\tau }(\mathbf{v})=\frac{\det \text{/}({\mathbf{v}}_1,\mathbf{v},{\mathbf{v}}_3)}{\det \text{/}({\mathbf{v}}_1,{\mathbf{v}}_2,{\mathbf{v}}_3)},b_{3,\tau }(\mathbf{v})=\frac{\det \text{/}({\mathbf{v}}_1,{\mathbf{v}}_2,\mathbf{v})}{\det \text{/}({\mathbf{v}}_1,{\mathbf{v}}_2,{\mathbf{v}}_3)},$$

$$\det \text{/}({\mathbf{v}}_1,{\mathbf{v}}_2,{\mathbf{v}}_3):=\det \text{/}\left(\begin{array}{ccc}{v}_1^x& v_2^x& v_3^x\\ v_1^y& v_2^y& v_3^y\\ v_1^z& v_2^z& v_3^z\end{array}\right).$$

$$B_{ijk}^{d,\tau }(\mathbf{v})=\frac{d!}{i!j!k!}{b}_{1,\tau }(\mathbf{v}{)}^i{b}_{2,\tau }(\mathbf{v}{)}^j{b}_{3,\tau }(\mathbf{v}{)}^k,i+j+k=d.$$

$$C(\mathbf{x},R)=\{\mathbf{y}\in \mathbb{S}:{\cos }^{-1}(\mathbf{x}\cdot \mathbf{y})\le R\}.$$

$$|\Delta |=\max \{|\tau |:\tau \in \Delta \}\text{and}{\rho }_{\Delta }=\min \{{\rho }_{\tau }:\tau \in \Delta \},$$

$$\left|\tau \right|\le \beta {\rho }_{\tau }\forall \tau \in \Delta$$

$$\left|\tau \right|\ge \gamma {\left|\Delta \right|}_{\text{texttimes}}\forall \tau \in \Delta .$$

$$h_{\tau }=\tan \left(\left|\tau \right|/2\right).$$

$${\rho }_{\tau }\le \left|\tau \right|\le 2h_{\tau }\forall \tau \in \Delta .$$

$$h_{\tau }\le {\beta }_1\tan \left(\frac{{\rho }_{\tau }}{2}\right)\text{or}{h}_{\tau }\le {\beta }_2{\rho }_{\tau }\forall \tau \in \Delta$$

$${\beta }_3{h}_{\tau }\le A_{\tau }^{1/2}\le {\beta }_4{h}_{\tau }\forall \tau \in \Delta ,$$

$$h_{\tau }\ge {\gamma }_1\left|\Delta \right|\forall \tau \in \Delta .$$

$$h_{\tau }\ge {\beta }_5\left|{\Omega }_{\tau }\right|\forall \tau \in \Delta ,$$

$$h_{\Delta }=\tan (\left|\Delta \right|/2).$$

$$\text{card}\left(T_{{\mathbf{v}}_i}^{\Delta }\right)\le L\forall {\mathbf{v}}_i\in V_{\Delta }.$$

$$H^{s}(\mathbb{S}):=\Big{\{}v\in{\mathcal{D}}^{\prime}(\mathbb{S}):\sum_{\ell=0}^{\infty}\,\sum_{m=-\ell}^{\ell}(\ell+1)^{2s}|\widehat{v}_{\ell,m}|^{2}<\infty\Big{\}},$$

$${\hat{v}}_{\ell ,m}={\int }_{\mathbb{S}}v(\mathbf{x})Y_{\ell ,m}(\mathbf{x})d{\sigma }_{\mathbf{x}}.$$

$$\Vert v{\Vert }_{H^s(\mathbb{S})}:={\left(\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }(\ell +1{)}^{2s}|{\hat{v}}_{\ell ,m}{|}^2\right)}^{1/2}$$

$${⟨v,w⟩}_{H^s(\mathbb{S})}:=\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }(\ell +1{)}^{2s}{\hat{v}}_{\ell ,m}{\hat{w}}_{\ell ,m}.$$

$$|{⟨v,w⟩}_{H^s(\mathbb{S})}|\le {\Vert v\Vert }_{H^s(\mathbb{S})}{\Vert w\Vert }_{H^s(\mathbb{S})}\forall v,w\in H^s(\mathbb{S}),\forall s\in \mathbb{R},$$

Full text

Keywords: Hypersingular integral equation; spherical spline; a posteriori error estimate; adaptivity.

AMS Subject Classification: 65N30, 65N38, 65N15, 65N50

1 Introduction

Hypersingular integral equations have many applications, for example in acoustics, fluid mechanics, elasticity and fracture mechanics [13]. These equations arise from the boundary-integral reformulation of the Neumann problem with the Laplacian in a bounded or unbounded domain, see e.g. [22, 40]. In this paper, we study the hypersingular integral equation on the unit sphere

[TABLE]

where $N$ is the hypersingular integral operator given by

[TABLE]

$\omega$ is some nonzero real constant, and $\mathbb{S}$ is the unit sphere in $\mathbb{R}^{3}$, that is, $\mathbb{S}=\{{\boldsymbol{x}}\in\mathbb{R}^{3}:\left|{\boldsymbol{x}}\right|=1\}$. Here $\partial/\partial\nu_{\boldsymbol{x}}$ is the normal derivative with respect to $\boldsymbol{x}$, and $\left|{\cdot}\right|$ denotes the Euclidean norm. The hypersingular integral equation on the unit sphere has applications in geophysics where people are solving Neumann problems in the interior or exterior of the surface of the Earth, see e.g. [18, 19, 34, 41, 42]. Efficient solutions to the hypersingular integral equation on the sphere become more demanding when given data are collected by satellites.

The equation (1.1) can be solved by using tensor products of univariate splines on regular grids which do not exist when the data is given by satellites. Spherical radial basis functions appear to be more suitable for solving problems with scattered data, see e.g. [29, 32, 37, 42] and references therein. However, the resulting matrix system from this approximation is very ill-conditioned. Even though overlapping additive Schwarz preconditioners can be designed for this problem, the condition number of the preconditioned system still depends on the number of subdomains and the angles between subspaces; see [43].

The space of spherical splines defined on a spherical triangulation seems particularly appropriate for use on the sphere [1, 2]. It consists of functions whose pieces are spherical homogeneous polynomials joined together with global smoothness, and thus has both the smoothness and high degree of flexibility [17]. That flexibility makes spherical splines become a powerful tool. These splines have been used successfully in interpolation and data approximation on spheres, see [3, 33]. In an attempt to use spherical splines in solving partial differential equations, Baramidze and Lai [5] use these functions to solve the Laplace–Beltrami equation on the unit sphere. Later, Pham et al. use spherical splines to solve pseudodifferential equations on the unit sphere [38]. The use of spherical splines has some significant advantages. One of them is the ability to write the approximate solutions of the equations in the form of linear combinations of Bernstein–Bézier polynomials which play an extremely important role in computer aided geometric design, data fitting and interpolation, computer vision and elsewhere; see e.g. [16, 21]. Another advantage is the ability to control the smoothness of a function and its derivatives across edges of the triangulations; see [1].

In this paper, the hypersingular integral equation (1.1) will be solved by using the Galerkin method with spherical splines. The linear system arising when solving this equation by using spherical splines is also ill-conditioned. However, preconditioners can be used to tackle this problem, see [36]. When solving the hypersingular integral equation (1.1) by using the Galerkin method with spherical splines associated with a regular and quasi-uniform spherical triangulation $\Delta$, an a priori error estimate is proved as follows

[TABLE]

see Theorem 5.1 in [38]. Here, $s$ is any real number satisfying $1/2\leq s\leq d+1$ where $d$ is the degree of spherical splines, and $C$ is a constant which is independent of the mesh size $h_{\Delta}$ and the exact (unknown) solution $u$. The a priori error estimate (1.3) reveals the rate of convergence in which the upper bound for the approximation error depends on the mesh size $h_{\Delta}$ and the unknown exact solution. However, the quasi-uniform condition on the mesh suggests that uniform refinements of all spherical triangles must be applied when one wish to improve approximation quality. This may lead to an unnecessary waste of computational efforts since contributions to the total error vary over different regions on the unit sphere.

A posteriori error estimates can provide numerical estimates of accuracy in terms of the source term and discrete solutions. In this paper, we shall prove two kinds of a posteriori upper bounds for the errors when solving the hypersingular integral equation on the unit sphere by using Galerkin method with spherical splines. Firstly, we shall prove an a posteriori residual estimate (see Theorem 3.6),

[TABLE]

where $s\in[0,1/2]$ and $C$ is a positive constant depending only on the smallest angle of $\Delta$. Here, the approximate solution $u_{\Delta}$ is found in the space $S_{d}^{r}(\Delta)$ of spherical splines of order $d$ and smoothness $r$ associated with $\Delta$ where $\Delta$ is a regular spherical triangulation. Secondly, when the approximate solution $u_{\Delta}$ is found in the space of continuous piecewise linear spherical splines, we shall prove another a posteriori error estimate (the hierarchical estimate),

[TABLE]

see Corollary 4.5. Here, $\Delta^{\prime}$ is a fictional refinement of $\Delta$ so that a saturation assumption is satisfied, $V_{\Delta^{\prime}}$ is the set of all vertices of $\Delta^{\prime}$, and $B_{\boldsymbol{v}_{i}}^{\prime}$ are nodal basis functions associated with vertices $\boldsymbol{v}_{i}$ of $\Delta^{\prime}$. Precise definitions of spherical triangulations, spherical splines and their basis functions, and Sobolev spaces defined on the unit sphere $\mathbb{S}$ will be presented in Section 2.

Based on these a posteriori error estimates, (1.4) and (1.5), we use adaptive mesh refinement techniques to create better approximation spaces. This results in a significant reduction in required degrees of freedom and computation time while preserving approximate accuracy.

This improvement is very important when we are solving geophysical problems which require considerably large numbers of data points. Furthermore, although all the results in this paper are established for problems on the unit sphere, they can be extended to more general (but related to the sphere) geometries, such as sphere-like geometries (see e.g. [3, 12, 23, 25]). This possible extension can broaden applications of our research.

The structure of the paper is as follows. In Section 2, we will review spherical splines, introduce the Sobolev spaces on the unit sphere to be used, present the quasi-interpolation operator and the hypersingular integral equation. The proof for an a posteriori residual upper bound for the error estimate is presented in Section 3. In Section 4, hierarchical basis techniques are used to prove a posteriori hierarchical error estimate when solving (1.1) by using continuous piecewise linear spherical splines. In Section 5, we discuss simple adaptive mesh refinement algorithms based on the a posteriori error estimates. The final section (Section 6) presents our numerical experiments which illustrate our theoretical results.

In this paper $C$ and $C_{i}$, for $i=1,\ldots,5$, denote generic constants which may take different values at different occurrences.

2 Preliminaries

In this section, we will first review spherical splines [1, 2, 3] and introduce our functional spaces on the unit sphere $\mathbb{S}\subset\mathbb{R}^{3}$. Then the quasi–interpolation operator and the hypersingular integral equation will be discussed.

2.1 Spherical splines

The trihedron $T$ generated by three linearly independent vectors $\{\boldsymbol{v}_{1},\boldsymbol{v}_{2},\boldsymbol{v}_{3}\}$ in ${\mathbb{R}}^{3}$ is defined by

[TABLE]

The intersection $\tau=T\cap\mathbb{S}$ is called a spherical triangle. Let $\Delta=\{\tau_{i}:i=1,\ldots,\mathcal{T}\}$ be a set of spherical triangles. Then $\Delta$ is called a spherical triangulation of the sphere $\mathbb{S}$ if there hold

• (i)

$\bigcup_{i=1}^{\mathcal{T}}\tau_{i}=\mathbb{S}$,

• (ii)

each pair of distinct triangles in $\Delta$ are either disjoint or share a common vertex or an edge.

Let $\Pi_{d}$ denote the space of trivariate homogeneous polynomials of degree $d$ in ${\mathbb{R}}^{3}$. The space of restrictions on the unit sphere $\mathbb{S}$ of all polynomials in $\Pi_{d}$ is denoted by $\Pi_{d}(\mathbb{S})$. Similarly, we also denote by $\mathcal{P}_{d}$ and $\mathcal{P}_{d}(\mathbb{S})$ the spaces of polynomials of degree $d$ in ${\mathbb{R}}^{3}$ and on $\mathbb{S}$, respectively. We define $S_{d}^{r}(\Delta)$ to be the space of piecewise homogeneous splines of degree $d$ and smoothness $r$ on a spherical triangulation $\Delta$, that is,

[TABLE]

Throughout this paper, we always assume that

[TABLE]

holds; see [1, 2, 3].

For a spherical triangle $\tau$ with vertices ${\boldsymbol{v}}_{1},{\boldsymbol{v}}_{2}$, and ${\boldsymbol{v}}_{3}$, let $b_{1,\tau}({\boldsymbol{v}}),b_{2,\tau}({\boldsymbol{v}})$, and $b_{3,\tau}({\boldsymbol{v}})$ denote the spherical barycentric coordinates as functions of ${\boldsymbol{v}}$ in $\tau$, i.e.,

[TABLE]

Suppose that $\boldsymbol{v}_{i}=(v_{i}^{x},v_{i}^{y},v_{i}^{z})$ for $i=1,2,3$ and $\boldsymbol{v}=(v^{x},v^{y},v^{z})$. Equation (2.2) defining the coordinates $b_{i,\tau}$, for $i=1,2,3$, can be written as a system of three linear equations

[TABLE]

Using Cramer’s rule, we have

[TABLE]

where

[TABLE]

We define the homogeneous Bernstein basis polynomials of degree $d$ relative to $\tau$ to be the polynomials

[TABLE]

As was shown in [1], we can use these polynomials as a basis for $\Pi_{d}$.

A spherical cap centred at ${\boldsymbol{x}}\in\mathbb{S}$ and having radius $R$ is defined by

[TABLE]

For any spherical triangle $\tau$, let $|\tau|$ denote the diameter of the smallest spherical cap containing $\tau$, and $\rho_{\tau}$ denote the diameter of the largest spherical cap contained in $\tau$. We define

[TABLE]

and refer to $|\Delta|$ as the mesh size. Our triangulations are said to be regular if for some given $\beta>1$, there holds

[TABLE]

and quasi-uniform if for some given positive number $\gamma<1$, there holds

[TABLE]

Roughly speaking, the regularity guarantees the smallest angles in our triangulations are sufficiently large so that there are no too narrow triangles and the quasi-uniformity guarantees that the sizes of triangles in a triangulation are not too much different.

To accompany the results used in [5, 33, 38] we also denote

[TABLE]

It is obvious that

[TABLE]

Noting (2.6) and (2.8), the regularity of a set of triangulations can also be written by

[TABLE]

for some positive numbers $\beta_{1}$ and $\beta_{2}$. For any $\tau\in\Delta$, we denote by $A_{\tau}$ the area of $\tau$. If $\Delta$ is regular, there holds

[TABLE]

for some positive constants $\beta_{3}$ and $\beta_{4}$. Similarly, the quasi-uniformity can be written as

[TABLE]

For any $\tau\in\Delta$, we denote $\Omega_{\tau}$ to be the union of all triangles in $\Delta$ which share with $\tau$ at least a common vertex or a common edge. If the triangulations $\Delta$ are regular, there holds

[TABLE]

for some $\beta_{5}>0$, see [24, Lemma 4.14]. We denote by $h_{\Delta}$ the mesh size of $\Delta$, i.e.,

[TABLE]

We denote by $V_{\Delta}$ the set of all vertices of the spherical triangulation $\Delta$. Let $\boldsymbol{v}_{i}\in V_{\Delta}$. We also denote by $T_{\boldsymbol{v}_{i}}^{\Delta}$ the set of triangles in $\Delta$ whose one of their vertices is $\boldsymbol{v}_{i}$. If $\Delta$ is regular, the smallest angle in $\Delta$ is bounded below. This suggests that the numbers of spherical triangles which share a common vertex is bounded, i.e., there is a positive integer $L$ (depending only on the smallest angle of $\Delta$) such that

[TABLE]

2.2 Sobolev spaces

For every $s\in{\mathbb{R}}$, the Sobolev space $H^{s}(\mathbb{S})$ defined on the whole unit sphere $\mathbb{S}$ can be defined by using Fourier expansion with spherical harmonics. A spherical harmonic of order $\ell$ on $\mathbb{S}$ is the restriction to $\mathbb{S}$ of a homogeneous harmonic polynomial of degree $\ell$ in ${\mathbb{R}}^{3}$. The space of all spherical harmonics of order $\ell$ is the eigenspace of the Laplace–Beltrami operator $\Delta_{\mathbb{S}}$ corresponding to the eigenvalue $\lambda_{\ell}=-\ell(\ell+1)$. The dimension of this space being $2\ell+1$, see e.g. [30], one may choose for it an orthonormal basis $\{Y_{\ell,m}\}_{m=-\ell}^{\ell}$. The collection of all the spherical harmonics $Y_{\ell,m}$, $m=-\ell,\ldots,\ell$ and $\ell=0,1,\ldots$, forms an orthonormal basis for $L_{2}(\mathbb{S})$. The Sobolev space $H^{s}(\mathbb{S})$ is defined as usual by

[TABLE]

where ${\mathcal{D}}^{\prime}(\mathbb{S})$ is the space of distributions on $\mathbb{S}$ and $\widehat{v}_{\ell,m}$ are the Fourier coefficients of $v$,

[TABLE]

The space $H^{s}(\mathbb{S})$ is equipped with the following norm and inner product:

[TABLE]

and

[TABLE]

When $s=0$ we write $\left\langle{\cdot},{\cdot}\right\rangle$ instead of $\left\langle{\cdot},{\cdot}\right\rangle_{H^{0}(\mathbb{S})}$; this is in fact the $L_{2}$-inner product. We note that

[TABLE]

and

[TABLE]

In particular, there holds

[TABLE]

In the case $k$ belongs to the set of nonnegative integers ${\mathbb{Z}}^{+}$, the Sobolev space $H^{k}(\Omega)$ on a subset $\Omega\subset\mathbb{S}$ can be defined by using an atlas for the unit sphere $\mathbb{S}$ [33]. Let $\{(\Gamma_{j},\phi_{j})\}_{j=1}^{J}$ be an atlas for $\Omega$, i.e, a finite collection of charts $(\Gamma_{j},\phi_{j})$, where $\Gamma_{j}$ are open subsets of $\Omega$, covering $\Omega$, and where $\phi_{j}:\Gamma_{j}\rightarrow B_{j}$ are infinitely differentiable mappings whose inverses $\phi_{j}^{-1}$ are also infinitely differentiable. Here $B_{j}$, $j=1,\ldots,J$, are open subsets in ${\mathbb{R}}^{2}$. Also, let $\{\psi_{j}\}_{j=1}^{J}$ be a partition of unity subordinate to the atlas $\{(\Gamma_{j},\phi_{j})\}_{j=1}^{J}$, i.e., a set of infinitely differentiable functions $\alpha_{j}$ on $\Omega$ vanishing outside the sets $\Gamma_{j}$, such that $\sum_{j=1}^{J}\psi_{j}=1$ on $\Omega$. For any $k\in{\mathbb{Z}}^{+}$, the Sobolev space $H^{k}(\Omega)$ on the unit sphere is defined as follows

[TABLE]

which is equipped with a norm defined by

[TABLE]

Here, $\left\|{\cdot}\right\|_{H^{k}(B_{j})}$ denotes the usual $H^{k}$-Sobolev norm defined on the subset $B_{j}$ of the plane ${\mathbb{R}}^{2}$. In the case $\Omega=\mathbb{S}$, this norm is equivalent to the norm defined in (2.17); see [26].

To accompany the results used in [5, 33, 38], we also present here a definition of Sobolev spaces defined on a subset of $\mathbb{S}$ by using homogeneous extensions of a function defined on $\mathbb{S}$. Let $\ell\in{\mathbb{N}}$ and let $v$ be a function defined on the unit sphere $\mathbb{S}$. We denote by $v_{\ell}$ the homogeneous extension of degree $\ell$ of $v$ to ${\mathbb{R}}^{3}$, i.e.,

[TABLE]

For every $v\in H^{k}(\Omega)$, we define Sobolev–type seminorms of $v$ by

[TABLE]

Here $\left\|{D^{\alpha}v_{\ell-1}}\right\|_{L_{2}(\Omega)}$ is understood as the $L_{2}$-norm of the restriction of the trivariate function $D^{\alpha}v_{\ell-1}$ to $\Omega$. When $\ell=0$ we define

[TABLE]

which can now be used together with (2.23) to define a norm in $H^{k}(\Omega)$:

[TABLE]

This norm is equivalent to the norm $\left\|{\cdot}\right\|_{H^{k}(\Omega)}^{*}$ defined by (2.22); see [33].

For every $s\in[0,1]$, the spaces $\widetilde{H}^{s}(\Omega)$ and $H^{s}(\Omega)$ are defined by Hilbert space interpolation [6] so that

[TABLE]

where $H_{0}^{1}(\Omega)=\left\{v\in H^{1}(\Omega):v=0\ \text{on }\partial\Omega\right\}$, and $[X_{0},X_{1}]_{s}$ denotes the $L_{2}$-interpolation of $X_{0}$ and $X_{1}$, see e.g. [6, 28]. Here, $H_{0}^{1}(\mathbb{S})$ is the space of all functions in $H^{1}(\mathbb{S})$ which vanish on the boundary $\partial\Omega$ of $\Omega$, i.e.,

[TABLE]

The spaces $H^{-s}(\Omega)$ and $\widetilde{H}^{-s}(\Omega)$ are defined as the dual spaces of $\widetilde{H}^{s}(\Omega)$ and $H^{s}(\Omega)$, respectively, with respect to the duality pairing which is the usual extension of the $L_{2}$-inner product on $\Omega$. In particular, the space $H^{-s}(\mathbb{S})$ is defined to be the dual space of $H^{s}(\mathbb{S})$. The $\left\|{\cdot}\right\|_{H^{s}(\mathbb{S})}$-norm defined by (2.17) turns out to be equivalent to the $\left\|{\cdot}\right\|_{H^{s}(\mathbb{S})}^{\prime}$-norm defined by (2.22), (2.25) and (2.20) when $\Omega=\mathbb{S}$ and $-1\leq s\leq 1$, i.e.,

[TABLE]

for some positive numbers $\gamma_{2}$ and $\gamma_{3}$, see e.g. [20, 26, 33, 34].

2.3 Quasi-Interpolation

We now briefly discuss the construction of a quasi-interpolation operator $Q:L_{2}(\mathbb{S})\rightarrow S_{d}^{r}(\Delta)$ which is defined in [33]. Firstly, we introduce the set of domain points of $\Delta$ to be

[TABLE]

Here, $\tau=\langle\boldsymbol{v}_{1},\boldsymbol{v}_{2},\boldsymbol{v}_{3}\rangle$ denotes the spherical triangle whose vertices are $\boldsymbol{v}_{1},\boldsymbol{v}_{2},\boldsymbol{v}_{3}$. We denote the domain points by $\xi_{1},\ldots,\xi_{D}$, where $D={{\rm dim}}\,S_{d}^{0}(\Delta)$. Let $\left\{B_{\ell}:\ell=1,\ldots,D\right\}$ be a basis for $S_{d}^{0}(\Delta)$ such that the restriction of $B_{\ell}$ on the triangle containing $\xi_{\ell}$ is $Bernstein$ polynomial of degree $d$ associated with this point, and that $B_{\ell}$ vanishes on other triangles.

A set $\mathcal{M}=\left\{\zeta_{\ell}\right\}_{\ell=1}^{M}\subset{\mathcal{D}}$ is called a minimal determining set for $S_{d}^{r}(\Delta)$ if, for every $s\in S_{d}^{r}(\Delta)$, all the coefficients $\nu_{\ell}(s)$ in the expression $s=\sum_{\ell=1}^{D}\nu_{\ell}(s)B_{\ell}$ are uniquely determined by the coefficients corresponding to the basis functions which are associated with points in $\mathcal{M}$. Given a minimal determining set, we construct a basis $\left\{B^{*}_{\ell}\right\}_{\ell=1}^{M}$ for $S_{d}^{r}(\Delta)$ by requiring

[TABLE]

By using Hahn-Banach theorem we extend the linear functions $\nu_{\ell}$, $\ell=1,\ldots,M$, to all of $L_{2}(\mathbb{S})$. We continue to use the same symbol for the extensions. The quasi-interpolation operator: $Q:L_{2}(\mathbb{S})\rightarrow S_{d}^{r}(\Delta)$ is now defined by

[TABLE]

2.4 The hypersingular integral equation

The hypersingular integral operator (1.2) arises from the boundary-integral reformulation of the Neumann problem with the Laplacian in the interior or the exterior of the sphere, see [41]. This operator (with minus sign) turns out to be a strongly elliptic pseudodifferential operator of order $1$, see e.g. [41, 38], i.e.

[TABLE]

In this paper, we solve the hypersingular integral equation (1.1),

[TABLE]

where $f\in H^{-1/2}(\mathbb{S})$. We denote by ${\cal N}^{*}:H^{s}(\mathbb{S})\rightarrow H^{s-1}(\mathbb{S})$ the operator which is given by

[TABLE]

Noting (2.17), (2.29) and (2.31), we have

[TABLE]

For every $\ell\geq 1$, there holds

[TABLE]

This together with (2.17) and (2.32) implies

[TABLE]

where

[TABLE]

To set up a weak formulation, we introduce the bilinear form

[TABLE]

where $\left\langle{v},{w}\right\rangle$ is the $H^{1/2}(\mathbb{S})$-duality pairing which coincides with the $L_{2}(\mathbb{S})$-inner product when $v$ and $w$ belong to $L_{2}(\mathbb{S})$. This bilinear form is clearly bounded and coercive, i.e.,

[TABLE]

and

[TABLE]

respectively. A natural weak formulation of equation (1.1) is: Find $u\in H^{1/2}(\mathbb{S})$ satisfying

[TABLE]

Let $\Delta$ be a spherical triangulation on $\mathbb{S}$. We denote by $u_{\Delta}\in S_{d}^{r}(\Delta)$ the Galerkin solution

[TABLE]

The unique existences of $u$ and $u_{\Delta}$ are guaranteed by the Lax–Milgram Theorem, noting the boundedness (2.36) and the coercivity (2.37) of the bilinear form $a(\cdot,\cdot)$. Furthermore, if $\Delta$ is a regular and quasi-uniform triangulation and if $u\in H^{s}(\mathbb{S})$ for some $1/2\leq s\leq d+1$, then there holds

[TABLE]

see [38]. Here, $h_{\Delta}$ is the mesh size of $\Delta$, see (2.14), and $\alpha_{3}$ is a positive constant depending only on $d$ and the smallest angle in $\Delta$. The a priori error estimate (2.40) reveals the convergence and stability of the Galerkin approximation (2.39). However, the upper bound of the error $\left\|{u-u_{\Delta}}\right\|_{H^{1/2}(\mathbb{S})}$ is given by the mesh size $h_{\Delta}$ and the norm $\left\|{u}\right\|_{H^{s}(\mathbb{S})}$ of the exact solution $u$ which is unknown. Furthermore, the quasi-uniform requirement means that one has to divide all spherical triangles in the current mesh whenever better accuracy is demanded. In the next section, we prove a residual upper bound for the error $\left\|{u-u_{\Delta}}\right\|_{H^{1/2}(\mathbb{S})}$ in terms of the given right hand side $f$ and the approximate solutions $u_{\Delta}$ of the corresponding discrete problems.

3 A posteriori residual error

estimate

In this section, the error $\left\|{u-u_{\Delta}}\right\|_{H^{s}(\mathbb{S})}$ will be bounded above by an a posteriori residual error estimator. We assume that $f\in L_{2}(\mathbb{S})$. Since $S_{d}^{r}(\Delta)\subset H^{r+1}(\mathbb{S})$ (see [38]), for each $u_{\Delta}\in S_{d}^{r}(\Delta)$, we have ${\cal N}^{*}u_{\Delta}\in H^{r}(\mathbb{S})\subset L_{2}(\mathbb{S})$. The residual $\mathcal{R}(u_{\Delta})\in L_{2}(\mathbb{S})$ is defined by

[TABLE]

This together with (2.35) and (2.38) gives

[TABLE]

It is obvious from (3.1) that the residual $\mathcal{R}(u_{\Delta})$ depends solely on the source term $f$ and the discrete solution $u_{\Delta}$. The following lemma states the equivalence of the error $\left\|{u-u_{\Delta}}\right\|_{H^{1/2}(\mathbb{S})}$ and the $H^{-1/2}(\mathbb{S})$-norm of the residual $\mathcal{R}(u_{\Delta})$.

Lemma 3.1**.**

Let $u$ and $u_{\Delta}$ be the weak and approximate solutions defined by (2.38) and (2.39), respectively. There holds

[TABLE]

where $\alpha_{1}$ and $\alpha_{2}$ are the coercivity and boundedness constants, see (2.37) and (2.36), respectively.

*Proof.  * Noting (2.20), we have

[TABLE]

It follows from the coercivity (2.37) of the bilinear form $a(\cdot,\cdot)$ and (3.2) that

[TABLE]

This together with (3.3) implies

[TABLE]

On the other hand, we derive

[TABLE]

noting (3.2) and the continuity (2.36) of the bilinear form $a(\cdot,\cdot)$. This implies

[TABLE]

finishing the proof of the lemma. $\Box$

For each $\tau\in\Delta$, we define the spherical triangle residual by

[TABLE]

and the local error estimator $\eta_{\Delta,s}(\tau)$ by

[TABLE]

where $0<s<1$. The residual estimators were used for solving the hypersingular integral equation with flat triangular elements, see [8]. In this paper, the local error estimators are defined on spherical triangles.

Note here that $f$ and ${\cal N}^{*}u_{\Delta}$ belong to $L_{2}(\mathbb{S})$. It follows from (3.1) and (3.4) that for any $v\in H^{1/2}(\mathbb{S})$, there holds

[TABLE]

The following lemma shows an approximation property of the quasi–interpolation operator $Q$ (defined in Subsection 2.3). This result extends Theorem 2 in [5] in which we relax on the quasi-uniform condition of $\Delta$.

Lemma 3.2**.**

Let $m$ be a positive integer satisfying

[TABLE]

Assume that $\Delta$ is a regular spherical triangulation such that $\left|\Omega_{\tau}\right|<1$ for all $\tau\in\Delta$. Recall that $Q:L_{2}(\mathbb{S})\rightarrow S_{d}^{r}(\Delta)$ is the quasi-interpolation operator defined by (2.28). For any $\tau\in\Delta$, if $v\in H^{m}(\Omega_{\tau})$, then there holds

[TABLE]

for all $k=0,\ldots,\min\left\{m-1,r+1\right\}$. Here, $\alpha_{4}$ is a positive constant depending only on $d$ and the smallest angle in $\Delta$.

*Proof.  * Note here that for any $m$ satisfying (3.7), we have $d-(m-1)$ is an even number, and thus $\left|{\boldsymbol{x}}\right|^{d-(m-1)}$, for $\boldsymbol{x}=(x_{1},x_{2},x_{3})$, is a homogeneous polynomial of degree $d-(m-1)$. Furthermore, for any $\boldsymbol{x}\in\mathbb{S}$, we have $\left|{\boldsymbol{x}}\right|^{d-(m-1)}=1$, and thus if $s\in\Pi_{m-1}(\mathbb{S})$, then

[TABLE]

By Theorem 4.2 in [33], for any $v\in H^{m}(\Omega_{\tau})$, there exists a spherical homogeneous polynomial $s\in\Pi_{m-1}(\mathbb{S})\subset\Pi_{d}(\mathbb{S})$ such that

[TABLE]

In particular, when $k=0$ we have

[TABLE]

Since $s$ is a spherical homogeneous polynomial of degree $d$ on $\mathbb{S}$, Lemma 9 in [33] assures that $s=Qs$ and

[TABLE]

This together with (3.10) implies

[TABLE]

Since $s=Qs$, by using the triangle inequality and noting (3.9), (3.12), we obtain

[TABLE]

Since $\Delta$ is regular, the inequality (3.8) is derived from (3.13) and noting (2.10) and (2.13).

$\Box$

The inequality (3.8) in Lemma 3.2 holds for any integer $m$ satisfying (3.7). In the following lemma, the inequality is proved when $k=0$ and $m$ is a real number between [math] and $1$.

Lemma 3.3**.**

Let $\Delta$ be a regular spherical triangulation such that $\left|\Omega_{\tau}\right|<1$ for all $\tau\in\Delta$, and let $Q:L_{2}(\mathbb{S})\rightarrow S_{d}^{r}(\Delta)$ be the quasi-interpolation operator defined by (2.28). For any $v\in H^{s}(\mathbb{S})$ where $0\leq s\leq 1$, there holds

[TABLE]

where $\alpha_{5}$ is a positive constant depending only on the smallest angle of triangles in $\Delta$.

*Proof.  * Using the result in [5, Lemma 9], we have

[TABLE]

where $C$ is a constant that depends only on the smallest angle of $\tau$. This together with the triangle inequality implies

[TABLE]

If $d$ is even, we apply Lemma 3.2 when $k=0$ and $m=1$ to obtain

[TABLE]

Noting (3.15), (3.16) and using [28, Theorem B.2] (for $\theta=s$ where $0\leq s\leq 1$), we obtain

[TABLE]

proving (3.14) when $d$ is even. We now prove (3.14) when $d$ is odd. Applying Lemma 3.2 when $k=0$ and $m=2$ we have

[TABLE]

Noting (3.15), (3.17) and applying [28, Theorem B.2] (for $\theta=s/2$ where $0\leq s\leq 2$) we obtain

[TABLE]

completing the proof of the lemma. $\Box$

Technical results in the following two lemmas will be used in the proof of Theorem 3.6.

Lemma 3.4**.**

Let $\Delta$ be a regular spherical triangulation on the unit sphere. There holds

[TABLE]

where $\alpha_{6}$ is a positive constant which depends only on the smallest angle of $\Delta$.

*Proof.  * Noting (2.15) there holds

[TABLE]

where $L$ is a positive constant depending only on the smallest angle of $\Delta$. If $\boldsymbol{v}_{1}$, $\boldsymbol{v}_{2}$ and $\boldsymbol{v}_{3}$ are the vertices of $\tau$ then $\Omega_{\tau}=\bigcup\left\{\widetilde{\tau}\in T_{\boldsymbol{v}_{i}}^{\Delta},i=1,2,3\right\}$ and thus

[TABLE]

Suppose that $\tau^{\prime}\in\Delta$ satisfies $\operatorname{{int\/}}\Omega_{\tau^{\prime}}\cap\operatorname{{int\/}}\Omega_{\tau}\not=\emptyset$. Then, there is a $\widetilde{\tau}\in\Delta$ such that $\widetilde{\tau}\subset\Omega_{\tau^{\prime}}\cap\Omega_{\tau}$. If $\widetilde{\tau}\subset\Omega_{\tau^{\prime}}$ then $\tau^{\prime}\subset\Omega_{\widetilde{\tau}}$. For every $\tau\in\Delta$, there are at most $3L$ options of choosing a $\widetilde{\tau}\subset\Omega_{\tau}$ by (3.19). On the other hand, for each $\widetilde{\tau}$ in $\Omega_{\tau}$, there are at most $3L$ options of choosing a $\tau^{\prime}\subset\Omega_{\widetilde{\tau}}$. Thus, there holds

[TABLE]

Denoting $\alpha_{6}={9L^{2}}$, we obtain the inequality (3.18), completing the proof of the lemma. $\Box$

Lemma 3.5**.**

Let $\Delta$ be a regular spherical triangulation and let $s\in[0,1]$. There exists a positive number $\alpha_{7}$ which depends only on the smallest angle of $\Delta$ such that

[TABLE]

*Proof.  * Since $\Delta$ is regular, by applying Lemma 3.4, there holds

[TABLE]

The set $\left\{\operatorname{{int\/}}\Omega_{\tau}:\tau\in\Delta\right\}$ is a set of overlapping subsets which covers the unit sphere $\mathbb{S}$. The coloring argument (see e.g. [9]) suggests that the set $\left\{\operatorname{{int\/}}\Omega_{\tau}:\tau\in\Delta\right\}$ can be divided into $C_{1}$ groups

[TABLE]

so that each group consists of mutually disjoint subsets. Here, the constant $C_{1}$ satisfies

[TABLE]

Since $\operatorname{{int\/}}\Omega_{\tau}\cap\operatorname{{int\/}}\Omega_{\tau^{\prime}}=\emptyset$ if $\tau$ and $\tau^{\prime}$ are two triangles that belong to the set $I_{k}$ and $\bigcup\left\{\Omega_{\tau}:\tau\in I_{k}\right\}\subset\mathbb{S}$, there holds

[TABLE]

see [8, 45]. We obtain

[TABLE]

noting (2.26). The inequality (3.20) can then be derived by denoting $\alpha_{7}=C_{1}\gamma_{3}^{2}$, completing the proof of the lemma. $\Box$

Recalling the local error estimator $\eta_{\Delta,s}(\tau)$ (see (3.5)), for a subset $\Omega\subset\mathbb{S}$, we define the error estimator $\eta_{\Delta,s}(\Omega)$ by

[TABLE]

In particular, we denote by $\eta_{\Delta,s}(\mathbb{S})$ the residual-type error estimator with respect to the mesh $\Delta$, i.e.,

[TABLE]

We are now ready to prove the main result of this section. The error $\left\|{u-u_{\Delta}}\right\|_{H^{s}(\mathbb{S})}$ will be bounded above by the residual error estimator $\eta_{\Delta,s}(\mathbb{S})$.

Theorem 3.6** **(A posteriori residual upper

bound).

Let $\Delta$ be a regular spherical triangulation such that $\left|\Omega_{\tau}\right|<1$ for all $\tau\in\Delta$. Let $u$ and $u_{\Delta}$ be the weak and approximate solutions defined by (2.38) and (2.39), respectively. There exists a positive constant $\alpha_{8}$ depending only on the smallest angle of $\Delta$ such that for all $0\leq s\leq 1/2$

[TABLE]

Here, $\alpha_{8}$ is a positive constant depending only on $d$ and the smallest angle of $\Delta$.

*Proof.  * Employing (3.2), (2.38) and (2.39), we derive

[TABLE]

Using the duality argument (2.19) and noting (3.24), we obtain

[TABLE]

Note that $v\in H^{1-s}(\mathbb{S})\subset L_{2}(\mathbb{S})$ for every $s\in[0,1/2]$, and $Qv\in S_{d}^{r}(\Delta)\subset H^{1}(\mathbb{S})$. By (3.6), we have

[TABLE]

where in the second step we apply Cauchy-Schwarz inequality. This together with the result in Lemma 3.3 gives

[TABLE]

By using Cauchy-Schwarz inequality and applying Lemma 3.5, we have

[TABLE]

Noting (3.1), (2.30) and (2.31), we have $\mathcal{R}(u_{\Delta})={\cal N}^{*}(u-u_{\Delta})$. Since $0\leq s\leq 1/2$, we have $u\in H^{1/2}(\mathbb{S})\subset H^{s}(\mathbb{S})$. Applying the inequality (2.33) and noting (3.22) and (3.5), we obtain

[TABLE]

finishing the proof of the theorem. $\Box$

4 An a posteriori hierarchical error

estimation

Hierarchical basis techniques have been used to prove a posteriori error estimates when solving hypersingular integral equation in two dimensions and linear elements, see e.g. [27, 31, 10, 15]. In this section, we discuss the use of these techniques to prove an a posteriori upper bound for the error $\left\|{u-u_{\Delta}}\right\|_{H^{1/2}(\mathbb{S})}$ when solving the hypersingular integral equation on the unit sphere, where the approximate solution $u_{\Delta}$ is found in the space $S_{1}^{0}(\Delta)$ and $\Delta$ is a spherical triangulation on $\mathbb{S}$. In the remainder of this paper, we use $S(\Delta)$ instead of $S_{1}^{0}(\Delta)$ for notational convenience. Suppose that the set $V_{\Delta}=\left\{\boldsymbol{v}_{1},\boldsymbol{v}_{2},\ldots,\boldsymbol{v}_{M}\right\}$ is the set of all vertices of $\Delta$. For each vertex $\boldsymbol{v}_{i}$, the associated basis function $B_{\boldsymbol{v}_{i}}$ is defined by

[TABLE]

where $b_{1,\tau}(\boldsymbol{x})$ is the first spherical barycentric coordinate of $\boldsymbol{x}$ with respect to $\tau$, see (2.2) and (2.3). We then have

[TABLE]

Recalling the definition of the quasi-interpolation operator with respect to the space $S_{d}^{r}(\Delta)$ in Subsection 2.3, the quasi-interpolation operator: $Q:L_{2}(\mathbb{S})\rightarrow S(\Delta)$ is given by

[TABLE]

Here, $\nu_{\boldsymbol{v}_{i}}(v)=v(\boldsymbol{v}_{i})$ for all $v\in S(\Delta)$. The quasi-interpolation operator is a projection onto $S(\Delta)$, i.e. $Q^{2}v=Qv$ for every $v\in L_{2}(\mathbb{S})$. Every $s\in S(\Delta)$ can uniquely be written as

[TABLE]

Lemma 4.1**.**

Let $\Delta$ be a regular spherical triangulation such that $h_{\tau}<1$ for all $\tau\in\Delta$. For every vertex $\boldsymbol{v}$ and any $\tau\in T_{\boldsymbol{v}}^{\Delta}$, the basis function $B_{\boldsymbol{v}}\in S(\Delta)$ associated with $\boldsymbol{v}$ (see (4.1)) satisfies

[TABLE]

where $\alpha_{9}$ is a constant which depends only on the smallest angle of $\Delta$.

*Proof.  * Since $\Delta$ is regular, the cardinality of $T_{\boldsymbol{v}}^{\Delta}$ is bounded, i.e., $\textrm{card}\left(T_{\boldsymbol{v}}^{\Delta}\right)\leq L$ for some positive integer $L$ depending only on the smallest angle of $\Delta$, see (2.15). If $\tau,\tau^{\prime}\in T_{\boldsymbol{v}}^{\Delta}$, then $\tau^{\prime}\subset\Omega_{\tau}$. This together with (2.13) implies

[TABLE]

for some positive constant $C_{1}$ depending only on the smallest angle in $\Delta$. We then have

[TABLE]

Statement (5) in [33, Proposition 5.1] and (2.11) give

[TABLE]

Since $\operatorname{{supp\/}}B_{\boldsymbol{v}}\subset\bigcup\left\{\tau^{\prime}:\tau^{\prime}\in T_{\boldsymbol{v}}^{\Delta}\right\}$ and noting (4.4), we obtain

[TABLE]

Similarly, the inequality (8) in [33, Proposition 5.1] together with (2.10) and (2.11) yields

[TABLE]

Since this is true for all $\tau^{\prime}\in T_{\boldsymbol{v}}^{\Delta}$ and $\operatorname{{supp\/}}B_{\boldsymbol{v}}=\bigcup\left\{\tau^{\prime}:\tau^{\prime}\in T_{\boldsymbol{v}}^{\Delta}\right\}$, we have

[TABLE]

On the other hand, the size $h_{\tau}$ is smaller than $1$ for every $\tau\in\Delta$. This together with (4.5) implies

[TABLE]

This together with (4.6) implies

[TABLE]

Noting (4.5), (4.7) and applying the interpolation inequality (see e.g. [28, Lemma B.1]), we derive

[TABLE]

This together with (2.26) yields (4.3) where $\alpha_{9}=\gamma_{3}\sqrt{C_{3}C_{6}}$, completing the proof of this lemma. $\Box$

A spherical triangulation $\Delta^{\prime}$ is said to be a refinement of another spherical triangulation $\Delta$ if every spherical triangle $\tau^{\prime}\in\Delta^{\prime}$ is a subtriangle of a triangle $\tau\in\Delta$. When $\Delta^{\prime}$ is a refinement of $\Delta$, we call $\Delta$ a coarser triangulation (coarser mesh) and $\Delta^{\prime}$ is a finer triangulation (finer mesh). In this case, the two spherical triangulations are said to be nested.

Suppose that $\Delta$ and $\Delta^{\prime}$ are two nested spherical triangulations, where $\Delta^{\prime}$ is the finer mesh. Then the space $S(\Delta)$ is a subspace of $S(\Delta^{\prime})$. We denote by $u_{\Delta}$ and $u_{\Delta^{\prime}}$ the Galerkin solutions to the hypersingular integral equation (2.30), i.e., $u_{\Delta}\in S(\Delta)$ and $u_{\Delta^{\prime}}\in S(\Delta^{\prime})$ satisfy

[TABLE]

and

[TABLE]

Following e.g. [15, 27, 31], we assume that the two triangulations $\Delta$ and $\Delta^{\prime}$ satisfy the saturation assumption:

[TABLE]

for some fixed $\eta\in(0,1)$. Here, the function $u$ in (4.10) is the weak solution to the hypersingular integral equation defined by (2.38). In our adaptive refinement strategy which will be discussed in Section 5, the approximate solution $u_{\Delta^{\prime}}$ is not computed and the finer mesh $\Delta^{\prime}$ only plays a role as a mean to evaluate only local error estimators which will then be used to conduct mesh refinement step and create better approximation spaces. In our numerical experiments (Section 6), $\Delta^{\prime}$ is created from $\Delta$ by joining midpoints of the three spherical edges in each spherical triangle of $\Delta$.

In this section, for each vertex $\boldsymbol{v}_{i}\in V_{\Delta}$, we denote by $B_{\boldsymbol{v}_{i}}$ the hat function in $S(\Delta)$ corresponding to the vertex $\boldsymbol{v}_{i}$, see (4.1). Since $V_{\Delta}\subset V_{\Delta^{\prime}}$, the vertex $\boldsymbol{v}_{i}$ is also a vertex in the spherical triangulation $\Delta^{\prime}$. If $\boldsymbol{v}_{i}\in V_{\Delta^{\prime}}$, we denote by $B_{\boldsymbol{v}_{i}}^{\prime}$ the hat function in $S(\Delta^{\prime})$ associated with the vertex $\boldsymbol{v}_{i}$. We recall here that $Q_{\Delta}$ and $Q_{\Delta^{\prime}}$ denote the quasi-interpolation operators associated with the spaces $S(\Delta)$ and $S(\Delta^{\prime})$, respectively. For each $\boldsymbol{v}_{i}\in V_{\Delta^{\prime}}$, we define a nodal estimator

[TABLE]

where $\mathcal{R}(u_{\Delta})=f-{\cal N}^{*}u_{\Delta}\in H^{-1/2}(\mathbb{S})$.

Lemma 4.2**.**

Let $\Delta$ and $\Delta^{\prime}$ be two nested spherical triangulations where $\Delta^{\prime}$ is the finer mesh. For every $v\in S(\Delta^{\prime})$, we denote

[TABLE]

Suppose that $v\in S(\Delta^{\prime})$ satisfies $I_{\Delta}v=0$, there holds

[TABLE]

Here, $\nu_{\boldsymbol{v}_{i}}^{\prime}$ is the linear functional which picks the coefficient associated with vertex $\boldsymbol{v}_{i}\in V_{\Delta^{\prime}}$.

*Proof.  * Since $v\in S(\Delta^{\prime})$, $v$ can uniquely be written as

[TABLE]

It follows that

[TABLE]

On the other hand, we have

[TABLE]

This together with (4.15) yields

[TABLE]

Since $I_{\Delta}v=0$, there holds

[TABLE]

This yields

[TABLE]

Equalities (4.14) and (4.16) imply (4.13), completing the proof of this lemma. $\Box$

Lemma 4.3**.**

Let $u_{\Delta}$ and $u_{\Delta^{\prime}}$ be Galerkin solutions defined by (4.8) and (4.9), respectively. Denote $e:=u_{\Delta^{\prime}}-u_{\Delta}$ and $w:=e-I_{\Delta}e$. There holds

[TABLE]

*Proof.  * Recall that $u_{\Delta}$ and $u_{\Delta^{\prime}}$ are the Galerkin solutions in the spaces $S(\Delta)$ and $S(\Delta^{\prime})$, respectively. Noting (4.8), (4.9) and (2.35), we have

[TABLE]

and

[TABLE]

Noting that $I_{\Delta}$ is a projection, i.e., $(I_{\Delta})^{2}=I_{\Delta}$, we obtain

[TABLE]

Noting that $S(\Delta)\subset S(\Delta^{\prime})$, we have $e=u_{\Delta^{\prime}}-u_{\Delta}\in S(\Delta^{\prime})$ and $w=e-I_{\Delta}e\in S(\Delta^{\prime})$. This together with (4.19) and the result in Lemma 4.2 implies

[TABLE]

By (4.17) and (4.18), we have

[TABLE]

noting that $I_{\Delta}e\in S(\Delta)\subset S(\Delta^{\prime})$. It follows from (4.21) and (4.20) that

[TABLE]

By the definition of $e$ and by using (4.18) (noting that $B_{\boldsymbol{v}_{i}}^{\prime}\in S(\Delta^{\prime})$), we obtain

[TABLE]

completing the proof of this lemma.

$\Box$

We are now ready to prove the main theorem of this section, an upper bound for the error $\left\|{u-u_{\Delta}}\right\|_{H^{1/2}(\mathbb{S})}$ in terms of the error estimators $\mu_{\boldsymbol{v}_{i}}$, see (4.11).

Theorem 4.4** (A posteriori hierarchical upper bound).**

Let $\Delta$ and $\Delta^{\prime}$ be two nested spherical triangulations (where $\Delta^{\prime}$ is the finer mesh) satisfying the saturation assumption (4.10). There exists a positive number $\alpha_{10}$ depending only on the smallest angle of the triangulations and the saturation assumption constant $\eta$ (see (4.10)) such that

[TABLE]

where $\mu_{\boldsymbol{v}_{i}}$ are the nodal estimators defined by (4.11).

*Proof.  * The triangle inequality and the saturation assumption (4.10) give

[TABLE]

It follows that

[TABLE]

Suppose that $e=u_{\Delta^{\prime}}-u_{\Delta}$ and $w=e-I_{\Delta}e$ as defined in Lemma 4.3. Then we have

[TABLE]

Applying Statement (4) in [33, Proposition 5.1] and (2.11), there exists a constant $C_{1}>0$ depending only on the smallest angle in $\Delta^{\prime}$ such that

[TABLE]

for every vertex $\boldsymbol{v}_{i}\in V_{\Delta^{\prime}}$ and for every $\tau_{i}\in\Delta^{\prime}$. Using (4.25) and the triangle inequality, we obtain

[TABLE]

noting that $Q_{\Delta}e=I_{\Delta}(Q_{\Delta}e)$. It follows from (4.12), (4.14) and the triangle inequality that

[TABLE]

Applying Proposition 5.1 (statement (4)) in [33] and (2.11) again, we have

[TABLE]

Statement (5) in [33, Proposition 5.1] and (2.11) give

[TABLE]

for some positive number $C_{2}>0$ depending only on the smallest angle of $\Delta^{\prime}$. The inequalities (4.26)–(4.29) and the result in Lemma 3.3 yield

[TABLE]

where $C_{3}=C_{1}\left(1+3C_{1}C_{2}\right)\alpha_{5}$. It follows from (4.24), the triangle inequality, (4.30) and the Cauchy–Schwarz inequality that

[TABLE]

Applying Lemma 4.1, we obtain

[TABLE]

We note that each $\tau_{i}$ can be chosen by at most three vertices (its vertices). Therefore, we have

[TABLE]

By applying the result in Lemma 3.5, we obtain

[TABLE]

It follows from (4.31)–(4.34) that

[TABLE]

This together with (2.35), (2.37) and (4.11) yields

[TABLE]

Noting that $e=u_{\Delta^{\prime}}-u_{\Delta}$ and (4.23) we obtain

[TABLE]

The desired inequality (4.22) can then be obtained by denoting

[TABLE]

completing the proof of the theorem. $\Box$

In Theorem 4.4, the error $\left\|{u-u_{\Delta}}\right\|_{H^{1/2}(\mathbb{S})}$ is bounded above by the sum of nodal estimators. For refinement purpose, the a posteriori error estimate can also be written in the the form of element estimators as in the following corollary.

Corollary 4.5**.**

Let all assumptions in Theorem 4.4 be satisfied. Then there holds

[TABLE]

where

[TABLE]

5 Mesh Refinement

In this section, we briefly discuss the mesh refinement technique that will be used to refine our spherical triangulations. The technique is based on the a posteriori error estimates proved in Theorems 3.6 and 4.4, and Corollary 4.5. Borrowing existing ideas in planar cases, see e.g. [4, 7, 8, 11, 35, 39, 44], our mesh refinement algorithms consist of two subroutines. One is constructing the indicators from the error estimators. The other is defining the rules that are used to divide the triangles. Here, indicator constructions are different for the two adaptive approaches which are based on the residual and the hierarchical estimates. Meanwhile, we use the same rule to divide the triangles for both adaptive procedures.

Residual adaptive approach: Starting with a spherical triangulation $\Delta_{k}$, we denote by $\widehat{\Delta}_{k}$ the subset of $\Delta_{k}$ containing all spherical triangles that will be refined. This can be achieved with the following marking strategy (see [14]):

Strategy: Given a parameter $0<\xi<1$, construct a minimal subset $\widehat{\Delta}_{k}$ of $\Delta_{k}$ such that

[TABLE]

*and mark all spherical triangles in $\widehat{\Delta}_{k}$ for refinement. Here, recall that $\eta_{\Delta_{k},1/2}(\tau)$ is defined by (3.5). *

Hierarchical adaptive approach: Starting with a spherical triangulation $\Delta_{k}$, we denote by $\Delta_{k}^{\prime}$ the finer mesh of $\Delta_{k}$ which is created by joining the midpoints of the three edges of all triangles in $\Delta_{k}$, see Figure 2. Note here that we only need the vertices of $\Delta_{k}^{\prime}$ in order to compute the nodal estimators

[TABLE]

see (4.11). The mesh $\Delta_{k}^{\prime}$ is not at all the finer mesh that we use to create approximation spaces. For each $\tau$ in $\Delta_{k}$, the local error estimator is computed by

[TABLE]

see (4.36). The subset $\widehat{\Delta}_{k}$ of spherical triangles in $\Delta_{k}$ which will be marked for refinement is determined by applying the above strategy:

Given a parameter $0<\xi<1$, construct a minimal subset $\widehat{\Delta}_{k}$ of $\Delta_{k}$ such that

[TABLE]

*and mark all spherical triangles in $\widehat{\Delta}_{k}$ for refinement. *

Once, the subset $\widehat{\Delta}_{k}$ of spherical triangles in $\Delta_{k}$ that are to be divided is obtained, mesh refinement techniques are then applied. When it comes to the mesh refinement, algorithms for cutting triangles in triangulations have been extensively discussed in [39]. These algorithms are based on the bisection of triangles by dividing the longest edges so that the following features are satisfied. Let $\Delta_{k}$ be a conforming triangulation, i.e. the intersection of two non-disjoint, nonidentical triangles is either a common vertex or common edge. With any refinement submesh $\widehat{\Delta}_{k}\in\Delta_{k}$, the algorithm produces a new conforming triangulation $\Delta_{k+1}$ with the following properties:

1. (i)

all elements of $\widehat{\Delta}_{k}$ are refined to create new elements in $\Delta_{k+1}$, 2. (ii)

$\Delta_{k+1}$ is nested in $\Delta_{k}$ in such a way that each refined triangle is embedded in one triangle of $\Delta_{k}$, 3. (iii)

$\Delta_{k+1}$ is non-degenerated, i.e. the interior angles of all triangles of $\Delta_{k+1}$ are guaranteed to be bounded away from 0, 4. (iv)

the transition between large and small triangles is not abrupt.

Following [44], the below steps are used to produce a totally refined and conforming triangulation $\Delta_{k+1}$ in the following way:

• Step 1:

Separate all $\tau$ in $\widehat{\Delta}_{k}$ into 4 pieces to obtain $\widetilde{\Delta_{k}}$, see Figure 1(a).

• Step 2:

Find all hanging nodes in $\widetilde{\Delta_{k}}$ and verify if each of these hanging nodes lies on the longest edge of a triangle or not.

• –

If the hanging node lies on the longest edge, join it with the opposite vertex to obtain 2 new triangles, see Figure 1(b).

• –

If the hanging node does not lie on the longest edge, join it with the middle point of the longest edge, together with joining the middle point of the longest edge with its opposite vertex to obtain 3 new triangles, see Figure 1(c).

6 Numerical Experiments

We consider the exterior Neumann problem

[TABLE]

where the boundary data $Z_{N}$ is one of the following functions

[TABLE]

and

[TABLE]

where $\boldsymbol{p}=(0,0,0.95)$ and $\boldsymbol{q}=(0,0,-0.95)$. Solving the problem (6.1) is equivalent to solving the hypersingular integral equation

[TABLE]

see e.g. [40, 42]. Here, the right hand side $f$ of (6.4) is given by

[TABLE]

for $k=1,2$, and the operator $D^{*}$ is defined by

[TABLE]

see [34, page 122]. The exact solution of the exterior Neumann problem (6.1) is

[TABLE]

and the exact solution to the hypersingular integral equation (6.4) is given by

[TABLE]

We solve (6.4) by using the Galerkin method with $S(\Delta)$, the space of continuous piecewise linear spherical splines. Here, the spherical triangulations $\Delta$ are obtained in three different ways: uniform, residual and hierarchical adaptive mesh refinements. For experimental purposes, we start with an initial triangulation of eight equal spherical triangles with six nodes (two at the poles and four on the equator). For the uniform meshes, every further refinement consists of partitioning every spherical triangle into four smaller spherical triangles by joining the midpoints of the edges, see Figure 2. This guarantees that all triangles in the spherical triangulations obtained after refinements are of a finite number of similarly distinct triangles. For the residual and hierarchical adaptive meshes, we apply the strategies in Section 5 to refine the meshes after estimating the element errors, $\eta_{\Delta,1/2}(\tau)$ and $\theta_{\Delta}(\tau)$, see (3.5) and (4.36), respectively.

Suppose that $V_{\Delta}=\left\{\boldsymbol{v}_{1},\ldots,\boldsymbol{v}_{M}\right\}$ is the set of all vertices in the spherical triangulation $\Delta$. We choose a basis for $S(\Delta)$ to be the set

[TABLE]

where $B_{\boldsymbol{v}_{i}}$ is the basis function associated with the vertex $\boldsymbol{v}_{i}$, see (4.1). We denote by $u_{\Delta}\in S(\Delta)$ the Galerkin solution to (6.4). Then $u_{\Delta}=\sum_{i=1}^{M}\nu_{i}B_{\boldsymbol{v}_{i}}$, where $\nu_{i}\in{\mathbb{R}}$ for $i=1,\ldots,M$, satisfies

[TABLE]

This results in the following matrix equation

[TABLE]

The entry $A_{ij}$, for $i,j=1,\ldots,M$, of the stiffness matrix $\boldsymbol{A}$ is computed by

[TABLE]

The first integral in (6.8) is computed by

[TABLE]

see [34, Theorem 3.3.2]. Here, $\overrightarrow{\rm{curl}}_{\mathbb{S}}v$ is the vectorial surface rotation defined by

[TABLE]

where $\overrightarrow{e_{\varphi}}$, $\overrightarrow{e_{\theta}}$ are the two unit vectors corresponding to the Euler angles.

Computation of the double integrals in (6.9) requires evaluation of integrals of the type

[TABLE]

where $\tau^{(1)}$ and $\tau^{(2)}$ are spherical triangles in $\Delta$ and the functions $f_{1}$ and $f_{2}$ are analytic for all $\boldsymbol{x}\in\tau^{(1)}$ and $\boldsymbol{y}\in\tau^{(2)}$. For more details about the above evaluation, please refer to [36, 38].

The right hand side $\boldsymbol{F}$ of the linear system (6.7) has entries given by

[TABLE]

for all $i=1,\ldots,M$. Once solving the matrix equation (6.7), we obtain the coefficient vector $\boldsymbol{\nu}=(\nu_{1},\ldots,\nu_{M})$ and thus the approximate solution $u_{\Delta}=\sum_{i=1}^{M}\nu_{i}B_{\boldsymbol{v}_{i}}$. The error $\left\|{u-u_{\Delta}}\right\|_{H^{1/2}(\mathbb{S})}$ is then computed by

[TABLE]

noting (2.36)–(2.39).

We solve (6.4) by using uniform, residual and hierarchical adaptive refinements for the right hand sides $f_{1}$ and $f_{2}$ being defined by (6.5). For both examples, we find approximate solutions, compute the errors, degrees of freedom and accumulating computation time, see Tables 1–4. We note here that the convergence rates of the uniform refinement method for both $f_{1}$ and $f_{2}$ are slightly smaller than theoretical results. The errors behave roughly $\mathcal{O}(M^{-1.24/2})$ instead of $\mathcal{O}(M^{-1.5/2})$ as suggested by (1.3). This may be due to the small number of uniform meshes that have have been used and the low number of elements in these meshes.

The numerical results suggest significant advantages of the two adaptive refinement approaches in terms of required degrees of freedom and accumulating computation time, see also Figures 3–6. For example, to obtain an accuracy of around 3.5% when solving (6.4) for $f_{1}$, while the uniform refinement approach requires $4098$ degrees of freedom (see Figure 7) and the corresponding computation time is almost $10.7$ hours, our residual and hierarchical adaptive refinement counterparts need only $211$ and $170$ vertices and it takes only more than $10$ minutes to complete the calculation, see Tables 1–2 and Figures 3–4. Similar advantages of the adaptive refinement approaches are also observed when solving (6.4) for $f_{2}$ given by (6.5) and (6.3), see Tables 3–4 and Figures 5–6. For example, to obtain an accuracy of $3.6\%$, uniform refinement method has to use the uniform mesh of $4098$ vertices and the calculation takes nearly $10$ hours to complete. Meanwhile, the residual adaptive method requires a mesh of $448$ nodes and the (accumulating) computation time is about $17.5$ minutes. The numbers for the hierarchical adaptive counterpart are $302$ nodes and $32.8$ minutes, respectively.

Figure 8 shows adaptive meshes obtained when we solve the equation (6.4) with the right hand side $f_{1}$ by using the residual and hierarchical refinement approaches. Denser areas of nodes surrounding the north pole are observed. The spherical triangulations shown in Figures 9 and 10 are the $448$-node and $302$-node meshes obtained when we solve (6.4) with the right hand side $f_{2}$ by using the two adaptive methods. In these two figures, we witness denser areas surrounding the north and south poles. These denser areas are due to the fact that their contributions to the total errors are higher than other regions on the unit sphere, and thus must be accordingly refined as discussed in Section 5.

Acknowledgement

This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 101.99–2016.13.