# A stable algorithm for divergence-free radial basis functions in the flat
limit

Kathryn P. Drake<EMAIL_ADDRESS>Grady B. Wright
<EMAIL_ADDRESS>Boise State University, Department of Mathematics,
1910 University Drive, Boise, Idaho 83725-1555, USA

###### Abstract

The direct method used for calculating smooth radial basis function (RBF)
interpolants in the flat limit becomes numerically unstable. The RBF-QR
algorithm bypasses this ill-conditioning using a clever change of basis
technique. We extend this method for computing interpolants involving matrix-
valued kernels, specifically surface divergence-free RBFs on the sphere, in
the flat limit. Results illustrating the effectiveness of this algorithm are
presented for a divergence-free vector field on the sphere from samples at
scattered points.

###### keywords:

RBF, shape parameter, ill-conditioning, RBF-QR, matrix-valued, divergence-
free, sphere

##  1 Introduction

Tangential vector fields to a sphere are important in many areas of
geophysical sciences, from the surface of the ocean to the ionosphere [1].
Often the values of these vector fields may only be known at “scattered”
locations, e.g., from measurement taken from rawinsondes, airplanes, buoys,
remote sensing devices, or from output from certain numerical models (cf. [2,
§4]), and values of the field must be approximated at other locations, e.g.,
on a grid or mesh. Additionally, these tangential fields may satisfy certain
physical constraints, such as being surface divergence-free (div-free) or
curl-free, that must be preserved in the approximation. A radial basis
function (RBF) technique was developed exactly for these applications in the
papers [3, 4]. The idea is to construct a positive definite kernel from a
radial basis function in such a way that shifts of the kernel can be linearly
combined to yield a div-free or curl-free interpolant of the underlying field.
This technique has the added benefits that it gives a well-posed interpolant
for scattered data, is devoid of any coordinate singularities, and naturally
allows a scalar potential for the field to be extracted [4].

When using a smooth RBF that features a shape parameter, $\varepsilon$, to
construct these div-free or curl-free kernels, one often finds that the best
accuracy of the interpolated field is achieved when $\varepsilon$ is small,
corresponding to a flat kernel, but that the direct way of computing the
interpolant (often called RBF-Direct) is prohibitively ill-conditioned. This
is exactly analogous to the standard RBF interpolation problem for scalar
functions. In the scalar setting, three distinctly different numerical
algorithms have been presented thus far in the literature to bypass this ill-
conditioning and open up the complete range of $\varepsilon$ that can be
considered. These are the RBF-RA method [5, 6], the RBF-QR method [7, 8, 9],
and the RBF-GA method [10]. The present short note focuses on the RBF-QR
method from [8], which is specific to RBF interpolation on the sphere. This
algorithm exploits the Mercer expansion of the scalar RBFs in terms of
spherical harmonics to change the interpolation basis so that the the ill-
conditioning associated with small $\varepsilon$ is analytically removed. It
also shows that scalar RBF interpolants converge to spherical harmonic
interpolants in the flat limit (i.e. $\varepsilon\rightarrow 0$). We show how
this algorithm, which we call the vector RBF-QR, can be extended to bypass the
ill-conditioning associated with surface div-free RBF interpolation for small
$\varepsilon$ and demonstrate that these interpolants converge to div-free
vector spherical harmonic interpolants as $\varepsilon\rightarrow 0$. The
algorithm also extends naturally to the surface curl-free case, but the sake
of brevity we leave out these details and refer the reader to [11].

The rest of the paper is organized as follows: In Section 2, we introduce some
notation and background on div-free RBF interpolation and scalar/vector
spherical harmonics. Section 3 contains the description of the new vector RBF-
QR algorithm and Section 4 contains a numerical example illustrating the
stability of the algorithm in the flat limit. We conclude with some brief
comments in Section 5.

##  2 Notation and preliminaries

###  2.1 Surface div-free vector fields

Any $C^{1}$ tangential velocity field $\mathbf{u}$ on $\mathbb{S}^{2}$ that is
surface div-free can be written as

$\displaystyle\mathbf{u}=\nabla^{\ast}\times(\psi\widehat{\mathbf{r}})=-\underbrace{\hat{\mathbf{r}}\times\nabla^{\ast}}_{\mathbf{L}^{\ast}}\psi,$

where $\nabla^{\ast}$ denotes the gradient in spherical coordinates,
$\hat{\mathbf{r}}$ is the unit radial vector in spherical coordinate basis,
and $\psi$ is some $C^{2}$ scalar-valued function on $\mathbb{S}^{2}$. The
function $\psi$ is called the stream function and is unique up to a constant
for a given $\mathbf{u}$ [12, Proposition 2.1].

The operator $\mathbf{L}^{\ast}$, which is sometimes called the surface-curl
operator or the _rot_ operator, can be written entirely in extrinsic
(Cartesian) coordinates as $\mathbf{L}=\hat{\mathbf{n}}\times\nabla$, where
$\nabla$ is the standard $\mathbb{R}^{3}$ gradient in the Cartesian basis, and
$\hat{\mathbf{n}}$ is the unit normal vector to $\mathbb{S}^{2}$ in the
Cartesian basis [3]. Here we have dropped the $\ast$ from $\mathbf{L}$ to
indicate it is defined in extrinsic coordinates. For scalar-valued functions
$\psi$ that can be extended smoothly from $\mathbb{S}^{2}$ to
$\mathbb{R}^{3}$, we can generate a tangential velocity field that is surface
div-free using $\mathbf{u}=\mathbf{L}\psi$. The operator $\mathbf{L}$ can be
simplified further by noting that $\hat{\mathbf{n}}$ at a point
$\mathbf{x}=(x,y,z)$ on $\mathbb{S}^{2}$ is just $\mathbf{x}$. In what
proceeds, we use the following notation for the operator $\mathbf{L}$:

$\displaystyle\mathbf{L}_{\mathbf{x}}=\underbrace{\begin{pmatrix}0&-z&y\\\
z&0&-x\\\ -y&x&0\end{pmatrix}}_{\displaystyle
Q(\mathbf{x})}\nabla_{\mathbf{x}},$ (1)

where $Q(\mathbf{x})$ applied to a vector in $\mathbb{R}^{3}$ gives the cross
product of $\mathbf{x}$ with that vector and the subscripts are used to
indicate what variable each operator is applied to. We note that one of the
main benefits of working in extrinsic coordinates is that artificial
coordinate singularities (e.g. pole singularities) can be avoided.

###  2.2 Surface div-free RBF interpolation

Surface div-free RBF interpolation is similar to scalar RBF interpolation in
the sense that one constructs interpolants from linear combinations of shifts
of a kernel at each of the given data sites; for a review of scalar RBF
interpolation, we refer the reader to [13]. The difference between the two
approaches is that surface div-free RBF interpolants use a matrix-valued
kernel $\Phi_{\rm div}$ whose columns are surface div-free. A detailed
discussion of the construction of $\Phi_{\rm div}$ is given in [3]. For the
sake of brevity, we do not repeat this derivation, but only present the final
result. Let $\phi$ be a scalar-valued, radial kernel on $\mathbb{R}^{3}$ with
at least two continuous derivatives and let
$\mathbf{x},\mathbf{y}\in\mathbb{S}^{2}$, then $\Phi_{\rm div}$ is constructed
using the operator $\mathbf{L}$ in (1) as

$\Phi_{\rm
div}(\mathbf{x},\mathbf{y})=\mathbf{L}_{\mathbf{x}}^{\phantom{T}}\mathbf{L}_{\mathbf{y}}^{T}\phi\left(\|\mathbf{x}-\mathbf{y}\|\right)=-Q(\mathbf{x})\left(\nabla_{\mathbf{x}}^{\phantom{T}}\nabla_{\mathbf{y}}^{T}\phi\left(\|\mathbf{x}-\mathbf{y}\|\right)\right)Q(\mathbf{y})=Q(\mathbf{x})\left(\nabla_{\mathbf{x}}^{\phantom{T}}\nabla_{\mathbf{x}}^{T}\phi\left(\|\mathbf{x}-\mathbf{y}\|\right)\right)Q(\mathbf{y}),$
(2)

where $\|\cdot\|$ denotes the vector two-norm and we have used the fact that
the $Q$ matrix in (1) is anti-symmetric and
$\nabla_{\mathbf{y}}^{T}\phi\left(\|\mathbf{x}-\mathbf{y}\|\right)=-\nabla_{\mathbf{x}}^{T}\phi\left(\|\mathbf{x}-\mathbf{y}\|\right)$.
It is straightforward to show that for any $\mathbf{c}\in\mathbb{R}^{3}$ that
is tangent to $\mathbb{S}^{2}$ at $\mathbf{y}$, the vector $\Phi_{\rm
div}(\mathbf{x},\mathbf{y})\mathbf{c}$ is surface div-free in $\mathbf{x}$ and
centered at the point $\mathbf{y}$; see Figure 1 for an illustration.

|
---|---
(a) | (b)
Figure 1: Illustration of the matrix-valued surface div-free kernel (2) using
the Hammer-Aitoff projection of the sphere [14]: (a) $\Phi_{\rm div}$ applied
to the zonal unit vector at $\mathbf{y}=(1,0,0)$, (a) $\Phi_{\rm div}$ applied
to the meridional unit vector at $\mathbf{y}=(1,0,0)$.

Let $\\{\mathbf{y}_{j}\\}_{j=1}^{n}$ be a distinct set of nodes on
$\mathbb{S}^{2}$ and $\\{\mathbf{u}_{j}\\}_{j=1}^{n}$ be samples of a surface
div-free tangent vector field sampled on these nodes. The surface div-free RBF
interpolant to this data takes the form

$\mathbf{s}(\mathbf{x})=\sum_{j=1}^{n}\Phi_{\rm
div}(\mathbf{x},\mathbf{y}_{j})\mathbf{c}_{j},$ (3)

where the coefficients $\mathbf{c}_{j}\in\mathbb{R}^{3}$ are tangent to
$\mathbb{S}^{2}$ at $\mathbf{y}_{j}$, which is necessary to make the
interpolation problem well-posed. A naïve approach to solving for the
$\\{\mathbf{c}_{j}\\}_{j=1}^{n}$ by imposing
$\mathbf{s}(\mathbf{y}_{j})=\mathbf{u}_{j}$, $j=1,\ldots,n$ will lead to a
singular system of equations since each $\mathbf{c}_{j}$ (and correspondingly
$\mathbf{u}_{j}$) can be expressed using only two degrees of freedom rather
than three. To see this, let
$\\{\mathbf{a}_{j},\mathbf{b}_{j},\mathbf{n}_{j}\\}$ be orthonormal vectors at
the node $\mathbf{y}_{j}$, where $\mathbf{n}_{j}$ is the outward normal to
$\mathbb{S}^{2}$, $\mathbf{b}_{j}$ is a unit tangent vector, and
$\mathbf{a}_{j}=\mathbf{n}_{j}\times\mathbf{b}_{j}$. Then we can write
$\mathbf{c}_{j}$ in this basis as
$\mathbf{c}_{j}=\alpha_{j}\mathbf{a}_{j}+\beta_{j}\mathbf{b}_{j}$, where
$\alpha_{j}=\mathbf{a}_{j}^{T}\mathbf{c}_{j}$ and
$\beta_{j}=\mathbf{b}_{j}^{T}\mathbf{c}_{j}$. Using this result, we can
express (3) as

$\displaystyle\mathbf{s}(\mathbf{x})=\sum_{j=1}^{n}\Phi_{\rm
div}(\mathbf{x},\mathbf{y}_{j})\underbrace{\left[\alpha_{j}\mathbf{a}_{j}+\beta_{j}\mathbf{b}_{j}\right]}_{\displaystyle\mathbf{c}_{j}},$
(4)

and write the interpolation conditions as
$\mathbf{a}_{i}^{T}\mathbf{s}(\mathbf{y}_{i})=\mathbf{a}_{i}^{T}\mathbf{u}_{i}=:\gamma_{i}$
and
$\mathbf{b}_{i}^{T}\mathbf{s}(\mathbf{y}_{i})=\mathbf{b}_{i}^{T}\mathbf{u}_{i}=:\delta_{i}$,
which leads to the $2n$-by-$2n$ system of equations

$\displaystyle\sum_{j=1}^{n}\underbrace{\left(\begin{bmatrix}\mathbf{a}_{i}^{T}\\\
\mathbf{b}_{i}^{T}\end{bmatrix}\Phi_{\rm
div}(\mathbf{y}_{i},\mathbf{y}_{j})\begin{bmatrix}\mathbf{a}_{j}&\mathbf{b}_{j}\end{bmatrix}\right)}_{\displaystyle
A_{i,j}}\begin{bmatrix}\alpha_{j}\\\
\beta_{j}\end{bmatrix}=\begin{bmatrix}\gamma_{i}\\\
\delta_{i}\end{bmatrix},\quad 1\leq i\leq n.$ (5)

The interpolation matrix that arises from this system (with entries given by
$A_{i,j}$) is positive definite if $\Phi_{\rm div}$ is constructed from an
appropriately chosen scalar-valued $\phi$ [3], such as any of those listed in
Table13. Simplifications of the entries of this matrix in terms of derivatives
of $\phi$ are given in [4].

In this work, we choose $\mathbf{a}_{j}$ and $\mathbf{b}_{j}$, for the point
$\mathbf{y}_{j}=(x_{j},y_{j},z_{j})$, to be the standard meridional and zonal
vectors, which can be expressed in Cartesian coordinates as

$\mathbf{a}_{j}=\frac{1}{\sqrt{1-z_{j}^{2}}}\begin{bmatrix}-z_{j}x_{j}\\\
-z_{j}y_{j}\\\
1-z_{j}^{2}\end{bmatrix},\quad\mathbf{b}_{j}=\frac{1}{\sqrt{1-z_{j}^{2}}}\begin{bmatrix}-y_{j}\\\
x_{j}\\\ 0\end{bmatrix}.$ (6)

If $\mathbf{y}_{j}=[0,0,\pm 1]$, then we can pick any orthogonal vectors in
the $xy$-plane.

We conclude by noting that once the coefficients $\mathbf{c}_{j}$ are
determined in (4), a stream-function $\psi(\mathbf{x})$ can obtained for the
interpolated field using (2)

$\displaystyle\mathbf{s}(\mathbf{x})=\underbrace{Q(\mathbf{x})\nabla_{\mathbf{x}}}_{\displaystyle\mathbf{L}_{\mathbf{x}}}\biggl{(}\underbrace{\sum_{j=1}^{n}\nabla_{\mathbf{x}}^{T}\phi\left(\|\mathbf{x}-\mathbf{y}_{j}\|\right)Q(\mathbf{y}_{j})\mathbf{c}_{j}}_{\displaystyle\psi(\mathbf{x})}\biggr{)}.$
(7)

Approximation results for $\psi$ in reconstructing the underlying stream
function for the field are given in [4].

Similar to scalar RBF interpolation, we refer to the method of computing the
div-free RBF interpolant by solving the system (5) as RBF-Direct. For
$\Phi_{\rm div}$ that are built from scalar kernels that depend on a shape
parameter $\varepsilon$, like those in Table 1, the RBF-Direct approach
becomes extremely ill-conditioned as $\varepsilon\rightarrow 0$. However, as
in the scalar RBF case, the interpolant remains well-conditioned and the RBF-
QR algorithm allows us to stably compute it.

###  2.3 Spherical harmonics

The RBF-QR algorithm presented in Section 3 relies on spherical harmonic
expansions, both scalar and vector ones, so we briefly review these expansions
here. We refer the reader to more detailed discussions in [15] and [16], for
the scalar and vector case, respectively.

Let $Y_{\mu}^{\nu}$ denote the scalar spherical harmonic of integer degree
$\mu\geq 0$ and integer order $-\mu\leq\nu\leq\mu$ on $\mathbb{S}^{2}$. These
functions are the eigenfunctions of the Laplace-Beltrami operator
$\Delta_{\mathbb{S}^{2}}$, i.e.,
$\Delta_{\mathbb{S}^{2}}Y_{\mu}^{\nu}=-\mu(\mu+1)Y_{\mu}^{\nu}$. We use the
real-form of the spherical harmonics functions in Cartesian coordinates, which
for $\mathbf{x}=(x,y,z)\in\mathbb{S}^{2}$ are given as

$\displaystyle
Y_{\mu}^{\nu}(\mathbf{x})=\begin{cases}\sqrt{\frac{2\mu+1}{4\pi}}\sqrt{\frac{(\mu-\nu)!}{(\mu+\nu)!}}P_{\mu}^{\nu}(z)\cos\left(\nu\tan^{-1}\left(\frac{y}{x}\right)\right),&\text{
}\nu=0,1,\dots,\mu,\\\
\sqrt{\frac{2\mu+1}{4\pi}}\sqrt{\frac{(\mu-\nu)!}{(\mu+\nu)!}}P_{\mu}^{\nu}(z)\sin\left(-\nu\tan^{-1}\left(\frac{y}{x}\right)\right),&\text{
}\nu=-\mu,\dots,-1\end{cases},$ (8)

where $P_{\mu}^{\nu}(z)$ are the associated Legendre functions of degree $\mu$
and order $\nu$. The spherical harmonics form a complete, orthonormal set of
basis functions for the space of square-integrable functions on
$\mathbb{S}^{2}$, which we denote by $L^{2}(\mathbb{S}^{2})$ [15]. Thus, any
function $f\in L^{2}(\mathbb{S}^{2})$ can be uniquely represented as

$\displaystyle
f(\mathbf{x})=\sum_{\mu=0}^{\infty}\sum_{\nu=-\mu}^{\mu}\hat{f}_{\mu,\nu}Y_{\mu}^{\nu}(\mathbf{x}),\;\text{where}\;\;\hat{f}_{\mu,\nu}=\int_{\mathbb{S}^{2}}f(\mathbf{x})Y_{\mu}^{\nu}(\mathbf{x})dS.$
(9)

Vector spherical harmonics are the vectorial analogue of scalar spherical
harmonics and can be used for representing vector-valued functions on
$\mathbb{S}^{2}$. There are three types of these vector harmonic functions:
one that is normal to $\mathbb{S}^{2}$, one that is tangent to
$\mathbb{S}^{2}$ and surface curl-free, and one that is tangent to
$\mathbb{S}^{2}$ and surface div-free [16]. We focus on the latter as they are
the ones used in this paper. Both tangent vector spherical harmonics are the
eigenfunctions of the vector Laplace-Beltrami operator.

The surface div-free vector spherical harmonics can be constructed in
Cartesian coordinates by applying operator $\mathbf{L}$ in (1) to the scalar
spherical harmonic functions (8) as follows [16]:

$\mathbf{w}_{\mu}^{\nu}(\mathbf{x})=\mathbf{L}_{\mathbf{x}}Y_{\mu}^{\nu}(\mathbf{x}),\;\mu>0,\;-\mu\leq\nu\leq\mu.$
(10)

These can be normalized as
$\overline{\mathbf{w}}_{\mu}^{\nu}=\frac{1}{\mu(\mu+1)}\mathbf{w}_{\mu}^{\nu}(vx)$
so they are orthonormal with respect to the $L^{2}(\mathbb{S}^{2})$-vector
inner product
$\langle\mathbf{f},\mathbf{g}\rangle=\int_{\mathbb{S}^{2}}\mathbf{f}^{T}\mathbf{g}dS$,
where $\mathbf{f}$ and $\mathbf{g}$ are tangent vector fields on
$\mathbb{S}^{2}$. The surface div-free vector spherical harmonics form a
complete orthonormal set of vector basis functions for the space of surface
div-free vector fields that are tangent to $\mathbb{S}^{2}$ and square
integrable. We denote this space as $L_{\rm div}(\mathbb{S}^{2})$. Any
function $\mathbf{u}\in L_{\rm div}^{2}(\mathbb{S}^{2})$ can be expanded in
these harmonics as

$\mathbf{u}(\mathbf{x})=\sum_{\mu=1}^{\infty}\sum_{\nu=-\mu}^{\mu}\hat{u}_{\mu,\nu}\overline{\mathbf{w}}_{\mu}^{\nu}(\mathbf{x}),\;\text{where}\;\;\hat{u}_{\mu,\nu}=\int_{\mathbb{S}^{2}}\mathbf{f(\mathbf{x})}^{T}\overline{\mathbf{w}}_{\mu}^{\nu}(\mathbf{x})dS.$
(11)

Note here that the sum in this expansion excludes $\mu=0$ since $Y_{0}^{0}$ is
annihilated by $\mathbf{L}$.

As discussed in Section 2.2, the interpolation problem on the sphere requires
that we utilize the tangent basis vectors at $\mathbf{x}\in\mathbb{S}^{2}$,
i.e., $\mathbf{a}_{\mathbf{x}}$ and $\mathbf{b}_{\mathbf{x}}$. Therefore, it
is relevant to introduce the following notation for the non-normalized surface
div-free and curl-free vector spherical harmonics in terms of these basis
vectors:

$\displaystyle\textnormal{meridional, div-free
:}\;G_{\mu}^{\nu}(\mathbf{x})=\mathbf{a}_{\mathbf{x}}^{T}\mathbf{L}_{\mathbf{x}}Y_{\mu}^{\nu}(\mathbf{x}),\qquad\textnormal{zonal,
div-free
:}\;H_{\mu}^{\nu}(\mathbf{x})=\mathbf{b}_{\mathbf{x}}^{T}\mathbf{L}_{\mathbf{x}}Y_{\mu}^{\nu}(\mathbf{x}).$
(12)

###  2.4 Spherical harmonic expansions of smooth radial kernels

The central idea of the scalar RBF-QR method is to replace the standard basis
consisting of shifts of radial kernels by an equivalent, but much better
conditioned basis that spans the same space in the case of small
$\varepsilon$. For the sphere, this is done by exploiting properties of the
Mercer expansion of smooth radial kernels on the sphere in terms of spherical
harmonics. As discussed in [8], these expansions are given by

$\phi(\|\mathbf{x}-\mathbf{y}_{j}\|)=\sum_{\mu=0}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{\nu=-\mu}^{\mu}\\{c_{\mu,\varepsilon}\varepsilon^{2\mu}Y_{\mu}^{\nu}(\mathbf{y}_{j})\\}Y_{\mu}^{\nu}(\mathbf{x}),$
(13)

where the symbol $\sum^{\prime}$ denotes that the $\nu=0$ term is halved. Note
that the spherical harmonic coefficients are independent of $\nu$, which
follows from the Funk-Hecke formula [17, 15] for zonal functions. Table 1
lists the coefficients $c_{\mu,\varepsilon}$ for many common radial kernels.
These were first computed by Hubbert and Baxter [18] for all the radial
kernels listed in this table except for the IQ, which was given in [8].

Radial kernel | Expression | Expansion coefficient, $c_{\mu,\varepsilon}$
---|---|---
Multiquadric (MQ) | $\phi(r)=(1+(\varepsilon r)^{2})^{\frac{1}{2}}$ | $\frac{-2\pi(2\varepsilon^{2}+1+(\mu+1/2)\sqrt{1+4\varepsilon^{2}})}{(\mu+3/2)(\mu+1/2)(\mu-1/2)}\left(\frac{2}{1+\sqrt{4\varepsilon^{2}+1}}\right)^{2\mu+1}$
Inverse multiquadric (IMQ) | $\phi(r)=(1+(\varepsilon r)^{2})^{-\frac{1}{2}}$ | $\frac{4\pi}{(\mu+1/2)}\left(\frac{2}{1+\sqrt{4\varepsilon^{2}+1}}\right)^{2\mu+1}$
Inverse quadratic (IQ) | $\phi(r)=(1+(\varepsilon r)^{2})^{-1}$ | $\frac{4\pi^{3/2}\mu!}{\Gamma(\mu+\frac{3}{2})(1+4\varepsilon^{2})^{\mu+1}}\prescript{}{2}{F}_{1}(\mu+1,\mu+1;2\mu+2;\frac{4\varepsilon^{2}}{1+4\varepsilon^{2}})$
Gaussian (GA) | $\phi(r)=\exp(-(\varepsilon r)^{2})$ | $\frac{4\pi^{3/2}}{\varepsilon^{2\mu+1}}e^{-2\varepsilon^{2}}I_{\mu+1/2}(2\varepsilon^{2})$
Table 1:  Coefficients in the spherical harmonic expansion (13) for various
smooth radial kernels on the sphere. For the IQ kernel,
$\prescript{}{2}{F}_{1}(\dots)$ denotes the hypergeometric function, and for
the GA kernel, $I_{\mu+1/2}$ denotes the Bessel function of the second kind.

It is important to note that the coefficients listed in Table 1 can be
calculated without the loss of any significant digits caused by numerical
cancellations, even as $\varepsilon\rightarrow 0$ (see [8] for a discussion).

The central idea behind the vector RBF-QR algorithm will also be to replace
the matrix-valued basis with a better basis built from vector spherical
harmonic expansions. We can use (2) and (13) to expand the matrix-valued
kernel $\Phi_{\rm div}$ in terms of vector spherical harmonics as follows:

$\displaystyle\Phi_{\rm
div}(\mathbf{x},\mathbf{y}_{j})=\sum_{\mu=1}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{\nu=-\mu}^{\mu}\\{\varepsilon^{2\mu}c_{\mu,\varepsilon}\mathbf{L}_{\mathbf{x}}Y_{\mu}^{\nu}(\mathbf{x})\\}\left(\left.\mathbf{L}_{\mathbf{y}}Y_{\mu}^{\nu}(\mathbf{y})\right|_{\mathbf{y}=\mathbf{y}_{j}}\right)^{T}=\sum_{\mu=1}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{\nu=-\mu}^{\mu}\\{\varepsilon^{2\mu}c_{\mu,\varepsilon}\mathbf{w}_{\mu}^{\nu}(\mathbf{x})\\}\left(\mathbf{w}_{\mu}^{\nu}(\mathbf{y}_{j})\right)^{T}.$
(14)

Here we have used the non-normalized div-free vector spherical harmonics
defined in (10) for simplicity.

The separation of the expansion (14) of a smooth matrix-valued kernel
$\Phi_{\rm div}$ in terms of increasing powers of $\varepsilon$ is essential
to the RBF-QR algorithm. It is these powers of $\varepsilon$ and not the
coefficients $c_{\mu,\varepsilon}$ that lead ill-conditioning in the matrix-
valued basis for small $\varepsilon$ in the RBF-Direct method. The RBF-QR
algorithm analytically factors out the effects of these powers of
$\varepsilon$ from the basis.

##  3 Vector RBF-QR algorithm

Recall from Section 2.2 that in order to interpolate div-free vector fields
tangent to the sphere with the matrix-valued div-free interpolant, we must
represent the coefficient vectors and target field samples in terms of the
orthonormal tangent basis vectors (e.g., using (6)). In (5) we saw that this
is equivalent to representing the kernel $\Phi_{\rm div}$ in terms of these
basis vectors. We will denote this kernel as $\widetilde{\Phi}_{\rm div}$:

$\widetilde{\Phi}_{\rm
div}(\mathbf{x},\mathbf{y}_{j})=\begin{bmatrix}\mathbf{a}_{\mathbf{x}}^{T}\\\
\mathbf{b}_{\mathbf{x}}^{T}\end{bmatrix}\Phi_{\rm
div}(\mathbf{x},\mathbf{y}_{j})\begin{bmatrix}\mathbf{a}_{j}&\mathbf{b}_{j}\end{bmatrix}.$
(15)

Using (14) on the right-hand side of (15) gives the expansion

$\widetilde{\Phi}_{\rm
div}(\mathbf{x},\mathbf{y}_{j})=\sum_{\mu=1}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{\nu=-\mu}^{\mu}\left(\varepsilon^{2\mu}c_{\mu,\varepsilon}\begin{bmatrix}\mathbf{a}_{\mathbf{x}}^{T}\\\
\mathbf{b}_{\mathbf{x}}^{T}\end{bmatrix}\mathbf{w}_{\mu}^{\nu}(\mathbf{x})\right)\left(\mathbf{w}_{\mu}^{\nu}(\mathbf{y}_{j})\right)^{T}\begin{bmatrix}\mathbf{a}_{j}&\mathbf{b}_{j}\end{bmatrix}.$
(16)

This is a $2$-by-$2$ matrix whose entries are in terms of the meridional and
zonal div-free vector spherical harmonics (12):

$\widetilde{\Phi}_{\rm
div}(\mathbf{x},\mathbf{y}_{j})=\begin{bmatrix}\begin{array}[]{c|c}(a)&(b)\\\
\hline\cr(c)&(d)\\\ \end{array}\end{bmatrix},\;\textnormal{where}$ (17)
$\displaystyle\text{(a)}\;\sum_{\mu=1}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{\nu=-\mu}^{\mu}\\{\varepsilon^{2\mu}c_{\mu,\varepsilon}G_{\mu}^{\nu}(\mathbf{x})\\}G_{\mu}^{\nu}(\mathbf{y}_{j})$
$\displaystyle\qquad\text{(b)}\;\sum_{\mu=1}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{\nu=-\mu}^{\mu}\\{\varepsilon^{2\mu}c_{\mu,\varepsilon}G_{\mu}^{\nu}(\mathbf{x})\\}H_{\mu}^{\nu}(\mathbf{y}_{j})$
$\displaystyle\text{(c)}\;\sum_{\mu=1}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{\nu=-\mu}^{\mu}\\{\varepsilon^{2\mu}c_{\mu,\varepsilon}H_{\mu}^{\nu}(\mathbf{x})\\}G_{\mu}^{\nu}(\mathbf{y}_{j})$
$\displaystyle\qquad\text{(d)}\;\sum_{\mu=1}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{\nu=-\mu}^{\mu}\\{\varepsilon^{2\mu}c_{\mu,\varepsilon}H_{\mu}^{\nu}(\mathbf{x})\\}H_{\mu}^{\nu}(\mathbf{y}_{j}).$

The goal of the vector RBF-QR algorithm is to express the space spanned by the
columns of the $2$-by-$2n$ array containing the shifts of the div-free matrix
valued kernel, $\begin{bmatrix}\widetilde{\Phi}_{\rm
div}(\mathbf{x},\mathbf{y}_{1})&\cdots&\widetilde{\Phi}_{\rm
div}(\mathbf{x},\mathbf{y}_{n})\end{bmatrix},$ using a basis that has the ill-
conditioning associated with small $\varepsilon$ removed. To this end, we
first use (17) to write this array (now in transposed form) as

$\begin{bmatrix}\widetilde{\Phi}_{\rm div}(\mathbf{x},\mathbf{y}_{1})^{T}\\\
\vdots\\\ \widetilde{\Phi}_{\rm
div}(\mathbf{x},\mathbf{y}_{n})^{T}\end{bmatrix}=\begin{bmatrix}\begin{bmatrix}\begin{array}[]{c|c}\displaystyle\sum_{\mu=1}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{\nu=-\mu}^{\mu}\\{\varepsilon^{2\mu}c_{\mu,\varepsilon}G_{\mu}^{\nu}(\mathbf{x})\\}G_{\mu}^{\nu}(\mathbf{y}_{1})&\displaystyle\sum_{\mu=1}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{\nu=-\mu}^{\mu}\\{\varepsilon^{2\mu}c_{\mu,\varepsilon}H_{\mu}^{\nu}(\mathbf{x})\\}G_{\mu}^{\nu}(\mathbf{y}_{1})\\\
\hline\cr\displaystyle\sum_{\mu=1}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{\nu=-\mu}^{\mu}\\{\varepsilon^{2\mu}c_{\mu,\varepsilon}G_{\mu}^{\nu}(\mathbf{x})\\}H_{\mu}^{\nu}(\mathbf{y}_{1})&\displaystyle\sum_{\mu=1}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{\nu=-\mu}^{\mu}\\{\varepsilon^{2\mu}c_{\mu,\varepsilon}H_{\mu}^{\nu}(\mathbf{x})\\}H_{\mu}^{\nu}(\mathbf{y}_{1})\\\
\end{array}\end{bmatrix}\\\ \vdots\\\
\begin{bmatrix}\begin{array}[]{c|c}\displaystyle\sum_{\mu=1}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{\nu=-\mu}^{\mu}\\{\varepsilon^{2\mu}c_{\mu,\varepsilon}G_{\mu}^{\nu}(\mathbf{x})\\}G_{\mu}^{\nu}(\mathbf{y}_{n})&\displaystyle\sum_{\mu=1}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{\nu=-\mu}^{\mu}\\{\varepsilon^{2\mu}c_{\mu,\varepsilon}H_{\mu}^{\nu}(\mathbf{x})\\}G_{\mu}^{\nu}(\mathbf{y}_{n})\\\
\hline\cr\displaystyle\sum_{\mu=1}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{\nu=-\mu}^{\mu}\\{\varepsilon^{2\mu}c_{\mu,\varepsilon}G_{\mu}^{\nu}(\mathbf{x})\\}H_{\mu}^{\nu}(\mathbf{y}_{n})&\displaystyle\sum_{\mu=1}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{\nu=-\mu}^{\mu}\\{\varepsilon^{2\mu}c_{\mu,\varepsilon}H_{\mu}^{\nu}(\mathbf{x})\\}H_{\mu}^{\nu}(\mathbf{y}_{n})\\\
\end{array}\end{bmatrix}\end{bmatrix}.$ (18)

We then rewrite this as the following infinite block matrix-matrix product,

$\displaystyle\underbrace{\begin{bmatrix}\displaystyle\widetilde{\Phi}_{\rm
div}(\mathbf{x},\mathbf{y}_{1})^{T}\\\ \vdots\\\
\displaystyle\widetilde{\Phi}_{\rm
div}(\mathbf{x},\mathbf{y}_{n})^{T}\end{bmatrix}}_{\displaystyle\widetilde{P}_{\rm
div}}=\underbrace{\begin{bmatrix}c_{1,\varepsilon}G_{1}^{-1}(\mathbf{y}_{1})&\frac{c_{1,\varepsilon}}{2}G_{1}^{0}(\mathbf{y}_{1})&c_{1,\varepsilon}G_{1}^{1}(\mathbf{y}_{1})\cdots\\\
c_{1,\varepsilon}H_{1}^{-1}(\mathbf{y}_{1})&\frac{c_{1,\varepsilon}}{2}H_{1}^{0}(\mathbf{y}_{1})&c_{1,\varepsilon}H_{1}^{1}(\mathbf{y}_{1})\cdots\\\
\vdots&\vdots&\vdots\\\
c_{1,\varepsilon}G_{1}^{-1}(\mathbf{y}_{n})&\frac{c_{1,\varepsilon}}{2}G_{1}^{0}(\mathbf{y}_{n})&c_{1,\varepsilon}G_{1}^{1}(\mathbf{y}_{n})\cdots\\\
c_{1,\varepsilon}H_{1}^{-1}(\mathbf{y}_{n})&\frac{c_{1,\varepsilon}}{2}H_{1}^{0}(\mathbf{y}_{n})&c_{1,\varepsilon}H_{1}^{1}(\mathbf{y}_{n})\cdots\\\
\end{bmatrix}}_{\displaystyle
B^{\infty}}\underbrace{\begin{bmatrix}\varepsilon^{2}&&\\\
&\varepsilon^{2}&\\\ &&\varepsilon^{2}\\\ &&&\ddots\vspace{-0.1in}\\\
&&&&\ddots\end{bmatrix}}_{\displaystyle
E^{\infty}}\underbrace{\begin{bmatrix}G_{1}^{-1}(\mathbf{x})&H_{1}^{-1}(\mathbf{x})\vspace{1.3ex}\\\
G_{1}^{0}(\mathbf{x})&H_{1}^{0}(\mathbf{x})\vspace{1.3ex}\\\
G_{1}^{1}(\mathbf{x})&H_{1}^{1}(\mathbf{x})\\\ \vdots&\vdots\\\
\end{bmatrix}}_{\displaystyle Y^{\infty}}.$ (19)

The next step is to truncate these infinite matrices to a vector spherical
harmonic degree value $\mu=\mu_{\text{trunc}}$. There are two stipulations for
this truncation degree. First, $\mu_{\text{trunc}}$ must be large enough that
the entries in $\widetilde{P}_{\rm div}$ are approximated to machine
precision; we refer to this value as $\mu_{\text{eps}}$. The second is that
$\mu_{\text{trunc}}$ must be large enough that the number of columns in the
truncated $B^{\infty}$ matrix is greater than or equal to $2n$; we refer to
this value as $\mu_{0}$. This condition on $\mu_{0}$ is given explicitly as
$\mu_{0}=\left\lceil\sqrt{2n+1}-1\right\rceil$ since the number of vector
spherical harmonic terms in a truncated expansion of $\mu_{0}$ is
$\mu_{0}(\mu_{0}+2)$. Putting the truncation requirements together gives the
condition $\mu_{\text{trunc}}\geq\max(\mu_{0},\mu_{\text{eps}})$. We denote
the truncated matrix-matrix product (19) as

$\displaystyle\widetilde{P}_{\rm div}\approx\begin{bmatrix}\\\ &&&B&&&\\\ \\\
\end{bmatrix}\begin{bmatrix}\\\ &&E&&\\\ \\\ \end{bmatrix}\begin{bmatrix}\\\
Y\\\ \\\ \end{bmatrix}.$ (20)

The matrix $B$ has $m=\mu_{\text{trunc}}(\mu_{\text{trunc}}+2)$ columns and
has the block form

$\displaystyle
B=\left[\begin{array}[]{c|c|c|c|c|c|c}B_{1}&B_{2}&\cdots&B_{\mu_{0}}&B_{\mu_{0}+1}&\cdots&B_{\mu_{\text{trunc}}}\end{array}\right]$
(22)

where $B_{\mu}$, $1\leq\mu\leq\mu_{\text{trunc}}$, are the block matrices of
size $2n$-by-$(2\mu+1)$ with block entries

$\left(B_{\mu}\right)_{i,j}=\left\\{\begin{matrix}\begin{bmatrix}c_{\mu,\varepsilon}G_{\mu}^{j-(\mu+1)}(\mathbf{y}_{i})\;\\\
c_{\mu,\varepsilon}H_{\mu}^{j-(\mu+1)}(\mathbf{y}_{i})\end{bmatrix}j\neq\mu+1,\vspace{2ex}\\\
\begin{bmatrix}\frac{c_{\mu,\varepsilon}}{2}G_{\mu}^{0}(\mathbf{y}_{i})\;\\\
\frac{c_{\mu,\varepsilon}}{2}H_{\mu}^{0}(\mathbf{y}_{i})\end{bmatrix}\;\qquad
j=\mu+1,\end{matrix}\right.\;\;j=1,\ldots,2\mu+1,\;i=1,\ldots,n.$

The diagonal $E$ matrix in (20) can be written as two square, diagonal blocks,
$E=\text{diag}([E_{1};E_{2}])$, where

$\displaystyle E_{1}=\begin{bmatrix}&\varepsilon^{2}I_{3}\\\
&&\varepsilon^{4}I_{5}\\\ &&&\ddots\\\
&&&&\varepsilon^{2\mu_{0}}I_{2\mu_{0}+1}\end{bmatrix}\;\text{and}\;E_{2}=\begin{bmatrix}\varepsilon^{2\mu_{0}+2}I_{2\mu_{0}+3}\\\
&\varepsilon^{2\mu_{0}+4}I_{2\mu_{0}+5}\\\ &&\ddots\\\
&&&\varepsilon^{2\mu_{\text{trunc}}}I_{2\mu_{\text{trunc}}+1}\end{bmatrix},$
(23)

and $I_{\tau}$ is the identity matrix of size $\tau$-by-$\tau$. Finally, the
$Y$ matrix in (20) is given in block form by

$\displaystyle
Y^{T}=\left[\begin{array}[]{c|c|c|c|c|c|c}Y_{1}^{T}&Y_{2}^{T}&\cdots&Y_{\mu_{0}}^{T}&Y_{\mu_{0}+1}^{T}&\cdots&Y_{\mu_{\text{trunc}}}^{T}\end{array}\right],$

where block $Y_{\mu}$ is given by

$\displaystyle\left(Y_{\mu}\right)_{j,1}=G_{\mu}^{j-(\mu+1)}(\mathbf{x}),\quad\left(Y_{\mu}\right)_{j,2}=H_{\mu}^{j-(\mu+1)}(\mathbf{x}),\quad
j=1,\ldots,2\mu+1.$

In the flat limit, the truncated basis is still highly ill-conditioned because
of the powers of $\varepsilon$ in (19) (recall that the expansion coefficients
$c_{\mu,\varepsilon}$ do not affect the conditioning). However, all of these
powers of $\varepsilon$ are confined to the $E$ matrix. To develop a better
conditioned basis, we need to factor out the ill-effects of the powers of
$\varepsilon$. To simplify the description of this step of the algorithm, we
make the assumption that $n=\mu_{0}(\mu_{0}+2)/2$ so that the sub-matrix block

$\displaystyle\left[\begin{array}[]{c|c|c|c}B_{1}&B_{2}&\cdots&B_{\mu_{0}}\end{array}\right]$
(25)

of $B$ in (22) and diagonal matrix $E_{1}$ in (23) are square and of size
$2n$-by-$2n$. We also assume that the interpolation node set
$\\{\mathbf{y}_{j}\\}_{j=1}^{n}$ is unisolvent with respect to the vector
spherical harmonics of degree $\mu_{0}$ so that (25) is invertible. This
restriction means that the interpolation matrix associated with the vector
spherical harmonic basis of degree $\mu_{0}$ is invertible. When using
“scattered” node sets, we have never encountered a situation where the point
set fails to be unisolvent.

The final step of the vector RBF-QR algorithm starts with a QR factorization
on $B$ in (20), which gives

$\widetilde{P}_{\rm div}\approx
Q\underbrace{\left[R_{1}\,|\,R_{2}\right]}_{\displaystyle
R}\left[\begin{array}[]{c|c}E_{1}&\\\ \hline\cr&E_{2}\\\ \end{array}\right]Y.$
(26)

Here we have partitioned $R$ into $R_{1}$ and $R_{2}$, where $R_{1}$ is
$2n$-by-$2n$ upper-triangular, and $R_{2}$ is an $2n$-by-$(m-2n)$ full matrix.
Since the sub-block (25) is invertible by our assumption on the node set, we
know $R_{1}$ is also invertible. This together with the block structure of $E$
allows us to re-write (26) as

$\displaystyle\widetilde{P}_{\rm div}\approx
QR_{1}\left[I_{2n}\,|\,R_{1}^{-1}R_{2}\right]\left[\begin{array}[]{c|c}E_{1}&\\\
\hline\cr&E_{2}\\\
\end{array}\right]Y=QR_{1}\left[E_{1}\,|\,R_{1}^{-1}R_{2}E_{2}\right]Y=QR_{1}E_{1}\underbrace{\left[I_{2n}\,|\,E_{1}^{-1}R_{1}^{-1}R_{2}E_{2}\right]}_{\displaystyle\widetilde{B}_{\rm
div}}Y.$ (29)

Since left matrix-multiplication is just a linear combination of the columns
of the matrices, it follows from this new expression that any element in the
span of the columns of $\widetilde{P}_{\rm div}^{T}$ (i.e., the span of the
original basis containing shifts of $\widetilde{\Phi}_{\rm div}$) can be
represented to machine precision by a linear combination of the columns of
$(\widetilde{B}_{\rm div}Y)^{T}$.

The form of $\widetilde{B}_{\rm div}Y$ in (29) is still not directly amenable
to computations with small $\varepsilon$ because it involves computing
$E_{1}^{-1}$. However, this matrix and $E_{2}$ are diagonal so we can
analytically remove the division by small $\varepsilon$ using [19, Lemma
5.1.2] to arrive at

$\widetilde{B}_{\rm
div}Y=\Bigl{[}I_{2n}\,|\,E_{1}^{-1}R_{1}^{-1}R_{2}E_{2}\Bigr{]}=\Bigl{[}I_{2n}\,|\,\left(R_{1}^{-1}R_{2}\right)\circ\underbrace{(E_{1}^{-1}J_{2n,m-2n}E_{2})}_{\displaystyle\widetilde{E}}\Bigr{]},$
(30)

where $J_{\sigma,\tau}$ is a $\sigma$-by-$\tau$ matrix of $1$’s, $\circ$
denotes the Hadamard product (or entry-wise multiplication), and the entries
of $\widetilde{E}$ can be determine explicitly as

$\widetilde{E}=\left[\begin{array}[]{c|c|cc|c}\varepsilon^{2\mu_{0}}J_{3,2\mu_{0}+3}&\varepsilon^{2\mu_{0}+2}J_{3,2\mu_{0}+5}&\cdots&\cdots&\varepsilon^{2\mu_{\text{trunc}}-2}J_{3,2\mu_{\text{trunc}}+1}\\\
\hline\cr\varepsilon^{2\mu_{0}-2}J_{5,2\mu_{0}+3}&\varepsilon^{2\mu_{0}}J_{5,2\mu_{0}+5}&\cdots&\cdots&\varepsilon^{2\mu_{\text{trunc}}-4}J_{5,2\mu_{\text{trunc}}+1}\\\
\hline\cr\vdots&\vdots&\ddots&\ddots&\vdots\\\
\hline\cr\varepsilon^{4}J_{2\mu_{0}-1,2\mu_{0}+3}&\varepsilon^{6}J_{2\mu_{0}-1,2\mu_{0}+5}&\cdots&\cdots&\varepsilon^{2\mu_{\text{trunc}}-2\mu_{0}+2}J_{2\mu_{0}-1,2\mu_{\text{trunc}}+1}\\\
\hline\cr\varepsilon^{2}J_{2\mu_{0}+1,2\mu_{0}+3}&\varepsilon^{4}J_{2\mu_{0}+1,2\mu_{0}+5}&\cdots&\cdots&\varepsilon^{2\mu_{\text{trunc}}-2\mu_{0}}J_{2\mu_{0}+1,2\mu_{\text{trunc}}+1}\end{array}\right].$
(31)

The columns $(\widetilde{B}_{\rm div}Y)^{T}$ in (30) can now be used as a
stable basis for the space spanned by $\\{\widetilde{{\Psi}}_{\rm
div}(\cdot,\mathbf{y}_{j})\\}_{j=1}^{n}$ for small $\varepsilon$. We note that
each element of this basis consists of a div-free vector spherical harmonic of
some degree $\leq\mu_{0}$ plus some $\mathcal{O}(\varepsilon^{2})$ combination
of div-free vector spherical harmonics of degree $>\mu_{0}$. This implies that
matrix-valued div-free RBF interpolant will converge to a div-free vector
spherical harmonic interpolant of degree $\mu_{0}$ in limit
$\varepsilon\rightarrow 0$, which is analogous to a scalar RBF interpolant
converging to a spherical harmonic interpolant in the flat limit [8].

######  Remark 3.1.

Note that just as with the scalar RBF-QR algorithm [8], it is possible to
include the spherical harmonic coefficients $c_{\mu,\varepsilon}$ in the
diagonal matrices $E_{1}$ and $E_{2}$ and generate a similar analytical
simplification for $\widetilde{E}$. This has the added advantage of removing
all $\varepsilon$ dependence in the actual QR numerical computation and allows
the central part of the vector RBF-QR to be performed independent of the
radial kernel used.

######  Remark 3.2.

When the condition $2n=\mu_{0}(\mu_{0}+2)$ is not met, the $R_{1}$ matrix in
(26) will not be square. The algorithm then needs to be modified to move
columns from $R_{1}$ to $R_{2}$ to make $R_{1}$ square. This requires a
similar move of the corresponding diagonals of $E_{1}$ to $E_{2}$. In this
case, the form of $\widetilde{E}$ in (31) also needs to change to include some
$\varepsilon^{0}$ terms.

######  Remark 3.3.

More sophisticated techniques for selecting the truncation level $\mu_{\rm
trunc}$ are discussed in [20] for the scalar case and may also be adopted
here. These methods can reduce the overall computational time without reducing
the approximation properties of the interpolants based on the stable basis.

##  4 Numerical Example

To test the vector RBF-QR algorithm, we use the target div-free field
displayed in Figure 2 (a). This field is generated from the stream function

$\displaystyle\begin{gathered}\psi(\mathbf{x})=-3z+2e^{-1.5((x-0.9)^{2}+(y+0.1)^{2})-8(z-0.2)^{2}}+3e^{-2((x+0.7)^{2}+(y-0.2)^{2})-8(z-0.25)^{2}}-\\\
2.5e^{-1.1((x+0.2)^{2}+(y-0.8)^{2})-8(z+0.19)^{2}}-2e^{-2.2((x+0.2)^{2}+(y+1)^{2})-8(z+0.21^{2})}\end{gathered}$
(34)

according to $\mathbf{u}=\mathbf{L}_{\mathbf{x}}\psi(\mathbf{x})$. For the
interpolation node set, we use the $n=924$ scattered nodes displayed in Figure
2 (b). These are an example of the Hammersley nodes, which give well-
distributed, but random sampling points for the sphere [21], and were obtained
from the SpherePts software package [22]. The number $n=924$ corresponds to a
truncated div-free vector spherical harmonic expansion of $\mu_{0}=42$, which
is commonly used in scalar spherical harmonic tests where it is denoted as
“T42” [8]. We only present results for the MQ kernel, but note that similar
results were obtained for other smooth kernels.

|
---|---
(a) Target field | (b) Interpolation nodes
Figure 2: Hammer-Aitoff projection [14] of the (a) target div-free field and
(b) $n=924$ node set used in the numerical example.

Using this target field and node set, we computed the surface div-free RBF
interpolants using the RBF-Direct and vector RBF-QR algorithms for various
values of $\varepsilon$. The relative max-norm error in the interpolant was
then computed by evaluating the difference between the interpolants and the
target field at a denser set of $4n$ points over the sphere. The results from
the experiment are displayed in Figure 3 (a). We see from the figure that the
error in the RBF-Direct approach decreases rapidly with decreasing
$\varepsilon$ until around $\varepsilon=1$ when it starts to increase
exponentially fast. This is where ill-conditioning in the standard div-free
RBF basis sets in. The vector RBF-QR algorithm on the other hand, shows no
issues with ill-conditioning for any values of $\varepsilon\leq 1$ and we can
use it to compute the resulting interpolant in a stable manner all the way to
the flat limit. Note that the error reaches a minimum at a non-zero value of
$\varepsilon$ and then starts to increase slightly. This is a feature that
comes from the target field, and not from instabilities in the algorithm.
Figure 3 (b) shows the max-norm errors in approximating the stream function
(34) as a function of $\varepsilon$. For this computation, we used (7) to
extract a stream function from the the interpolants and then adjusted it to
have the same mean as the target stream function (34). We see that there is a
similar trend in these results, which is to be expected since the stream
function comes from the interpolant.

|
---|---
(a) Field $\mathbf{u}$ | (b) Stream function $\psi$
Figure 3: Comparison of the relative max-norm errors in reconstructing the
target (a) div-free field and (b) stream function (34) as a function of
$\varepsilon$ using the vector RBF-QR algorithm and RBF-Direct for the MQ
kernel.

##  5 Concluding remarks

The vector RBF-QR algorithm presented here can be used to bypass the ill-
conditioning associated with surface div-free RBF interpolation on a sphere in
the $\varepsilon\rightarrow 0$ limit. This allows for more comprehensive
studies of how $\varepsilon$ affects the accuracy of the interpolants. The
derivation of the algorithm additionally demonstrates the connection between
these interpolants and div-free vector spherical harmonics in the flat limit.
The algorithm can also be applied straightforwardly to surface curl-free RBF
interpolation, since in this case, one simply has to replace the matrix-valued
kernel $\Phi_{\rm div}$ with the kernel $\Phi_{\rm
curl}=(P({\mathbf{x}})\nabla_{\mathbf{x}})(P({\mathbf{y}})\nabla_{\mathbf{y}})^{T}\phi(\|\mathbf{x}-\mathbf{y}\|)$,
where $P(\boldsymbol{\xi})=(I-\boldsymbol{\xi}\boldsymbol{\xi}^{T})$ [11].
Div-free and curl-free RBFs are also available for interpolation in
$\mathbb{R}^{d}$ [23]. For this problem it may be possible to extend the RBF-
QR algorithms [7, 9] developed for scalar RBF interpolation in
$\mathbb{R}^{d}$. However, the algorithm in [7] is limited to the Gaussian
kernel and the algorithm in [9] is limited to kernels with known Mercer series
expansions.

## Acknowledgments

KPD’s work was partially supported by the SMART Scholarship, which is funded
by The Under Secretary of Defense-Research and Engineering, National Defense
Education Program / BA-1, Basic Research. GBW’s work was partially supported
by National Science Foundation grant CCF 1717556.

## References

  * [1] M. Fan, D. Paul, T. C. M. Lee, T. Matsuo, Modeling tangential vector fields on a sphere, Journal of the American Statistical Association 113 (524) (2018) 1625–1636.
  * [2] D. L. Williamson, The evolution of dynamical cores for global atmospheric models, JMSJ 85B (2007) 241–269.
  * [3] F. J. Narcowich, J. D. Ward, G. B. Wright, Divergence-free RBFs on surfaces, Fourier Anal. Appl. 13 (2007) 643–663.
  * [4] E. J. Fuselier, G. B. Wright, Stability and error estimates for vector field interpolation and decomposition on the sphere with RBFs, SIAM J. Num. Anal. 47 (2009) 3213–3239.
  * [5] B. Fornberg, G. B. Wright, Stable computation of multiquadric interpolants for all values of the shape parameter, Comp. Math. Appl. 48 (2004) 853–867.
  * [6] G. B. Wright, B. Fornberg, Stable computation with flat radial basis functions using vector-valued rational approximations, J. Comput. Phys. 331 (2017) 137–156.
  * [7] B. Fornberg, E. Larsson, N. Flyer, Stable computations with Gaussian radial basis functions, SIAM J. Sci. Comput. 33 (2011) 869–892.
  * [8] B. Fornberg, C. Piret, A stable algorithm for flat radial basis functions on a sphere, SIAM J. Sci. Comput. 30 (2007) 60–80.
  * [9] G. E. Fasshauer, M. J. McCourt, Stable evaluation of Gaussian radial basis function interpolants, SIAM J. Sci. Comput. 34 (2012) A737–A762.
  * [10] B. Fornberg, E. Lehto, C. Powell, Stable calculation of Gaussian-based RBF-FD stencils, Comp. Math. Applic. 65 (2013) 627–637.
  * [11] K. P. Drake, A stable algorithm for divergence-free and curl-free radial basis functions in the flat limit, Master’s thesis, Boise State University, USA, 1290 (2017).  doi:https://doi.org/10.18122/B2G70C.
  * [12] E. J. Fuselier, F. J. Narcowich, J. D. Ward, G. B. Wright, Error and stability estimates for divergence-free RBF interpolants on the sphere, Math. Comp. 78 (2009) 2157–2186.
  * [13] G. E. Fasshauer, Meshfree Approximation Methods with MATLAB, Interdisciplinary Mathematical Sciences - Vol. 6, World Scientific Publishers, Singapore, 2007.
  * [14] J. P. Snyder, Flattening the Earth: Two Thousand Years of Map Projections, University of Chicago Press, Chicago, 1993.
  * [15] K. Atkinson, W. Han, Spherical Harmonics and Approximations on the Unit Sphere: An Introduction, Springer Berlin Heidelberg, 2012.
  * [16] P. N. Swarztrauber, The vector harmonic transform method for solving partial differential equations in spherical geometry, Mon. Wea. Rev. 121 (1993) 3415–3437.
  * [17] C. Müller, Spherical Harmonics: Lecture Notes in Mathematics, Vol. 17, Springer Berlin Heidelberg, New York, 1966.
  * [18] S. Hubbert, B. J. C. Baxter, Radial basis functions for the sphere, in: W. Haussmann, K. Jetter, M. Reimer (Eds.), Recent Progress in Multivariate Approximation, Proc. of the 4th Intern. Conf., Witten-Bommerholz, Germany, Vol. 137 of International Series of Numerical Mathematics, Birkhäuser, Basel, 2001, pp. 33–47.
  * [19] R. Horn, C. R. Johnson, Topics in matrix analysis, Cambridge University Press, Cambridge, 1991.
  * [20] K. Kormann, C. Lasser, A. Yurova, Stable interpolation with isotropic and anisotropic Gaussians using Hermite generating function, SIAM J. Sci. Comput. 41 (6) (2019) A3839–A3859.  doi:10.1137/19M1262449.
  * [21] J. Cui, W. Freeden, Equidistribution on the sphere, SIAM J. Sci. Comput. 18 (2) (1997) 595–609.
  * [22] G. B. Wright, SpherePts, https://github.com/gradywright/spherepts/ (2016).
  * [23] E. Fuselier, Sobolev-type approximation rates for divergence-free and curl-free RBF interpolants, Math. Comp. 77 (2008) 1407–1423.