Radial Basis Functions Theory

Radial Basis Functions (RBF) use only a distance (typically Euclidean) when constructing the basis. For example, an interpolator for $u(\mathbf{x})$ is given by the linear combination of RBFs

$$u(\mathbf{x})=\sum _{i=1}^N{\alpha }_i\varphi (|\mathbf{x}-{\mathbf{x}}_i|)$$

where $\mid \cdot \mid$ is a norm (we will use Euclidean from here on) and so $|\mathbf{x}-{\mathbf{x}}_i|=r$ is the Euclidean distance (although any norm can be used) and $N$ is the number of data points.

There are several types of RBFs to choose from, some with a tunable shape parameter, $\epsilon$. Here are some popular ones:

Type	Function
Polyharmonic Spline	$\varphi (r)=r^n$ where $n=1,3,5,7,\cdots$
Multiquadric	$\varphi (r)=\sqrt{(r\epsilon {)}^2+1}$
Inverse Multiquadric	$\varphi (r)=1/\sqrt{(r\epsilon {)}^2+1}$
Gaussian	$\varphi (r)=e^{-(r\epsilon {)}^2}$

Augmenting with Monomials

The interpolant may be augmented with a polynomial:

$$u(\mathbf{x})=\sum _{i=1}^N{\alpha }_i\varphi (|\mathbf{x}-{\mathbf{x}}_i|)+\sum _{i=1}^{N_p}{\gamma }_i{p}_i(\mathbf{x})$$

where $N_p=\left(\begin{array}{c}m+d\\ m\end{array}\right)$ is the number of monomials ($m$ is the monomial order and $d$ is the dimension of $\mathbf{x}$) and $p_i(\mathbf{x})$ is the monomial term.

For instance, in 2D with $m=2$, we have $N_p=6$ and

$$p_i(\mathbf{x})\in \{1,x,y,x^2,xy,y^2\}$$

When we require the interpolation to be exact on a set of data points $\{{\mathbf{x}}_i,u({\mathbf{x}}_i)\}$ where $i=1,\cdots ,N$, we obtain the following linear system for the interpolation coefficients:

$$\left[\begin{array}{cc}\mathbf{A}& \mathbf{P}\\ {\mathbf{P}}^{\mathrm{T}}& 0\end{array}\right]\left[\begin{array}{c}\mathit{\alpha }\\ \mathit{\gamma }\end{array}\right]=\left[\begin{array}{c}\mathbf{u}\\ 0\end{array}\right]$$

where

$$\mathbf{A}=\left[\begin{array}{ccc}\varphi (|{\mathbf{x}}_1-{\mathbf{x}}_1|)& \cdots & \varphi (|{\mathbf{x}}_1-{\mathbf{x}}_N|)\\ ⋮& & ⋮\\ \varphi (|{\mathbf{x}}_N-{\mathbf{x}}_1|)& \cdots & \varphi (|{\mathbf{x}}_N-{\mathbf{x}}_N|)\end{array}\right]\mathbf{P}=\left[\begin{array}{ccc}{p}_1({\mathbf{x}}_1)& \cdots & p_N({\mathbf{x}}_1)\\ ⋮& & ⋮\\ p_1({\mathbf{x}}_N)& \cdots & p_N({\mathbf{x}}_N)\end{array}\right]$$

and $\mathbf{u}$ is the vector of dependent data points

$$\mathbf{u}=\left[\begin{array}{c}u({\mathbf{x}}_1)\\ ⋮\\ u({\mathbf{x}}_N)\end{array}\right]$$

and $\mathit{\alpha }$ and $\mathit{\gamma }$ are the interpolation coefficients. Note that the equations relating to ${\mathbf{P}}^{\mathrm{T}}$ are included to ensure optimal interpolation and unique solvability, given that conditionally positive radial functions are used and the nodes in the subdomain form a unisolvent set. See (Fasshauer, et al. - Meshfree Approximation Methods with Matlab) and (Wendland, et al. - Scattered Data Approximation).

Polynomial augmentation of the system has two benefits:

1. Increases accuracy, especially for polynomial fields and near boundaries.

2. Ensures the linear system has a unique solution for certain types of RBFs (conditionally positive definite).

See (Flyer, et al. - On the role of polynomials in RBF-FD approximations: I. Interpolation and accuracy) for more information on this.

Local Collocation

The traditional Kansa approach used in most RBF methods is based on constructing a unique interpolant for all the nodes in the domain. This involves coupling all nodes in the domain simultaneously and therefore makes it a global method. Such a global approach, while theoretically exact, scales poorly: the resulting dense system becomes prohibitively expensive and increasingly ill-conditioned as the number of nodes grows, particularly in 3D, due to the curse of dimensionality. Instead, RadialBasisFunctions.jl employs a local approach, where each node is influenced only by its $k$ nearest neighbors.

Boundary Conditions

Radial Basis Functions can also be used to solve PDEs with various types of boundary conditions (Dirichlet, Neumann, Robin, etc.), when this is done in a meshless context, special considerations must be made to ensure that all local systems, including those near boundaries, remain well-posed. This package enables different approaches for handling boundary conditions, among them, the Hermite approach is suggested when the PDE is to be solved only at interior nodes. Furthermore, it can also be used to interpolate data near the boundary while including boundary conditions where these are known.

Hermite Approach for Boundary Stencils

When a stencil is centered around an internal node but includes boundary nodes, standard RBF collocation can lead to ill-conditioning and singularity issues. This occurs because applying boundary operators $\mathcal{B}$ (such as normal derivatives for Neumann conditions) to the interpolation conditions breaks the symmetry of the local matrix $\mathbf{A}$.

The RBF-HFD (Hermite Finite Difference) method resolves this issue by modifying the basis functions only for stencils near boundaries. Instead of keeping the same basis regardless of boundary conditions and applying the operator to the interpolation conditions (which creates asymmetry), the Hermite approach modifies the basis itself.

For a stencil with $m_I$ internal nodes and $m_B$ boundary nodes ($m=m_I+m_B$), the approximate solution $u^h$ (the RBF interpolant) takes the form:

$$u^h({\mathbf{x}}_c)=\sum _{j=1}^{m_I}{\alpha }_j\varphi (|{\mathbf{x}}_c-{\mathbf{x}}_j|)+\sum _{j=m_I+1}^m{\alpha }_j{\mathcal{B}}_2\varphi (|{\mathbf{x}}_c-{\mathbf{x}}_j|)+\sum _{k=1}^{N_p}{\gamma }_k{p}_k({\mathbf{x}}_c)$$

where ${\mathbf{x}}_c$ is the stencil center (evaluation point), ${\alpha }_j$ are RBF coefficients, ${\gamma }_k$ are polynomial coefficients, and ${\mathcal{B}}_2$ denotes the boundary operator applied to the second argument of the kernel (i.e., to ${\mathbf{x}}_j$), and ${\mathcal{B}}_1$ denotes application to the first argument (i.e., to ${\mathbf{x}}_c$). The key insight is that the basis function is changed from $\varphi (\cdot ,{\mathbf{x}}_j)$ to ${\mathcal{B}}_2\varphi (\cdot ,{\mathbf{x}}_j)$ for boundary nodes.

The local system becomes:

$$\left[\begin{array}{ccc}{\mathbf{A}}_{I,I}& {\mathcal{B}}_2{\mathbf{A}}_{I,B}& {\mathbf{P}}_I\\ {\mathcal{B}}_1{\mathbf{A}}_{B,I}& {\mathcal{B}}_1{\mathcal{B}}_2{\mathbf{A}}_{B,B}& \mathcal{B}{\mathbf{P}}_B\\ {\mathbf{P}}_I^{\mathrm{T}}& (\mathcal{B}{\mathbf{P}}_B{)}^{\mathrm{T}}& 0\end{array}\right]\left[\begin{array}{c}{\mathit{\alpha }}_I\\ {\mathit{\alpha }}_B\\ \mathit{\gamma }\end{array}\right]=\left[\begin{array}{c}{\mathbf{u}}_I\\ \mathbf{g}\\ 0\end{array}\right]$$

where subscripts $I$ and $B$ denote internal and boundary quantities, respectively. The matrix blocks ${\mathbf{A}}_{I,I}$, ${\mathbf{A}}_{I,B}$, ${\mathbf{A}}_{B,I}$, and ${\mathbf{A}}_{B,B}$ represent RBF evaluations between internal-internal, internal-boundary, boundary-internal, and boundary-boundary nodes. The vectors ${\mathit{\alpha }}_I$ and ${\mathit{\alpha }}_B$ are the RBF coefficients for internal and boundary nodes, $\mathit{\gamma }$ are the polynomial coefficients, ${\mathbf{u}}_I$ contains function values at internal nodes, and $\mathbf{g}$ contains boundary condition values. This system is now symmetric and positive definite (for appropriate RBF kernels), ensuring unique solvability regardless of the boundary condition type.

We remark that the Hermite approach is only applied to stencils that include boundary nodes. For internal stencils far from boundaries, the standard RBF formulation remains unchanged, maintaining computational efficiency where boundary effects are not present.

Constructing an Operator

In the Radial Basis Function - Finite Difference method (RBF-FD), a stencil is built to approximate derivatives using the same neighborhoods/subdomains of $N$ points. <!– This is used in the [[MeshlessMultiphysics.jl]] package. –> For example, if $\mathcal{L}$ represents a linear differential operator, one can express the differentiation of the field variable $u$ at the center of the subdomain ${\mathbf{x}}_c$ in terms of some weights $\mathbf{w}$ and the field variable values on all the nodes within the subdomain as

$$\mathcal{L}u({\mathbf{x}}_c)=\sum _{i=1}^N{w}_{i}u({\mathbf{x}}_i)$$

We can find $\mathbf{w}$ by satisfying

$$\sum _{i=1}^N{w}_i{\varphi }_j({\mathbf{x}}_i)=\mathcal{L}{\varphi }_j({\mathbf{x}}_c)$$

for each basis function ${\varphi }_j$ (where ${\varphi }_j({\mathbf{x}}_i)=\varphi (|{\mathbf{x}}_i-{\mathbf{x}}_j|)$) and $j=1,\cdots ,N$, and if you wish to augment with monomials, we also must satisfy

$$\sum _{i=1}^{N_p}{\lambda }_i{p}_j({\mathbf{x}}_i)=\mathcal{L}{p}_j({\mathbf{x}}_c)$$

which leads to an overdetermined problem

$$\min (\frac{1}{2}\mathbf{w}{\mathbf{A}}^{⊺}\mathbf{w}-{\mathbf{w}}^{⊺}\mathcal{L}\varphi ),\text{ subject to }{\mathbf{P}}^{⊺}\mathbf{w}=\mathcal{L}\mathbf{p}$$

which is practically solved as a linear system for the weights $\mathbf{w}$ as

$$\left[\begin{array}{cc}\mathbf{A}& \mathbf{P}\\ {\mathbf{P}}^{\mathrm{T}}& 0\end{array}\right]\left[\begin{array}{c}\mathbf{w}\\ \mathit{\lambda }\end{array}\right]=\left[\begin{array}{c}\mathcal{L}\mathit{\varphi }\\ \mathcal{L}\mathbf{p}\end{array}\right]$$

where $\mathit{\lambda }$ are treated as Lagrange multipliers and are discarded after solving the linear system. The vectors are defined as

$$\mathcal{L}\mathit{\varphi }=\left[\begin{array}{c}\mathcal{L}\mathit{\varphi }(|{\mathbf{x}}_1-{\mathbf{x}}_c|)\\ ⋮\\ \mathcal{L}\mathit{\varphi }(|{\mathbf{x}}_N-{\mathbf{x}}_c|)\end{array}\right]\mathcal{L}\mathbf{p}=\left[\begin{array}{c}\mathcal{L}{p}_1({\mathbf{x}}_c)\\ ⋮\\ \mathcal{L}{p}_{N_p}({\mathbf{x}}_c)\end{array}\right]$$

where $\mathcal{L}\mathit{\varphi }$ is the vector of the operator applied to each RBF basis function evaluated at the stencil nodes, and $\mathcal{L}\mathbf{p}$ is the vector of the operator applied to each polynomial basis function.

Hermite Operator Construction

When constructing operators for stencils near boundaries using the Hermite approach, the system is modified to:

$$\left[\begin{array}{ccc}{\mathbf{A}}_{I,I}& {\mathcal{B}}_2{\mathbf{A}}_{I,B}& {\mathbf{P}}_I\\ {\mathcal{B}}_1{\mathbf{A}}_{B,I}& {\mathcal{B}}_1{\mathcal{B}}_2{\mathbf{A}}_{B,B}& \mathcal{B}{\mathbf{P}}_B\\ {\mathbf{P}}_I^{\mathrm{T}}& (\mathcal{B}{\mathbf{P}}_B{)}^{\mathrm{T}}& 0\end{array}\right]\left[\begin{array}{c}{\mathbf{w}}_I\\ {\mathbf{w}}_B\\ \mathit{\lambda }\end{array}\right]=\left[\begin{array}{c}{\mathcal{L}}_1\mathit{\varphi }({\mathbf{x}}_c,{\mathcal{X}}_I)\\ {\mathcal{L}}_1{\mathcal{B}}_2\mathit{\varphi }({\mathbf{x}}_c,{\mathcal{X}}_B)\\ \mathcal{L}\mathbf{p}({\mathbf{x}}_c)\end{array}\right]$$

where ${\mathcal{L}}_1$ denotes the differential operator applied to the first argument of the kernel, ${\mathcal{X}}_I$ is the set of internal nodes in the stencil, and ${\mathcal{X}}_B$ is the set of boundary nodes in the stencil. This yields weight vectors ${\mathbf{w}}_I$ and ${\mathbf{w}}_B$ that properly account for boundary conditions while maintaining symmetry. The global system assembly then proceeds as:

$$\mathcal{L}{u}^h({\mathbf{x}}_i)=\sum _{j\in {\mathcal{X}}_{i,I}}{w}_j({\mathbf{x}}_i)u({\mathbf{x}}_j)+\sum _{j\in {\mathcal{X}}_{i,B}}{w}_j({\mathbf{x}}_i)g({\mathbf{x}}_j)$$

where $g({\mathbf{x}}_j)$ represents the boundary condition values at boundary nodes.

Constructing an Operator Treating Boundary Nodes as Unknowns

In some applications, particularly multi-region or coupled problems, it may be advantageous to solve the governing equation at boundary nodes as well, treating all nodes (interior and boundary) as unknowns in the global system. When this strategy is adopted, an alternative implementation to the Hermite approach becomes available and the usual interpolation scheme can be preserved.

Rather than modifying the basis functions for boundary nodes, this approach maintains the standard RBF basis $\{\varphi (|\mathbf{x}-{\mathbf{x}}_j|)\}$ for all nodes regardless of their position. When treating boundary nodes as unknowns in the global system, the operator construction simplifies significantly. The key distinctions are:

• When the stencil includes boundary nodes but the evaluation point is interior: Apply the standard RBF-FD method unchanged. Boundary neighbors contribute as regular unknowns with no special treatment.

• When the evaluation point itself is on the boundary: Instead of modifying basis functions, modify the right-hand side of the local system. For a boundary point ${\mathbf{x}}_c$ (the stencil center) with operator $\mathcal{B}$, construct the RHS as $\mathcal{B}\mathit{\varphi }({\mathbf{x}}_c)$ and $\mathcal{B}\mathbf{p}({\mathbf{x}}_c)$ rather than using the standard differential operator.

This means the collocation matrix $\mathbf{A}$ always uses the standard kernel evaluation:

$$[\mathbf{A}{]}_{ij}=\varphi (|{\mathbf{x}}_i-{\mathbf{x}}_j|)$$

maintaining symmetry trivially. The local system for a boundary evaluation point becomes:

$$\left[\begin{array}{cc}\mathbf{A}& \mathbf{P}\\ {\mathbf{P}}^{\mathrm{T}}& 0\end{array}\right]\left[\begin{array}{c}\mathbf{w}\\ \mathit{\lambda }\end{array}\right]=\left[\begin{array}{c}\mathcal{B}\mathit{\varphi }({\mathbf{x}}_c)\\ \mathcal{B}\mathbf{p}({\mathbf{x}}_c)\end{array}\right]$$

This approach is significantly simpler than the Hermite method and stencil classification depends only on the evaluation point type, the same RBF basis is used everywhere, and interior stencils with boundary neighbors require no special treatment.

References

• Fasshauer, G. E., & McCourt, M. (2015). Kernel-based Approximation Methods using MATLAB. World Scientific. https://doi.org/10.1142/9335

• Flyer, N., Fornberg, B., Bayona, V., & Barnett, G. A. (2016). On the role of polynomials in RBF-FD approximations: I. Interpolation and accuracy. Journal of Computational Physics, 321, 21-38. https://doi.org/10.1016/j.jcp.2016.05.026

• Shankar, V., Wright, G. B., Kirby, R. M., & Fogelson, A. L. (2015). A radial basis function (RBF)-finite difference (FD) method for diffusion and reaction-diffusion equations on surfaces. Journal of Scientific Computing, 63(3), 745-768. https://doi.org/10.1007/s10915-014-9914-1

• Wendland, H. (2004). Scattered Data Approximation. Cambridge University Press. https://doi.org/10.1017/CBO9780511617539

• Wright, G. B., & Fornberg, B. (2006). Scattered node compact finite difference-type formulas generated from radial basis functions. Journal of Computational Physics, 212(1), 99-123. https://doi.org/10.1016/j.jcp.2005.05.030