# Accurate Stabilization Techniques for RBF-FD Meshless Discretizations with
Neumann Boundary Conditions

Riccardo Zamolo<EMAIL_ADDRESS>Davide Miotti  Enrico Nobile

###### Abstract

A major obstacle to the application of the standard Radial Basis Function-
generated Finite Difference (RBF-FD) meshless method is constituted by its
inability to accurately and consistently solve boundary value problems
involving Neumann boundary conditions (BCs). This is also due to ill-
conditioning issues affecting the interpolation matrix when boundary
derivatives are imposed in strong form. In this paper these ill-conditioning
issues and subsequent instabilities affecting the application of the RBF-FD
method in presence of Neumann BCs are analyzed both theoretically and
numerically. The theoretical motivations for the onset of such issues are
derived by highlighting the dependence of the determinant of the local
interpolation matrix upon the boundary normals. Qualitative investigations are
also carried out numerically by studying a reference stencil and looking for
correlations between its geometry and the properties of the associated
interpolation matrix. Based on the previous analyses, two approaches are
derived to overcome the initial problem. The corresponding stabilization
properties are finally assessed by succesfully applying such approaches to the
stabilization of the Helmholtz-Hodge decomposition.

###### keywords:

RBF-FD , Meshless , Neumann boundary conditions , Stabilization

###### MSC:

65N06 , 65N12

††journal: Elsevier

[inst1] organization=Department of Engineering and Architecture, University of
Trieste, via Alfonso Valerio 10, 34127, city=Trieste, country=Italy

##  1 Introduction

Some of the main issues with traditional CAE (Computer Aided Engineering)
solvers, i.e. those based on Finite Volume (FVM) [1] and Finite Elements (FEM)
methods [2], are caused by the limited flexibility of the mesh data structure
and the complexity of its generation [3]. In order to leverage the increased
CPU capabilities and allow the user to focus on other aspects of the
simulation process, meshless methods [4, 5, 6, 7, 8, 9] might be adopted as a
replacement for the more traditional mesh based ones.

Before they can be considered as an improvement over the traditional ones,
however, meshless methods must guarantee reliable results also in presence of
generic geometries, e.g., arbitrarily complex or affected by deformations, and
boundary conditions, e.g., Neumann boundary conditions in particular. This
would allow the user to skip the mesh generation process altogether without
sacrificing the accuracy or the stability of the solution procedure.

In the present paper the Radial Basis Function - Finite Difference (RBF-FD)
collocation method augmented with polynomials is employed, a brief discussion
on the main features is presented in the following sections, however an
interested reader is encouraged to also see [10].

The discretization of the domain is usually achieved by placing a number of
nodes both on the boundary and on the interior, when the RBF-FD collocation
method is applied there is no need to store connectivity information and
partial differential equations can be solved in strong form directly [3].

In practical cases the basic RBF-FD method has to be enhanced in order to
ensure good convergence close to the boundaries [11]. Among others, the most
popular improvement is polynomial augmentation which allows to reach higher
convergence rates and to improve interpolation near boundaries as long as
Dirichlet boundary conditions (BCs) are considered [12, 13, 14, 15]. The
accuracy of this approximation can therefore be controlled by properly
increasing the polynomial order and the number of nodes.

Neumann BCs can be enforced both at the stencil level, i.e., when the local
interpolants exactly satisfy the BCs at the boundary nodes [16, 17, 18, 10,
19], or at the assembly level, i.e., when the equations for the BCs at the
boundary nodes appear explicitly in the final sparse matrix [20, 21]. In both
cases the RBF-FD collocation method, however, may lead to large errors and/or
ill-conditioning problems. These issues constitute a major obstacle for the
robust and accurate solution of boundary value problems of practical interest.
Some effective solutions for this last problem have been suggested, however no
one has been universally adopted yet.

In [22, 23, 24] the stabilization is achieved through different least squares
procedures. The one in [24] seems especially capable of providing both stable
and accurate results by using two different node sets.

Another approach to the stabilization is the one based on the so called ghost
nodes. In [25] one layer of ghost nodes is placed outside the domain at a
distance roughly continuing the pattern of the interior nodes in order to
address the issues related to one-sided boundary stencils. In [26] the authors
propose instead the use of fictitious nodes that extend far outside the
domain, resembling the approach from the method of fundamental solutions [27].
In another implementation [28] the authors add a set of nodes, which can lie
inside or outside of the domain, adjacent to the boundary and,
correspondingly, add an additional set of equations obtained via collocation
of the PDE on the boundary. The strategy of ghost nodes is successfully
employed also in [29], for the stable RBF-FD solution of elliptic PDEs
problems on complex 3D geometries and even in higher dimension in [21].

Alternatively to the ghost nodes, in [30] additional nodes are placed close to
the boundary, but within the domain, in order to reduce the error in those
areas in the context of the overlapped RBF-FD method. The overlapped variant
differs from the traditional RBF-FD in the selection of the interpolation
stencils and can be intended as a generalization of the latter.

Yet another approach is adopted in [31, 16], where every time a Neumann
boundary node has to be included in the stencil, its normal is checked
according to different geometric criteria.

The aim of the present paper is to present some remarks on the roots of those
issues related to the enforcement of Neumann BCs at the stencil level and to
suggest new improvements of the traditional interpolation scheme. It is
important to remark that the final implementation must be capable of dealing
with any possible shape of the boundary, as a consequence it is not possible
to make any a priori assumption on the direction of the normals, i.e., the
outward facing vectors defined at each boundary node and orthogonal to the
local tangent plane.

Building on top of the simplest RBF-FD scheme, the main source of errors and
stability issues is found in the poor conditioning of the local interpolation
matrix in presence of Neumann BCs. We found that, given a certain node
distribution, it is always possible to have normals that make the local
interpolation matrix singular. Such singular directions can be calculated as
the ones that send the determinant of the said matrix to zero, therefore a
control can be implemented in order not to include the corresponding boundary
nodes in the stencil. This approach, somewhat similar to the one adopted in
[16], is more expensive but totally robust and it is tested with good results.
At the best of the authors’ knowledge it represents the only possible solution
that carries on the traditional RBF-FD method totally unchanged and doesn’t
require any intervention on the node generation. The presented mathematical
formulation leaves great freedom for further developments aimed at improving
both accuracy and stability of the RBF-FD scheme furthermore.

Alternatively, we found that it is possible to properly move the nodes on the
boundary in order to rule out the possibility of singular interpolation
matrices. This operations appears to be the most computationally efficient and
the simplest to implement on smooth and regular domains, although requiring
non-trivial geometrical operations on generic domains. However, a robust
implementation of this approach on complex-shaped domains is not simple and
would require a proper modification of the node generation algorithm.

The proposed strategies for the prevention of ill-conditioning problems due to
Neumann BCs, which doesn’t require human intervention, paves the way for the
implementation of robust RBF-FD meshless methods in combination with fully
automated node generation.

##  2 RBF-FD overview

Follows a brief description of the RBF-FD method which constitutes the
starting point from which all the variants discussed above are developed.

The actual solution of a boundary value problem can be split in $3$ steps:

  1. 1.

node generation;

  2. 2.

discretization through the RBF-FD method;

  3. 3.

solution procedure for the resulting sparse linear system.

###  2.1 Node generation

The employed node generation procedure consists of two phases:

  1. 1.

Generation of a node distribution within the domain $\Omega$ that satisfies a
defined spacing function $s(\boldsymbol{x})$ on average. The boundary
$\partial\Omega$ is not taken into account yet.

  2. 2.

Iterative refinement of the initial node distribution. At this stage the
already existing nodes can be moved and projected onto $\partial\Omega$
appropriately.

Implementation details fall outside the topic of this paper, more details can
be found in [32]. Throughout this paper we will assume that nodes are evenly
distributed, this can be attained with an adequate number of refinement
iterations. This allows us to assume a negligible difference between fill
distance and separation distance [33, 34]. As a consequence such a difference
can reasonably be ignored as a source of interpolation instabilities.

###  2.2 RBF-FD method

Smooth | Type |  Definition $\Phi(r)$ | Positive definite, order
---|---|---|---
Infinitely | Gaussian (GA) | $e^{-(\varepsilon\,r)^{2}}$ | strictly
| Multiquadric (MQ) | $\sqrt{1+(\varepsilon\,r)^{2}}$ | conditionally, 1
| Inverse Multiquadric (IMQ) | $1/\sqrt{1+(\varepsilon\,r)^{2}}$ | strictly
| Inverse Quadratic (IQ) | $1/(1+(\varepsilon\,r)^{2})$ | strictly
Piecewise | Monomial PHS |  $r^{2k+1}$, $k\in\mathbb{N}$ |  conditionally, $k+1$
| Thin Plate Spline PHS |  $r^{2k}\log r$, $k\in\mathbb{N}$ |  conditionally, $k+1$
Table 1: most common Radial Basis Functions

In the simplest version of this method, any real valued function can be
interpolated at a generic point $\boldsymbol{x}$ as a linear combination of
$m$ Radial Basis Functions (RBFs)
$\Phi(||\boldsymbol{x}-\boldsymbol{x}_{i}||_{2})$, and the set
$\\{\boldsymbol{x}_{i}\\}_{1}^{m}$ is the interpolation stencil, composed by
distinct nodes. The most commonly used RBFs are reported in Table 1, where
they are defined in terms of $r=||\boldsymbol{x}-\boldsymbol{x}_{i}||_{2}$,
for $\boldsymbol{x}_{i}$ fixed.

By using the notation
$\varphi_{i}(\boldsymbol{x})=\Phi(||\boldsymbol{x}-\boldsymbol{x}_{i}||_{2})$,
the interpolation $u^{h}(\boldsymbol{x})$ of a function $u(\boldsymbol{x})$
thus takes the following form:

$u^{h}(\boldsymbol{x})=\sum_{i=1}^{m}\alpha_{i}\varphi_{i}(\boldsymbol{x})$
(1)

In practical cases the previous interpolation has to be enhanced in order to
ensure appropriate convergence properties [11]. A very popular improvement is
polynomial augmentation, where the bare RBF interpolation scheme in equation
(1) is modified by including the polynomial basis $\\{p_{j}\\}_{1}^{q}$ of the
multivariate polynomial space $\Pi_{P}^{d}$ of total degree $P$ in $d$
dimensions [12] as follows:

$u^{h}(\boldsymbol{x})=\sum_{i=1}^{m}\alpha_{i}\varphi_{i}(\boldsymbol{x})+\sum_{j=1}^{q}\beta_{j}p_{j}(\boldsymbol{x})$
(2)

We assume the nodes in the interpolation stencil are listed in the following
order: the first $m_{I}$ are internal nodes, i.e., nodes lying within
$\Omega$, while the following $m_{B}$ are boundary nodes, i.e., those lying on
$\partial\Omega$, and $m_{I}+m_{B}=m$.

The following interpolation conditions at the internal nodes hold:

$u^{h}(\boldsymbol{x}_{i})=u(\boldsymbol{x}_{i})\,,\quad i=1,\dots,m_{I}$ (3)

Since we will focus only on Neumann BCs, i.e., $\partial
u/\partial\boldsymbol{n}=g$, we will use the following notation for the normal
derivative of a generic function $F$ at the boundary node
$\boldsymbol{x}_{k}$:

$\partial F(\boldsymbol{x}_{k}):=\frac{\partial
F}{\partial\boldsymbol{n}_{k}}(\boldsymbol{x}_{k})$ (4)

where $\boldsymbol{n}_{k}$ is the unit normal at $\boldsymbol{x}_{k}$.

The interpolation conditions (3) are then supplemented with Neumann BCs at the
boundary nodes as follows:

$\partial u^{h}(\boldsymbol{x}_{k})=\partial
u(\boldsymbol{x}_{k})=g(\boldsymbol{x}_{k})\,,\quad k=m_{J},\dots,m$ (5)

where $m_{J}=m_{I}+1$.

In the case of polynomial augmentation, the following orthogonality conditions
are enforced [35]:

$\sum_{i=1}^{m}\alpha_{i}p_{j}(\boldsymbol{x}_{i})=0\,,\quad j=1,\dots,q$ (6)

By collecting equations (3) and (5)-(6) in matrix form, the following linear
system is attained:

$\underbrace{\begin{bmatrix}\boldsymbol{\varphi}_{BC}&\mathbf{P}_{BC}\\\
\mathbf{P}^{T}&\boldsymbol{0}\end{bmatrix}}_{\displaystyle\mathbf{M}}\begin{Bmatrix}\boldsymbol{\alpha}\\\
\boldsymbol{\beta}\end{Bmatrix}=\begin{Bmatrix}\boldsymbol{u}\\\
\boldsymbol{g}\\\ \boldsymbol{0}\end{Bmatrix}$ (7)

where $\boldsymbol{\alpha}=\\{\alpha_{i}\\}_{1}^{m}$,
$\boldsymbol{\beta}=\\{\beta_{j}\\}_{1}^{q}$,
$\boldsymbol{u}=\\{u(\boldsymbol{x}_{i})\\}_{1}^{m_{I}}$ and
$\boldsymbol{g}=\\{g(\boldsymbol{x}_{k})\\}_{m_{J}}^{m}$ are column vectors.
The matrix $\mathbf{M}$ in equation (7) will be called the interpolation
matrix and its blocks are expanded as follows:

$\boldsymbol{\varphi}_{BC}=\begin{bmatrix}\varphi_{1}(\boldsymbol{x}_{1})\\!\\!\\!&\cdots\\!\\!\\!&\varphi_{m}(\boldsymbol{x}_{1})\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
\varphi_{1}(\boldsymbol{x}_{m_{I}})\\!\\!\\!&\cdots\\!\\!\\!&\varphi_{m}(\boldsymbol{x}_{m_{I}})\\\
\partial\varphi_{1}(\boldsymbol{x}_{m_{J}})\\!\\!\\!&\cdots\\!\\!\\!&\partial\varphi_{m}(\boldsymbol{x}_{m_{J}})\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
\partial\varphi_{1}(\boldsymbol{x}_{m})\\!\\!\\!&\cdots\\!\\!\\!&\partial\varphi_{m}(\boldsymbol{x}_{m})\end{bmatrix}$
(8)

$\mathbf{P}_{BC}=\begin{bmatrix}p_{1}(\boldsymbol{x}_{1})\\!\\!\\!&\cdots\\!\\!\\!&p_{q}(\boldsymbol{x}_{1})\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
p_{1}(\boldsymbol{x}_{m_{I}})\\!\\!\\!&\cdots\\!\\!\\!&p_{q}(\boldsymbol{x}_{m_{I}})\\\
\partial p_{1}(\boldsymbol{x}_{m_{J}})\\!\\!\\!&\cdots\\!\\!\\!&\partial
p_{q}(\boldsymbol{x}_{m_{J}})\\\ \vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
\partial p_{1}(\boldsymbol{x}_{m})\\!\\!\\!&\cdots\\!\\!\\!&\partial
p_{q}(\boldsymbol{x}_{m})\end{bmatrix}$ (9)

$\mathbf{P}^{T}=\begin{bmatrix}p_{1}(\boldsymbol{x}_{1})\\!\\!\\!&\cdots\\!\\!\\!&p_{1}(\boldsymbol{x}_{m})\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
p_{q}(\boldsymbol{x}_{1})\\!\\!\\!&\cdots\\!\\!\\!&p_{q}(\boldsymbol{x}_{m})\end{bmatrix}$
(10)

The normal derivatives of RBFs in equation (8) can be expressed as follows:

$\partial\varphi_{i}(\boldsymbol{x}_{k})=\frac{\partial\varphi_{i}}{\partial\boldsymbol{n}_{k}}(\boldsymbol{x}_{k})=\Phi^{\prime}(r_{i,k})\,\boldsymbol{e}_{i,k}\cdot\boldsymbol{n}_{k}$
(11)

where $r_{i,k}=||\boldsymbol{x}_{k}-\boldsymbol{x}_{i}||_{2}$ and
$\boldsymbol{e}_{i,k}=(\boldsymbol{x}_{k}-\boldsymbol{x}_{i})/r_{i,k}$ is the
unit direction from node $i$ towards node $k$. For the problem of interest,
i.e., solution of PDEs, we also require the RBFs to satisfy equation (12) for
$\varphi_{i}$ to be differentiable also at the nodes:

$\Phi^{\prime}(0)=0$ (12)

therefore, the diagonal entries $\partial\varphi_{i}(\boldsymbol{x}_{i})$ of
the last $m_{B}$ rows of $\boldsymbol{\varphi}_{BC}$ are $0$.

Given a linear partial differential equation $\mathfrak{L}(u)=f$ in the
unknown field $u$, the Weighted Residual Method is applied in combination with
the Collocation technique as explained in [3]. The linear operator
$\mathfrak{L}$ is therefore applied to $u^{h}$ and the following equation is
made valid at each internal node $\boldsymbol{x}$:

$\mathfrak{L}u^{h}(\boldsymbol{x})=f(\boldsymbol{x})$ (13)

When we substitute the definition of $u^{h}$ given in equation (2) into
equation (13), the linear operator can be applied to the basis functions as
follows:

$\sum_{i=1}^{m}\alpha_{i}\,\mathfrak{L}\varphi_{i}(\boldsymbol{x})+\sum_{j=1}^{q}\beta_{j}\,\mathfrak{L}p_{j}(\boldsymbol{x})=f(\boldsymbol{x})$
(14)

Equation (14) can be recast in vector form as follows, where
$\boldsymbol{\Psi}(\boldsymbol{x})=\\{\mathfrak{L}\varphi_{i}(\boldsymbol{x})\\}_{1}^{m}$
and
$\boldsymbol{\Pi}(\boldsymbol{x})=\\{\mathfrak{L}p_{j}(\boldsymbol{x})\\}_{1}^{q}$
are column vectors:

$\begin{Bmatrix}\boldsymbol{\Psi}(\boldsymbol{x})\\\
\boldsymbol{\Pi}(\boldsymbol{x})\end{Bmatrix}^{T}\begin{Bmatrix}\boldsymbol{\alpha}\\\
\boldsymbol{\beta}\end{Bmatrix}=f(\boldsymbol{x})$ (15)

By substituting vectors $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ from
equation (7) into equation (15) we obtain:

$\underbrace{\begin{Bmatrix}\boldsymbol{\Psi}(\boldsymbol{x})\\\
\boldsymbol{\Pi}(\boldsymbol{x})\end{Bmatrix}^{T}\mathbf{M}^{-1}}_{\displaystyle\boldsymbol{c}^{T}}\begin{Bmatrix}\boldsymbol{u}\\\
\boldsymbol{g}\\\ \boldsymbol{0}\end{Bmatrix}=f(\boldsymbol{x})$ (16)

We point out that, instead of calculating the inverse matrix
$\mathbf{M}^{-1}$, we can solve the following adjoint system for the
coefficients $\boldsymbol{c}$:

$\mathbf{M}^{T}\boldsymbol{c}=\begin{Bmatrix}\boldsymbol{\Psi}(\boldsymbol{x})\\\
\boldsymbol{\Pi}(\boldsymbol{x})\end{Bmatrix}$ (17)

Given an interpolation stencil built around point $\boldsymbol{x}$, the vector
$\boldsymbol{c}$ represents therefore the stencil coefficients for the RBF-FD
approximation of the linear operator $\mathfrak{L}$ at point $\boldsymbol{x}$.
Only the first $m$ elements of $\boldsymbol{c}$ are kept, i.e., those
multiplying nonzero elements in equation (16). These elements of
$\boldsymbol{c}$ are then plugged in the corresponding row of a large
$N_{I}\times N_{I}$ sparse matrix $\mathbf{L}$, $N_{I}$ being the total number
of internal nodes, leading to the following linear system:

$\mathbf{L}\boldsymbol{u}^{h}=\boldsymbol{q}+\boldsymbol{f}$ (18)

where $\boldsymbol{q}$ comes from the boundary conditions, i.e., from those
terms of $\boldsymbol{c}$ which are associated to boundary nodes in equation
(16), while $\boldsymbol{u}^{h}$ and $\boldsymbol{f}$ are the solution vector
and the vector of the known term $f$, respectively, at the internal nodes.
Equation (18), which is the discrete version of the PDE $\mathfrak{L}u=f$ we
started with, can now be solved for $\boldsymbol{u}^{h}\in\mathbb{R}^{N_{I}}$.

We point out that, if $\mathbf{M}$ is ill conditioned in equation (7), so is
$\mathbf{M}^{T}$ of equation (17). We will show that the onset of
instabilities at the boundary is due to ill conditioning of the interpolation
matrix $\mathbf{M}$, and that such ill conditioning is essentially influenced
by the direction of the boundary normals appearing in the corresponding rows,
see equations (8) and (9).

The final sparse linear system of equation (18) can now be solved using a
proper iterative method. Because of the lack of connectivity information,
local interpolation systems $\eqref{eq:adjoint}$ can be solved independently
from one another.

###  2.3 Cardinal functions and Lebesgue constant

The interpolant $u^{h}$ can also be obtained by considering the identity
operator $\mathfrak{L}=\mathcal{I}$ in the left side of equation (13). Then,
by considering the left side of equation (16), we obtain explicitly:

$\mathcal{I}u^{h}(\boldsymbol{x})=u^{h}(\boldsymbol{x})=\sum_{i=1}^{m_{I}}u(\boldsymbol{x}_{i})\psi_{i}(\boldsymbol{x})+\\!\sum_{k=m_{J}}^{m}\\!\partial
u(\boldsymbol{x}_{k})\psi_{k}(\boldsymbol{x})$ (19)

where the cardinal functions $\\{\psi_{i}(\boldsymbol{x})\\}_{1}^{m}$ are
therefore given by the first $m$ elements of $\boldsymbol{c}$ in equation (17)
in the case $\mathfrak{L}=\mathcal{I}$, i.e., when the right side of the same
equation is the column vector of the bare basis functions:

$\mathbf{M}^{T}\begin{Bmatrix}\\{\psi_{i}(\boldsymbol{x})\\}_{1}^{m}\\\
\boldsymbol{0}\end{Bmatrix}=\begin{Bmatrix}\\{\varphi_{i}(\boldsymbol{x})\\}_{1}^{m}\\\
\\{p_{j}(\boldsymbol{x})\\}_{1}^{q}\\\ \end{Bmatrix}$ (20)

These cardinal functions extend the definition given in [33, 36] to the case
with Neumann BCs. The cardinal functions associated to the internal nodes,
i.e., $\\{\psi_{i}(\boldsymbol{x})\\}_{1}^{m_{I}}$, satisfy the following
properties:

$\begin{cases}\begin{aligned}
\psi_{i}(\boldsymbol{x}_{j})&=\delta_{ij}&&\text{if $j\leq m_{I}$}\\\
\partial\psi_{i}(\boldsymbol{x}_{j})&=0&&\text{otherwise}\end{aligned}\end{cases}$
(21)

while the ones associated to the boundary nodes, i.e.,
$\\{\psi_{k}(\boldsymbol{x})\\}_{m_{J}}^{m}$, satisfy the following
properties:

$\begin{cases}\begin{aligned} \psi_{k}(\boldsymbol{x}_{j})&=0&&\text{if $j\leq
m_{I}$}\\\
\partial\psi_{k}(\boldsymbol{x}_{j})&=\delta_{kj}&&\text{otherwise}\end{aligned}\end{cases}$
(22)

where $\delta_{ij}$ is the Kronecker delta.

In order to decouple the contribution of internal and boundary nodes
explicitly, we define two Lebesgue functions:

$\begin{split}\lambda_{I}(\boldsymbol{x})&=\sum_{i=1}^{m_{I}}|\psi_{i}(\boldsymbol{x})|\\\
\lambda_{B}(\boldsymbol{x})&=\sum_{k=m_{J}}^{m}|\psi_{k}(\boldsymbol{x})|\end{split}$
(23)

where $\lambda_{I}$ is associated to the internal nodes, whereas $\lambda_{B}$
is associated to the boundary nodes. This can be seen as an extension of the
traditional definition of a single Lebesgue function [33].

When the interpolation stencil is restricted to lie in a compact set
$K\subset\mathbb{R}^{d}$, then the Lebesgue constant is defined as the maximum
of the Lebesgue function [33, 36]. In our case we will consider $K$ to be the
minimal convex set containing all of the stencil nodes (including those lying
on the boundary). In order to keep the distinction between boundary nodes and
inner nodes, we define 2 constants:

$\begin{split}\Lambda_{I}&=\max_{\boldsymbol{x}\in
K}\lambda_{I}(\boldsymbol{x})\\\ \Lambda_{B}&=\max_{\boldsymbol{x}\in
K}\lambda_{B}(\boldsymbol{x})\end{split}$ (24)

The Lebesgue constant is especially useful since it provides an upper bound
for the interpolation error, which in turn is closely related to the
conditioning of matrix $\mathbf{M}^{T}$. From equation (19) $u^{h}$ can be
seen as the image of a function $u$ with respect to the linear interpolation
operator $u\mapsto u^{h}$. We can find an upper bound for
$||u^{h}||_{{\infty}}$ on $K$ as follows:

$\begin{split}||u^{h}||_{{\infty}}&=\max_{x\in
K}\left|\sum_{i=1}^{m_{I}}u(\boldsymbol{x}_{i})\psi_{i}(\boldsymbol{x})+\sum_{k=m_{J}}^{m}\partial
u(\boldsymbol{x}_{k})\psi_{k}(\boldsymbol{x})\right|\\\
&\leq\max_{\boldsymbol{x}\in
K}\left|\sum_{i=1}^{m_{I}}u(\boldsymbol{x}_{i})\psi_{i}(\boldsymbol{x})\right|+\max_{\boldsymbol{x}\in
K}\left|\sum_{k=m_{J}}^{m}\partial
u(\boldsymbol{x}_{k})\psi_{k}(\boldsymbol{x})\right|\\\ &\leq\left(\max_{i\leq
m_{I}}|u(\boldsymbol{x}_{i})|\right)\Lambda_{I}+\left(\max_{m_{J}\leq k\leq
m}|\partial u(\boldsymbol{x}_{k})|\right)\Lambda_{B}\\\
&\leq||u||_{{\infty}}\Lambda_{I}+\underbrace{\max_{\boldsymbol{x}\in
K}||\nabla u(\boldsymbol{x})||_{2}}_{\displaystyle
C_{\infty}(u)}\Lambda_{B}\end{split}$ (25)

where in the last inequality the Euclidean length of the gradient $\nabla u$
is employed as the upper bound for the absolute value of the normal derivative
and $C_{\infty}(u)$ is the $L^{\infty}$-norm on $K$ of the Euclidean length of
$\nabla u$. When the stencil does not include any boundary node,
$\Lambda_{B}=0$ and equation (25) takes the usual form [33, 36]:

$||u^{h}||_{{\infty}}\leq\Lambda_{I}||u||_{{\infty}}$ (26)

where $\Lambda_{I}$ can therefore be seen as an upper bound for the operator
norm of the RBF interpolation operator.

In the general case, from equation (25) we can derive an estimate for the
interpolation error; assume that $u^{*}$ is the best RBF approximation
possible of $u$, while $u^{h}$ is a generic approximation:

$\begin{split}||u-u^{h}||_{{\infty}}&\leq||u-u^{*}||_{{\infty}}+||u^{*}-u^{h}||_{{\infty}}\\\
&=||u-u^{*}||_{{\infty}}+||(u-u^{*})^{h}||_{{\infty}}\\\
&\leq(1+\Lambda_{I})||u-u^{*}||_{{\infty}}+\Lambda_{B}\underbrace{\max_{\boldsymbol{x}\in
K}||\nabla(u-u^{*})(\boldsymbol{x})||_{2}}_{\displaystyle
C_{\infty}(u-u^{*})}\end{split}$ (27)

Once again, in the case without boundary nodes we have $\Lambda_{B}=0$.

In the attempt of reaching a more compact form, similar to the one in (26), we
now assume that equation (12) holds. Using multi-index notation and the
definition of the $||\cdot||_{W^{1,\infty}}$ norm given in [37] we can write:

$\begin{split}||u^{h}||_{W^{1,\infty}}&=\sum_{|\alpha|\leq
1}\max_{\boldsymbol{x}\in
K}\left|\sum_{i=1}^{m_{I}}u(\boldsymbol{x}_{i})\psi_{i}^{(\alpha)}(\boldsymbol{x})+\sum_{k=m_{J}}^{m}\partial
u(\boldsymbol{x}_{k})\psi_{k}^{(\alpha)}(\boldsymbol{x})\right|\\\
&\leq\sum_{|\alpha|\leq 1}\max_{\boldsymbol{x}\in
K}\left(\left|\sum_{i=1}^{m_{I}}u(\boldsymbol{x}_{i})\psi_{i}^{(\alpha)}(\boldsymbol{x})\right|+\left|\sum_{k=m_{J}}^{m}\partial
u(\boldsymbol{x}_{k})\psi_{k}^{(\alpha)}(\boldsymbol{x})\right|\right)\\\
&\leq\sum_{|\alpha|\leq 1}\left(\max_{i\leq
m_{I}}|u(\boldsymbol{x}_{i})|\right)\Lambda_{I}^{(\alpha)}+\sum_{|\alpha|\leq
1}\left(\max_{m_{J}\leq k\leq m}|\partial
u(\boldsymbol{x}_{k})|\right)\Lambda_{B}^{(\alpha)}\\\ &=\left(\max_{i\leq
m_{I}}|u(\boldsymbol{x}_{i})|\right)\sum_{|\alpha|\leq
1}\Lambda_{I}^{(\alpha)}+\left(\max_{m_{J}\leq k\leq m}|\partial
u(\boldsymbol{x}_{k})|\right)\sum_{|\alpha|\leq 1}\Lambda_{B}^{(\alpha)}\\\
&\leq||u||_{W^{1,\infty}}\underbrace{\sum_{|\alpha|\leq
1}(\Lambda_{I}^{(\alpha)}+\Lambda_{B}^{(\alpha)})}_{\displaystyle\Lambda_{m}}\end{split}$
(28)

where $\Lambda_{I}^{(\alpha)}$ and $\Lambda_{B}^{(\alpha)}$ are now defined by
equations (23)-(24) with $\psi_{i}^{(\alpha)}$ instead of $\psi_{i}$. From
equation (28) we see that $\Lambda_{m}$ can be used as an upper bound for an
operator norm of the RBF interpolation operator when it is defined on the
Sobolev space $W^{1,\infty}$ [37]. If $u\in W^{1,\infty}$ and equation (12)
holds, we can derive another estimate for the interpolation error which
follows from equation (28):

$\begin{split}||u-u^{h}||_{W^{1,\infty}}&\leq||u-u^{*}||_{W^{1,\infty}}+||u^{*}-u^{h}||_{W^{1,\infty}}\\\
&=||u-u^{*}||_{W^{1,\infty}}+||(u-u^{*})^{h}||_{W^{1,\infty}}\\\
&\leq(1+\Lambda_{m})||u-u^{*}||_{W^{1,\infty}}\end{split}$ (29)

In conclusion, we remark that any $\psi_{i}^{(\alpha)}(\boldsymbol{x})$ is
obtained as the $i^{th}$ component of the solution of a linear system similar
to the one in equation (20), where the the left hand side is kept the same and
$\\{\varphi_{i}(\boldsymbol{x})\\}_{1}^{m}$ is changed with
$\\{\varphi_{i}^{(\alpha)}(\boldsymbol{x})\\}_{1}^{m}$ and
$\\{p_{j}(\boldsymbol{x})\\}_{1}^{q}$ with
$\\{p_{j}^{(\alpha)}(\boldsymbol{x})\\}_{1}^{q}$. We now see that, whenever
the matrix $\mathbf{M}^{T}$ is ill-conditioned, the terms
$\Lambda_{I}^{(\alpha)}$ and $\Lambda_{B}^{(\alpha)}$ are unbounded.

##  3 Neumann Boundary Conditions

###  3.1 Preliminary considerations

Figure 1: reference stencil with $m_{I}=15$ internal nodes (filled dots) and
$m_{B}=7$ boundary nodes (empty circles).

Figure 2: condition number of the interpolation matrix (left), and Lebesgue
constant $\Lambda_{I}$ (right), for the RBF interpolation on the reference
stencil of Figure 1.

$\alpha=-\displaystyle\frac{\pi}{4}-0.02$

$\alpha=\displaystyle\frac{\pi}{4}$

Figure 3: 3D visualization of the Lebesgue functions
$\lambda_{I}(\boldsymbol{x})$ and $\lambda_{B}(\boldsymbol{x})$ for the MQ RBF
interpolation ($\varepsilon s=0.5$ and no polynomial augmentation) on the
stencil of Figure 1. Note the large difference in the $z$-axis scaling for the
two values of $\alpha$.

Figure 1 represents a stand-alone simplified model of the typical situation we
encounter when the interpolation system of equation (17) is solved for a
stencil close to the boundary. Nodes are arranged in a perfectly hexagonal
layout and each unit normal is arranged along the line joining the
corresponding boundary node with point $G$, and pointing outward. The
direction of the normals can thus be controlled by moving the point $G$ along
the axis of symmetry, this way an angle $\alpha\in[-\pi/2\,,\,\pi/2]$ remains
defined by the normal located further to the right. $\alpha>0$ if $G$ is above
the line of the boundary nodes, $\alpha<0$ otherwise.

In Figure 2 we report the condition number $\kappa(\mathbf{M})$ of the
interpolation matrix and the constant $\Lambda_{I}$, defined in equation (24),
against the angle $\alpha$ in the case of the reference stencil of Figure 1.
Whenever one of these quantities grows unbounded, we have a very ill-
conditioned system in equation (17), which translates into large errors and
most likely stability issues. All plots of Figures 2 and 3 are attained using
MQ RBF with shape parameter $\varepsilon$ satisfying $\varepsilon s=0.5$ (see
Table 1), where $s$ is the distance between the nodes.

It is possible to see that the negative values of $\alpha$ are the most
dangerous, they are associated with locations of the point $G$ placed below
the boundary in Figure 1, i.e., with converging normals. From Figure 2 we see
that multiple singularities appear for certain negative values of $\alpha$,
near which both $\kappa(\mathbf{M})$ and $\Lambda_{I}$ assume very large
values. Furthermore such singular configurations exist for modest values of
$\alpha$ regardless of the degree of the polynomial augmentation and this
remains true even with larger stencils and polynomials of higher degrees.
These remarks seem to tell us that no matter the node density or the
polynomial degree, even a slight local concavity of the boundary might cause
severe ill-conditioning issues for the considered RBF interpolations.

In Figure 3 it is possible to qualitatively appreciate the functions
$\lambda_{I}(\boldsymbol{x})$ and $\lambda_{B}(\boldsymbol{x})$ for two
particular values of the angle $\alpha$, $\lambda_{I}$ and $\lambda_{B}$ being
defined in equation (23). These plots are again attained using the same MQ RBF
scheme with no polynomial augmentation, while the inclusion of the polynomial
augmentation would have increased even more the values of $\lambda_{I}$ and
$\lambda_{B}$ near the edges of the stencil.

Once again, we can see that for certain negative values of $\alpha$, i.e., for
critical directions of the normals, large interpolation errors can occur. We
can therefore see that Neumann boundary conditions might compromise the
accuracy and the well-posedness of the interpolation scheme of equation (17)
and, through this, compromise the final solution of system (18). In the
following sections we aim at providing a more precise description for this
phenomenon, along with some strategies for avoiding singular configurations of
the normals.

###  3.2 Bare RBF with one boundary node

####  3.2.1 Singular direction

Consider for simplicity a stencil with $m_{B}=1$ boundary node
$\boldsymbol{x}_{m}$, i.e., $m_{I}=m-1$, and an interpolant with no polynomial
augmentation as in equation (1). The $m\times m$ interpolation matrix is
therefore:

$\mathbf{M}=\boldsymbol{\varphi}_{BC}=\begin{bmatrix}\varphi_{1}(\boldsymbol{x}_{1})\\!\\!\\!&\cdots\\!\\!\\!&\varphi_{m}(\boldsymbol{x}_{1})\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
\varphi_{1}(\boldsymbol{x}_{m_{I}})\\!\\!\\!&\cdots\\!\\!\\!&\varphi_{m}(\boldsymbol{x}_{m_{I}})\\\
\partial\varphi_{1}(\boldsymbol{x}_{m})\\!\\!\\!&\cdots\\!\\!\\!&\partial\varphi_{m}(\boldsymbol{x}_{m})\end{bmatrix}$
(30)

The determinant of $\mathbf{M}$ can be expressed through a cofactor expansion
along the last row:

$\det(\mathbf{M})=\sum_{j=1}^{m_{I}}C_{m,j}\partial\varphi_{j}(\boldsymbol{x}_{m})$
(31)

In equation (31), $C_{m,j}=(-1)^{m+j}D_{m,j}$ is the cofactor of entry
$M_{m,j}=\partial\varphi_{j}(\boldsymbol{x}_{m})$ in $\mathbf{M}$, where
$D_{m,j}$ is the minor of the corresponding entry, i.e., the determinant of
the submatrix obtained by deleting row $m$ and column $j$ from $\mathbf{M}$.
We point out that the expansion in equation (31) has $m_{I}=m-1$ terms because
the diagonal entry of the last row is
$\partial\varphi_{m}(\boldsymbol{x}_{m})=0$, as outlined in section 2.2.

Recalling the form of the normal derivative and the notation of equation (11),
the determinant can be expressed as:

$\det(\mathbf{M})=\sum_{j=1}^{m_{I}}C_{m,j}\Phi^{\prime}(r_{j,m})\,\boldsymbol{e}_{j,m}\cdot\boldsymbol{n}=\boldsymbol{v}\cdot\boldsymbol{n}$
(32)

where $\boldsymbol{n}=\boldsymbol{n}_{m}$ is the unit normal at the boundary
node $\boldsymbol{x}_{m}$. Thus $\boldsymbol{v}$ turns out to be a linear
combination of $\boldsymbol{e}_{j,m}$, i.e., the unit directions pointing
towards the boundary node $\boldsymbol{x}_{m}$, as defined in section 2.2.
Denoting with $w_{j,m}=C_{m,j}\Phi^{\prime}(r_{j,m})$ the coefficient of such
linear combination, it takes the form:

$\boldsymbol{v}=\sum_{j=1}^{m_{I}}w_{j,m}\boldsymbol{e}_{j,m}$ (33)

When the location of the boundary node $\boldsymbol{x}_{m}$ is fixed and the
unit normal $\boldsymbol{n}$ can vary, the direction of $\boldsymbol{v}$ can
be interpreted as the optimal direction for $\boldsymbol{n}$ which maximizes
$\det(\mathbf{M})$ (or minimizes it in the opposite direction). On the other
hand, equation (32) states that if $\boldsymbol{n}$ is orthogonal to
$\boldsymbol{v}$, then $\mathbf{M}$ becomes singular.

The existence of such singular directions can also be directly deduced from
the following fact. If we denote the interpolation matrix in equation (30) by
$\mathbf{M}(\boldsymbol{n})$, i.e., as a function of the unit normal
$\boldsymbol{n}$, then $\det(\mathbf{M}(\boldsymbol{n}))$ is a continuous
function of $\boldsymbol{n}$ and we observe that
$\det(\mathbf{M}(-\boldsymbol{n}))=-\det(\mathbf{M}(\boldsymbol{n}))$ for any
$\boldsymbol{n}$, since the reversal of $\boldsymbol{n}$ results in the change
of sign of the last row of $\mathbf{M}(\boldsymbol{n})$. Since
$\boldsymbol{n}$ can be changed continuously to reach $-\boldsymbol{n}$, there
exist at least one normal for which the determinant vanishes.

####  3.2.2 Possible remedies

The aforementioned negative result can be addressed by using different
strategies, as outlined in section 1. In this paper we will focus on the
following two approaches:

Approach 1. |  Discard a boundary node if the corresponding dot product $|\boldsymbol{v}\cdot\boldsymbol{n}|/\|\boldsymbol{v}\|_{2}$ is too small.
---|---
Approach 2. |  Change the boundary node location according to the corresponding normal.

Both approaches are inspired by the negative result on scattered data
interpolation in multiple dimensions expressed by the Curtis-Mairhuber Theorem
[35]. In the case of interpolation with Neumann BCs this translates into the
following: the basis functions should depend not only on the data locations,
but also on the normals.

Approach 2 will be discussed in section 4.2. Approach 1, on the other hand,
requires either a calculation of $\boldsymbol{v}$ or a suitable approximation
of it. Given the definition of $\boldsymbol{v}$ in equation (33), the easiest
approximation is to consider the direction pointing from the boundary node
towards some reference point, e.g., the stencil centroid. The more constant
the coefficients $w_{j,m}$ in equation (33), the more accurate this
approximation would be. Unfortunately, this condition never occurs, as shown
in Figure 4. In this figure the magnitude and the sign of the coefficients
$w_{j,m}$ are shown in the case of the MQ RBF for a boundary stencil, i.e.,
one-sided stencil close to a straight boundary, with hexagonal node
arrangement and for two different positions of the boundary node, shown as a
red cross. Three different values for the shape parameter $\varepsilon$ are
considered, where the distances are scaled with respect to the constant nodal
spacing $s$, i.e., $\varepsilon s=constant$ (the case $\varepsilon s=0.1$ has
been computed with MATLAB variable precision arithmetic due to the extreme
cancellations occuring for small $\varepsilon$ [38]).

$\varepsilon s=0.1$

$\varepsilon s=1$

$\varepsilon s=10$

Figure 4: coefficients $w_{j,m}$ in equation (32) for the optimal direction,
shown as a blue arrow, associated to the boundary node, shown as a red cross.
Positive and negative coefficients are shown as filled dots and empty circles,
respectively, while the size of the markers is proportional to the magnitude
of the coefficients. Top row: symmetric position of the boundary node; bottom
row: asymmetric position of the boundary node.

Figure 4 shows that moving away from the boundary node, the coefficients
$w_{j,m}$ have alternating signs and a decay rate growing with $\varepsilon$,
which is in perfect accordance with theoretical observations for RBF
interpolation on equispaced infinite lattices [39]. More specifically, by
qualitative arguments based on the adaptation to a finite stencil, the
cardinal expansion coefficients $\lambda_{k}$ defined in [39] correspond to
the cofactors $C_{m,j}$ in equation (31) up to a multiplicative constant,
where the subscript $k$ is the nondimensional distance from the boundary node.

The qualitative behaviour of the magnitude of the coefficients
$w_{j,m}=C_{m,j}\Phi^{\prime}(r_{j,m})$ is therefore given by the following
equation, where the second proportionality can be found in [39] for the
exponential regime:

$|w_{j,m}|\propto|\lambda_{k}\Phi^{\prime}(ks)|\propto e^{-\mu
k}|\Phi^{\prime}(ks)|$ (34)

In the case of MQ RBF, equation (34) can be expressed in the form:

$|w_{j,m}|\propto\frac{ke^{-\mu k}}{\sqrt{1+(\varepsilon
sk)^{2}}}=\tilde{w}(k)$ (35)

When $\varepsilon s=0.1$, then $\mu=0.15$ [39] and $\tilde{w}(k)$ has a
maximum for $k\approx 5$, which in Figure 4 corresponds to the large magnitude
coefficients $w_{j,m}$ far from the boundary node. When $\varepsilon s=1$,
then $\mu=1.04$ [39] and $\tilde{w}(k)$ has a maximum for $k\approx 1$, which
in Figure 4 corresponds to the large magnitude coefficients $w_{j,m}$
immediately close to the boundary node, followed by a rapid decay. When
$\varepsilon s=10$ the coefficients with largest magnitude are again
associated to the nodes immediately close to the boundary node, followed by a
stronger decay.

Figure 4 also shows the optimal direction, i.e., the direction of
$\boldsymbol{v}$ defined in equation (33), as a blue vector. In the case of
the symmetric position of the boundary node (top row in Figure 4), the optimal
direction is vertical, regardless of $\varepsilon s$, as expected for
symmetry. In the case of the asymmetric position of the boundary node (bottom
row in Figure 4), the optimal direction changes from almost vertical for
$\varepsilon s=0.1$ to a certain inclination for $\varepsilon s=10$, in
accordance to the previous observations. Indeed, the larger the shape
parameter $\varepsilon$, the larger the coefficients $w_{j,m}$ associated to
the closest nodes and therefore the more the optimal direction points towards
these closest nodes.

Figure 5: optimal vectors $\boldsymbol{v}$ (solid black lines) and their
envelopes (solid red curves) for two stencils with $m_{I}=3$ (left) and
$m_{I}=15$ (right) internal nodes as the boundary node moves along the blue
dashed curve. This last curve represents the points having constant distance
$s$, i.e., the nodal spacing, from the nearest internal node.

Figure 5 shows the direction and the magnitude of the optimal vector
$\boldsymbol{v}$, depicted with solid black lines, for two hexagonal stencils,
with $m_{I}=3$ (left) and $m_{I}=15$ (right) internal nodes, as the boundary
node moves along the blue dashed curve. This blue dashed curve is formed by
those points that have constant distance $s$ from the nearest internal node.
MQ RBF with $\varepsilon s=0.5$ has been employed. The envelopes of the
optimal directions of $\boldsymbol{v}$ are shown as solid red curves. In the
case $m_{I}=3$ the optimal vectors seem to point qualitatively toward some
reference mid point related to the internal nodes, although their envelope
shows there is no such single point, or region, which can be heuristically
assumed as a good reference point: the envelope indeed falls far out of the
region delimited by the three internal nodes. We also note that the magnitude
of $\boldsymbol{v}$, i.e., the length of the solid black lines, is slightly
larger when the boundary node is located further to the right of the blue
dashed curve. In such cases the optimal direction is pointing toward the
closest internal node on the right and the envelope for these locations gets
indeed closer to this node. This means that, given the role of
$\boldsymbol{v}$ in equation (32), these are the locations for the boundary
nodes on the blue dashed curve that maximize $\det(\mathbf{M})$ when
$\boldsymbol{n}=\boldsymbol{v}$.

In the case $m_{I}=3$, although a first look might suggests that all optimal
directions converge at the same point, the fact that the envelope falls far
out of the region delimited by the three internal nodes refutes this initial
impression. Therefore, it is not possible to naively approximate the optimal
directions as vectors radiating from a precise region.

The length of the solid black lines in Figure 5 are scaled according to the
value of $||\boldsymbol{v}||_{2}$. We see that whenever the boundary node is
placed at the right margin, then the optimal direction, as well as
$||\boldsymbol{v}||_{2}$, suggest that a better conditioning is achieved when
the normal is parallel to the vector joining the boundary node with its
closest internal neighbor.

The case $m_{I}=15$ is very similar to the previous one, with an envelope
falling further away from the region delimited by the internal nodes, thus
confirming that there is no such reasonable mid point for a simple geometric
approximation of the optimal directions. This holds even in the case of very
regular stencils with hexagonal arrangement like those employed here. The
length of the solid black lines confirms that also in this case it is
convenient to have a boundary normal pointing toward the nearest internal node
at the right margin. By heuristic arguments, the enforcement of a Neumann BC
at a boundary node should be “supported” by a close internal node along the
direction of the normal.

We conclude that an estimate for the optimal direction $\boldsymbol{v}$ can
not be obtained through simple geometric heuristics and a more robust
calculation is therefore mandatory in the general case. Such calculation
follows here below.

####  3.2.3 Practical calculation of $\boldsymbol{v}$

The formulation given in equation (33) for $\boldsymbol{v}$ is helpful to
highlight the previous observations but is not of practical interest. A more
convenient approach can be derived by considering the following splitting of
the interpolation matrix $\mathbf{M}$ given in equation (30):

$\mathbf{M}=\begin{bmatrix}\begin{array}[]{ccc|c}\varphi_{1}(\boldsymbol{x}_{1})\\!\\!\\!&\cdots\\!\\!\\!&\varphi_{m_{I}}(\boldsymbol{x}_{1})&\varphi_{m}(\boldsymbol{x}_{1})\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots&\vdots\\\
\varphi_{1}(\boldsymbol{x}_{m_{I}})\\!\\!\\!&\cdots\\!\\!\\!&\varphi_{m_{I}}(\boldsymbol{x}_{m_{I}})&\varphi_{m}(\boldsymbol{x}_{m_{I}})\\\
\hline\cr\partial\varphi_{1}(\boldsymbol{x}_{m})\\!\\!\\!&\cdots\\!\\!\\!&\partial\varphi_{m_{I}}(\boldsymbol{x}_{m})&0\end{array}\end{bmatrix}=\begin{bmatrix}\begin{array}[]{c|c}\boldsymbol{\varphi}_{II}&\boldsymbol{\varphi}_{IB}\\\
\hline\cr\partial\boldsymbol{\varphi}_{BI}&0\end{array}\end{bmatrix}$ (36)

where the bottom-right diagonal entry is
$\partial\varphi_{m}(\boldsymbol{x}_{m})=0$. $\boldsymbol{\varphi}_{II}$ is
the symmetric $m_{I}\times m_{I}$ RBF interpolation matrix obtained when
considering the $m_{I}$ internal nodes without the boundary node
$\boldsymbol{x}_{m}$, $\boldsymbol{\varphi}_{IB}$ is the column vector of the
values of function $\varphi_{m}$, i.e., the RBF centered at
$\boldsymbol{x}_{m}$, evaluated at the $m_{I}$ internal nodes, while
$\partial\boldsymbol{\varphi}_{BI}$ is the row vector of the normal
derivatives of the RBFs centered at the $m_{I}$ internal nodes and evaluated
at $\boldsymbol{x}_{m}$.

By taking the Schur complement $S_{BB}$ of the block
$\boldsymbol{\varphi}_{II}$ in equation (36), we obtain the following scalar
Schur complement:

$S_{BB}=-\partial\boldsymbol{\varphi}_{BI}\boldsymbol{\varphi}_{II}^{-1}\boldsymbol{\varphi}_{IB}=-\boldsymbol{\varphi}_{IB}^{T}\boldsymbol{\varphi}_{II}^{-1}\partial\boldsymbol{\varphi}_{BI}^{T}$
(37)

where the second equality derives from a simple transposition since
$\boldsymbol{\varphi}_{II}$ is symmetric. The following equality also holds:

$\det(\mathbf{M})=\det(\boldsymbol{\varphi}_{II})S_{BB}$ (38)

where $\det(\boldsymbol{\varphi}_{II})\neq 0$ when using strictly positive
definite or conditionally positive definite RBFs $\Phi(r)$ of order 1.

Because of the symmetry of the RBFs, i.e.,
$\varphi_{m}(\boldsymbol{x}_{j})=\varphi_{j}(\boldsymbol{x}_{m})$, the column
vector $\boldsymbol{\varphi}_{IB}$ can also be interpreted as the values of
the RBFs centered at the $m_{I}$ internal nodes, evaluated at the boundary
node $\boldsymbol{x}_{m}$. Therefore, the column vector
$\boldsymbol{\varphi}_{II}^{-1}\boldsymbol{\varphi}_{IB}$ in the first matrix
expression in equation (37) can be interpreted as the values of the cardinal
functions $\bar{\psi}_{j}$ of the RBF interpolation based on the internal
nodes only, evaluated at $\boldsymbol{x}_{m}$:

$\bar{\boldsymbol{\psi}}=\boldsymbol{\varphi}_{II}^{-1}\boldsymbol{\varphi}_{IB}$
(39)

where
$\bar{\boldsymbol{\psi}}=\\{\bar{\psi}_{j}(\boldsymbol{x}_{m})\\}_{j=1}^{m_{I}}$.

In the end, by comparing equation (31) with equation (38) and by using the
first matrix expression in equation (37), we can see that, up to the
multiplicative constant $-\det(\boldsymbol{\varphi}_{II})$, the cofactors
$C_{m,j}$ in equation (31) correspond to the aforementioned cardinal functions
$\bar{\psi}_{j}$:

$C_{m,j}=-\det(\boldsymbol{\varphi}_{II})\bar{\psi}_{j}(\boldsymbol{x}_{m})\,,\quad
j=1,\dots,m_{I}$ (40)

Therefore, $\boldsymbol{v}$ can be computed from equation (33) recalling that
$w_{j,m}=C_{m,j}\Phi^{\prime}(r_{j,m})$:

$\boldsymbol{v}=\sum_{j=1}^{m_{I}}C_{m,j}\Phi^{\prime}(r_{j,m})\boldsymbol{e}_{j,m}=-\det(\boldsymbol{\varphi}_{II})\sum_{j=1}^{m_{I}}\bar{\psi}_{j}(\boldsymbol{x}_{m})\Phi^{\prime}(r_{j,m})\boldsymbol{e}_{j,m}$
(41)

where $\det(\boldsymbol{\varphi}_{II})$ does not need to be evaluated since it
is a multiplicative constant and the values
$\bar{\psi}_{j}(\boldsymbol{x}_{m})$ can be computed from equation (39),
intended as the solution of the linear system
$\boldsymbol{\varphi}_{II}\bar{\boldsymbol{\psi}}=\boldsymbol{\varphi}_{IB}$.

####  3.2.4 Qualitative interpretation of $S_{BB}$

The second matrix expression in equation (37), i.e.,
$S_{BB}=-\boldsymbol{\varphi}_{IB}^{T}\boldsymbol{\varphi}_{II}^{-1}\partial\boldsymbol{\varphi}_{BI}^{T}$,
provides the following explicit interpretation of $S_{BB}$.

The entries of $\partial\boldsymbol{\varphi}_{BI}$ are antisymmetric, i.e.,
$\partial\varphi_{j}(\boldsymbol{x}_{m})=-\partial\varphi_{m}(\boldsymbol{x}_{j})$,
from equation (11), since we always refer to the same normal $\boldsymbol{n}$.
Therefore, the vector $\partial\boldsymbol{\varphi}_{BI}$ can be interpreted
as the values of the normal derivative of the function $-\varphi_{m}$
evaluated at the $m_{I}$ internal nodes. The matrix-vector operation
$\boldsymbol{\varphi}_{II}^{-1}\partial\boldsymbol{\varphi}_{BI}^{T}$ hence
gives the column vector of the RBF expansion coefficients
$\\{\alpha_{i}\\}_{1}^{m_{I}}$ of the RBF interpolation, based on the $m_{I}$
internal nodes only, of the function $-\partial\varphi_{m}$. The product
$-\boldsymbol{\varphi}_{IB}^{T}(\boldsymbol{\varphi}_{II}^{-1}\partial\boldsymbol{\varphi}_{BI}^{T})$
then gives the extrapolation of the function $\partial\varphi_{m}$ at
$\boldsymbol{x}_{m}$.

Since $\partial\varphi_{m}(\boldsymbol{x}_{m})=0$, $S_{BB}$ can now be
interpreted as the difference between the value of the function
$\partial\varphi_{m}$ at $\boldsymbol{x}_{m}$, i.e., 0, and the extrapolation
of the same function at that same location $\boldsymbol{x}_{m}$ when using a
RBF scheme based on the internal nodes only. This interpretation, which will
continue to hold in the case of multiple boundary nodes presented in section
3.3.1, tells us that the interpolation matrix $\mathbf{M}$ in equation (30)
will be ill-conditioned if the internal nodes are sufficient for a good RBF
extrapolation of the function $\partial\varphi_{m}$ at the location of the
boundary node. Indeed, if the RBF interpolation based on the internal nodes is
not capable of that, the additional information provided by the Neumann BC,
i.e., the last row of $\mathbf{M}$, is not redundant, thus leading to a well-
conditioned interpolation matrix.

####  3.2.5 Additional remarks

$m_{I}=3$

$m_{I}=5$

$m_{I}=12$

Figure 6: norm of $\boldsymbol{v}$ for different node arrangements as the
position of the boundary node changes. The $m_{I}$ internal nodes are shown as
black empty circles and the red dots represent the positions of the boundary
node for which $\boldsymbol{v}=\boldsymbol{0}$. Red dots represent therefore
forbidden positions for the boundary node. Top row: hexagonal arrangements;
bottom row: hexagonal arrangements perturbed by random displacements.
$\|\boldsymbol{v}_{max}\|$ represents the maximum norm of $\boldsymbol{v}$ for
each subfigure.

A final remark before we conclude the discussion about the case with a single
boundary node concerns the following questions. Given a node arrangement for
the $m_{I}$ internal nodes, are there particular locations of the boundary
node for which $\det(\mathbf{M})=0$ regardless of the direction of the normal
$\boldsymbol{n}$, i.e., for which $\boldsymbol{v}=\boldsymbol{0}$ in equation
(32) or, equivalently, $S_{BB}=0$ in equation (38), provided that
$\det(\boldsymbol{\varphi}_{II})\neq 0$? If so, is it possible that a Neumann
boundary condition ends up being defined in one of such singular locations?

Fortunately, this seems not to be the case: Figure 6 shows these singular
locations as red dots for different node arrangements in the case of MQ RBF
with $\varepsilon s=0.5$, together with a representation of
$\|\boldsymbol{v}\|$. From this figure it can be observed that all singular
locations fall “inside” the stencil and dangerous locations, i.e., locations
of the boundary node for which $\|\boldsymbol{v}\|$ is close to 0, are also in
the neighbourhood of the stencil. In other words, the further the boundary
node is from the internal nodes, the further $\mathbf{M}$ is from being
singular, if the employed normal is not orthogonal to the optimal direction.
Whenever we have a boundary node “outside” the region delimited by the
internal nodes, we do not need therefore to worry about the singular
locations. This principle has been verified numerically for each type of RBF
presented in Table 1 which is at least conditionally positive definite of
order 1, i.e., RBFs that can be used for RBF interpolation without the need of
any polynomial augmentation. For instance, the use of $\Phi(r)=r^{2}\log r$,
which is conditionally positive definite of order 2 and therefore would
require an augmentation with a linear polynomial, leads to the existence of
singular locations far away from internal nodes.

###  3.3 Bare RBF with multiple boundary nodes

####  3.3.1 Preliminary results

Consider now a stencil with $m_{B}$ boundary nodes and, again, an interpolant
with no polynomial augmentation. The $m\times m$ interpolation matrix is
therefore:

$\mathbf{M}=\boldsymbol{\varphi}_{BC}=\begin{bmatrix}\varphi_{1}(\boldsymbol{x}_{1})\\!\\!\\!&\cdots\\!\\!\\!&\varphi_{m}(\boldsymbol{x}_{1})\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
\varphi_{1}(\boldsymbol{x}_{m_{I}})\\!\\!\\!&\cdots\\!\\!\\!&\varphi_{m}(\boldsymbol{x}_{m_{I}})\\\
\partial\varphi_{1}(\boldsymbol{x}_{m_{J}})\\!\\!\\!&\cdots\\!\\!\\!&\partial\varphi_{m}(\boldsymbol{x}_{m_{J}})\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
\partial\varphi_{1}(\boldsymbol{x}_{m})\\!\\!\\!&\cdots\\!\\!\\!&\partial\varphi_{m}(\boldsymbol{x}_{m})\end{bmatrix}$
(42)

where $m_{J}=m_{I}+1$. Similarly to what was done in section 3.2.3, consider
the following splitting of $\mathbf{M}$:

$\begin{split}\mathbf{M}=&\begin{bmatrix}\begin{array}[]{ccc|ccc}\varphi_{1}(\boldsymbol{x}_{1})\\!\\!\\!&\cdots\\!\\!\\!&\varphi_{m_{I}}(\boldsymbol{x}_{1})&\varphi_{m_{J}}(\boldsymbol{x}_{1})\\!\\!\\!&\cdots\\!\\!\\!&\varphi_{m}(\boldsymbol{x}_{1})\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots&\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
\varphi_{1}(\boldsymbol{x}_{m_{I}})\\!\\!\\!&\cdots\\!\\!\\!&\varphi_{m_{I}}(\boldsymbol{x}_{m_{I}})&\varphi_{m_{J}}(\boldsymbol{x}_{m_{I}})\\!\\!\\!&\cdots\\!\\!\\!&\varphi_{m}(\boldsymbol{x}_{m_{I}})\\\
\hline\cr\partial\varphi_{1}(\boldsymbol{x}_{m_{J}})\\!\\!\\!&\cdots\\!\\!\\!&\partial\varphi_{m_{I}}(\boldsymbol{x}_{m_{J}})&\partial\varphi_{m_{J}}(\boldsymbol{x}_{m_{J}})\\!\\!\\!&\cdots\\!\\!\\!&\partial\varphi_{m}(\boldsymbol{x}_{m_{J}})\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots&\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
\partial\varphi_{1}(\boldsymbol{x}_{m})\\!\\!\\!&\cdots\\!\\!\\!&\partial\varphi_{m_{I}}(\boldsymbol{x}_{m})&\partial\varphi_{m_{J}}(\boldsymbol{x}_{m})\\!\\!\\!&\cdots\\!\\!\\!&\partial\varphi_{m}(\boldsymbol{x}_{m})\end{array}\end{bmatrix}\\\
=&\begin{bmatrix}\begin{array}[]{c|c}\boldsymbol{\varphi}_{II}&\boldsymbol{\varphi}_{IB}\\\
\hline\cr\partial\boldsymbol{\varphi}_{BI}&\partial\boldsymbol{\varphi}_{BB}\end{array}\end{bmatrix}\end{split}$
(43)

where the $m_{B}$ diagonal entries of $\partial\boldsymbol{\varphi}_{BB}$ are
$\partial\varphi_{i}(\boldsymbol{x}_{i})=0$, $i=m_{J},\ldots,m$, as previously
stated. $\boldsymbol{\varphi}_{II}$ is again the symmetric $m_{I}\times m_{I}$
RBF interpolation matrix obtained when considering the $m_{I}$ internal nodes
only. $\boldsymbol{\varphi}_{IB}$ is the $m_{I}\times m_{B}$ matrix whose
columns are the RBFs centered at the $m_{B}$ boundary nodes and evaluated at
the $m_{I}$ internal nodes, while $\partial\boldsymbol{\varphi}_{BI}$ is the
$m_{B}\times m_{I}$ matrix whose rows are the normal derivatives of the RBFs
centered at the $m_{I}$ internal nodes and evaluated at the $m_{B}$ boundary
nodes.

By taking the Schur complement $\mathbf{S}_{BB}$ of the
$\boldsymbol{\varphi}_{II}$ block of $\mathbf{M}$ in equation (43), we obtain
the following $m_{B}\times m_{B}$ Schur complement matrix:

$\mathbf{S}_{BB}=\partial\boldsymbol{\varphi}_{BB}-\partial\boldsymbol{\varphi}_{BI}\boldsymbol{\varphi}_{II}^{-1}\boldsymbol{\varphi}_{IB}$
(44)

The following equality also holds:

$\det(\mathbf{M})=\det(\boldsymbol{\varphi}_{II})\det(\mathbf{S}_{BB})$ (45)

where, again, $\det(\boldsymbol{\varphi}_{II})\neq 0$ when using strictly
positive definite or conditionally positive definite functions $\Phi(r)$ or
order 1.

####  3.3.2 Qualitative interpretation of $\mathbf{S}_{BB}$

Equation (44) provides again an explicit interpretation of $\mathbf{S}_{BB}$
as follows. For RBFs $\varphi_{k}$ centered at the $m_{B}$ boundary nodes
$\boldsymbol{x}_{k}$ we define:

$\partial^{*}\\!\varphi_{k}(\boldsymbol{x}_{i}):=\frac{\partial\varphi_{k}}{\partial\boldsymbol{n}_{k}}(\boldsymbol{x}_{i})$
(46)

Because of the symmetry of the RBFs, from equation (11) we have:

$\partial\varphi_{i}(\boldsymbol{x}_{k})=\frac{\partial\varphi_{i}}{\partial\boldsymbol{n}_{k}}(\boldsymbol{x}_{k})=-\frac{\partial\varphi_{k}}{\partial\boldsymbol{n}_{k}}(\boldsymbol{x}_{i})=-\partial^{*}\\!\varphi_{k}(\boldsymbol{x}_{i})$
(47)

and then from equation (44) we have:

$\mathbf{S}_{BB}=-\big{(}\partial^{*}\\!\boldsymbol{\varphi}_{BB}-\partial^{*}\\!\boldsymbol{\varphi}_{BI}\boldsymbol{\varphi}_{II}^{-1}\boldsymbol{\varphi}_{IB}\big{)}$
(48)

where $\partial^{*}\\!\boldsymbol{\varphi}_{BB}$ and
$\partial^{*}\\!\boldsymbol{\varphi}_{BI}$ are obtained from
$\partial\boldsymbol{\varphi}_{BB}$ and $\partial\boldsymbol{\varphi}_{BI}$,
respectively, by replacing $\partial$ with $\partial^{*}$ and reversing the
indexes of RBF centers with the indexes of evaluation points, as follows:

$\partial^{*}\\!\boldsymbol{\varphi}_{BB}=\begin{bmatrix}\partial^{*}\\!\varphi_{m_{J}}(\boldsymbol{x}_{m_{J}})\\!\\!\\!&\cdots\\!\\!\\!&\partial^{*}\\!\varphi_{m_{J}}(\boldsymbol{x}_{m})\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
\partial^{*}\\!\varphi_{m}(\boldsymbol{x}_{m_{J}})\\!\\!\\!&\cdots\\!\\!\\!&\partial^{*}\\!\varphi_{m}(\boldsymbol{x}_{m})\\\
\end{bmatrix}$ (49)

$\partial^{*}\\!\boldsymbol{\varphi}_{BI}=\begin{bmatrix}\partial^{*}\\!\varphi_{m_{J}}(\boldsymbol{x}_{1})\\!\\!\\!&\cdots\\!\\!\\!&\partial^{*}\\!\varphi_{m_{J}}(\boldsymbol{x}_{m_{I}})\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
\partial^{*}\\!\varphi_{m}(\boldsymbol{x}_{1})\\!\\!\\!&\cdots\\!\\!\\!&\partial^{*}\\!\varphi_{m}(\boldsymbol{x}_{m_{I}})\end{bmatrix}$
(50)

By following the same arguments presented in section 3.2.3, we have that
column $j$ of $\boldsymbol{\varphi}_{II}^{-1}\boldsymbol{\varphi}_{IB}$ in
equation (48) contains the values of the $m_{I}$ cardinal functions of the RBF
interpolation, based on the $m_{I}$ internal nodes only, evaluated at the
$j^{th}$ boundary node. Row $i$ of $\partial^{*}\\!\boldsymbol{\varphi}_{BI}$
in equation (50), on the other hand, contains the values of the function
$\partial^{*}\\!\varphi_{m_{I}+i}$, thus centered at the $i^{th}$ boundary
node, and evaluated at the $m_{I}$ internal nodes. Therefore, entry $(i,j)$ of
matrix
$\partial^{*}\\!\boldsymbol{\varphi}_{BI}\boldsymbol{\varphi}_{II}^{-1}\boldsymbol{\varphi}_{IB}$
in equation (48) is the extrapolation at the $j^{th}$ boundary node of the
function $\partial^{*}\\!\varphi_{m_{I}+i}$ based on its values at the
internal nodes only. On the other hand, entry $(i,j)$ of matrix
$\partial^{*}\\!\boldsymbol{\varphi}_{BB}$ in equation (49) is the actual
value of $\partial^{*}\\!\varphi_{m_{I}+i}$ at the $j^{th}$ boundary node.
Given the definition of $\partial^{*}$, we conclude that entry $(i,j)$ of
$\mathbf{S}_{BB}$ in equation (48) can be interpreted as the difference
between two values: the value of the normal derivative, intended with respect
to the normal at the $i^{th}$ boundary node, of the RBF centered at the same
$i^{th}$ boundary node and the extrapolated value of the same normal
derivative, obtained through a RBF interpolation based on the $m_{I}$ internal
nodes only, both evaluated at the $j^{th}$ boundary node. Similarly to the
case with one boundary node, the previous interpretation tells us that the
interpolation matrix $\mathbf{M}$ in equation (42) will be qualitatively well-
conditioned if the additional information provided by the Neumann BCs, i.e.,
the last $m_{B}$ rows in $\mathbf{M}$, are not redundant with the ability of
the RBF interpolation (based on the internal nodes only) to provide a good
extrapolation at the boundary nodes of the normal derivatives of the RBFs
centered at the same boundary nodes.

####  3.3.3 Dependance upon the normals

Equation (44) can be written as follows:

$\mathbf{S}_{BB}=\partial\boldsymbol{\varphi}_{BB}-\partial\boldsymbol{\varphi}_{BI}\bar{\boldsymbol{\psi}}$
(51)

where
$\bar{\boldsymbol{\psi}}=\boldsymbol{\varphi}_{II}^{-1}\boldsymbol{\varphi}_{IB}$,
already defined in equation (39), is now a $m_{I}\times m_{B}$ matrix in the
case of $m_{B}$ boundary nodes.

Despite the matrices in equation (51) could be handled by standard linear
algebra operations, it is anyhow convenient to use the notation introduced in
section A for simplicity as follows. By denoting with $\boldsymbol{d}_{i,k}$
the following term appearing in equation (11):

$\boldsymbol{d}_{i,k}=\Phi^{\prime}(r_{i,k})\,\boldsymbol{e}_{i,k}$ (52)

the normal derivatives in $\partial\boldsymbol{\varphi}_{BB}$ and in
$\partial\boldsymbol{\varphi}_{BI}$ of equation (43) can be written as simple
dot products:

$\partial\varphi_{i}(\boldsymbol{x}_{k})=\boldsymbol{d}_{i,k}\cdot\boldsymbol{n}_{k}$
(53)

and therefore $\partial\boldsymbol{\varphi}_{BB}$ and
$\partial\boldsymbol{\varphi}_{BI}$ can be written using the operator $H$, as
defined in equation (80), as follows:

$\begin{split}\partial\boldsymbol{\varphi}_{BB}&=H(\mathcal{D}_{BB},\mathcal{N})\\\
\partial\boldsymbol{\varphi}_{BI}&=H(\mathcal{D}_{BI},\mathcal{N})\end{split}$
(54)

In the previous equations $\mathcal{D}_{BB}\in\mathbb{R}^{m_{B}\times
m_{B},d}$ and $\mathcal{D}_{BI}\in\mathbb{R}^{m_{B}\times m_{I},d}$ are
$d$-matrices storing purely geometric information:

$\mathcal{D}_{BB}=\begin{bmatrix}\boldsymbol{d}_{m_{J},m_{J}}\\!\\!\\!&\cdots\\!\\!\\!&\boldsymbol{d}_{m,m_{J}}\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
\boldsymbol{d}_{m_{J},m}\\!\\!\\!&\cdots\\!\\!\\!&\boldsymbol{d}_{m,m}\\\
\end{bmatrix}\,,\quad\mathcal{D}_{BI}=\begin{bmatrix}\boldsymbol{d}_{1,m_{J}}\\!\\!\\!&\cdots\\!\\!\\!&\boldsymbol{d}_{m_{I},m_{J}}\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
\boldsymbol{d}_{1,m}\\!\\!\\!&\cdots\\!\\!\\!&\boldsymbol{d}_{m_{I},m}\\\
\end{bmatrix}$ (55)

and
$\mathcal{N}=\\{\boldsymbol{n}_{m_{J}},\ldots,\boldsymbol{n}_{m}\\}\in\mathbb{R}^{m_{B},d}$
is the $d$-matrix storing the unit normals at the $m_{B}$ boundary nodes. By
considering the vector entries of $\mathcal{D}_{BB}$, we see that
$\boldsymbol{d}_{i,k}=-\boldsymbol{d}_{k,i}$ for $k\neq i$ because of equation
(52), while the diagonal entries are null vectors because of equation (12),
i.e., $\mathcal{D}_{BB}$ is an antisymmetric $d$-matrix.

By using the previous notations and properties (81)-(82), equation (51)
becomes:

$\mathbf{S}_{BB}=H(\mathcal{D}_{BB},\mathcal{N})-H(\mathcal{D}_{BI},\mathcal{N})\bar{\boldsymbol{\psi}}=H(\underbrace{\mathcal{D}_{BB}-\mathcal{D}_{BI}\bar{\boldsymbol{\psi}}}_{\displaystyle\mathcal{G}_{BB}},\mathcal{N})$
(56)

where $\mathcal{G}_{BB}=(\boldsymbol{g}_{ij})\in\mathbb{R}^{m_{B}\times
m_{B},d}$. Given the interpretation of $\mathbf{S}_{BB}$ given in section
3.3.2, vector $\boldsymbol{g}_{ij}$ of $\mathcal{G}_{BB}$ represents therefore
the difference between the gradient of the RBF centered at the $i^{th}$
boundary node and the extrapolated value of the same gradient, obtained
through a RBF interpolation based on the $m_{I}$ internal nodes only, both
evaluated at the $j^{th}$ boundary node. When considering increasingly flat
RBFs [38], i.e., $\varepsilon\to 0$ for the RBFs reported in Table 1, the
previous RBF extrapolation gets more and more accurate and therefore
$\boldsymbol{g}_{ij}\to\boldsymbol{0}$ if the internal stencil is unisolvent
for the lowest-order interpolating polynomial. Hence, in the last case, given
the definition of $\mathcal{G}_{BB}$ in equation (56), we have
$\mathcal{D}_{BI}\bar{\boldsymbol{\psi}}\to\mathcal{D}_{BB}$ in the limit
$\varepsilon\to 0$, i.e., $\mathcal{D}_{BI}\bar{\boldsymbol{\psi}}$ tends to
an antisymmetric $d$-matrix as well.

Equation (56) thus accounts for the explicit influence of the boundary normals
on $\mathbf{S}_{BB}$ and, by equation (45), $\mathbf{S}_{BB}$ in turn
influences the conditioning of the whole interpolation matrix.

From equation (56), the determinant of $\mathbf{S}_{BB}$ can be explicitly
written as follows:

$\det(\mathbf{S}_{BB})=\det\\!\big{(}H(\mathcal{G}_{BB},\mathcal{N})\big{)}=\begin{vmatrix}\boldsymbol{g}_{1,1}\cdot\bar{\boldsymbol{n}}_{1}&\cdots&\boldsymbol{g}_{1,m_{B}}\cdot\bar{\boldsymbol{n}}_{1}\\\
\vdots&\ddots&\vdots\\\
\boldsymbol{g}_{m_{B},1}\cdot\bar{\boldsymbol{n}}_{m_{B}}&\cdots&\boldsymbol{g}_{m_{B},m_{B}}\cdot\bar{\boldsymbol{n}}_{m_{B}}\end{vmatrix}$
(57)

where $\bar{\boldsymbol{n}}_{i}=\boldsymbol{n}_{m_{I}+i}$, i.e., the unit
normal at the $i^{th}$ boundary node, is introduced for simplicity of
notation. We note that in the case with $m_{B}=1$ boundary nodes, equation
(57) becomes:

$\det(\mathbf{S}_{BB})=S_{BB}=\boldsymbol{g}_{1,1}\cdot\bar{\boldsymbol{n}}_{1}$
(58)

which corresponds to equation (32) with
$\boldsymbol{v}=\det(\boldsymbol{\varphi}_{II})\boldsymbol{g}_{1,1}$ because
of equation (38).

Finally we note that even in the case with $m_{B}>1$ boundary nodes, the
determinant of the interpolation matrix can still be expressed in the form of
a dot product with respect to a specific normal $\bar{\boldsymbol{n}}_{i}$ as
in equation (32), i.e.,
$\det(\mathbf{M})=\boldsymbol{v}_{i}\cdot\bar{\boldsymbol{n}}_{i}$ where
$\boldsymbol{v}_{i}$ therefore depends upon the remaining $m_{B}-1$ normals
$\bar{\boldsymbol{n}}_{j\neq i}$. This fact can be derived by directly
expanding the determinant in both equations (43) and (57) with respect to the
row corresponding to the $i^{th}$ boundary node.

####  3.3.4 Optimal directions

Given the formulation of equation (57), one may ask which are the $m_{B}$ unit
vectors $\bar{\boldsymbol{n}}_{i}$ that maximize $\det(\mathbf{S}_{BB})$,
i.e., the optimal directions of the $m_{B}$ normals leading to a well-
conditioned interpolation matrix. Such optimal directions can therefore be
obtained by solving the following constrained maximization problem:

$\begin{array}[]{lll}\text{maximize}&&\det\\!\big{(}H(\mathcal{G}_{BB},\mathcal{N})\big{)}\\\\[2.5pt]
\text{subject to}&&\|\bar{\boldsymbol{n}}_{i}\|_{2}^{2}=1,\quad
i=1,\ldots,m_{B}\\\ \end{array}$ (59)

where
$\mathcal{N}=\\{\bar{\boldsymbol{n}}_{1},\ldots,\bar{\boldsymbol{n}}_{m_{B}}\\}\in\mathbb{R}^{m_{B},d}$.

Problem (59) can be solved by the Lagrange multipliers method as presented in
B. Such optimal directions are shown in Figure 7 for the reference stencil of
Figure 1 when using MQ RBF with different values of the shape factor
$\varepsilon$. In order to highlight the effect of the relative position of
the nodes, optimal directions are also shown when the reference stencil is
perturbed with random nodal displacements.

$\varepsilon s=0.1$

$\varepsilon s=1$

$\varepsilon s=10$

Figure 7: optimal directions for the reference stencil of Figure 1 (top row)
and for the same stencil with random nodal displacements (bottom row).

In the case of the reference stencil, the optimal directions are approximately
pointing towards some reference mid point of the stencil, which however
depends upon the shape factor $\varepsilon$, analogously to the case with one
boundary node depicted in Figure 4. When $\varepsilon s=0.1$, the long range
influence of the internal nodes on the optimal normals seems to dominate: as a
result, the optimal directions have a radial pattern originating from the
geometrical center of the internal nodes. For larger shape factors, i.e.,
$\varepsilon s=1$ and especially $\varepsilon s=10$, conversely, the short
range mutual influence of the boundary nodes seems to dominate, leading to a
more flat pattern, again in accordance with the qualitative arguments
presented in section 3.2.2 related to the coefficients $w_{j,m}$ shown in
Figure 4. In the case of the perturbed stencil, however, the optimal
directions are quite unpredictable, as there are boundary nodes close to each
other with quite different optimal directions, confirming again that a simple
geometric approximation of the optimal directions is not feasible.

###  3.4 Influence of the polynomial part

In order to highlight the effect of the polynomial part of the interpolant (2)
on the previous results, let us consider again the interpolation matrix
$\mathbf{M}$ in equation (7). By using the same matrix splitting employed in
equation (43), submatrix $\boldsymbol{\varphi}_{BC}$ initially defined in
equation (8) can be expressed as:

$\boldsymbol{\varphi}_{BC}=\begin{bmatrix}\begin{array}[]{c|c}\boldsymbol{\varphi}_{II}&\boldsymbol{\varphi}_{IB}\\\
\hline\cr\partial\boldsymbol{\varphi}_{BI}&\partial\boldsymbol{\varphi}_{BB}\end{array}\end{bmatrix}$
(60)

while matrices $\mathbf{P}_{BC}$ and $\mathbf{P}^{T}$, defined in equations
(9) and (10), can be expressed as follows:

$\mathbf{P}_{BC}=\begin{bmatrix}\begin{array}[]{ccc}p_{1}(\boldsymbol{x}_{1})\\!\\!\\!&\cdots\\!\\!\\!&p_{q}(\boldsymbol{x}_{1})\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
p_{1}(\boldsymbol{x}_{m_{I}})\\!\\!\\!&\cdots\\!\\!\\!&p_{q}(\boldsymbol{x}_{m_{I}})\\\
\hline\cr\partial
p_{1}(\boldsymbol{x}_{m_{J}})\\!\\!\\!&\cdots\\!\\!\\!&\partial
p_{q}(\boldsymbol{x}_{m_{J}})\\\ \vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
\partial p_{1}(\boldsymbol{x}_{m})\\!\\!\\!&\cdots\\!\\!\\!&\partial
p_{q}(\boldsymbol{x}_{m})\end{array}\end{bmatrix}=\begin{bmatrix}\begin{array}[]{c}\mathbf{P}_{I}\\\
\hline\cr\partial\mathbf{P}_{B}\end{array}\end{bmatrix}$ (61)

$\mathbf{P}^{T}=\begin{bmatrix}\begin{array}[]{ccc|ccc}p_{1}(\boldsymbol{x}_{1})\\!\\!\\!&\cdots\\!\\!\\!&p_{1}(\boldsymbol{x}_{m_{I}})&p_{1}(\boldsymbol{x}_{m_{J}})\\!\\!\\!&\cdots\\!\\!\\!&p_{1}(\boldsymbol{x}_{m})\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots&\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
p_{q}(\boldsymbol{x}_{1})\\!\\!\\!&\cdots\\!\\!\\!&p_{q}(\boldsymbol{x}_{m_{I}})&p_{q}(\boldsymbol{x}_{m_{J}})\\!\\!\\!&\cdots\\!\\!\\!&p_{q}(\boldsymbol{x}_{m})\end{array}\end{bmatrix}=\begin{bmatrix}\begin{array}[]{c|c}\mathbf{P}^{T}_{I}&\mathbf{P}^{T}_{B}\end{array}\end{bmatrix}$
(62)

By rearranging the rows and the columns of $\mathbf{M}$ according to internal
nodes, polynomial terms and boundary nodes, respectively, the interpolation
matrix in block form becomes:

$\mathbf{M}=\begin{bmatrix}\begin{array}[]{cc|c}\boldsymbol{\varphi}_{II}&\mathbf{P}_{I}&\boldsymbol{\varphi}_{IB}\\\\[2.5pt]
\mathbf{P}^{T}_{I}&\boldsymbol{0}&\mathbf{P}^{T}_{B}\\\
\hline\cr\partial\boldsymbol{\varphi}_{BI}&\partial\mathbf{P}_{B}&\partial\boldsymbol{\varphi}_{BB}\end{array}\end{bmatrix}=\begin{bmatrix}\begin{array}[]{c|c}\mathbf{M}_{II}&\mathbf{Q}_{IB}\\\
\hline\cr\partial\mathbf{Q}_{BI}&\partial\boldsymbol{\varphi}_{BB}\end{array}\end{bmatrix}$
(63)

Similarly to the approach employed in section 3.3.1, we consider the Schur
complement $\mathbf{S}_{BB}$ of the $\mathbf{M}_{II}$ block of $\mathbf{M}$:

$\mathbf{S}_{BB}=\partial\boldsymbol{\varphi}_{BB}-\underbrace{\partial\mathbf{Q}_{BI}\mathbf{M}_{II}^{-1}\mathbf{Q}_{IB}}_{\displaystyle\mathbf{W}_{BB}}$
(64)

where $\mathbf{M}_{II}$ is not singular as long as the internal stencil is
unisolvent for the chosen degree of the polynomial augmentation.

By comparing equation (64) with the corresponding equation (44), obtained
without polynomial augmentation, we see that the only difference is in the
double matrix-product term $\mathbf{W}_{BB}$, that in the case of polynomial
augmentation can be written as follows:

$\mathbf{W}_{BB}=\begin{bmatrix}\partial\boldsymbol{\varphi}_{BI}&\partial\mathbf{P}_{B}\end{bmatrix}\begin{bmatrix}\boldsymbol{\varphi}_{II}&\mathbf{P}_{I}\\\\[2.5pt]
\mathbf{P}^{T}_{I}&\boldsymbol{0}\end{bmatrix}^{-1}\begin{bmatrix}\boldsymbol{\varphi}_{IB}\\\\[2.5pt]
\mathbf{P}^{T}_{B}\end{bmatrix}=\partial\boldsymbol{\varphi}_{BI}\mathbf{K}_{1}+\partial\mathbf{P}_{B}\mathbf{K}_{2}$
(65)

where $\mathbf{K}_{1}$ and $\mathbf{K}_{2}$ can be obtained by simple
algebraic operations as follows:

$\begin{split}\mathbf{K}_{1}&=\boldsymbol{\varphi}_{II}^{-1}\big{(}\boldsymbol{\varphi}_{IB}-\mathbf{P}_{I}\mathbf{K}_{2}\big{)}\\\
\mathbf{K}_{2}&=\underbrace{\big{(}\mathbf{P}_{I}^{T}\boldsymbol{\varphi}_{II}^{-1}\mathbf{P}_{I}\big{)}^{-1}}_{\displaystyle\mathbf{S}}\big{(}\boldsymbol{\varphi}_{IB}^{T}\boldsymbol{\varphi}_{II}^{-1}\mathbf{P}_{I}-\mathbf{P}_{B}\big{)}^{T}\end{split}$
(66)

Further algebraic manipulations lead to the following form:

$\mathbf{W}_{BB}=\partial\boldsymbol{\varphi}_{BI}\boldsymbol{\varphi}_{II}^{-1}\boldsymbol{\varphi}_{IB}-\underbrace{\big{(}\partial\boldsymbol{\varphi}_{BI}\boldsymbol{\varphi}_{II}^{-1}\mathbf{P}_{I}-\partial\mathbf{P}_{B}\big{)}\mathbf{K}_{2}}_{\displaystyle\Delta\mathbf{W}}$
(67)

which, compared to the corresponding term
$\partial\boldsymbol{\varphi}_{BI}\boldsymbol{\varphi}_{II}^{-1}\boldsymbol{\varphi}_{IB}$
in equation (44), reveals that the difference in the Schur complement of the
interpolation matrix, with and without polynomial augmentation, is solely due
to the term $\Delta\mathbf{W}$:

$\Delta\mathbf{W}=\underbrace{\big{(}\partial\boldsymbol{\varphi}_{BI}\boldsymbol{\varphi}_{II}^{-1}\mathbf{P}_{I}-\partial\mathbf{P}_{B}\big{)}}_{\displaystyle\partial\mathbf{E}}\mathbf{S}\big{(}\underbrace{\boldsymbol{\varphi}_{IB}^{T}\boldsymbol{\varphi}_{II}^{-1}\mathbf{P}_{I}-\mathbf{P}_{B}}_{\displaystyle\mathbf{E}}\big{)}^{T}$
(68)

Similarly to the qualitative interpretations given in section 3.3.2, we have
that entry $(i,j)$ of matrix $\partial\mathbf{E}$ can be interpreted as the
difference between two values: the value of the normal derivative, intended
with respect to the normal at the $i^{th}$ boundary node, of the $j^{th}$
function of the polynomial base, i.e., $p_{j}$, and the extrapolated value of
the same normal derivative, obtained through a RBF interpolation based on the
$m_{I}$ internal nodes only, both evaluated at the $i^{th}$ boundary node. On
the other hand, entry $(i,j)$ of matrix $\mathbf{E}$ can again be interpreted
as the difference between two values: the value of $p_{j}$ and the
extrapolated value of the same function $p_{j}$, obtained through a RBF
interpolation based on the $m_{I}$ internal nodes only, both evaluated at the
$i^{th}$ boundary node.

Therefore, if the internal stencil is unisolvent for the chosen degree of the
polynomial augmentation, each entry of matrices $\partial\mathbf{E}$ and
$\mathbf{E}$ vanishes as $\varepsilon\to 0$ and the same holds for the entries
of matrix $\mathbf{S}$ except for the diagonal entry corresponding to the
constant monomial of the polynomial base, which tends to a constant value. In
the end, matrix $\Delta\mathbf{W}$ vanishes rapidly as $\varepsilon\to 0$, and
therefore in this limit the polynomial augmentation does not affect the
formulation presented in the previous sections to study the effect of the
normals on the interpolation matrix.

$\varepsilon s=0.1$

$\varepsilon s=1$

$\varepsilon s=10$

Figure 8: optimal directions for the reference stencil of Figure 1 without
polynomial augmentation (black solid vectors) and with polynomial augmentation
with degree $P=2$ (red dashed vectors).

The influence of the polynomial augmentation, with degree $P=2$, on the
optimal directions for the reference stencil of Figure 1 is shown in Figure 8
in the case of MQ RBF for different values of the shape parameter
$\varepsilon$. In the case $\varepsilon s=0.1$ there is no noticeable
difference in the optimal directions with and without polynomial augmentation,
confirming the previous theoretical observations that the influence of the
polynomial augmentation is negligible for sufficiently small values of
$\varepsilon$. For larger values of $\varepsilon$, i.e., in the cases
$\varepsilon s=1$ and $\varepsilon s=10$, the differences in the optimal
direction start to be visually noticeable, although very small and almost
negligible from a practical point of view.

In conclusion we remark that, since in practical cases $\varepsilon
s\in[0.2,0.6]$ with MQ RBF [19] and $m>2q$ [14] in equation (2), the
calculation of the optimal directions can be reasonably carried out without
taking into account the polynomial terms.

##  4 Techniques for an improved interpolation

###  4.1 Approach 1: boundary node selection based on optimal directions

Once the optimal directions $\hat{\boldsymbol{n}}_{i}$ have been calculated
according to equation (59), it is possible to extend Approach 1 in section
3.2.2 to the case with multiple boundary nodes, i.e., discarding from the
stencil those boundary nodes whose normal vectors are too different from the
corresponding optimal normals. This strategy allows to rule out the
possibility of a singular interpolation matrix $\mathbf{M}$ by dropping both
the rows and columns associated with dangerous nodes.

The boundary node selection can be implemented as an iterative process (see
B.2) that allows to discard the $i^{th}$ boundary node whenever the dot
product between the associated normal $\bar{\boldsymbol{n}}_{i}$ and the
corresponding optimal direction $\hat{\boldsymbol{n}}_{i}$ is less than a
certain threshold value $d_{min}$:

$\text{discard boundary node $i$ if
}|\bar{\boldsymbol{n}}_{i}\cdot\hat{\boldsymbol{n}}_{i}|<d_{min}$ (69)

where the absolute value comes from the fact that the reversal of any normal
has no effect on the interpolation.

The condition number $\kappa(\mathbf{M})$, the Lebesgue constant $\Lambda_{I}$
and the number of removed nodes $N_{rem}$ for the reference stencil of Figure
1 are shown in Figure 9 as functions of the angle $\alpha$. From a comparison
with Figure 2 we can see that no more singularities appear for $\alpha<0$,
being now $\kappa(\mathbf{M})$ and $\Lambda_{I}$ bounded, thus leading to a
more stable interpolation. Figure 9 is obtained with MQ RBF ($\varepsilon
s=0.5$) and with a threshold value $d_{min}=0.6$ in equation (69). Lower
values are not able to bound the Lebesgue constant $\Lambda_{I}$ when $\alpha$
approaches the limit values of $\pi/2$ and $-\pi/2$. When 6 nodes are removed,
only the central boundary node associated to a vertical normal is included in
the stencil. $N_{rem}$ is always even due to the symmetry of the reference
stencil.

Figure 9: condition number $\kappa(\mathbf{M})$ of the interpolation matrix
(top left), and Lebesgue constant $\Lambda_{I}$ (top right) for the RBF
interpolation on the reference stencil of Figure 1 when the selection of
boundary nodes is based on the optimal directions. Bottom row: number of
removed nodes $N_{rem}$ (the same plot is reported twice in order to simplify
comparisons).

###  4.2 Approach 2: optimal placement for boundary nodes

As hinted in section 3.2.2, another strategy for avoiding a singular
interpolation matrix consists in moving the stencil nodes in an appropriate
way. It follows that any implementation of this strategy involves a
modification of the node placement which also depends on the direction of the
normals.

In section C the optimization of the boundary node position is discussed, it
emerges that the interpolation can be improved by projecting on the boundary
along the normals those inner nodes which are in the immediate neighborhood.
In in Figure 16 the projected nodes are depicted as asterisks, when the
curvature of the boundary is limited they provide a very good approximation
for the optimal locations for boundary nodes. In practice, the adoption of the
projected boundary nodes can be implemented with minimal computational cost,
here we discuss the attained improvement on the interpolation properties when
this is done for the reference stencil of Figure 1.

Figure 10: condition number of the interpolation matrix (left) and Lebesgue
constant $\Lambda_{I}$ (right) for the RBF interpolation on the reference
stencil of Figure 1 with projected boundary nodes.

In Figure 10 the condition number $\kappa(\mathbf{M})$ and the Lebesgue
constant $\Lambda_{I}$ are plotted against the angle $\alpha$ when the
boundary nodes of the reference stencil are replaced by the projected ones.
The angle $\alpha$ varies in a smaller interval, $\alpha\in[-\pi/3\,,\,\pi/3]$
since for large values the projected nodes are placed too close together or
too far from one another. From Figure 10 emerges that, for moderate angles,
i.e. approx $\alpha\geq-\pi/8$, the node projection alone is able to avoid the
appearance of ill-conditioning problems. About this, we remark that the
reference stencil of Figure 1 can be used to adequately model the case of a
curved boundary only for small angles.

Finally, since the use of the projected boundary nodes doesn’t impact the
solution process it can be applied in conjunction with the selection based on
the optimal normals explained above, the two strategies are not mutually
exclusive.

Figure 11: interpolation error
$|u(\boldsymbol{x}_{eval})-u^{h}(\boldsymbol{x}_{eval})|$ at a given location
$\boldsymbol{x}_{eval}=(s/2,s/2)$ for the reference stencil of Figure 1. No
improvement strategy (top left), node selection based on the optimal
directions (top right) with the corresponding number $N_{rem}$ of removed
boundary nodes (bottom right), boundary node projection (bottom left). The
interpolated function is $u=\exp{(x+2y)}$.

In Figure 11 the interpolation error obtained with the different techniques
described so far is displayed against the angle $\alpha$. The interpolation
error $|u(\boldsymbol{x}_{eval})-u^{h}(\boldsymbol{x}_{eval})|$ is evaluated
at the point $\boldsymbol{x}_{eval}=(s/2,s/2)$, willingly chosen not to be on
the symmetry axis, where $s$ is the distance between any two inner nodes of
the reference stencil of Figure 1 and the analytic function is defined as
$u=\exp{(x+2y)}$. For the selection of boundary nodes based on the optimal
directions (right column) the same settings were used as in Figure 9.

We can see once again that both stabilization techniques succeed in improving
the interpolation for small variations of the angle $\alpha$. For the
projection of the boundary nodes, $\alpha$ varies within the interval
$\alpha\in[-\pi/4\,,\,\pi/4]$ for the reasons explained above. We can see that
the spikes in the errors on the left column, corresponding to singular
configurations, are located as in Figures 2 and 10, this confirms the
reliability of the Lebesgue constant $\Lambda_{I}$ as a measure for accuracy
and stability of the interpolation scheme.

##  5 Applications

###  5.1 Stability of the Helmholtz-Hodge Decomposition

In this section the boundary node selection and the node projection
strategies, discussed above, are applied to the stabilization of the
Helmholtz-Hodge decomposition (HHD) [40]. The HHD plays a pivotal role in many
theoretical and practical applications [41], among others, it is at the base
of the projection methods for the numerical solution of the incompressible
Navier-Stokes equations [42, 43]. According to the HHD, under certain
hypothesis, any vector field $\boldsymbol{u}^{*}$ can be expressed as a sum of
the gradient of a scalar potential $\phi$ and a divergence-free vector field
$\boldsymbol{u}$:

$\boldsymbol{u}^{*}=\nabla\phi+\boldsymbol{u}$ (70)

When solving incompressible Navier-Stokes equations through projection
methods, the computation of the unknown velocity field $\boldsymbol{u}$ is
carried out iteratively or by marching in time. In any case, each iteration or
time step is subdivided into two substeps: first an intermediate velocity
$\boldsymbol{u}^{*}$ is obtained from the momentum equations, then
$\boldsymbol{u}$ is obtained from equation (70) as follows:

$\boldsymbol{u}=\boldsymbol{u}^{*}-\nabla\phi$ (71)

where the unknown scalar potential $\phi$ is calculated from equation (72),
that results from taking the divergence of equation (70), recalling
$\nabla\cdot\boldsymbol{u}=0$:

$\nabla^{2}\phi=\nabla\cdot\boldsymbol{u}^{*}$ (72)

In the case of known velocity at the boundary, e.g., in a cavity, homogeneous
Neumann BCs must be enforced in order to solve Poisson equation (72):

$\partial\phi/\partial\boldsymbol{n}=0$ (73)

Figure 12: test domain with $N_{I}\approx 3,300$ internal nodes.

Since the HHD expressed by equations (71)-(72) is repeated many times, once
for every iteration (or time step), it is critical that the discretization of
such equations is stable, i.e., it does not introduce any spurious mode which
grows indefinitely, thus compromizing the whole simulation. In order to assess
the stability properties of the strategies presented in this work, it was
decided to apply multiple times the projection scheme described in equations
(71)-(72) starting from a given vector field $\boldsymbol{w}$, defined on the
complex-shaped domain depicted in Figure 12. The boundary is characterized by
locally convex and concave features with varying curvature, designed to
trigger ill-conditioning issues described in section 3. A Lagrange multiplier
is employed to solve the linear systems arising from the RBF-FD discretization
of equation (72) since Neumann BCs are prescribed on the whole boundary [44].

The pseudocode for the stability test is presented in Algorithm 1. It is worth
mentioning that any conservative discretization scheme, e.g., FVM, implicitly
satisfies the previous stability test since the value of the discretized
divergence operator vanishes after the first iteration. On the other hand, the
basic formulation of the RBF-FD discretization, as presented in section 2.2,
is not conservative, therefore the value of the discretized divergence
operator never vanishes at the nodes. Algorithm 1 might therefore diverge if
the Poisson equation is solved with the basic RBF-FD method without a proper
treatment of the Neumann BCs.

$P=2$

$P=3$

$P=4$

Figure 13: stability range for the $d_{min}$ parameter in the case of approach
1, i.e., boundary node selection based on optimal directions. Upper and lower
bounds for $d_{min}$ values are displayed with red and blue curves,
respectively. Algorithm 1 Pseudocode for the stability test of repeated HHD

1:Initialize: $\boldsymbol{u}^{(1)}\leftarrow\boldsymbol{w}$

2:for $i=1,2,\ldots$ do

3: Solve $\nabla^{2}\phi=\nabla\cdot\boldsymbol{u}^{(i)}$ for $\phi$ with
homogeneous Neumann BCs

4: Update: $\boldsymbol{u}^{(i+1)}\leftarrow\boldsymbol{u}^{(i)}-\nabla\phi$

5:end for

Since the equations involved in Algorithm 1 are linear, the choice of the
initial vector field $\boldsymbol{w}$ does not affect the stability of the
iterative process. The employed initial vector field $\boldsymbol{w}$ is
chosen to be an irrotational vortex for simplicity: the tangential velocity is
$w_{\theta}=r^{-1}$ where $r$ is the distance from the origin.

The results for approach 1 presented in section 4.1, i.e., boundary node
selection, are shown in Figure 13 in the case of MQ RBF with polynomial
degrees $P=2,3,4$ and for different values of the shape parameter
$\varepsilon$, such that $\varepsilon s\in[0.2,0.9]$. The total number of
nodes inside the domain is $N_{I}\approx 15,000$. The number of internal nodes
for each stencil is chosen to follow the rule $m_{I}=2q$ [14], see equation
(2): $m_{I}=20$ for $P=3$ and $m_{I}=30$ for $P=4$. In the case $P=2$ it was
decided to use slightly larger stencils with $m_{I}=15$ internal nodes, since
this enhanced stability. In Figure 13 the values of $d_{min}$ allowing a
stable computation were identified as those lying between the two curves. As
expected, too small values of $d_{min}$, e.g., $d_{min}<0.2$, lead to
instability due to stencils with boundary nodes whose actual normal is too
different from the corresponding optimal direction, resulting in ill-
conditioned local interpolants. On the other hand, too large values of
$d_{min}$, e.g., $d_{min}>0.95$, lead to instability due to stencils with too
few boundary nodes enforcing the prescribed Neumann BCs (73). From the same
figure it is possible to observe that the stability range is slightly reduced
when $P$ is increased from $P=2$ to $P=4$, while a general and reasonably good
choice is $d_{min}=0.7$.

No stabilization

Approach 1

Approach 2

Figure 14: condition number of the local interpolation matrix. Top row, from
left to right: no stabilization, boundary selection based on optimal
directions, boundary projection. Bottom row: enlarged visualization of a
portion of the plots of the top row.

In the case of approach 2 presented in section 4.2, i.e., by employing
projected boundary nodes, stability is always attained for $P=2,3,4$ in the
range $\varepsilon s\in[0.2,0.9]$ previously considered, highlighting the
remarkable stabilization effect of this approach. Nevertheless, stability
issues are still encountered for shape factors outside this range.

Figure 14 shows the effect of the proposed stabilization approaches on the
condition number of the local interpolation matrix, in the case of MQ RBF with
$\varepsilon s=0.5$ and $P=3$ (results with different shape parameters
$\varepsilon$ and polynomial degrees $P$ are qualitatively similar). As
expected from the theoretical observations presented in section 3, stencils
close to locally concave boundaries can be very ill-conditioned if no
stabilization technique is employed (red nodes in Figure 14, left), leading to
large errors and stability issues.

By using the stabilization approach 1 with $d_{min}=0.7$, this issue is
properly addressed and ill-conditioned stencil are no longer present (Figure
14, center). The same result is attained by employing the stabilization
approach 2 (Figure 14, right). In the latter case we observe an increase in
the condition number for the stencils close to locally convex boundaries: this
is due to the employed projection technique that is performed from the
boundary inward. Therefore, internal nodes projected from locally convex
boundaries are very close to each other, leading to an increased condition
number. This geometrical issue is common to any simple geometrical projection
technique, regardless of the direction and starting locations employed for the
projection. With reference to arbitrary geometries, such problem can only be
addressed by taking into account the need for projected boundary nodes already
at the node generation phase, thus requiring a proper modification of the node
generation algorithm.

###  5.2 Accuracy

Besides the stability considerations presented in the previous section, we
also provide some insight about the influence of the proposed stabilization
approaches on the accuracy of the RBF-FD discretization. In order to do so, we
consider the following Poisson equation:

$\nabla^{2}u=f$ (74)

defined over the previously employed domain depicted in Figure 12. Neumann BCs
are prescribed on the whole boundary:

$\frac{\partial u}{\partial\boldsymbol{n}}=g$ (75)

and the analytic solution is chosen to be $u=r^{-1}$, where $r$ is the
distance from the origin. The node distributions and the parameters for the
RBF-FD discretization are the same as those used in the previous section 5.1,
and a Lagrange multiplier is again employed to solve the corresponding
discretized linear system.

Figure 15: Errors for the Poisson problem (74) when employing stabilization
approach 1: initial node distribution (left), projected boundary nodes
(right).

The normalized RMS errors are shown in Figure 15 when approach 1 is applied to
both node distributions without and with projected boundary nodes. The latter
case is therefore a combination of both the proposed stabilization techniques.
In the former case, and for all polynomial degrees $P=2,3,4$, the error is
lower when $d_{min}$ ranges from 0.4 to 0.8. As expected, the error grows when
$d_{min}$ is decreased below 0.4 because of the increasing number of ill-
conditioned stencils near the boundary. The limit case $d_{min}=0$ corresponds
indeed to no stabilization action at all. On the other hand, the error also
grows when $d_{min}$ is increased beyond 0.8 because of the increasing number
of boundary stencils with an insufficient number of boundary nodes where BCs
are enforced.

When both the stabilization approaches are employed, Figure 15 right, the
increase in the error for $d_{min}<0.4$ is no longer present since the ill-
conditioning issues are already addressed by the projected nodes on the
boundary, while the increase for $d_{min}>0.8$ is obviously still present,
i.e., when too few boundary nodes are included in the boundary stencils. In
this case we also observe that the boundary node selection based on the
optimal directions does not affect the accuracy of the discretization at all
if $d_{min}<0.8$, as expected.

##  6 Conclusions

In this paper the problem of ill-conditioning and the subsequent instabilities
affecting the application of RBF-FD method in presence of Neumann boundary
condition have been analyzed both theoretically and numerically. The
theoretical motivations for the onset of such issues are derived by
highlighting the dependency of the determinant of the local interpolation
matrix upon the boundary normals. Qualitative investigations are also carried
out numerically by studying a reference stencil and looking for correlations
between its geometry and the properties of the associated interpolation
matrix.

Once the problem had been described with sufficient detail, two approaches are
proposed, namely approach 1, i.e., the boundary node selection based on the
optimal normals, and approach 2, i.e., the optimal placement of boundary
nodes. Many evidences indicate that both can be applied to the solution of
numerical problems of practical interest.

In order to provide a better demonstration of the onset of instabilities due
to ill-conditioning of the interpolation matrix, the stabilization of the
Helmholtz-Hodge decomposition is discussed in the last section. This
demonstration helped to put into context the role of the approaches in the
stabilization of the solution process of any boundary value problem involving
Neumann boundary conditions.

Summing up, the main pros and cons of the two stabilization strategies are the
following. Approach 1 is the most robust and it does not require any control
over the node placement. However, the price to pay for its application is a
higher computational effort, required to estimate the optimal directions for
all the boundary stencils. We remark that this cost can be significantly
reduced if the optimal directions are approximated with lower accuracy. Future
developments might significantly improve the computational efficiency of
approach 1 through the application of machine learning techniques, for
example. Furthermore, the mathematical formulation behind this approach leaves
great freedom for further developments. Approach 2, on the other side, does
not involve any additional computation and is even more effective at ensuring
the stabilization of the RBF-FD method. Unfortunately, it requires the
implementation of a more sophisticated node generation algorithm, indeed, more
complex geometries might present high curvatures at the boundary thus making
it very difficult to develop a generic node projection technique. At last we
remark that both approaches can be applied without any modification to the 3D
cases, where all of the considerations made in this paper remain valid.

The application of the stabilization techniques described above does not
require the user to intervene during the simulation process and does not
require any parameter tuning. Even the $d_{min}$ parameter can be set once for
all to be $d_{min}=0.7$ without any contraindication. This fact is especially
important since the major advantage of the RBF-FD method lies in its ability
to free the user from all of the considerations involved in the mesh
generation, which in turn require a specific expertise. Indeed, for this
method to be competitive with traditional ones, it must succeed in making the
whole simulation process easier to handle, without sacrificing the reliability
of the attained solutions.

Furthermore, the presented approaches can be suitably combined to the
appropriate extent to provide even more robust RBF-FD meshless discretizations
of boundary value problems.

In conclusion, we believe that the approaches presented in this work are
successful not only at solving the stability issues of the RBF-FD method, but
also at preserving its potential advantages.

##  Appendix A $d$-matrices

######  Definition ($d$-matrix).

An $m\times n$ $d$-matrix is defined as the following $m\times n$ matrix
$\mathcal{A}=(\boldsymbol{a}_{ij})$ of $d$-dimensional vectors
$\boldsymbol{a}_{ij}\in\mathbb{R}^{d}$:

$\mathcal{A}=\begin{bmatrix}\boldsymbol{a}_{11}\\!\\!\\!&\cdots\\!\\!\\!&\boldsymbol{a}_{1n}\\\
\vdots\\!\\!\\!&\ddots\\!\\!\\!&\vdots\\\
\boldsymbol{a}_{m1}\\!\\!\\!&\cdots\\!\\!\\!&\boldsymbol{a}_{mn}\\\
\end{bmatrix}$ (76)

and the set of all $m\times n$ $d$-matrices will be denoted by
$\mathbb{R}^{m\times n,d}$. For simplicity $\mathbb{R}^{m\times 1,d}$ will be
called $\mathbb{R}^{m,d}$.

The usual operations are defined as follows:

  * 1.

product of $\mathcal{A}=(\boldsymbol{a}_{ij})\in\mathbb{R}^{m\times n,d}$ with
a scalar $\lambda\in\mathbb{R}$:

$\lambda\mathcal{A}=(\lambda\boldsymbol{a}_{ij})\in\mathbb{R}^{m\times n,d}$
(77)

  * 2.

sum of $\mathcal{A}=(\boldsymbol{a}_{ij})\in\mathbb{R}^{m\times n,d}$ and
$\mathcal{B}=(\boldsymbol{b}_{ij})\in\mathbb{R}^{m\times n,d}$:

$\mathcal{A}+\mathcal{B}=(\boldsymbol{a}_{ij}+\boldsymbol{b}_{ij})\in\mathbb{R}^{m\times
n,d}$ (78)

  * 3.

product of $\mathcal{A}=(\boldsymbol{a}_{ij})\in\mathbb{R}^{m\times n,d}$ with
$\mathbf{Q}=(q_{ij})\in\mathbb{R}^{n\times p}$:

$\mathcal{A}\mathbf{Q}=\left(\sum_{k=1}^{n}\boldsymbol{a}_{ik}q_{kj}\right)\in\mathbb{R}^{m\times
p,d}$ (79)

For $\mathcal{A}=(\boldsymbol{a}_{ij})\in\mathbb{R}^{m\times n,d}$ and
$\mathcal{V}=(\boldsymbol{v}_{i})\in\mathbb{R}^{m,d}$, by using the dot
product $\cdot$ in $\mathbb{R}^{d}$ we can define the following operator H:

$H(\mathcal{A},\mathcal{V})=(\boldsymbol{a}_{ij}\cdot\boldsymbol{v}_{i})\in\mathbb{R}^{m\times
n}$ (80)

Given $\mathcal{A},\mathcal{B}\in\mathbb{R}^{m\times n,d}$,
$\mathcal{V}=(\boldsymbol{v}_{i})\in\mathbb{R}^{m,d}$ and
$\mathbf{Q}\in\mathbb{R}^{n\times p}$, the following equalities hold:

$\displaystyle H(\mathcal{A},\mathcal{V})+H(\mathcal{B},\mathcal{V})$
$\displaystyle=H(\mathcal{A}+\mathcal{B},\mathcal{V})$ (81) $\displaystyle
H(\mathcal{A},\mathcal{V})\mathbf{Q}$
$\displaystyle=H(\mathcal{A}\mathbf{Q},\mathcal{V})$ (82)
$\displaystyle\frac{\partial}{\partial
v_{i,\eta}}\det\\!\big{(}H(\mathcal{A},\mathcal{V})\big{)}$
$\displaystyle=\det\\!\big{(}H(\mathcal{A},\mathcal{V}_{i,\eta})\big{)}$ (83)

for $i=1,\ldots,m$ and $\eta=1,\ldots,d$. By considering the canonical basis
$\\{\boldsymbol{e}_{\eta}\\}_{1}^{d}$ of $\mathbb{R}^{d}$ for simplicity, in
the last equation $v_{i,\eta}$ is the $\eta^{th}$ component of the $i^{th}$
vector $\boldsymbol{v}_{i}$ in $\mathcal{V}$, while $\mathcal{V}_{i,\eta}$ is
obtained from $\mathcal{V}$ by replacing its $i^{th}$ vector
$\boldsymbol{v}_{i}$ with $\boldsymbol{e}_{\eta}$. Equations (81)-(82) follow
from the simple distributive property of the dot product, while equation (83)
follows from the multilinearity of the determinant and from the definition of
the operator $H$ in equation (80).

##  Appendix B Optimal directions

###  B.1 Generic properties

By defining the Lagrangian function $\mathcal{L}$ with multipliers
$\boldsymbol{\nu}=\\{\nu_{1},\ldots,\nu_{m_{B}}\\}$ as follows:

$\mathcal{L}(\mathcal{N},\boldsymbol{\nu})=\det\\!\big{(}H(\mathcal{G}_{BB},\mathcal{N})\big{)}-\sum_{i=1}^{m_{B}}\nu_{i}\left(\|\bar{\boldsymbol{n}}_{i}\|_{2}^{2}-1\right)$
(84)

the solution of problem (59) is given by the proper solution of the following
system of equations:

$\begin{cases}\displaystyle\frac{\partial}{\partial\hat{n}_{i,\eta}}\mathcal{L}(\hat{\mathcal{N}},\boldsymbol{\nu})=\frac{\partial}{\partial\hat{n}_{i,\eta}}\det\\!\big{(}H(\mathcal{G}_{BB},\hat{\mathcal{N}})\big{)}-2\nu_{i}\hat{n}_{i,\eta}=0\\\\[10.00002pt]
\|\hat{\boldsymbol{n}}_{i}\|_{2}^{2}=1\end{cases}$ (85)

for $i=1,\ldots,m_{B}$ and $\eta=1,\ldots,d$. By looking at the first equation
of system (85), we see that the Lagrange multipliers $\nu_{i}$ act as
multiplicative factors allowing the constraints in problem (59) to be
fulfilled. Therefore, the solution
$\hat{\mathcal{N}}=\\{\hat{\boldsymbol{n}}_{1},\ldots,\hat{\boldsymbol{n}}_{m_{B}}\\}\in\mathbb{R}^{m_{B},d}$
can be explicitly written through the $m_{B}$ auxiliary vectors
$\boldsymbol{t}_{i}=\\{t_{i,1},\ldots,t_{i,d}\\}\in\mathbb{R}^{d}$ as follows:

$\displaystyle t_{i,\eta}$
$\displaystyle=\frac{\partial}{\partial\hat{n}_{i,\eta}}\det\\!\big{(}H(\mathcal{G}_{BB},\hat{\mathcal{N}})\big{)}=\det\\!\big{(}H(\mathcal{G}_{BB},\hat{\mathcal{N}}_{i,\eta})\big{)}$
(86) $\displaystyle\hat{\boldsymbol{n}}_{i}$
$\displaystyle=\pm\frac{\boldsymbol{t}_{i}}{\|\boldsymbol{t}_{i}\|_{2}}$ (87)

for $i=1,\ldots,m_{B}$ and $\eta=1,\ldots,d$. Property (83) is used in
equation (86), where $\hat{\mathcal{N}}_{i,\eta}$ is obtained from
$\hat{\mathcal{N}}$ as explained in appendix A. We note that the term
$\hat{\mathcal{N}}_{i,\eta}$ depends on $\hat{\boldsymbol{n}}_{j\neq i}$,
i.e., it does not depend on $\hat{\boldsymbol{n}}_{i}$ itself. Equation (87)
then acts as a normalization of $\boldsymbol{t}_{i}$ through the Lagrange
multiplier $\nu_{i}$ as previously noted, i.e.,
$2\nu_{i}=\pm\|\boldsymbol{t}_{i}\|_{2}$. Equations (86) and (87) simply state
that each solution vector $\hat{\boldsymbol{n}}_{i}$ is parallel to the
gradient of the target function with respect to that normal, in accordance
with the principles of the constrained optimization.

From equation (57) we see that $\det(\mathbf{S}_{BB})$ is an odd function with
respect to each normal $\bar{\boldsymbol{n}}_{i}$, i.e., the reversal of any
$\bar{\boldsymbol{n}}_{i}$ results in a change of sign of
$\det(\mathbf{S}_{BB})$ because of the properties of the determinant.
Therefore the solution $\hat{\mathcal{N}}$ of the maximization problem (59) is
also a solution when an even number of normals are reversed. Furthermore,
analogously to the case with one boundary node, it is also a solution of the
corresponding minimization problem when an odd number of normals are reversed,
in accordance with the previous interpretation of the Lagrange multipliers. In
other words, the target function in equation (59) has $2^{m_{B}}$ local
extrema which are equivalent, hence the symbol $\pm$ in equation (87) for each
of the $m_{B}$ normals. This is obviously expected: the reversal in the
directions of the normals in the Neumann BCs (5) has no effect on the
properties of the interpolant, indeed we are interested in the normals that
maximize the magnitude of the target function $\det(\mathbf{S}_{BB})$
regardless of its sign. When $m_{B}>1$ the target function has saddle points
which therefore are also solutions of equations (86)-(87), and it can also
have multiple distinct local extrema, i.e., considering the reversal of any
normal as an equivalent solution, as shown as follows.

Let us consider the case $m_{B}=2$ as an example. The target function
$\det(\mathbf{S}_{BB})$ defined by equation (57) can be expressed as follows:

$\begin{split}\det(\mathbf{S}_{BB})&=\begin{vmatrix}\boldsymbol{g}_{1,1}\cdot\bar{\boldsymbol{n}}_{1}&\boldsymbol{g}_{1,2}\cdot\bar{\boldsymbol{n}}_{1}\\\
\boldsymbol{g}_{2,1}\cdot\bar{\boldsymbol{n}}_{2}&\boldsymbol{g}_{2,2}\cdot\bar{\boldsymbol{n}}_{2}\end{vmatrix}=\\\
&=(\boldsymbol{g}_{1,1}\cdot\bar{\boldsymbol{n}}_{1})(\boldsymbol{g}_{2,2}\cdot\bar{\boldsymbol{n}}_{2})-(\boldsymbol{g}_{1,2}\cdot\bar{\boldsymbol{n}}_{1})(\boldsymbol{g}_{2,1}\cdot\bar{\boldsymbol{n}}_{2})=\bar{\boldsymbol{n}}_{1}\cdot\mathbf{G}\bar{\boldsymbol{n}}_{2}\end{split}$
(88)

where:

$\mathbf{G}=\boldsymbol{g}_{1,1}\otimes\boldsymbol{g}_{2,2}-\boldsymbol{g}_{1,2}\otimes\boldsymbol{g}_{2,1}$
(89)

is thus a $d\times d$ matrix having exactly $d-2$ null eigenvalues in the
general case where $\boldsymbol{g}_{1,1}$ and $\boldsymbol{g}_{2,2}$ are not
parallel to $\boldsymbol{g}_{1,2}$ and $\boldsymbol{g}_{2,1}$, respectively,
and all $\boldsymbol{g}_{i,j}$ are non-zero vectors. The explicit solution
expressed by equations (86)-(87) can therefore be simplified as follows:

$\begin{split}\hat{\boldsymbol{n}}_{1}&=\pm\frac{\mathbf{G}\hat{\boldsymbol{n}}_{2}}{\|\mathbf{G}\hat{\boldsymbol{n}}_{2}\|_{2}}\\\
\hat{\boldsymbol{n}}_{2}&=\pm\frac{\mathbf{G}^{T}\hat{\boldsymbol{n}}_{1}}{\|\mathbf{G}^{T}\hat{\boldsymbol{n}}_{1}\|_{2}}\end{split}$
(90)

which can be combined to yield the following decoupled equations:

$\begin{split}\hat{\boldsymbol{n}}_{1}&=\frac{\mathbf{Q}_{1}\hat{\boldsymbol{n}}_{1}}{\|\mathbf{Q}_{1}\hat{\boldsymbol{n}}_{1}\|_{2}}\\\
\hat{\boldsymbol{n}}_{2}&=\frac{\mathbf{Q}_{2}\hat{\boldsymbol{n}}_{2}}{\|\mathbf{Q}_{2}\hat{\boldsymbol{n}}_{2}\|_{2}}\end{split}$
(91)

where $\mathbf{Q}_{1}=\mathbf{G}\mathbf{G}^{T}$ and
$\mathbf{Q}_{2}=\mathbf{G}^{T}\mathbf{G}$. Equations (91) represent therefore
two classic eigenvalue problems with matrix $\mathbf{Q}_{1}$ for
$\hat{\boldsymbol{n}}_{1}$ and matrix $\mathbf{Q}_{2}$ for
$\hat{\boldsymbol{n}}_{2}$. $\mathbf{Q}_{1}$ and $\mathbf{Q}_{2}$ have the
same real non-negative eigenvalues $\\{\lambda_{i}\\}_{1}^{d}$ since these
matrices are symmetric and because their definitions imply
$\mathbf{G}^{T}\mathbf{Q}_{1}=\mathbf{Q}_{2}\mathbf{G}^{T}$ and
$\mathbf{G}\mathbf{Q}_{2}=\mathbf{Q}_{1}\mathbf{G}$.

In $d=2$ dimensions the eigenvalues are $0<\lambda_{1}<\lambda_{2}$ in the
general case, therefore the possible distinct solutions of equations (90) are
four: two for each of the two normals. It can be shown that the solutions
corresponding to complementary eigenvalues, i.e., $(\lambda_{1},\lambda_{2})$
or $(\lambda_{2},\lambda_{1})$ for $\hat{\boldsymbol{n}}_{1}$ and
$\hat{\boldsymbol{n}}_{2}$, although solutions of equations (91), are not
solutions of equations (90). Therefore the two remaining solutions are those
associated to the same eigenvalue for both normals, i.e.,
$(\lambda_{1},\lambda_{1})$ or $(\lambda_{2},\lambda_{2})$. By combining the
last expression of equation (88) with the first equation (90) with positive
sign, we obtain:

$\det(\mathbf{S}_{BB})\Big{|}_{\begin{subarray}{c}\bar{\boldsymbol{n}}_{1}=\hat{\boldsymbol{n}}_{1}\\\
\bar{\boldsymbol{n}}_{2}=\hat{\boldsymbol{n}}_{2}\end{subarray}}=\|\mathbf{G}\hat{\boldsymbol{n}}_{2}\|_{2}=\sqrt{\lambda_{i}}$
(92)

where the last equality comes from the explicit form of the second eigenvalue
problem (91), i.e.,
$\mathbf{G}^{T}\mathbf{G}\hat{\boldsymbol{n}}_{2}=\lambda_{i}\hat{\boldsymbol{n}}_{2}$.
From equation (92) we can see that the solution associated to the largest
eigenvalue $\lambda_{2}$ is the sought solution of problem (59), while the
solution associated to the smallest eigenvalue $\lambda_{1}$ is therefore a
saddle point.

In $d=3$ dimensions there is a null eigenvalue since $\mathbf{G}$ also has one
null eigenvalue. This implies that in the 3D case with $m_{B}=2$ boundary
nodes, for each boundary node there is a normal which leads to a singular
interpolation matrix regardless of the other normal. Nonetheless, the normals
associated to null eigenvalues can not be solutions of equations (90),
therefore even in the 3D case there are only two solutions, namely a local
extremum and a saddle point.

Unfortunately, when $m_{B}=3$ in $d=2$ dimensions, there can be two distinct
local extrema in addition to saddle points. Indeed, if we consider the
following particular symmetric choice for the vectors $\boldsymbol{g}_{i,j}$:

$\begin{split}\boldsymbol{g}_{1,1}&=\boldsymbol{g}_{2,3}=\boldsymbol{g}_{3,2}=\\{1,0\\}^{T}\\\
\boldsymbol{g}_{2,2}&=\boldsymbol{g}_{1,3}=\boldsymbol{g}_{3,1}=\\{0,0\\}^{T}\\\
\boldsymbol{g}_{3,3}&=\boldsymbol{g}_{1,2}=\boldsymbol{g}_{2,1}=\\{0,1\\}^{T}\end{split}$
(93)

then the determinant in equation (57) becomes:

$\begin{split}\det(\mathbf{S}_{BB})&=\begin{vmatrix}\boldsymbol{g}_{1,1}\cdot\bar{\boldsymbol{n}}_{1}&\boldsymbol{g}_{1,2}\cdot\bar{\boldsymbol{n}}_{1}&\boldsymbol{g}_{1,3}\cdot\bar{\boldsymbol{n}}_{1}\\\
\boldsymbol{g}_{2,1}\cdot\bar{\boldsymbol{n}}_{2}&\boldsymbol{g}_{2,2}\cdot\bar{\boldsymbol{n}}_{2}&\boldsymbol{g}_{2,3}\cdot\bar{\boldsymbol{n}}_{2}\\\
\boldsymbol{g}_{3,1}\cdot\bar{\boldsymbol{n}}_{3}&\boldsymbol{g}_{3,2}\cdot\bar{\boldsymbol{n}}_{3}&\boldsymbol{g}_{3,3}\cdot\bar{\boldsymbol{n}}_{3}\end{vmatrix}=\\\
&=\begin{vmatrix}\cos\alpha_{1}&\sin\alpha_{1}&0\\\
\sin\alpha_{2}&0&\cos\alpha_{2}\\\
0&\cos\alpha_{3}&\sin\alpha_{3}\end{vmatrix}=\\\
&=-\big{(}\cos\alpha_{1}\cos\alpha_{2}\cos\alpha_{3}+\sin\alpha_{1}\sin\alpha_{2}\sin\alpha_{3}\big{)}\end{split}$
(94)

where the 2D normals are parametrized as usual as
$\bar{\boldsymbol{n}}_{i}=\\{\cos\alpha_{i},\sin\alpha_{i}\\}^{T}$.

From the last trigonometric expression in equation (94) it can be seen that
the extremum $\alpha_{1}=\alpha_{2}=\alpha_{3}=0$, including the equivalent
solutions obtained by reversing any of the normals, i.e., $\alpha_{i}=\pi$, is
complemented by another distinct extremum
$\alpha_{1}=\alpha_{2}=\alpha_{3}=\pi/2$, including the corresponding
equivalent solutions.

It is therefore to be expected that multiple distinct extrema can also occur
with $m_{B}>3$ boundary nodes, in both 2D and 3D cases. Nonetheless, we note
that $\boldsymbol{g}_{i,j}$ are not arbitrary vectors but rather they are
defined in equation (56) according to the RBF interpolant, i.e., according to
the relative positions of the stencil nodes. From a practical point of view,
we think that a sufficiently regular distribution of nodes in the stencil can
guarantee a unique extremum, as confirmed by numerical experiments.

###  B.2 Computation of the optimal directions

Equations (86)-(87) represent a nonlinear system of coupled algebraic
equations in the $m_{B}$ unknown vectors $\hat{\boldsymbol{n}}_{i}$, which
therefore requires an iterative solution. Such iterative process can be
directly obtained from equation (86) as follows:

$t_{i,\eta}^{(k+1)}=\det\\!\big{(}H(\mathcal{G}_{BB},\hat{\mathcal{N}}_{i,\eta}^{(k)})\big{)}$
(95)

where the superscript $(k)$ refers to iteration $k$. The values
$t_{i,\eta}^{(k+1)}$ for $i=1,\ldots,m_{B}$ and $\eta=1,\ldots,d$ are
therefore computed explicitly by using equation (95) where
$\hat{\mathcal{N}}_{i,\eta}^{(k)}$ contains the normals at iteration $k$. Then
the actual normals $\hat{\boldsymbol{n}}_{i}^{(k+1)}$ at iteration $k+1$ are
obtained once again through the normalization expressed by equation (87).

It is convenient to compute the determinant in equation (95) as follows. Let
us consider the explicit form of equation (86) where both row $i$ and column
$i$ of matrix $H(\mathcal{G}_{BB},\hat{\mathcal{N}}_{i,\eta})$ are moved to
the last row and to the last column, respectively:

$\begin{split}t_{i,\eta}&=\det\\!\big{(}H(\mathcal{G}_{BB},\hat{\mathcal{N}}_{i,\eta})\big{)}=\\\
&=\begin{vmatrix}\begin{array}[]{ccc|c}\boldsymbol{g}_{1,1}\cdot\hat{\boldsymbol{n}}_{1}&\cdots&\boldsymbol{g}_{1,m_{B}}\cdot\hat{\boldsymbol{n}}_{1}&\boldsymbol{g}_{1,i}\cdot\hat{\boldsymbol{n}}_{1}\\\
\vdots&\ddots&\vdots&\vdots\\\
\boldsymbol{g}_{m_{B},1}\cdot\hat{\boldsymbol{n}}_{m_{B}}&\cdots&\boldsymbol{g}_{m_{B},m_{B}}\cdot\hat{\boldsymbol{n}}_{m_{B}}&\boldsymbol{g}_{m_{B},i}\cdot\hat{\boldsymbol{n}}_{m_{B}}\\\\[2.5pt]
\hline\cr\boldsymbol{g}_{i,1}\cdot\boldsymbol{e}_{\eta}&\cdots&\boldsymbol{g}_{i,m_{B}}\cdot\boldsymbol{e}_{\eta}&\boldsymbol{g}_{i,i}\cdot\boldsymbol{e}_{\eta}\end{array}\end{vmatrix}=\begin{vmatrix}\begin{array}[]{c|c}\mathbf{S}_{i}&\boldsymbol{c}_{i}\\\
\hline\cr\boldsymbol{r}_{i,\eta}&S_{ii,\eta}\end{array}\end{vmatrix}\end{split}$
(96)

In other words: starting from matrix
$H(\mathcal{G}_{BB},\hat{\mathcal{N}}_{i,\eta})$, matrix $\mathbf{S}_{i}$ is
obtained by removing row $i$ and column $i$, column vector
$\boldsymbol{c}_{i}$ is obtained by removing entry $i$ of column $i$, and row
vector $\boldsymbol{r}_{i,\eta}$ is obtained by removing entry $i$ of row $i$.

By the properties of the Schur complement of the block $\mathbf{S}_{i}$ of the
last matrix in equation (96), we have:

$t_{i,\eta}=\begin{vmatrix}\begin{array}[]{c|c}\mathbf{S}_{i}&\boldsymbol{c}_{i}\\\
\hline\cr\boldsymbol{r}_{i,\eta}&S_{ii,\eta}\end{array}\end{vmatrix}=\det(\mathbf{S}_{i})\big{(}S_{ii,\eta}-\boldsymbol{r}_{i,\eta}\mathbf{S}_{i}^{-1}\boldsymbol{c}_{i}\big{)}$
(97)

Since equation (97) holds for each of the $d$ components of the vector
$\boldsymbol{t}_{i}=\\{t_{i,\eta}\\}_{\eta=1}^{d}$, it can be written
explicitly as follows:

$\boldsymbol{t}_{i}=\det(\mathbf{S}_{i})\Big{(}\boldsymbol{g}_{i,i}-\sum_{j=1}^{m_{B}-1}w_{j}\boldsymbol{g}_{i,j^{\prime}}\Big{)}$
(98)

where $w_{j}$ are the $m_{B}-1$ entries of the column vector
$\boldsymbol{w}=\mathbf{S}_{i}^{-1}\boldsymbol{c}_{i}$ and $j^{\prime}$ is
used to index vectors $\boldsymbol{g}_{i,j}$ excluding the case $j=i$, i.e.,
$j^{\prime}=j$ if $j<i$, $j^{\prime}=j+1$ otherwise. We finally note that
$\det(\mathbf{S}_{i})$ in equation (98) does not need to be computed because
of the subsequent normalization of equation (87).

Moving back to the required iterative scheme formerly expressed by equation
(95), it can finally be expressed by equation (98) where the dependence of the
right side upon the iteration index $k$ arises only in the terms
$w_{j}=w_{j}^{(k)}$. We observe that the vector $\boldsymbol{w}$ carrying the
values $w_{j}$ does not depend upon any specific component $\eta$, and
therefore it can be computed only once at each iteration for each of the
$m_{B}$ vectors $\boldsymbol{t}_{i}$.

Since the system of equations (86)-(87) is nonlinear, it is important to
choose a proper set of initial vectors $\hat{\boldsymbol{n}}_{i}^{(0)}$ for
the iterative scheme (95) to converge to the correct solution. Numerical
experiments showed that a good starting point can be obtained from the
diagonal vectors of $d$-matrix $\mathcal{G}_{BB}$ defined in equation (56):

$\hat{\boldsymbol{n}}_{i}^{(0)}=\frac{\boldsymbol{g}_{ii}}{\|\boldsymbol{g}_{ii}\|_{2}}$
(99)

for $i=1,\ldots,m_{B}$. We note that this choice corresponds exactly to
equation (98) with $w_{j}=0$. Equation (98) can therefore be employed to
provide an improved starting point if any suitable approximation of $w_{j}$ is
available. Nonetheless, since in general the norm of $\boldsymbol{g}_{ii}$ and
$\boldsymbol{g}_{ij^{\prime}}$ are in the same order of magnitude, naive
approximations of $w_{j}\neq 0$ can lead to large errors.

Another good starting point can be obtained from the following geometric
heuristic:

$\hat{\boldsymbol{n}}_{i}^{(0)}=\frac{\boldsymbol{x}_{m_{I}+i}-\bar{\boldsymbol{x}}}{\|\boldsymbol{x}_{m_{I}+i}-\bar{\boldsymbol{x}}\|_{2}}$
(100)

where $\bar{\boldsymbol{x}}$ is some geometrical reference point for the
stencil, e.g., the centroid of the internal nodes.

###  B.3 Computational costs

The computation of the optimal directions consists of two phases: the
computation of the $d$-matrix $\mathcal{G}_{BB}$, equation (56), and the
iterative solution of the associated system, equation (95). The computation of
$\mathcal{G}_{BB}$ requires the computation of the $m_{I}\times m_{B}$ matrix
$\bar{\boldsymbol{\psi}}$ from equation (39), whose cost is therefore
$\mathcal{O}(m_{I}^{3})+\mathcal{O}(m_{B}\cdot m_{I}^{2})$, i.e.,
factorization of $\boldsymbol{\varphi}_{II}$ \+ solution of $m_{B}$ associated
linear systems. Then, each iteration updates the $m_{B}$ optimal directions
through equation (98), each of which requires the $(m_{B}-1)\times(m_{B}-1)$
linear system $\mathbf{S}_{i}\boldsymbol{w}=\boldsymbol{c}_{i}$ to be solved
for the weights $\boldsymbol{w}$, thus leading to a cost
$\mathcal{O}(N_{it}\cdot m_{B}^{4})$, where $N_{it}$ is the number of
iterations which depends upon the desired level of convergence. We note that
the optimal directions do not need to be computed with high accuracy since
they are only used in equation (69).

From a practical point of view, our numerical experiments showed that most of
the computational time is due to the iterative phase, even if the number of
boundary nodes $m_{B}$ is usually much smaller than the number of internal
nodes $m_{I}$ of a single stencil. Nonetheless, we trivially note that the
computation of the optimal directions is required only for those stencils
having boundary nodes, and the number of these boundary stencils is usually
much smaller than the total number of stencils because of obvious geometrical
reasons.

We finally note that the iterative strategy presented as approach 1 in section
4.1 requires a single calculation of the $d$-matrix $\mathcal{G}_{BB}$, since
the successive elimination of boundary nodes corresponds to the removal of the
corresponding rows and columns of $\mathcal{G}_{BB}$.

##  Appendix C Optimal position for boundary nodes

###  C.1 Optimization process

We are interested in improving the interpolation properties of the reference
stencil of Figure 1. First of all, let’s consider that moving the inner nodes
would probably affect the interpolation properties of the neighboring
stencils. Indeed, if we decided to move the inner nodes of one stencil, in
order to preserve the local uniformity of the node distribution, which ensures
that the difference between fill distance and separation distance remains
negligible, we would need to adjust the position of many other nodes
accordingly. On top of that, if we take into account the matrix $\mathbf{M}$
in the general case of equation (43) and consider that
$\det(\boldsymbol{\varphi}_{II})\neq 0$ in the cases of interest, we are once
again prompted to focus our attention on the boundary nodes.

Since the the interpolation becomes unstable as we approach those
configurations where $\det(\mathbf{S}_{BB})=0$, we might consider the problem
of maximizing $\det(\mathbf{M})$ or $\det(\mathbf{S}_{BB})$ by controlling the
locations of boundary nodes. Unfortunately, we would find that, because of the
radial basis functions, the further apart we spread the boundary nodes, the
larger the determinant becomes.

In section 2.3 the Lebesgue constants were pointed as an upper bound for the
interpolation error, it was also shown that an ill conditioned interpolation
matrix $\mathbf{M}$ would lead to very high Lebesgue constants. Starting from
these remarks it seems natural to minimize the Lebesgue functions (as defined
in equation (23)) by controlling the position of boundary nodes. Here follows
a brief discussion of the method adopted for this purpose.

When no polynomial augmentation is used, cardinal functions
$\psi_{i}(\boldsymbol{x})$ can be found as the solution of the following
linear system:

$\mathbf{M}^{T}\\{\psi_{i}(\boldsymbol{x})\\}_{1}^{m}=\\{\varphi_{i}(\boldsymbol{x})\\}_{1}^{m}$
(101)

where the same notation of equation (20) was used.

In some cases $\lambda_{I}$ and $\lambda_{B}$ exhibit large values at the
boundary when Neumann boundary conditions are employed, see for example Figure
3.

Willing to find the configuration of boundary nodes that minimizes the
Lebesgue constants $\Lambda_{I}$ and $\Lambda_{B}$, as defined in equation
(24), the cost function is defined as follows:

$\mathcal{F}(\boldsymbol{x}_{m_{J}},\dots,\boldsymbol{x}_{m})=\sum_{k=m_{J}}^{m}\lambda_{I}(\boldsymbol{x}_{k})=\sum_{k=m_{J}}^{m}\sum_{i=1}^{m_{I}}|\psi_{i}(\boldsymbol{x}_{k})|$
(102)

where the stencil nodes are enumerated as in section 2.2. In equation (102) we
considered only internal Lebesgue functions $\lambda_{I}$, however, it was
observed that optimizing with respect to the boundary Lebesgue functions
$\lambda_{B}$ leads to very similar optimal configurations.

By using a more compact vector notation, the sensitivities of the cost
function with respect to the position of a boundary node
$\boldsymbol{x}_{j}=\\{x_{j,1}\,,\,x_{j,2}\\}$ in 2D are:

$\frac{\partial\mathcal{F}}{\partial
x_{j,\eta}}=\sum_{k=m_{J}}^{m}\operatorname{sgn}\big{(}\boldsymbol{\psi}_{I}(\boldsymbol{x}_{k})\big{)}^{T}\frac{\partial\boldsymbol{\psi}(\boldsymbol{x}_{k})}{\partial
x_{j,\eta}}\\\ $ (103)

where $\boldsymbol{\psi}(\boldsymbol{x}_{k})$ is the column vector
$\\{\psi_{i}(\boldsymbol{x}_{k})\\}_{i=1}^{m}$ and
$\operatorname{sgn}(\boldsymbol{\psi}_{I}(\boldsymbol{x}_{k}))$ is the column
vector having the signs of $\boldsymbol{\psi}(\boldsymbol{x}_{k})$ in the
first $m_{I}$ entries, and zeros for the last $m_{B}$ indices corresponding to
the boundary nodes.

The term $\frac{\partial\boldsymbol{\psi}(\boldsymbol{x}_{k})}{\partial
x_{j,\eta}}$ in equation (103) can be obtained by deriving equation (101) as
follows:

$\frac{\partial\boldsymbol{\psi}(\boldsymbol{x}_{k})}{\partial
x_{j,\eta}}=(\mathbf{M}^{T})^{-1}\left(\frac{\partial\boldsymbol{\varphi}(\boldsymbol{x}_{k})}{\partial
x_{j,\eta}}-\frac{\partial\mathbf{M}^{T}}{\partial
x_{j,\eta}}\,\boldsymbol{\psi}(\boldsymbol{x}_{k})\right)$ (104)

where $\frac{\partial\mathbf{M}^{T}}{\partial x_{j,\eta}}$ is obtained by
taking the derivative of all the matrix entries. If we substitute equation
(104) into equation (103), the following term appears:

$\operatorname{sgn}\big{(}\boldsymbol{\psi}_{I}(\boldsymbol{x}_{k})\big{)}^{T}(\mathbf{M}^{T})^{-1}=\Big{(}\mathbf{M}^{-1}\operatorname{sgn}\big{(}\boldsymbol{\psi}_{I}(\boldsymbol{x}_{k})\big{)}\Big{)}^{T}=\boldsymbol{c}^{T}$
(105)

where $\boldsymbol{c}$ is therefore the column vector given as the solution of
the following adjoint equation:

$\mathbf{M}\,\boldsymbol{c}=\operatorname{sgn}\big{(}\boldsymbol{\psi}_{I}(\boldsymbol{x}_{k})\big{)}$
(106)

The computational cost for the calculation of the sensitivities of (103) then
becomes independent of the dimension of the problem and requires the solution
of $m_{B}$ linear systems like equation (106), one for each boundary node.

The sensitivities can then be used to perform a gradient based optimization of
the boundary nodes positions, this is done by taking the components orthogonal
to the normals.

Whenever two nodes ended up overlapping during the optimization, one of them
was removed, this is why in Figure 16 only $6$ boundary nodes appear.

###  C.2 Results

We have performed the optimization process described above by starting from
the reference stencil of Figure 1 and using MQ RBF with shape parameter
$\varepsilon$ satisfying $\varepsilon s=0.5$, where $s$ is the distance
between the nodes. Once the convergence is reached, we are left with the
situation depicted in Figure 16. We can see that different values of the angle
$\alpha$ correspond to different optimal placements of the boundary nodes.

$\alpha=-\pi/12$

$\alpha=0$

$\alpha=\pi/12$

Figure 16: optimal placement of boundary nodes obtained starting from the
reference stencil of Figure 1. Black solid lines are the actual normals, red
dotted lines are the optimal directions, asterisks are the position of
boundary nodes obtained by projecting the first layer of inner nodes along the
normals.

We remark that, even if the procedure described so far gives the best
placement for the boundary nodes, this comes at a relatively high
computational cost. Indeed, it is required to solve $m_{B}$ linear systems
like (101) in order to calculate the cardinal functions, in addition to
$m_{B}$ linear systems like (106) in order to estimate the sensitivities at
each step of the optimization. On top of that, a control on the normals, like
the one explained in section 4.1, is still needed when the boundary has high
curvature, see the angle between actual and optimal normals for the case
$\alpha=-\pi/12$ in Figure 16.

The results obtained here suggest that the interpolation can be improved by
placing the boundary nodes on the locations identified by projecting along the
normals the inner nodes of the first layer. Replacing standard boundary nodes
with the projected ones can be done after the node generation phase and
requires minimal additional computational cost on smooth and regular domains.

## References

  * [1] H. K. Versteeg, W. Malalasekera, An introduction to computational fluid dynamics: the finite volume method, Pearson education, 2007.
URL https://bit.ly/3OhalEr

  * [2] O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu, The finite element method: its basis and fundamentals, Elsevier, 2005.
URL https://bit.ly/3xtSC5S

  * [3] G.-R. Liu, Y.-T. Gu, An introduction to meshfree methods and their programming, Springer Science & Business Media, 2005.  doi:10.1007/1-4020-3468-7.
  * [4] H. Li, S. S. Mulay, Meshless methods and their numerical properties, CRC press, 2013\.  doi:10.1201/b14492.
  * [5] V. P. Nguyen, T. Rabczuk, S. Bordas, M. Duflot, Meshless methods: a review and computer implementation aspects, Math. Comput. Simulat. 79 (3) (2008) 763–813.  doi:10.1016/j.matcom.2008.01.003.
  * [6] V. G. Patel, N. V. Rachchh, Meshless method–Review on recent developments, Mater. Today-Proc. 26 (2020) 1598–1603.  doi:10.1016/j.matpr.2020.02.328.
  * [7] G.-R. Liu, Meshfree methods: moving beyond the finite element method, CRC press, 2009.  doi:10.1201/9781420082104.
  * [8] C. A. Duarte, J. Oden, A review of some meshless methods to solve partial differential equations, Texas Institute for Computational and Applied Mathematics Austin, TX, 1995.  doi:10.13140/RG.2.2.21361.48489.
  * [9] J. Slak, G. Kosec, Medusa: a C++ library for solving PDEs using strong form mesh-free methods, ACM T. Math. Software 47 (3) (2021) 1–25.  doi:10.1145/3450966.
  * [10] D. Miotti, R. Zamolo, E. Nobile, A fully meshless approach to the numerical simulation of heat conduction problems over arbitrary 3D geometries, Energies 14 (5) (2021) 1351.  doi:10.3390/en14051351.
  * [11] B. Fornberg, T. A. Driscoll, G. Wright, R. Charles, Observations on the behavior of radial basis function approximations near boundaries, Comput. Math. Appl. 43 (3-5) (2002) 473–490.  doi:10.1016/S0898-1221(01)00299-1.
  * [12] V. Bayona, An insight into RBF-FD approximations augmented with polynomials, Comput. Math. Appl. 77 (9) (2019) 2337–2353.  doi:10.1016/j.camwa.2018.12.029.
  * [13] N. Flyer, B. Fornberg, V. Bayona, G. A. Barnett, On the role of polynomials in RBF-FD approximations: I. Interpolation and accuracy, J. Comput. Phys. 321 (2016) 21–38.  doi:10.1016/j.jcp.2016.05.026.
  * [14] V. Bayona, N. Flyer, B. Fornberg, G. A. Barnett, On the role of polynomials in RBF-FD approximations: II. Numerical solution of elliptic PDEs, J. Comput. Phys. 332 (2017) 257–273.  doi:10.1016/j.jcp.2016.12.008.
  * [15] V. Bayona, N. Flyer, B. Fornberg, On the role of polynomials in RBF-FD approximations: III. Behavior near domain boundaries, J. Comput. Phys. 380 (2019) 378–399.  doi:10.1016/j.jcp.2018.12.013.
  * [16] B. Mavrič, Meshless modeling of thermo-mechanics of low-frequency electromagnetic direct chill casting, Ph.D. thesis, University of Nova Gorica, Graduate School (2017).
URL
https://repozitorij.ung.si/IzpisGradiva.php?id=3108&lang=eng&prip=rul:10911498:d1

  * [17] B. Mavrič, B. Šarler, Application of the RBF collocation method to transient coupled thermoelasticity, int. J. Numer. Method. H (2017).  doi:10.1108/HFF-03-2016-0110.
  * [18] B. Mavrič, B. Šarler, Local radial basis function collocation method for linear thermoelasticity in two dimensions, int. J. Numer. Method. H (2015).  doi:10.1108/HFF-11-2014-0359.
  * [19] R. Zamolo, L. Parussini, Analysis of geometric uncertainties in CFD problems solved by RBF-FD meshless method, J. Comput. Phys. 421 (2020) 109730.  doi:10.1016/j.jcp.2020.109730.
  * [20] J. Slak, Adaptive RBF-FD method, Ph.D. thesis, University of Ljubljana, Faculty of Mathematics and Physics (2020).
URL https://e6.ijs.si/~jslak/files/phd.pdf

  * [21] M. Jančič, J. Slak, G. Kosec, Monomial augmentation guidelines for RBF-FD from accuracy versus computational time perspective, J. Sci. Comput. 87 (4 2021).  doi:10.1007/s10915-020-01401-y.
  * [22] G. Liu, B. B. Kee, L. Chun, A stabilized least-squares radial point collocation method (LS-RPCM) for adaptive analysis, Comput. Method. Appl. M. 195 (37-40) (2006) 4843–4861.  doi:10.1016/j.cma.2005.11.015.
  * [23] I. Tominec, E. Larsson, A. Heryudono, A least squares radial basis function finite difference method with improved stability properties, SIAM J. Sci. Comput. 43 (2) (2021) A1441–A1471.  doi:10.1137/20M1320079.
  * [24] I. Tominec, E. Breznik, An unfitted RBF-FD method in a least-squares setting for elliptic PDEs on complex geometries, J. Comput. Phys. 436 (2021) 110283\.  doi:10.1016/j.jcp.2021.110283.