Energy-stable global radial basis function methods on summation-by-parts form

Radial basis function methods are powerful tools in numerical analysis and have demonstrated good properties in many different simulations. However, for time-dependent partial differential equations, only a few stability results are known. In particular, if boundary conditions are included, stability issues frequently occur. The question we address in this paper is how provable stability for RBF methods can be obtained. We develop and construct energy-stable radial basis function methods using the general framework of summation-by-parts operators often used in the Finite Difference and Finite Element communities.



page 1

page 2

page 3

page 4


Towards Stable Radial Basis Function Methods for Linear Advection Problems

In this work, we investigate (energy) stability of global radial basis f...

Stability estimates for radial basis function methods applied to time-dependent hyperbolic PDEs

We derive stability estimates for three commonly used radial basis funct...

Summation-by-parts operators for general function spaces

Summation-by-parts (SBP) operators are popular building blocks for syste...

Towards stability of radial basis function based cubature formulas

Cubature formulas (CFs) based on radial basis functions (RBFs) have beco...

An unfitted radial basis function generated finite difference method applied to thoracic diaphragm simulations

The thoracic diaphragm is the muscle that drives the respiratory cycle o...

Stabilizing Radial Basis Function Methods for Conservation Laws Using Weakly Enforced Boundary Conditions

It is well understood that boundary conditions (BCs) may cause global ra...

Conditions for Digit Stability in Iterative Methods Using the Redundant Number Representation

Iterative methods play an important role in science and engineering appl...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

We investigate energy stability of global radial basis function (RBF) methods for time-dependent partial differential equations (PDEs). Unlike finite differences (FD) or finite element (FE) methods, RBF schemes are mesh-free, making them flexible with respect to the geometry of the computational domain since the only used geometrical property is the pairwise distance between two centers. Further, they are suitable for problems with scattered data like in climate [flyer2012guide, lazzaro2002radial] or stock market [cuomo2020greeks, pettersson2008improved] simulations. Finally, for smooth solutions, one can reach spectral convergence [flyer2016role, fornberg2005accuracy]. In addition, they have recently become more and more popular for solving time-dependent problems in quantum mechanics, fluid dynamics, etc. [dehghan2017numerical, hesthaven2020rbf, iske2020ten, tominec2021residual]. One distinguishes between global RBF methods (Kansa’s methods) [kansa1990multiquadrics] and local RBF methods, such as the RBF generated finite difference (RBF-FD) [tolstykh2000using] and RBF partition of unity (RBF-PUM) [wendland2002fast] method. See the monograph [fornberg2015primer] and references therein.

Even though their efficiency and good performance have been demonstrated for various problems, only a few stability results are known for advection-dominated problems. For example, an eigenvalue analysis was performed for a linear advection equation in

[platte2006eigenvalue], and it was found that RBF discretizations often produced eigenvalues lending to an exponential increase of the norm when boundary conditions were introduced. To illustrate this, consider the following example (also found in [glaubitz2021stabilizing, Section 6.1]):


with , , and where periodic boundary conditions are applied. In this example, a bump is traveling to the right, leaving the domain and coming back to the left.

(a) Numerical solution at
(b) Energy profile developing in time
Figure 1: Gaussian kernel with points (equidistant points) after 10 periods

In Figure 1, we plot the numerical solution and its energy up to using a global RBF method with a Gaussian kernel and points. An increase of the size of the bump and of the energy can be seen. For longer times, the computation breaks down. The discrete setting does not reflect the continuous one with zero energy growth and demonstrates the stability problems.
To overcome those, it was shown in [glaubitz2021stabilizing, glaubitz2021towards] that a weak formulation could result in a stable method. Recently, estimates were obtained using an oversampling technique [tominec2021stability]. Both these efforts use special techniques, and the question we address in this paper is how to stabilize RBF methods in a general way.
Classical summation-by-parts (SBP) operators were introduced during the 1970s in the context of FD schemes and they allow for a systematic development of energy-stable semi-discretizations of well-posed initial-boundary-value problems (IBVPs) [fernandez2014review, svard2014review]. The SBP property is a discrete analog to integration by parts, and proofs from the continuous setting carry over directly to the discrete framework [nordstrom2017roadmap] if proper boundary procedures are added [svard2014review]. First based on polynomial approximations, the SBP theory has recently been extended to general function spaces developing so-called FSBP operators in [glaubitz2022nonpolnyomial]. Here, we investigate stability of global RBF methods through the lens of the FSBP theory.
We demonstrate that many existing RBF discretizations do not satisfy the FSBP property, which opens up for instabilities in these methods. Based on these findings, we show how RBF discretizations can be modified to obtain an SBP property. This then allows for a systematic development of energy-stable RBF methods. We give a couple of concrete examples including the most frequently used RBFs, where estimates are derived using an oversampling technique. For simplicity, we focus on the univariate setting for developing an SBP theory in the context of global RBF methods. That said, RBF methods and SBP operators can easily be extended to the multivariate setting, which is also demonstrated in our numerical tests.
The rest of this work is organized as follows. In Section 2, we provide some preliminaries on energy-stability of IBVPs and global RBF methods. Next, the concept of FSBP operators is shortly revisited in Section 3. We adapt the FSBP theory to RBF function spaces in Section 4. Here, it is also demonstrated that many existing RBF methods do not satisfy the SBP property and how to construct RBF operators in SBP form (RBFSBP). In Section 5, we give a couple of concrete examples of RBFSBP operators resulting in energy-stable methods. Finally, we provide numerical tests in Section 6 and concluding thoughts in Section 7.

2 Preliminaries

We now provide a few preliminaries on IBVPs and RBF methods.

2.1 Well-posedness and Energy Stability

Following [gustafsson1995time, nordstrom2017roadmap, svard2014review], we consider


where is the solution and is a differential operator with smooth coefficients. Further, and are operators defining the boundary conditions, is a forcing function, is the initial data, and denote the boundary data. Examples of Eq. 2 include the advection equation


with constant , the diffusion equation


with depending on , as well as combinations of Eqs. 4 and 3. Let us now formalize what we mean by the IBVP Eq. 2 being well-posed.

Definition 1

The IBVP Eq. 2 with is well-posed, if for every that vanishes in a neighborhood of , Eq. 2 has a unique smooth solution that satisfies


where are constants independent of . Moreover, the IBVP Eq. 2 is strongly well-posed, if it is well-posed and


holds, where the function is bounded for finite and independent of , and .

Switching to the discrete framework, our numerical approximation of Eq. 2 should be constructed in such a way that similar estimates to Eq. 5 and Eq. 6 are obtained. We denote our grid quantity (a measure of the grid size) by . In the context of RBF methods, denotes the maximum distance between two neighboring points. We henceforth denote by a discrete version of the -norm and represents a discrete boundary norm. Then, we define stability of the numerical solution as follows.

Definition 2

Let and be an adequate projection of the initial data which vanishes at the boundaries. The approximation is stable if


holds for all sufficiently small , where and are constants independent of . The approximated solution is called strongly energy stable if it is stable and


holds for all sufficiently small . The function is bounded for finite and independent of , and .

2.2 Discretization

To discretize the IBVP Eq. 2

, we apply the method of lines. The space discretization is done using a global RBF method resulting in a system of ordinary differential equations (ODEs):



denotes the vector of coefficients and

represents the spatial operator. We used the explicit strong stability preserving (SSP) Runge–Kutta (RK) method of third-order with three stages (SSPRK(3,3)) [shu1988total] for all subsequent numerical tests.

2.2.1 Radial Basis Function Interpolation

RBFs are powerful tools for interpolation and approximation

[wendland2004scattered, fasshauer2007meshfree, fornberg2015primer]. In the context of the present work, we are especially interested in RBF interpolants. Let be a scalar valued function and a set of interpolation points, referred to as centers. The RBF interpolant of is


Here, is the RBF (also called kernel) and is a basis for the space of polynomials up to degree , denoted by . Furthermore, the RBF interpolant Eq. 10 is uniquely determined by the conditions


Note that Eq. 11 and Eq. 12 can be reformulated as a linear system for the coefficient vectors and :


where and


Incorporating polynomial terms of degree up to in the RBF interpolant Eq. 10 is important for several reasons:

  1. The RBF interpolant Eq. 10 becomes exact for polynomials of degree up to , i. e., for .

  2. For some (conditionally positive) kernels , the RBF interpolant Eq. 10 only exists uniquely when polynomials up to a certain degree are incorporated.

In addition, we will show that (i) is needed for the RBF method to be conservative [glaubitz2021stabilizing, glaubitz2022nonpolnyomial]. The property (ii) is explained in more detail in Appendix A as well as in [fasshauer2007meshfree, Chapter 7] and [glaubitz2020shock, Chapter 3.1]. For simplicity and clarity, we will focus on the choices of RBFs listed in Table 1. More types of RBFs and their properties can be found in the monographs [wendland2004scattered, fasshauer2007meshfree, fornberg2015primer].

RBF parameter order
Gaussian 0

Polyharmonic splines (odd)

Polyharmonic splines (even)
Table 1: Some frequently used RBFs

Note that the set of all RBF interpolants Eq. 10 forms a -dimensional linear space, denoted by . This space is spanned by the cardinal functions


which are uniquely determined by the cardinal property


and condition Eq. 12. They also provide us with the following (nodal) representation of the RBF interpolant:


2.2.2 Radial Basis Function Methods

We outline the standard global RBF method for the IBVP Eq. 2. The domain on which we solve (2) is discretized using two point sets:

  • The nodal point set (centers) used for constructing the cardinal basis functions Eq. 15.

  • The grid (evaluation) point set for describing the IBVP (2), where .

By selecting , we get a collocation method, and with , a method using oversampling. The numerical solution is defined by the values of at and the operator by using the spatial derivative of the RBF interpolant , also at . The RBF discretization can be summarized in the following three steps:

  1. Determine the RBF interpolant .

  2. Define in the semidiscrete equation by inserting Eq. 17 into the continuous spatial operator. This yields

  3. Use a classical time integration scheme to evolve Eq. 9.

Global RBF methods come with several free parameters. These include the center and evaluation points and , the kernel , the degree of the polynomial term included in the RBF interpolant Eq. 10. The kernel might come with additional free parameters such as the shape parameter . Finally, we note that also the basis of the RBF approximation space , that one uses for numerically computing the RBF approximation and its derivatives, can influence how well-conditioned the RBF method is in practice. Discussions of appropriate choices for these parameters are filling whole books [wendland2004scattered, fasshauer2007meshfree, fornberg2011stable, fornberg2015primer] and are avoided here. In this work, we have a different point in mind and focus on the basic stability conditions of RBF methods.

3 Summation-by-parts Operators on General Function Spaces

SBP operators were developed to mimic the behavior of integration by parts in the continuous setting and provide a systematic way to build energy-stable semi-discrete approximations. First, constructed for an underlying polynomial approximation in space, the theory was recently extended to general function spaces in [glaubitz2022nonpolnyomial]. For completeness, we shortly review the extended framework of FSBP operators and repeat their basic properties. We consider the FSBP concept on the interval where the boundary points are included in the evaluation points . Using this framework, we give the following definition originally found in [glaubitz2022nonpolnyomial]:

Definition 3 (FSBP operators)

Let be a finite-dimensional function space. An operator is an -based SBP (FSBP) operator if

  1. for all ,

  2. is a symmetric positive definite matrix, and

  3. .

Here, and respectively denote the vector of the function values of and its derivative at the evaluation points .
Further, denotes the differentiation matrix and is a matrix defining a discrete norm. In order to produce an energy estimate, must be positive definite and symmetric. In this manuscript and in [glaubitz2022nonpolnyomial], we focus for stability reasons on diagonal norm FSBP operators [linders2020properties, gassner2016split, ranocha2017extended]. The matrix

is nearly skew-symmetric and can be seen as the stiffness matrix in context of FE. With these operators, integration-by-parts is mimicked discretely as:


for all .

3.1 Properties of FSBP Operators

In [glaubitz2022nonpolnyomial], the authors proved that the FSBP-SAT semi-discretization of the linear advection equation yields an energy stable semi-discretization. The so-called SAT term imposes the boundary condition weakly. Moreover, the underyling function space should contain constants in order to ensure conservations.

In context of RBF methods, constants have to be included in the RBF interpolants Eq. 10, also for the reasons discussed above.
We will extend the previous investigation to the linear advection-diffusion equation.


where is a constant and can depend on and . The problem (20) is strongly well-posed, as can be seen by the energy rate


with To translate this estimate to the discrete setting, we discretize (20). The most straightforward FSBP-SAT discretization reads


with and


We can prove the following result using the FSBP Definition 3.

Theorem 1

The scheme (22) is strongly stable with and .


We use the energy method together with the FSBP property to obtain


with . This is similar to the continuous estimate (21). Note that and have to be diagonal to ensure that defines a norm.

Clearly, the FSBP operators automatically reproduce the results from the continuous setting, similar to the classical SBP operators based on polynomial approximations [svard2014review]. Note that no details are assumed on the specific function space, grid or the underlying methods. The only factors of importance is that the FSBP property is fulfilled and that well posed boundary condition are used. In what follows, we will adapt the FSBP theory to radial basis functions.

4 SBP operators for RBFs

First, we adapt the FSBP theory in Section 2.2 to the RBF framework. Next, we investigate classical RBF methods concerning the FSBP property, and demonstrate that standard global RBF schemes does not fulfill this property. Finally, we describe how RBFSBP operators can be constructed that lead to stability.

4.1 RBF-based SBP operators

The function space for RBF methods is defined by the description in Subsection 2.2. Consider a set of points, . The set of all RBF interpolants Eq. 10 forms a -dimensional approximation space, which we denote by . Let be a basis in . Further, we have the grid points which include the boundaries. They are used to define the RBFSBP operators.

Definition 4 (RBF Summation-by-Parts Operators)

An operator is an RBFSBP operator on the grid points if

  1. for and ,

  2. is a symmetric positive definite matrix, and

  3. .

In the classical RBF discretizations, the exactness of the derivatives of the cardinal functions is the only condition which is imposed. However, to construct energy stable RBF methods, the existence of an adequate norm is as important as the condition on the derivative matrix. Hence it is often necessary to use a higher number of grid points than centers to ensure the existence of a positive quadrature formula to guarantee the conditions in Definition 4.
The norm matrix in Definition 4 has only been assumed to be symmetric positive definite. However, as mentioned above for the remainder of this work, we restrict ourselves to diagonal norm matrices where is the associated quadrature weight because Diagonal-norm operators are

  1. required for certain splitting techniques [gassner2016split, nordstrom2006conservative, offner2019error], and variable coefficients, see for example (24).

  2. better suited to conserve nonquadratic quantities for nonlinear stability [kitson2003skew],

  3. easier to extend to, for instance, curvilinear coordinates [chan2019efficient, ranocha2017extended, svard2004coordinate].

Remark 1

In Definition 4, we have two sets of points, the interpolation points and the grid points . The derivative matrix is constructed with respect to the exactness of the cardinal functions related to the interpolation points . However, all operators are constructed with respect to the grid points , i.e. . This is in particular essential when ensuring the existence of suitable norm matrix . This means that the size of the SBP operator is determine by the quadrature formula. So, the number of grid points and their placing highly effects the size of the operators and so the efficiency of the underlying method itself. In the future, this will be investigated in more detail.

4.2 Existing Collocation RBF Methods and the FSBP Property

In the classical collocation RBF approach, the centers intersect with the grid points, i.e. . It was shown in [glaubitz2022nonpolnyomial] that a diagonal-norm -exact SBP operator exists on the grid if and only if a positive and -exact quadrature formula exists on the same grid (the same requirement as for classical SBP operators). The differentiation matrix of a collocation RBF method can thus only satisfy the FSBP property if there exists a positive and -exact quadrature formula on the grid . The weights of such a quadrature formula would have to satisfy


with the coefficient matrix

and vector of moments

given by


In (26), is a basis of the function space . In many cases, the dimension of is larger than the dimension of . In this case, and the linear system in Eq. 25 is overdetermined and has no solution. This is demonstrated in Table 2, which reports on the residual and smallest element of the least squares solution (solution with minimal -error) of Eq. 25 for different cases. In all of our considered tests, the residuals were always larger than zero indicating that the operator is not in SBP form. Similar results are obtained for non-diagonal norm matrices , which is outlined in Appendix B.

Equidistant points

Halton points

Random points
Table 2: Residual and smallest elements for the cubic PHS-RBF on equidistant, Halton, and random points.

4.3 Existence and Construction of RBFSBP Operators

Translating the main result from [glaubitz2022nonpolnyomial], we need quadrature formulas to ensure the exact integration of . For RBF spaces, we use least-squares formulas, which can be used on almost arbitrary sets of grid points and to any degree of exactness. The least squares ansatz always leads to a positive and -exact quadrature formula as long a sufficiently large number of data points is used.

Remark 2

Existing results on positivity and exactness of least squares quadrature formulas usually assume that the function space contains constants [glaubitz2020stableQF, glaubitz2021constructing]. Translating this to our setting, we need this property to be fulfilled for . Therefore, should contain constants and linear functions. However, this assumption is primarily made for technical reasons and can be relaxed. Indeed, even when only contained constants, we were still able to construct positive and -exact least squares quadrature formulas in all our examples. Future work will provide a theoretical justification for this.

Due to the least-square ansatz, we may always assume that we have a positive and
-exact quadrature formula. With that ensured, we summarize the algorithm to construct a diagonal norm RBFSBP operators in the following steps:

  1. Build by setting the quadrature weights on the diagonal.

  2. Split into its known symmetric and unknown anti-symmetric part .

  3. Calculate by using

    and is defined analogous to where is a basis of the K-dimensional function space.

  4. Use in to calculate .

  5. gives the RBFSBP operator.

In the RBF context, one can always use cardinal functions as the basis. However, for simplicity reason is can be wise to use another basis representation, derived from the cardinal functions.

5 RBFSBP Operators

Next, we construct RBFSBP operators for a few frequently used kernels111The matlab code to replicate the results is provided in the corresponding repository We consider a set of points, , and assume that these include the boundaries and . Henceforth, we will consider the kernels listed in Table 1 and augment them with constants. The set of all RBF interpolants including constants Eq. 10 forms a -dimensional approximation space, which we denote by . This space is spanned by the cardinal functions which are uniquely determined by (16). The matching constraint is then simply That is,


with the approximation space having dimension .
The product space and its derivative space are respectively given by


Note that the right-hand sides of Eq. 28 and Eq. 29 both use elements to span the product space and its derivative space . However, these elements are not linearly independent and the dimensions of and are smaller than . Indeed, we can observe that and the dimension of Eq. 28 is therefore bounded from above by


Finally, we point out that in the calculation of the operators and below, we will round the numbers to the second decimal place.

5.1 RBFSBP Operators using Polyharmonic Splines

In the first test, we work with cubic polyharmonic splines, . On and for the centers , the three-dimensional cubic RBF approximation space Eq. 27 is given by with cardinal functions


and alternative basis functions222This basis can be constructed using a simple Gauss elimination method.

We make the transformation to the basis representation to simplify the determination of . In this alternative basis representation, the product space and its derivative space are respectively given by


Next, we have to find an -exact quadrature formula with positive weights. For the chosen equidistant grid points, the least-squares quadrature formula has positive weights and is -exact. The points and weights are and . The corresponding matrices and of the RBFSBP operator obtained from the construction procedure described before are


This example was presented with less details in [glaubitz2022nonpolnyomial].

5.2 RBFSBP Operators using Gaussian Kernels

Next, we consider the Gaussian kernel on for the centers . The three-dimensional Gaussian RBF approximation space Eq. 27 is given by with cardinal functions


Again for equidistant grid points in the least square quadrature formula, we obtain exactness and positive weights. They are and . The corresponding matrices and of the RBFSBP operator obtained from the construction procedure described before are