Guidelines for RBF-FD Discretization: Numerical Experiments on the Interplay of a Multitude of Parameter Choices

There exist several discretization techniques for the numerical solution of partial differential equations. In addition to classical finite difference, finite element and finite volume techniques, a more recent approach employs radial basis functions to generate differentiation stencils on unstructured point sets. This approach, abbreviated by RBF-FD (radial basis function-finite difference), has gained in popularity since it enjoys several advantages: It is (relatively) straightforward, does not require a mesh and generalizes easily to higher spatial dimensions. However, its application is not quite as blackbox as it may appear at first sight. The computed solution might suffer severely from various sources of errors if RBF-FD parameters are not selected carefully. Through comprehensive numerical experiments, we study the influence of several of these parameters on the condition numbers of intermediate (local) weight matrices, on the condition number of the resulting (global) stiffness matrix and ultimately on the approximation error of the computed discrete solution to the partial differential equation. The parameters of investigation include the type of RBF (and its shape or other parameters if applicable), the degree of polynomial augmentation, the discretization stencil size, the underlying type of point set (structured/unstructured), and the total number of (interior and boundary) points to discretize the PDE, here chosen as a three-dimensional Poisson’s problem with Dirichlet boundary conditions. Numerical tests on a sphere as well as tests for the convection-diffusion equation are included in a supplement and demonstrate that the results obtained for the Laplace problem on a cube generalize to wider problem classes. The purpose of this paper is to provide a comprehensive survey on the various components of the basic algorithms for RBF-FD discretization and steer away from potential pitfalls such as computationally more expensive setups which not always lead to more accurate numerical solutions. We guide toward a compatible selection of the multitude of RBF-FD parameters in the basic version of RBF-FD. For many of its components we refer to the literature for more advanced versions.


Introduction
(Linear) partial differential equations can be solved numerically by a suitable discretization followed by the solution of a linear system of equations. Among possible discretization techniques, the radial basis function-finite difference (RBF-FD) method can be viewed as a generalization of the finite difference (FD) method to unstructured point sets. The earlier papers on RBF-FD include [72,73]. Since then, it has gained in popularity and entered into several applications including geosciences [25,26], heat flow [44], the financial sector [45], linear elasticity [66], fluid mechanics [11,36], and neuronal dynamics [52].
As in classical finite difference methods, RBF-FD replaces (linear) differential operators by differentiation stencils which encode formulas for weighted sums of function values at neighboring points. The weights are chosen such that the formulas become exact for a set of test functions which typically include radial as well as polynomial functions.
The discretization error in the numerical solution of a boundary value problem is affected by several aspects of the solution process including the underlying point set (the number of interior/boundary nodes and their distribution), the number (and distribution) of nodes in the discretization stencil as well as the RBF (including its shape or other parameters) and the degree of polynomial augmentation used to generate the stencil weights.
The application of RBF-FD discretization can lead to (auxiliary) highly ill-conditioned linear systems of equations to be solved for the weights, resulting in instability which prevents convergence of the numerical to the exact solution [24,25,39]. One possibility to address the ill-conditioning lies in a scaling of the shape parameter which, however, leads to a stagnation effect [1]. Another possibility lies in the use of hybrid kernels [46,47], another in the transition to stable versions of RBFs, e. g. RBF-GA [27,41,42], RBF-RA [31,80] or RBF-QR [39]. These, however, typically incur increased computational costs. A computationally even more expensive possibility to address the ill-conditioning is provided by (sufficiently high) extended precision [23]. Further possibilities to address the ill-conditioning are spatially variable shape parameters [8,29,56], Tikhonov regularization [69] or smoothing splines [20,82]. More recent papers often suggest to use polyharmonic splines (PHS) in combination with high order polynomials to generate the stencil weights. In [23,25], the Gaussian (GA) and PHS generating functions are studied and compared with respect to polynomial augmentation, extended precision and stable algorithms. In [58], PHS with polynomial augmentation (for a fixed stencil size and a numerically determined optimal degree of polynomial augmentation) is compared to the stable RBF-GA method with a small shape parameter (i. e., representing the limit case ε → 0). The advantages and disadvantages of the choice of parameters and the performance of the PHS generating function with polynomial augmentation are further studied in [2-4, 22, 35, 62].
In this paper, we will review some of these approaches with a focus on the most commonly used and (to us) most promising ones. We will discuss the multitude of involved parameter choices, illustrate their influence through numerical tests, and draw comparisons between the different approaches. There already exist quite a few numerical studies for the performance of RBF interpolation, approximation [7] and the solution of partial differential equations, but we address yet another angle that we have not found in the literature yet. Namely, the novel contributions of this paper consist of • A comprehensive view and comparison of the basic but often still competitive RBF-FD methods from the literature, in particular comprehensive numerical tests to illustrate the condition numbers of weight and stiffness matrices side-by-side with the discretization errors for a broad set of parameter choices across the landscape of RBF-FD discretization; the intention is to illustrate the ill-conditioning of weight matrices as a cause for numerical instabilities and the scaling of the shape parameter which is introduced to control this ill-conditioning as a cause for stagnation, show both the potential and limitations of individual RBF classes and offer a comparison to alternative RBF-FD setups (see Sects While in this paper all linear systems are solved directly, our focus on condition numbers of the resulting stiffness matrices will also be of interest for the subsequent development of iterative solvers for the systems.
The remainder of this paper is organized as follows. In Sect. 2, we review RBFs, RBF-FD approximation of linear differential operators and its application to the discretization of linear partial differential equations. The details pertaining to the node generation, the stencil formation and computation of the stencil weights, the choice of the RBF and the role of polynomial augmentation are discussed in Sect. 3. Section 4 is devoted to comprehensive numerical tests for a multitude of parameter choices in RBF-FD discretization applied to our model problem, the three-dimensional Poisson's equation with Dirichlet boundary conditions. In order to illustrate that qualitative results also apply to a wider class of model problems, we provide a supplement with numerical results on a different domain, i. e., a sphere, as well as for a different partial differential equation, i. e., the convection-diffusion equation with a recirculating convection. Finally, a conclusion and an outlook are presented in Sect. 5.

Definition 2.1
For a spatial dimension d ∈ N, an open domain ⊆ R d , a shape parameter ε > 0, a set of pairwise distinct nodes X := {x 1 , . . . , x n } ⊂ with n ∈ N, and a generating function φ ε : [0, ∞) → R, we define radial basis functions (RBFs) centered at a node x i as Sometimes, the generating function φ ε is referred to as a kernel function or even as a radial basis function itself. Ideally, the functions ε,x i are linearly independent and hence form a basis of their span so that the terminology basis function is justified. Furthermore, in order to later be able to guarantee the existence of unique solutions of (discrete) interpolation problems using radial basis functions, we will require the kernel function to be (conditionally) positive definite, see Definition 2.2. Typical examples for such kernel functions are listed in Table 1.
We list several reasons for the popularity of RBFs for multivariate interpolation and approximation [78]: • RBFs can be used in any spatial dimension.
• RBFs are applicable to arbitrarily scattered data. In particular, no mesh is required.
• Interpolants/approximants have a simple representation and inherit the smoothness of the radial basis functions.
RBFs can be used to approximate the action of a linear differential operator as follows. Let L denote a linear differential operator, u : → R be a sufficiently smooth function, and let be sets of distinct interior and boundary nodes, resp., with N I , N ∈ N, N I ≤ N (see Sect. 3.1 for more details about possibilities for the generation of these node sets). Such node sets X := X ∪ X ∂ also occur in finite difference discretizations, but now the nodes may be scattered and need not be connected by edges of a (structured tensor product) grid. For a node x j ∈ X , j ∈ {1, . . . , N I }, and stencil size n ≤ N , let X j ⊆ X be a stencil associated to the (stencil) center x j (see Sect. 3.2 for more information on the stencil computation). We use I j ⊆ {1, . . . , N } to denote the index set of nodes included in the stencil X j , i. e., The goal is to determine weights {w j 1 , . . . , w j n } ⊂ R such that the differential operator applied to the function u and evaluated at x j can be expressed (approximately) as a weighted sum of function values at the nodes in the stencil X j , i. e., (2. 2) The stencil weights are determined such that the weighted sum yields the correct result for the interpolant w. r. t. the stencil nodes, i. e., equality is enforced for (basis) functions k 1 , . . . , k n : → R. If these are cardinal basis functions w. r. t. the stencil nodes in X j , i. e., if there holds k (x s j i ) = δ i with the Kronecker delta, then p u : is the interpolant of u, and the application of a differential operator L at some x c ∈ is given by We will next illustrate the computation of these weights for typical RBF spaces with (or without) polynomial augmentation. Let { ε,x i | i ∈ I j } denote the n radial basis functions associated with the nodes of the stencil X j , and let denote the space of d-variate polynomials of degree at most ∈ N 0 . There holds and we denote a basis of by { p 1 , . . . , p M }. The n-dimensional constrained function space of RBFs associated with the stencil X j (2.1) and augmented by polynomials up to degree is defined by It is straightforward that (the weights of) an interpolant s ∈ R j of a function f : → R can be computed as the solution of the linear system with coefficient and right hand side vectors and system matrix (2.6) The interpolation problem has a unique solution if this matrix A j is regular. This can be guaranteed if the kernel function and node set X satisfy the properties given in Definitions 2.2 and 2.3, resp.

Definition 2.2 The generating function
for any (finite) set X = {x 1 , . . . , x N } ⊂ R n of distinct nodes and coefficients The generating function φ is called positive definite if (2.7) holds without any restrictions on λ j imposed through polynomials.
More details concerning (conditionally) positive definite functions are included in most textbooks on radial basis functions, see e. g. [9,77]. In terms of matrices, -unisolvency guarantees that the submatrix P j in (2.6) has full column rank, and if in addition the generating function is conditionally positive definite (of order + 1 or less), then A j is regular.
The cardinal functions of the space R j may now be represented by where e i ∈ R n denotes the i'th unit vector. Application of the differential operator and evaluation at the stencil center x c := x j yields for a single weight and for all weights with identity matrix I n ∈ R n×n and zero matrix 0 M×n ∈ R M×n . Augmentation of the matrix on the right to an (n + M) × (n + M) identity matrix and introduction of dummy variables w j i (which may later be discarded since they originated from the arbitrary augmentation to the identity matrix) leads to the linear system of equations for the desired stencil weights w j i in the approximation (2.2) [25]. Without polynomial augmentation, the saddle point system (2.8) simplifies to A j w j = L j . However, there exist several reasons to add polynomials in RBF-FD interpolation and approximation [23,25], including the following.
• For some RBFs, the matrix A j in (2.8) is guaranteed to be positive definite. For other RBFs, it can be shown that A j is positive definite on the subspace that satisfies the constraints for an appropriately chosen polynomial degree (this is e. g. the case for conditionally positive definite generating functions and a unisolvent stencil). • Exact/accurate interpolation of low order polynomials (in particular constants) instead of an oscillatory representation which is often obtained without polynomial augmentation. • The accuracy of the approximation in (2.2) may be significantly improved.
• Stagnation errors may be reduced.

RBF-FD for Partial Differential Equations
In the previous subsection, we derived an approximation of the application of a linear differential operator to a function and its evaluation at a fixed node by a sum of weighted function values at stencil nodes. We next extend this approach to the discretization of a linear partial differential equation where f : → R is a sufficiently smooth function and g : ∂ → R is a function to specify the Dirichlet boundary conditions. For each interior node x j ∈ X , j ∈ {1, . . . , N I }, of the point set X = X ∪ X ∂ ⊂ , we compute a stencil X j ⊂ X of size n (2.1) and stencil weights {w  (2.10) Given the Dirichlet boundary data from (2.9b), we set u s j substitute this boundary information into (2.10) and obtain i∈{1,...,n} Hence, we define the right hand side vector f := ( f 1 , . . . , f N I ) T ∈ R N I and the global stiffness matrix B ∈ R N I ×N I , containing row-wise the weights of the interior nodes of the j'th stencil, : else, yields the RBF-FD approximation u j ≈ u(x j ) to the solution of the partial differential equation (2.9a), (2.9b) at the interior nodes x j for j ∈ {1, . . . , N I }.

Sources for Errors in the RBF-FD Method
There are several sources that contribute to (numerical) approximation errors |u j − u(x j )|, j ∈ {1, . . . , N I } in the RBF-FD method where u ∈ R N I denotes the computed solution of (2.12). First of all, there occurs a discretization error when approximating the differential operator by a weighted sum of function values at the stencil nodes (2.2). The computation of the weights requires the solution of (small, possibly ill-conditioned) linear systems of equations subject to floating point arithmetic [33]. The subsequent computation of the RBF-FD approximate solution u requires the solution of another (now large and sparse) linear system of equations whose matrix coefficients consist of those weights, i. e., solutions of ill-conditioned systems. Errors in u can hence be attributed to errors in the stiffness matrix B as well as errors in the right hand side f in (2.12) which might result in amplified relative errors in u if the stiffness matrix B is ill-conditioned.
To illustrate the propagation of errors, we recall the well-known result from numerical linear algebra bounding the relative error between the solutionx of a perturbed linear system of equationsÂx =b and the solution x of Ax = b with respect to the relative errors in the matrix entries, A , and right hand side, b , if A · κ(A) < 1. If we assume for the weight matrices and right hand sides in (2.8) relative errors of order A j , b j ∈ O(10 −16 ) (e. g. rounding errors when computing in double precision accuracy), then the relative error in the computed weights is bounded by κ(A j ) · 10 −16 . Hence, we might lose up to log 10 κ(A j ) digits of accuracy in the weights which in turn define via (2.11) the matrix entries of the stiffness matrix B in (2.12), leading to a possibly large relative error B in these matrix entries. Hence, both the condition numbers of the weight matrices A j in (2.8) and the stiffness matrix B in (2.12) give insight into the numerical approximation errors observed in the RBF-FD solution of a partial differential equation.

Node and Stencil Generation, Choice of RBFs and Polynomial Augmentation
In this section, we discuss some methods for the node generation (Sect. 3.1), the stencil formation (Sect. 3.2), the choice of the RBF (Sect. 3.3) and the role of polynomial augmentation in RBF-FD (Sect. 3.4).

Generation of Node Distributions
The underlying application may already provide the set of (scattered) nodes. However, we may have the luxury to generate our own node set for a given domain and its boundary which is typically the case for the numerical solution of partial differential equations. We discuss several options to do so and refer to the literature for further sophisticated methods that may be necessary for more complex applications. The different ways to generate a node distribution can be roughly classified as regular/structured (e. g., Cartesian, Chebyshev or hexagonal) [46] or irregular/unstructured/scattered (e. g., coordinates generated by some type of random number generator). There are also variants that are neither structured nor (completely) unstructured such as node sets generated by a sequence of quasi-random numbers, or initially regular node distributions that have been subsequently perturbed [55]. Sophisticated algorithms to create node distributions with variable node density or strategies for local modification and refinement can be found in [19,25,48,51,63,67,75], and an overview of generation methods for node distributions, focusing on meshless finite difference methods including RBF-FD, is given in [70].
Two important quantities which are relevant for the RBF approximation error are the fill distance and separation distance of a point set. Definition 3.1 Given a node set X in a domain , the fill distance h X , (radius of largest circle not intersecting with X ) and the separation distance q X (shortest distance of two disjoint nodes) are defined by x − x j 2 and q X := 1 2 min In general, the (global) RBF approximation error depends on the fill distance (the smaller the better) as well as the separation distance (the larger the better for stability reasons) of the node set [67,74,78], hence these two aspects should be kept in mind for its construction.
We will use the following four prototypes for node distributions in our numerical tests in Sect. 4: regular Cartesian node distributions, a recently introduced node placing algorithm based on Poisson disk sampling which rejects nodes that are too close to each other (according to a user-defined minimal spacing) abbreviated as PNP [67] (created by [65]), Halton node distributions [21,32] (created by [10]) and irregular random node distributions (created by [53]).
The construction of these four types of node sets in the interior of our model domain (0, 1) d is straightforward. The extension of the Cartesian node set to the boundary ∂ is also straighforward whereas different variants exist in the case of PNP, Halton or random nodes, both with respect to the number and location of boundary nodes. One could e. g. use PNP/Halton/random nodes in a lower spatial dimension. We will use the same number of boundary nodes (as determined by the Cartesian grid) for all four types of node sets. The PNP and Halton interior node sets are complemented with boundary nodes of the Cartesian grid whereas the random interior node set is complemented with d − 1 dimensional random boundary nodes. Figure 1 illustrates these four types of node distributions for d = 2 spatial dimensions (even though all our numerical experiments will be performed for d = 3 spatial dimensions).

Generation of Stencils
There are several aspects to take into consideration when selecting the number and location of nodes in a stencil. If a partial differential equation is to be solved, then the stencil will be used to approximate partial/directional derivatives or more general differential operators. Since these describe local characteristics of a function, it is natural to include nodes that are close to the stencil center, e. g. the n nearest neighbors or all neighbors within a certain distance to the center.
Finite difference stencils often have a symmetric structure, i. e., node x i belongs to stencil X j if and only if node x j belongs to stencil X i (possibly even with identical weights w j i = w i j , e. g. for self-adjoint differential operators). This symmetry property typically no longer holds if upwinding is used for convection-type operators [64] or if a node lies near the boundary. In this case, one may use an extrapolation approach, introducing new nodes outside the closure , or asymmetric stencils [76]. For scattered node sets, we cannot expect any symmetry for RBF-FD stencils.
In the finite difference method, the stencil size n is typically fixed for all nodes. In the RBF-FD method, it is easily possible (and sometimes advantageous) to utilize different stencil sizes n j for different stencils X j , e. g., by selecting a fixed radius and choosing stencils that consist of all nodes within this radius with respect to the stencil center [46]. Further propositions and algorithms to form stencils can be found in [14,15,51,60] but are not pursued in more detail here. Figure 1 illustrates nearest neighbor stencils (marked in red) for stencil centers (marked in black) for the four types of node distributions. The remaining nodes are marked green if they lie on the boundary or blue if they lie in the interior.
If stencil nodes are very close to each other, the small separation distance often results in high condition numbers of the weight matrices. The criteria for good node distributions for (global) RBF approximations as mentioned in Sect. 3.1 carry over to the local RBF-FD approach, i. e., a good stencil X j should consist of a small fill distance (w. r. t. a suitable domain j containing x j ) and a large separation distance (see Definition 3.1) but does not necessarily have to be structured.
In the case of a structured grid, the weights and resulting truncation errors of RBF-FD may be computed analytically [5,6]. Stencils may be optimized with respect to these truncation errors, e. g. for the discretization of the Laplace operator on Cartesian (tensor product) nodes, we observed that so-called sparse stencils along the coordinate axes can be advantageous compared to nearest neighbor stencils. However, these advantages typically disappear (or are hard to exploit) on unstructured node sets. Stencils on structured grids more likely run into problems with unisolvency (Definition 2.3) than stencils on unstructured grids. A typical example of a node set that is not 2 -unisolvent in R 2 consists of any three nodes that are collinear, i. e., lying on a line. In [77, Theorem 2.7], a criterion is provided that guarantees for certain (subsets of Cartesian) node sets in R 2 to be -unisolvent, a generalization to higher spatial dimensions along with additional information on polynomial interpolation in several variables is given in [30,43]. Here we will pursue only the nearest neighbor approach (with a fixed stencil size for all nodes). For more advanced stencil selection algorithms (which may lead to different conclusions for the choice of RBF-FD parameters) and a machine learning approach to classify the quality of stencils we refer to [13,14,17,51,57].
If not only the nodes but an underlying, possibly unstructured mesh is available, then a nearest neighbor search could be performed on the graph with the Euclidean distance between nodes being replaced by the shortest path in the graph.
In the case of an unstructured node distribution without a mesh, nearest neighbor stencils may be determined efficiently using a kd-tree [81]. However, especially in the case of randomly distributed nodes, problems may occur if, as illustrated in Fig. 1 (right), the nearest neighbors are located in a somewhat imbalanced way around the stencil center, possibly with relatively large fill distance and small separation distance. Lop-sided stencils may also occur on structured grids for stencil centers close to the boundary ∂ . Table 1 shows some of the most commonly used generating functions which can be roughly classified as
Additional generating functions can be found in [9,25], including compactly supported generating functions which lead to a sparse matrix A j (2.8).
For the first class of generating functions, a major challenge in interpolation and approximation is the selection of a suitable shape parameter ε. Choosing it too large results in poor approximation accuracies since the RBFs either form peaks (GA, IMQ, IQ) or approximate piecewise polynomials (GMQ(ν) with ν > 0) [18]. On the other side, as the shape parameter ε goes to zero, the RBFs converge to a constant function and the matrix A j in (2.8) becomes increasingly ill-conditioned, rendering a solution in finite floating point arithmetic numerically infeasible. However, the underlying interpolants/approximants may still exist in the limit and be determined analytically in which case they are related to multivariate polynomials [28,59] (which also explains why polynomial augmentation is not needed for small shape parameters [39]). E. g., for d = 1, the Lagrange minimal-degree interpolating polynomial can be obtained for the underlying interpolant/approximant in the limit ε → 0. Hence, the FD method can be seen as the limit ε → 0 in the RBF-FD method [5,6,18]. On the other hand, this explains the occurrence of the Runge phenomenon (i. e., large boundary oscillations in interpolations with higher degree polynomials) in (global) RBF interpolation with small shape parameters. In the context of RBFs, the Runge phenomenon can also be observed if the separation distance becomes very small [29]. For each node set X and stencil X j , there typically exists a shape parameter ε j > 0 that is optimal for the numerical computation of the corresponding stencil weights (2.8). Such an ε j depends on several factors such as the function to be approximated, the generating function, the stencil nodes and even the machine precision [38]. Hence, using spatially varying (i. e., stencil dependent) shape parameters ε j can be (theoretically) advantageous [8,25,29], especially for functions with gradient discontinuities [56], but it is not pursued any further here since it adds even more parameters/choices into the setting of RBF-FD. On Cartesian grids, optimal shape parameters can be derived analytically [5,6].
Another way to address the problem of ill-conditioning is using the hybrid kernel HYB(γ ) [46,47]. The idea is to achieve numerical stability (i. e., lower condition numbers of the weight and stiffness matrices as well as no large spurious eigenvalues of the stiffness matrix for smaller shape parameters) in a direct approach, i. e., with lower computational costs than a stable RBF algorithm. HYB(γ ) typically leads to smaller errors and larger condition numbers by decreasing γ and vice versa (i. e., a too small as well as a too large γ value can be problematic in terms of stability or accuracy, respectively). Hence, often is utilized [46,52]. Yet another (computationally more expensive) possibility to address the ill-conditioning is to switch to stable versions of RBFs, e. g., RBF-GA [27,41,42], RBF-RA [31,80] or RBF-QR [39]. These stable algorithms typically increase the computational costs by at least a factor of 10 [61]. Using (sufficiently high) extended precision is another approach, but in general not advisable for large scale problems since it is computationally even more expensive [23,39].

Polynomial Augmentation in RBF-FD
Polynomial interpolation can be viewed as a special case of RBF-FD with polynomial augmentation in which the dimension of the polynomial space M equals the stencil size n, but no radial basis functions are used, i. e., the stencil weights are computed as the solution of P T j w j = Lp j (2.8). As mentioned in Sect. 2.1, polynomial interpolation can lead to a singular P j if min{d, n} > 1. If polynomial augmentation is used in RBF-FD, necessary (but not sufficient) criteria for the weight matrix A j in (2.8) to be regular are given by n ≥ M and rank(P j ) = M (full column rank).
Often, for a numerically stable computation, a lower bound n min = n min ( ) is suggested for the stencil size n that depends on the degree of polynomial augmentation (or, equivalently, the dimension M of the respective polynomial space (2.4)). In [23], the lower bound n min = M is suggested for RBF-FD interpolation and approximation. For the solution of elliptic and hyperbolic PDEs, n min ≈ 2M has been proposed [4,22,35]. A slightly different bound, n min ≈ 2M + ln(2M) is proposed in [62] for the advection-diffusion equation. In [22,58], it is observed that unstructured node distributions allow for smaller n min than Cartesian node distributions with an increasing difference as increases (e. g., n min ≈ 2M is too small for Cartesian nodes with > 4 whereas unstructured nodes can achieve lowest errors for n ≈ M + ln(2M) [58]. In general, the optimal relation between polynomial degree and stencil size n depends on the underlying problem [62].
There are different effects of polynomial augmentation in the RBF-FD method which depend on the type of generating function that they complement. Polynomial augmentation (with close too maximal permissible degree ) is studied in [23,25] for RBF-FD interpolation and approximation for the GA generating function (with a shape parameter ε, leading to infinitely smooth RBFs) and the PHS(k) generating function (without a shape parameter, leading to not infinitely smooth RBFs). In both cases, stagnation errors disappear and higher accuracies (for the same convergence order) compared to polynomial least squares approximation with the same are observed (for sufficiently large ). The convergence order for if the approximated function is sufficiently smooth [63,74] (i. e., D = 0 for function interpolation/approximation and D = 2 for the Laplace operator L = or the convection-diffusion operator L = + b T ∇) and (approximate) internodal spacing N −1/d [2,23,39]. Hence, the convergence order depends on (which is connected to n as a larger requires a larger n). While n, ε and k do not enter into the above order of convergence, varying these parameters can still lead to significant changes in the accuracy [23], i. e., affect the constants hidden in the O notation. For GA, the accuracy can be significantly increased by polynomial augmentation for shape parameters larger than the optimal one [27] (while possibly increasing the condition number of A j in (2.8)) and the optimal shape parameter is (similar to the case without polynomial augmentation) mostly invariant with respect to refinement (i. e., the number of nodes N ).
For PHS (or other only conditionally positive definite RBFs), polynomial augmentation is typically needed for convergence [35,46]. The condition number of the matrix A j in (2.8) may increase (also for PHS) when increasing the polynomial degree . Furthermore, when increasing the stencil size n, the stencil weights become larger (by absolute value) near the stencil center and smaller elsewhere, and if the stencil size n is large enough, then there is no Runge phenomenon [2,3].
Infinitely smooth RBFs have the potential for spectral accuracy in global RBF approaches [29], but this benefit is typically sacrificed in the local RBF-FD approach [5] since the stencil size n is significantly smaller than the number of nodes N . Additionally, GA and PHS with polynomial augmentation (and n only slightly larger than M) can achieve (especially for unstructured nodes since they allow for smaller n min than Cartesian nodes [58]) comparable accuracies with lower computational costs than stable RBF variants. Using the GA generating function (with a stable algorithm or close to maximal permissible degree of polynomial augmentation) leads to the challenge of choosing an appropriate shape parameter, and using extended precision is computationally too expensive for RBF-FD approximation. Hence, the PHS generating function with close to maximal permissible degree of polynomial augmentation is recommended (taking into consideration the approximation accuracy and computational cost) [23,25].
Many results for RBF-FD approximation using PHS with polynomial augmentation carry over to the solution of elliptic and hyperbolic PDEs, i. e., stagnation errors disappear, the convergence order is again given by (3.2), hence depends on the polynomial degree (which requires a minimal stencil size n min ). Increasing the PHS degree determined by k and stencil size n may lead to slight improvements of accuracy. The stencil weights become larger (by absolute value) near the stencil center and smaller elsewhere by increasing the stencil size n. If the stencil size n is large enough, then there is no Runge phenomenon [4,22].
For the Poisson problem (involving a differential operator L of order D = 2 (3.2)), a convergence order of O(N − /d ) is numerically observed for ∈ {2, 4, 6, 8} (and no convergence for = 0) [35], i. e., slightly higher than the (theoretical) convergence order for RBF-FD interpolation and approximation (3.2). Additionally, it is observed in [4,22] that using PHS in the form r 2k−1 or r 2k log(r ) leads to similar results (thus the form r 2k−1 is preferred), quasi-uniform node distributions lead to similar (often even higher) accuracies compared to Cartesian node distributions, and if the stencil size n is sufficiently large, then no ghost nodes (i. e., additional nodes outside the domain to increase stability and accuracy by avoiding lop-sided stencils near the boundary ∂ ) are needed. It has been observed in [74] that ghost nodes can increase stability and accuracy especially in the presence of Neumann boundary conditions. Moreover, whereas a larger PHS degree can lead to a higher accuracy (accompanied by the potential for more numerical instabilities such as larger spurious eigenvalues of the stiffness matrix and larger absolute values of the stencil weights), a PHS degree that is chosen too small can also be problematic. Hence, given the degree of augmentation ∈ N, in [62] it is suggested to select k in the PHS r 2k−1 according to if is even.
While the optimal choice for the relation between PHS exponent k and polynomial degree could depend on the type of the node distribution, the relation k ≤ k max := + 1 ( 3 . 4 ) has to be satisfied to prove unisolvency in this setting [62]. Furthermore, there has to hold for the order D of the differential operator (3.2) to ensure that the right hand side L ε,x if is even (3.6) and propose to select k ∈ [k min , k max ] depending on the model problem and number of points N (see also (4.1) for our choice in the numerical experiments). A different idea lies in the use of variable degrees of polynomial augmentation j and variable sizes n j for each stencil X j [49]. While it can reduce the computational costs and increase the accuracy, it adds even more parameters/choices into the setup of RBF-FD, hence this approach is not pursued any further here.

Numerical Results
This section contains numerical tests for the RBF-FD discretization of Poisson's equation in d = 3 spatial dimensions on a unit cube, i. e., in (2.9a), (2.9b), we set = (0, 1) 3 and L = . We refer to the supplement for numerical tests for further model problems. We choose the right hand side functions f , g such that the analytical solution u is given by the three-dimensional version of Franke's function F : R 3 → R [40,54] with We will study the influence of the following parameters on the performance of the RBF-FD method to illustrate causes for potential pitfalls and hopefully steer away from them: • The type of the node distribution (i. e., Cartesian, PNP, Halton or random); • The stencil size n; • The total number of nodes N ; • The generating function φ ε and its intrinsic parameters (shape parameter ε, the parameter γ in HYB(γ ), the PHS degree specified by k); • the degree of polynomial augmentation .
We will illustrate the resulting approximation errors of the discrete approximate solution side-by-side with the condition numbers of the involved weight and stiffness matrices. To this end, all of the following Figs

General Information
In order to compare results for the different types of node distributions (see Sect. 3.1), we generate PNP, Halton and random node distributions with the same number of interior and boundary nodes as in a typical Cartesian grid, i. e., N I = ν d interior nodes and N = (ν + 2) d total nodes for some ν ∈ N. While the construction of Halton or random nodes with a prescribed number N I of interior nodes is straighforward, it is slightly more complicated for PNP nodes since their construction is based on a local distance function h : → (0, ∞) that specifies the minimal nodal spacing at point x ∈ via h(x). A constant nodal spacing function h ≡ (ν +1) −1 (i. e., the same spacing as for the Cartesian nodes) leads to fewer nodes than N I . Hence, we use a scaled version h ≡ α(ν + 1) −1 with α ∈ [0.9256, 0.9284] experimentally determined such that exactly N I nodes are generated. (The tests with h ≡ (ν + 1) −1 lead to qualitatively comparable observations as h ≡ α(ν + 1) −1 with slightly increased errors and decreased condition numbers due to the typically larger fill and separation distances).
When showing results on Cartesian grids, for comparison we also include test results (for the error (4.2) and the condition number of the stiffness matrix) obtained for finite difference discretizations with a standard central 7-point stencil (i. e., second-order approximation) as well as for a compact/Hermite 19-point stencil (i. e., the 19 stencil nodes consist of the 3 × 3 × 3 nodes around the center without the 8 corner nodes, and not only function values but also derivatives are included in the approximation in (2.2), leading to a fourth-order approximation, see [68,79]).
For the computation of stencils (see Sect. 3.2), we use a library [50] that utilizes a kd-tree for the nearest neighbor search with respect to the norm · μ for some μ ∈ N. In our experiments, we use the Euclidean norm by setting μ = 2.
We focus on the generating functions in Table 1, i. e., we consider the generating functions GA, MQ, IMQ, IQ, HYB(γ ) and PHS(k) of the form r 2k−1 with k ≥ 2 determined by for polynomial augmentation ≥ 2 (which is required for convergence, see (3.2)). Even though some generating functions (e. g., GA, IMQ, IQ [25]) guarantee the resulting matrices A j in (2.8) to be symmetric positive definite, they might be highly ill-conditioned so that a direct solver based on a Cholesky factorization may break down in practice. Hence, for solving linear systems, we will always use direct solvers based on LU factorizations with pivoting. For stencils that are not unisolvent, the matrix block P j in (2.8) may no longer have full column rank, hence A j is no longer regular but (2.8) may still have (no longer unique) solutions. These can be computed by a nullspace approach [12,15] which may also be beneficial in the unisolvent case. For the computation of (lower bounds of) condition numbers κ 1 of the arising matrices, we use the LAPACK [37] function dgecon with the 1norm (we also computed κ ∞ which led to qualitatively similar observations for the stiffness matrix and made no difference for the symmetric weight matrices, hence these results are not reported here).
In order to measure the approximation accuracy of the RBF-FD computed discrete solution, we determine the relative 2-norm error between the exact solution of the PDE We have also performed numerical experiments and measured the error in the (relative) infinity as well as the (relative) 1-norm but do not report these results here since they are qualitatively very similar (see also [19,35] for results in different error norms).
Since the error e 2 in (2.12) does not provide insight into the component-wise absolute or relative errors |(u s (·) − u(·))| or |(u s (·) − u(·))/u s (·)|, resp., we include Fig. 2 which illustrates these pointwise errors with respect to their distance to the center (0.5, 0.5, 0.5) T of the domain [0, 1] 3 for three different types of generating functions on a Halton grid. While we observe that the pointwise errors possess a significant deviation from their mean and tend to become smaller closer to the (Dirichlet) boundary, we argue that the error e 2 (added as a horizontal line to all plots, scaled by 1/ √ N I for the bottom row) is a reasonable measure to assess the overall RBF-FD discretization error.

Numerical Results for (Infinitely Smooth) Gaussian RBFs
In this subsection, we provide numerical tests using the GA generating function which leads to infinitely smooth RBFs. We have also performed the same tests for generalized multiquadric GMQ(ν) for ν ∈ {1, −1, −2} (i. e., for MQ, IMQ and IQ) generating functions but do not report on these results since they were qualitatively similar (with adjusted shape parameters and constants c in (4.3)).
Similarly, since random nodes lead to qualitatively comparable observations as Halton nodes (with slightly increased errors and condition numbers due to the typically larger fill distance and smaller separation distance), we also do not include our test results for random nodes here.
We will later (in Sect. 4.4) provide numerical results comparing several generating functions, including MQ, IMQ and IQ, on all four types of grids (Cartesian, PNP, Halton, random).
We begin with an overview of the numerical tests before a more detailed discussion of the results. In all of these tests, we show plots for the RBF-FD discretization error (4.2), the • Figure 3: Influence of the shape parameter ε on different grid types (Cartesian, PNP, Halton) for different stencil sizes n, fixed problem size N . • Figure 4: Influence of the stencil size n on different grid types (Cartesian, PNP, Halton) for different shape parameters ε, fixed problem size N , polynomial augmentation of degree = 0 (a constant term). • Figure 5: Influence of the problem size N on different grid types (Cartesian, PNP, Halton) for different stencil sizes n, a problem size dependent shape parameter ε (4.3), polynomial augmentation of degree = 0.
In the three left plots, we see that there exists an optimal shape parameter ε ∈ (1, 4) that increases with increasing stencil size, especially on the Cartesian and PNP grids. For (much) larger shape parameters, the discretization error increases since the basis functions turn more and more into spikes so that the spanned function space does not offer accurate approximations. For (much) smaller shape parameters, numerical instabilities occur which become apparent in the middle and right plots showing condition numbers for both the weight and stiffness matrices. Larger stencil sizes incur larger condition numbers in the weight matrices but also have the potential to reach smaller discretization accuracies for suitably chosen shape parameters. In particular, the (best possible) RBF-FD discretization accuracies obtained on the PNP and Halton grids are comparable to the ones on the Cartesian grid (or even better) as long as the stencil size is large enough.
In Fig. 4, the setting is similar to the one in Fig. 3, only that now the stencil size n ∈ {7, . . . , 200} is shown along the x-axis and each color/marker represents a different shape parameter ε ∈ {1, 2, 3, 4, 5}. Furthermore, polynomial augmentation of degree = 0 (i. e., a constant term) has been added. This choice typically increases the accuracy slightly with only marginal increase in the computational cost (see Sect. 3.4 for more details). This figure once more illustrates the numerical instabilities that occur for (too) small shape parameters. On the Cartesian grid, one can observe a noticeable jump that occurs when increasing the stencil size from n ≤ 32 to n ≥ 33 which was also observed in [5] and coincides with the following observation: The stencil nodes that are most valuable for the approximation of the Laplace operator on the Cartesian grid appear to be those along the coordinate axes. The first 27 nodes in a nearest neighbor search are those in the 3 × 3 × 3 grid around the stencil center, the next 6 are the extensions to both sides along the respective three coordinate axes. When all 6 are added and a stencil of 33 nodes is formed, one sees the jumping improvement in the approximation accuracy. On the Cartesian grid, a further increase in the stencil size leads to hardly any improvement. On the Halton grid, stagnation sets in for a stencil size around n = 100 whereas the error continues to (slowly) decrease up to a stencil size around n = 150.
In Fig. 5, we show plots with respect to 3 √ N where N ∈ {(ν + 2) 3 : ν ∈ N, 6 ≤ ν ≤ 25} denotes the number of total nodes with N I = ν 3 ∈ [216, 15, 625] interior nodes. As before, we use the GA generating function with polynomial augmentation of degree = 0. For the shape parameter, we follow the strategy to increase it along with the problem size as in for a constant c to avoid increasing ill-conditioning in the weight matrices [22,24,46,71]. The constant c influences the shape parameter, a smaller c leads to larger condition numbers of the weight (stiffness) matrices, a larger c leads to "spiky" basis functions that no longer provide good approximations on finer grids. A similar approach using two constants c 1 , c 2 with ε = c 1 N 1/d − c 2 > 0 (and the determination of these constants) has been studied in [24]. We have performed several tests for constants c ∈ {0.05, 0.1, 0.15, 0.2, 0.25} (not reported here) and found c = 0.15 to yield the best results for a stencil size n = 50 (and good results for a wider range of stencil sizes). In the plots in the middle of Fig. 5, using (4.3) to determine the shape parameter, it is confirmed that with such a choice of shape parameter, the condition numbers of the weight matrices remain (close to) constant (with different constants for different c) as N increases while the condition numbers of the stiffness matrix still increase with N .
On a Cartesian grid, the RBF-FD error for a stencil of size n = 7 is comparable to the 7-point-FD error (even better up to 3 √ N ≈ 25 and worse for larger N since c = 0.15 is too large for this n), and stencil sizes n > 33 appear to lead to no further advantage compared to n = 33. On PNP and Halton grids, RBF-FD with stencil size n = 7 performs worse than on a Cartesian grid. For larger stencils, the error decreases both with increasing stencil and problem sizes n, N , reaching accuracies even better than the finite difference 19-point-stencil.

Numerical Results for Polyharmonic Spline (PHS) RBFs
An attractive feature of the PHS generating function with polynomial augmentation lies in their simplicity due to the lack of a shape parameter ε to be finetuned and their scalability. In fact, it can be shown [16,34] that the weight vector w j in the solution of the system (2.8) can be recovered from the solution w j h of a scaled (and typically much better conditioned) system where the scaling parameter h ∈ R is selected such that the maximum distance from the center x j to another stencil node is one. Since the new center node after the shift is the origin, the evaluation of Lp h 0 greatly simplifies in the case of a monomial basis of : For the Laplace operator, the entries corresponding to x 2 and y 2 are 2, all others are zero. Since the high condition numbers in (2.8) hence can be reduced by a (row and column) scaling, they are harmless in actual computations and have no negative effect on the computed interpolant. This has also been observed in Sect. 5.2 of [23]. We will illustrate the (lack of) difference between the numerical solutions with and without shift-and-scale in Fig. 7 (last two rows of plots) and use (4.4) instead of (2.8) in all other numerical tests involving PHS.
However, even without a shape parameter, there still remain several parameters to choose such as the stencil size n, the exponent in the polyharmonic spline determined by k and the degree of polynomial augmentation. Similarly as in the previous subsection, PNP leads to results somewhat between those for Cartesian and Halton nodes while random nodes lead to worse, but qualitatively comparable observations as Halton nodes. Hence, we do not include test results for PNP and random nodes (Sect. 4.4 will include results for PHS on PNP and random nodes in comparison to other RBFs). In this subsection, we provide the following numerical results to illustrate the interplay of these quantities and their resulting discretization errors.
• Figure 6: Influence of the stencil size n on a Halton grid for the PHS(k) generating function and different degrees of polynomial augmentation where k is determined by via (3.6) or (4.1). • Figure 7: Influence of the problem size N on different grid types (Cartesian, Halton) for different k in PHS(k) and degrees of polynomial augmentation which in turn determine the stencil size n via (4.5).  Table 2).
The comparison between plots in the two rows demonstrates that larger PHS degrees can decrease the error but may eventually lead to numerical instabilities. Larger degrees than given by (3.6) seem to be especially beneficial for smaller polynomial degrees ≤ 4 while (4.1) can lead to numerical instabilities for ≥ 5.
We do not show plots for our results on Cartesian grids but only briefly summarize our findings: For Cartesian nodes and ∈ {2, 3, 5}, there occured breakdowns in the solver (i. e., in the LU factorization of the weight matrices) for some of the smaller stencil sizes. The nullspace approach of [12,15] could have been used instead of an LU factorization to compute weights but was not pursued here. The second column in Table 2 summarizes the minimal stencil sizes n min which were necessary in these setups.
As discussed in Sect. 3.4, using a stencil size of approximately n ≈ 2M seems reasonable for Cartesian nodes with ≤ 3 or unstructured nodes. Further tests on unstructured grids (not included here) have shown an only slight advantage when using n = 2M + ln(2M) compared to n = 2M for the error (4.2) and the condition numbers of the weight and stiffness matrices, especially in the case of larger degrees of polynomial augmentation (or on random nodes). Since the difference between n = 2M and n = 2M + ln(2M) (or some similar stencil size) is typically rather small (as long as is sufficiently large (3.2)), we recommend (and will use ourselves in subsequent tests) the following stencil sizes, A larger stencil size not only increases the computational cost to compute the stencil weights but also leads to more non-zero entries in the stiffness matrix. We hence conclude that RBF-FD using PHS and polynomial augmentation with a large degree is not suitable/efficient for Cartesian nodes since the stencil size should increase with the degree of polynomial augmentation as for example suggested in (4.5). Compared to the search for a best shape parameter in infinitely smooth RBFs, the sensitivity of the results when moving away from the optimal setting is much reduced here.
In Fig. 7, we show results for the PHS(k) generating functions with stencil sizes given by (4.5), polynomial augmentation of degree = 3 with k ∈ {2, . . . , 5}, = 5 with k ∈ {2, . . . , 7} and ∈ {6, 8} with k ∈ {2, . . . , 8} ( ≥ 5 only for Halton nodes since the stencil size for Cartesian nodes according to (4.5) would be rather large). In particular, in these tests k and are no longer connected via (3.6) or (4.1). However, the polynomial degree now determines the stencil size via (4.5). These tests show that violating the unisolvency condition (3.4) (by using (k, ) combinations (5, 3), (7,5), (8,6)) significantly increases the error (4.2) and the condition numbers of the weight and stiffness matrices. The results for = 8 in combination with k = 2 show that additionally (as mentioned in [62]) it can be problematic if k is too small. Furthermore, it can be observed that using higher values of k than those given by (3.6) can further decrease the error. The equation (3.2) for the convergence order can be interpreted in the way that for increasing N and , the other parameters n and k become less influential. Our numerical tests indicate that for our settings (4.1) leads to smaller approximation errors compared to (3.6). However, [62] points out that using a larger k than given by (3.6) can lead to numerical instabilities (especially for a large degree ).
The plots in the last row show results that have been computed without scaling and shifting, i. e., using (2.8) instead of (4.4). Comparing with the plots directly above, it can be seen that shifting and scaling has an effect on the condition numbers of the weight matrices but leads to hardly noticible changes in the approximation errors (except for k = 2).
In Fig. 8, we show results with respect to the total number of nodes N for the generating function PHS(k). Each color/marker represents a different degree of polynomial augmentation ∈ {2, 3, 4, 5, 6, 7, 8}. This degree determines k by (4.1) and n by (4.5). Further tests (not included here) on unstructured nodes with stencil sizes n determined via n ∈ {2M, 3M} instead of (4.5) do not lead to significant changes in the error (4.2) or the condition numbers of the weight or stiffness matrices whereas n = M + ln(2M) turned out to be too small (in most cases), leading to significantly larger errors and condition numbers of the weight and stiffness matrices. For Cartesian nodes, we show results only for ≤ 5 since larger degrees require very large stencil sizes, rendering the RBF-FD approach computationally very costly and hence not competitive (e. g., n = 420 for = 6 on a Cartesian grid versus n = 335 for = 8 on a Halton grid).
We find that for ∈ {2, 3, 4, 5}, we obtain comparable errors on Cartesian and Halton node sets (for significantly larger stencil sizes on Cartesian grids for ≥ 4 (4.5)). Larger degrees of polynomial augmentation (and with it larger stencils) lead to smaller discretization errors but larger condition numbers in the weight matrices.
The tests in this subsection illustrated that the PHS(k) generating function can be used successfully as long as the parameter k, the degree of polynomial augmentation and the stencil size n are adjusted correspondingly.

Numerical Results for Hybrid (HYB) RBFs and Comparison with Other RBFs
In this subsection, we show numerical results for hybrid generating functions HYB(γ ) and comparisons to GA, IMQ, IQ, MQ and PHS (see Table 1). We will mostly use γ = 10 −5 (which fulfills (3.1)) but also include some numerical experiments with γ = 10 −9 or γ = 10 −1 . The results of this subsection are presented through the following Figures.
• Figure 9   The results of Fig. 9 may be compared to those in Fig. 3 since we have only exchanged the Gaussian generating function and now use HYB(10 −5 ). A noticeable difference can be observed as the shape parameter ε goes to zero: the condition numbers of weight and stiffness matrices now remain bounded, and while there still is an optimal shape parameter ε ∈ (1, 4), the errors/condition numbers no longer blow up for ε → 0.  Fig. 4 in the same sense that only the Gaussian generating function has been replaced by HYB(10 −5 ). On the Cartesian grid, we observe the same jump in accuracy at the stencil size n = 33. When increasing the stencil size n, the condition numbers now remain bounded (and significantly smaller than for the Gaussian) so that the discretization error no longer blows up due to ill-conditioning for large n. Figure 11 shows test results for the same setting as in Fig. 5, again for HYB(10 −5 ) instead of GA. It shows that HYB(10 −5 ) can lead to similar errors as the GA generating function with the additional benefit of numerical stability for smaller shape parameters, larger stencil sizes and finer discretizations.
In Fig. 12, we use stencil sizes n = 50 (top four rows) and n = 131 (bottom two rows) to compare results for the generating functions GA, IMQ, IQ, MQ, HYB(10 −9 ), HYB(10 −5 ), HYB(10 −1 ) on Cartesian, PNP, Halton and random nodes. For each combination of generating function, grid type and stencil size, we beforehand numerically determined the constant c ∈ {0.05, 0.1, 0.15, 0.2, 0.25} in (4.3) leading to the smallest approximation errors over a wide range of tested problem sizes N , see Table 3. We observe that on all except the random node sets, the RBF-FD variants lead to more accurate approximations than the 7-point FD stencil (except for very small problem sizes N ). On PNP and Halton grids, all except HYB(10 −1 ) lead to approximation errors comparable or even better than the 19-point FD stencil on the Cartesian grid. Comparing the three hybrid versions, we note that HYB(10 −1 ) (compared to HYB(10 −5 )) reduces the condition numbers but increases the error, while HYB(10 −9 ) leads to condition numbers and approximation errors comparable to those for the infinitely smooth RBFs.
For n = 131, we include only results for PNP and Halton nodes since there is no improvement on Cartesian nodes compared to n = 50 (or, in fact, n = 33), and the results for random nodes are qualitatively similar to those for Halton nodes. Increasing from n = 50 to n = 131 on PNP or Halton nodes leads to a significant reduction in the error (accompanied by an increase in the condition numbers of the weight matrices) for all generating functions except for HYB(10 −1 ). (HYB(10 −1 ) behaves similar to PHS, i. e., the error is mainly decreased by increasing the polynomial degree and not the stencil size n.) In Fig. 13, we add polynomial augmentation to the HYB(10 −5 ) generating function. We also add two horizontal (dashed gray/green) lines for the (shape parameter independent) PHS(2) results using the two polynomial augmentation degrees ∈ {3, 5}, resp. (and otherwise identical setting). On the Cartesian grid, HYB(10 −5 ) with = 5 failed for ε ∈ {0.5, 1} (with a breakdown in the LU factorization of the weight matrices) which we indicate by the dummy values 10 39 and 10 22 for the condition numbers of the weight and stiffness matrices, resp. We recall that a minimal stencil size n depending on the polynomial augmentation degree is necessary for the weight matrices to be non-singular. For = 6, the polynomial space 6 (2.4) for d = 3 has dimension 84 (2.5), hence the fixed stencil size n = 90 satisfies  Table 3, polynomial augmentation of degree = 0, stencil size n = 50 for Cartesian, PNP, Halton and random nodes and stencil size n = 131 for PNP and Halton nodes The best results in terms of the error (4.2) and condition numbers of the weight and stiffness matrices are obtained for Cartesian nodes for = 2 and for Halton nodes for = 4. Larger degrees do not appear advisable since they increase the computational complexity without reducing the approximation error. The magnitude of the optimal shape parameter increases (and the range of suitable shape parameters widens) as the degree of polynomial augmentation increases (especially on the Halton grid), and an appropriately chosen shape parameter leads to significantly better results than PHS(2) by itself. While ε → 0 leads to blow-up (and no improvement though polynomial augmentation) for the Gaussian kernel, HYB(γ ) with polynomial augmentation behaves similar to PHS with polynomial augmentation as ε → 0 (since the Gaussian approximates a constant contained in the polynomial augmentation).
This figure as well as many additional tests (not included here) with a setting as in Fig. 6 but for HYB(γ ) for γ ∈ {10 −9 , 10 −5 , 10 −1 } have shown that some results for PHS regarding polynomial augmentation carry over to HYB (e. g., Cartesian nodes need larger stencil sizes than unstructured nodes). Hence, for the remaining tests, we use (4.5) to determine the stencil size n for HYB when used in combination with larger degrees .
The setup in Fig. 14 is similar to Fig. 12 with PHS instead of the generating functions GA, IMQ, IQ, MQ, and now larger degrees . Each color/marker represents a different generating function (HYB(10 −9 ), HYB(10 −5 ), HYB(10 −1 ) and PHS(k)). The top four rows of plots show results for = 3 on Cartesian, PNP, Halton and random nodes, resp. The polynomial degree determines k and n via (4.1) and (4.5), leading to ( , n, k) = (3,43,4). In the bottom two rows, we show results for the two parameter combinations ( , n, k) ∈ {(5, 116, 5), (8,335,8)}. For each combination of γ value in HYB(γ ), grid type and degree of polynomial augmentation, we beforehand numerically determined the constant c ∈ {0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4} for the problem size dependent shape parameter in (4.3) leading to the smallest approximation errors over a wide range of tested problem sizes N , see Table 4. The PHS generating function needs polynomial augmentation of degree ≥ D (3.2) in order to converge and we use respective degrees for the HYB(γ ) generating function as well. The γ value seems to be mainly important for the numerical stability and not so much for the error in the case of larger degrees . Tests with = 5 (not included here) showed that Cartesian nodes are not an advisable option for RBF-FD with HYB in combination with larger degrees of polynomial augmentation since the required stencil sizes are significantly larger than for unstructured nodes without leading to smaller errors, similar to the PHS case studied in Sect. 4.3. Tests for ∈ {5, 8} with k according to (4.1) show a decrease in the approximation error as increases with hardly any noticible difference between HYB(10 −9 ) and PHS. However, for large , PHS seems to be the better option in terms of numerical stability.

Summary
In Sect. 4.2, we focused on the main difficulty of the (infinitely smooth) GA generating function: the selection of a shape parameter.
We concentrated in Sect. 4.3 on the PHS generating function. Tests on Cartesian nodes (not included here) exhibited problems for larger (and required excessively large stencil sizes). For non-Cartesian grids, determining the stencil size via (4.5) and the PHS degree via (4.1) (or in case of numerical instabilities via the smaller (3.6)) seems reasonable in terms of accuracy and numerical stability.
In Sect. 4.4, we showed tests for the HYB generating function which combines GA with PHS(2) and compared results with those obtained for several infinitely smooth RBFs and PHS. On one hand, we repeated the tests in Sect. 4.2 with HYB(γ ) instead of GA. Our observations (in the setup with no or only small degrees of polynomial augmentation) agree with [46] that typically smaller errors and larger condition numbers are obtained by decreasing γ . Additionally, we showed that HYB(γ ), GA, IQ, IMQ and MQ generating functions can lead to comparable errors (for carefully selected shape parameters, see Table 3, Table 4 and (4.3) for general guidelines) and if γ is chosen according to (3.1).
We also focused on the role of polynomial augmentation. For GA, IQ, IMQ and MQ generating functions, typically a small degree (e. g., ≤ 2) is suggested. For HYB(γ ), with a sufficiently large γ , we observed that increasing can lead to smaller errors than small (always paired with suitable shape parameters). Therefore, we compared PHS and HYB generating functions with larger degrees of polynomial augmentation. These tests showed that > 3 on Cartesian nodes is also not advisable  Table 4, for the HYB generating functions on Halton nodes (bottom two rows) for HYB. For ≤ 5 on PNP and Halton nodes, HYB can lead to smaller errors than PHS. For sufficiently large , HYB converges to PHS(2) as ε → 0 which is reflected in the illustrated approximation errors (the condition numbers of the weight and stiffness matrices depend on γ ). An advantage of using degrees ≥ 3 for HYB is the decrease of sensitivity with respect to the choice of γ or the shape parameter.
When it comes to the underlying nodes, one observes "PNP ≤ Halton ≤ random" where "≤" relates the approximations errors (as well as condition numbers of the weight and stiffness matrices) obtained for those node sets. However, for PHS and HYB generating functions with larger degrees of polynomial augmentation, the differences are rather small, making HYB and PHS good options for adaptive node refinement and stencil selection as long as stencil sizes according to (4.5) (or larger) are used (for nearest neighbor stencils).
Hence, our recommendations for a good RBF-FD discretization setup using basic algorithms as described in this paper may be summarized as follows: • If Cartesian nodes are used, then a standard finite difference discretization or HYB with a small degree ∈ {0, 1, 2} of polynomial augmentation and γ according to (3.1) should be used. (Alternatives such as stabilized RBF-GA [58] may also be viable.) • On unstructured nodes, PHS with polynomial augmentation offers a shape-parameterindependent approach. The required parameters , n, k may be determined via (3.2), (4.5) and (4.1) (or (3.6) if numerical instabilities occur), resp.. Best results are observed for large exponents and polynomial degrees k = = 8 which require rather large stencils. • On unstructured nodes, GA, MQ, IMQ, IQ or HYB with carefully finetuned parameters and constant polynomial augmentation provide options to achieve similar accuracies as PHS with smaller stencils and hence may decrease the computational time or memory compared to the PHS setup. • (Unstructured) PNP node sets often lead to the best results (smallest approximation error) among the various node sets included in our tests and are a promising option for more complicated domains , adaptive node refinement or variable density node sets.
These recommendations almost certainly will have to be adjusted when more advanced algorithms for node generation or stencil selection, stable versions of radial basis functions or computation of weights for deficient sets, etc., are used which are beyond the scope of this paper.

Conclusion and Outlook
RBF-FD provides a meshfree and relatively easy to implement approach for the discretization of partial differential equations. Its realization, however, involves a multitude of parameter choices, many of which (but by far not all) have been tested and illustrated in this work. Our conclusions for these have been summarized at the end of the previous section. It turns out that several methods work well if parameters are chosen accordingly, and there is no clear winner (at least not among the tested variants). However, if parameters are not selected carefully, there are many "losers" leading to stagnation/divergence of the approximation error. A straightforward extension of this work would lie in the consideration of further RBF-FD variants, including stabilized versions, variable shape parameters, other than nearest neighbor stencils, further types of generating functions or their combinations. It would furthermore be of interest to extend our study to other types of partial differential equations.
We have focused on the discretization process and resulting approximation errors of RBF-FD discretization and mostly left out considerations of computational complexity. The computation of stencil weights may be executed in parallel, hence an increase in stencil size (or polynomial augmentation degree) may not be a severe drawback with respect to computational time for the setup phase. However, it will lead to a less sparse and worse conditioned stiffness matrix and hence increase the computational time required for solving it (directly or iteratively). We plan to design and analyse iterative solvers for linear systems arising in RBF-FD discretization next.
Author Contributions All authors contributed to the work presented in the manuscript. All authors read and approved the final manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL. The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.
Data Availability Data sharing is not applicable to this article as no datasets were generated or analysed during the current study. (Scattered data sets as a basis for RBF-FD algorithms were generated with well documented and easily reproducible standard algorithms.)

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.