The D-RBF-PU method for solving surface PDEs

In this paper a new direct RBF partition of unity (D-RBF-PU) method is developed for numerical solution of partial differential equations deﬁned on smooth orientable surfaces which are discretized with sets of scattered nodes and with approximations to normal vectors at each of the nodes. The accuracy, stability and eﬃciency of the new method are studied through some theoretical and experimental results. This method is a localized RBF based technique, results in a perfectly sparse ﬁnal linear system, uses only scattered nodes on the surface rather than a connected mesh, and is applicable for a large class of PDEs on manifolds. Applications to some biological and chemical reaction-diffusion models are also given. Results show that the new method outperforms other comparable techniques for surface PDEs. © 2023 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons .org /licenses /by /4 .0/).


Introduction
Radial basis function (RBF) approximation is a powerful tool for reconstructing an unknown function from its given scattered data values [7,56].This approximation has become increasingly popular for solving partial differential equations (PDEs) and other applications due to its easy implementation, high-order convergence rates and its ability to handle scattered node layouts on arbitrary domains [16,18].In this paper we focus on solving surface PDEs using a new technique based on this interesting approach.
Solving PDEs on surfaces have received growing interest due to some applications in meteorology, geometry, biology, image processing, etc.The finite difference (FD) method is a standard technique for solving PDEs in Euclidean domains where the one dimensional FD formulas can be combined to create FD formulas in higher dimensions for partial derivatives.The first limitation of this FD approach for surface PDEs is that the extension of univariate formulas for surface derivatives is not straightforward.The generalized multivariate FD formulas can be instead organized on stencils with scattered nodes, but this approach also raises the question of how the stencil weights should be computed for example in cases in which a polynomial interpolation on a surface is singular.For the special case of spheres the Euclidean polynomials can be replaced by spherical harmonics and a theoretical solvability is guaranteed.However, if the space of spherical harmonics is restricted to spherical caps (for local stencils) then the linear system for computing the stencil weights becomes highly ill-conditioned.On the other hand, an alternative to spherical harmonics for other manifolds is not available, directly.These drawbacks can be bypassed if RBF interpolation is used instead to generate the stencil weights in the FD formulas (e.g.[17,18,28,36,46,49,50]).This approach is called RBF-FD and results in a sparse and well-conditioned final system contrary to the global RBF approximation which produces full and ill-conditioned matrices.We note that RBF-FD has been frequently used for solving PDE problems in Euclidean spaces as well; see for example [18] and references therein.
Another approach for localization and evolving RBF into a fast technique applicable on large-scale problems is the RBF partition of unity (RBF-PU) method.The first combination of a PU method with RBF interpolation goes back to [55] and the analysis is given in [56].See also [16].A fast algorithm based on PU for two-dimensional interpolation of large scattered data sets was presented in [9].The RBF-PU collocation method for solving transport equations on the unit sphere was developed in [1].A least squares RBF-PU method as an alternative to the standard collocation method was studied in [26].RBF-PU has also been shown to be successful for a number of other PDE problems (e.g.[10,12,13,40,47,48]).For applications to rational interpolation see [11,15] and for applications to parametric structural and shape optimization see [23,24].
In [30] the new D-RBF-PU method was presented to simplify the implementation and increase the efficiency of the standard RBF-PU method.The direct approximation of PDE operators avoids differentiation against PU weight functions and allows the use of some useful discontinuous weights functions.This method is also related to the RBF-FD method but is much faster for setting up the stiffness and mass matrices.The efficiency of this approach for surface PDEs is even more notable because the computation and implementation of surface derivatives of PU weight functions is a challenging problem.In this paper, we develop the D-RBF-PU method for solving PDEs defined on smooth embedded submanifolds of Euclidean spaces.As a basis function, we use the restriction of Euclidean positive definite kernels on surfaces [20].On the PU part we apply both smooth and discontinuous weights for blending local RBF approximants and forming a global solution for the PDE problem.
The remainder of this paper is organized as follows.In sections 2 and 3 the continuous differential operators on embedded submanifolds in R d and their discrete versions via a specific RBF approximation are reviewed, respectively.In section 4 the well-known RBF-FD method is briefly discussed, and in section 5 the idea of PU approximation on surfaces and the standard approach for solving PDEs are presented.In section 6, the D-RBF-PU method is developed for numerical solution of surface PDEs.Then its advantages over other techniques, its connection to RBF-FD and its error analysis are given.In section 7, some numerical experiments, comparisons with RBF-FD and applications to some reaction-diffusion models are presented.

Continuous differential operators on surfaces
Throughout this paper, by terms surface or manifold we mean a smooth embedded submanifold of codimension 1 in R d with no boundary.The specific case d = 3 is considered here but straightforward extensions are possible for manifolds of higher codimension, or manifolds of codimension 1 embedded in higher or lower dimensional spaces.We denote a generic manifold by .
The standard way for defining the differential operators on manifolds is through the use of intrinsic coordinates.We instead follow [21,28,46] to formulate these operators entirely in Cartesian (or extrinsic) coordinates to avoid singularities associated with intrinsic coordinate systems.In this methodology the surface gradient of a differentiable function u can be simply formulated as follows.We apply the Euclidean gradient ∇ to u at a point x = (x, y, z) ∈ and then project the resulting vector into the tangent space of at x.More precisely where n = n(x) is the normal vector on at the point x.This shows that the surface gradient operator can be written entirely in Cartesian coordinates as ∇ := ∇ where = I − nn T is the projection operator.In an extensive form for d = 3 the surface gradient is written as where (n x , n y , n z ) =: n.The surface divergence of a smooth vector field f = ( The Laplace-Beltrami operator at x ∈ is then defined by Note that all expressions above are in Cartesian coordinates.

Discrete differential operators on surfaces
In this section, we present a RBF based method for discretizing the surface gradient and the Laplace-Beltrami operators.First, the RBF interpolation in Euclidean spaces is briefly reviewed and then approximation via basis functions constructed by restricting the RBF on a manifold is discussed.Finally two approaches for approximating the Laplace-Beltrami operator are given.The main part of the section follows the methodology given in [21].

RBF interpolation
We start with RBF interpolation which is essential for the method we aim to develop in next sections.Assume that ⊂ R d and : × → R is a kernel with property (x, y) = φ( x − y ) for a continuous and even univariate function φ : R → R.Here • stands for Euclidean norm in R d .Such kernels are called radial functions.We assume that φ is positive definite in the following sense.

Definition 3.1. A continuous and even function
where P m−1 (R d ) is the polynomial space of degree m − 1 (order m), then φ is called conditionally positive definite of order m.A positive definite function is conditionally positive definite of order m = 0 and any order m > 0. A radial and (conditionally) positive definite function φ is called a radial basis function (RBF).
For a given set X = {x 1 , . . ., x N } ⊂ and a continuous function u : → R sampled at points in X we consider the RBF interpolant of u on X of the form for a positive definite or first order conditionally positive definite function φ.The coefficients c j are obtained by imposing the interpolation conditions I φ u(x k ) = u(x k ), 1 k N and the side condition c 1 + • • • + c N = 0, or equivalently by solving the linear system positive definite or conditionally positive definite of order 1 then the above system is uniquely solvable [56].The inverse matrix can be obtained as where Z ∈ R N×(N−1) is any matrix whose columns form a basis for ker(e T ), ẽ = 1 N (I − S X A X )e and a = − 1 N e T A X ẽ.For more details about the properties and the inverse of a saddle point matrix of the above form (and even in a more general form) see [15].We note that Z T A X Z is necessarily positive definite on R N−1 .Thus, the interpolation coefficients are obtained as c = S X u, c N+1 = ẽT u.
Here, we are interested in RBF interpolations on manifolds but we still use the standard Euclidean distance • in R d instead of geodesic distances intrinsic to surfaces.In fact, we use the restriction of kernel on submanifold .Since is a subset of R d , the restriction | × inherits the positive definiteness from the original kernel which makes it suited for interpolation problems on .

Discretizations of surface gradient and Laplacian
Let the RBF φ have at least one continuous derivative, x = (x, y, z) ∈ , x j = (x j , y j , z j ) ∈ , and r j (x) = x− x j .Referring to (2.1), the x-component of the surface gradient of translates of φ by using the chain rule can be written as where denotes differentiation with respect to r.Note that the singularity in the above expression will analytically cancel, because φ is assumed to be an even function with at least one continuous derivative.Other components of the surface gradient of φ can be obtained similarly by

x) .
Note that, in the above expressions the normal vector components n x = n x (x), n y = n y (x) and n z = n z (x) depend on variable x although this dependence is ignored in the notation for brevity.The surface gradient of a function u : → R is then approximated by the surface gradient of its RBF interpolation on , i.e.
where the symmetric matrix S X is defined in (3.1).If we assume for 1 k N then the approximate surface gradient at points in X is obtained as X S X are differential matrices associated to continuous operators G x , G y and G z , respectively.
To discretize the Laplace-Beltrami operator the standard approach uses where in the form of (2.2), the gradient operators G x , G y and G z need to be applied twice requiring the first partial derivatives of normal vector n = n(x).This makes the approximation rather complicated especially when explicit representation of normal vector is not available or when the manifold is defined by point clouds rather than an explicit representation.On the other hand, second derivatives of basis functions are required.To bypass these difficulties, a repeated interpolation and differentiation approach is proposed in [19] as If we imitate the continuous form (2.2) by replacing the differentiation matrices G x X , G y X and G z X by the operators G x , G y and G z , respectively, then the approximation is obtained for the Laplace-Beltrami of u at nodes in X .The N-by-N differentiation matrix L X approximates the operator on the set X .

RBF-FD method on surfaces
The RBF-FD method arises naturally as a generalization of standard FD scheme where the one dimensional formulas are combined to create FD formulas in higher dimensions for partial derivatives.In a generalized FD method the stencil weights are generated by imposing exactness on (multivariate) polynomials.Due to drawbacks outlined in the introduction section, the RBF-generated formulas generate stencil weights by assuming exactness on the space of shifts of the RBF to scattered points in the stencil.
For a linear surface operator L, the value (Lu)(x 0 ) can be recovered from values of u at a stencil X 0 ⊂ X of points nearing x 0 .Indeed, we obtain the weights ξ 0, j such that where J 0 is the set of indices of points in X 0 .The weight vector ξ 0 will be computed by assuming exactness of the formula for basis functions u(x) = φ( x − x j ), j ∈ J 0 .Usually, exactness for constant function u(x) ≡ 1 is also imposed.This leads to system where Lφ 0 = (Lφ( x 0 − x j )) j∈ J 0 .Comparing with (3.3) we observe that ξ 0 = Lψ(x 0 ) where ψ is the vector of Lagrange basis functions on set X .If φ is conditionally positive definite of order 0 or 1 then the above system is guaranteed to be solvable.For solving a linear PDE problem by RBF-FD, we consider two sets X = {x 1 , . . ., x N } (trial set) and Y = {y 1 , . . ., y M } (test set) of scattered points on .
Then we loop over test points y k and give the role of x 0 to all y k successively.The stencil weights ξ k on points X k in neighborhood of y k , when interspersed with zeros corresponding to points outside the stencils, form the rows of a then gives the approximate nodal values at the trial set X .For postprocessing calculations, the RBF-FD formula can be used again to find the approximate values of u and its derivatives at any evaluation point x ∈ .
To tackle the problem of solvability we may use a large enough set Y to have a full rank overdetermined matrix A L and look for a unique least squares solution u instead.See [42][43][44] for more details.However, in practice the degeneracy happens in rare situations and usually the case Y = X gives an invertible matrix A L .For some applications of RBF-FD to surface PDEs see [28,46].This ensures that every point x ∈ is necessarily covered by at least one patch .See Fig. 1 for an example of covering on the unit sphere and on a torus.A family of nonnegative functions {w } N c =1 is called a partition of unity (PU) with respect to where supp(w) = {x ∈ M : w(x) = 0}.The title 'partition of unity' is induced from the second property.
We start with an overlapping covering { } N c =1 of .If we assume V is an approximation space on is a global interpolation of u on which is formed by joining the local approximants s via PU weights w .A PU approximation provides a collection of local approximants (small problems) and blends them via a set of PU weights to form a global approximation.This idea is numerically efficient because of avoiding a large (and possibly unstable) system in favor of several small (and possibly stable) local systems.
The smoothness of the global solution depends on the smoothness of both local approximants and PU weights.A possible choice for w are Shepard's weights where are nonnegative, nonvanishing and compactly supported functions on .Assume that the linear PDE problem (4.2) is posed on a manifold .A standard PU method for discretizing this PDE finds the global PU approximation s of u and then approximates Lu by Ls: In this scenario, the surface operator L should act on products w s which may result in a complicated computation as it requires a kind of Leibniz's rule for surface derivatives of weight functions.The rational structure of Shepard's weights makes it more challenging especially when higher derivatives are demanded.This may be the reason why the PU approach has been rarely employed for surface PDEs.We refer the reader to [1] as an available implementation on the sphere.In the present paper we suggest an alternative approach which avoids the differentiation against weight functions and allows the use of some discontinuous weights to develop faster algorithms.

D-RBF-PU method
Assume that X = {x 1 , . . ., x N } is a set of well distributed scattered points on surface and { } is a covering for .Besides, let X = X ∩ .When the local approximants s are local RBF interpolations on X the resulting PU method is called RBF-PU method.This method is an attractive alternative to global RBFs as it performs better in terms of accuracy per computational cost because its resulting differentiation matrix is sparse.On the other side, RBF-PU inherits the advantages of global RBF methods in terms of handling irregular geometries and scattered node layouts.However, RBF-PU when applying on surface PDEs suffer from some disadvantages which were outlined at the end of previous section.To overcome this we follow the D-RBF-PU approach introduced in [30] for solving problems on bounded Euclidean domains and extend it for surface PDEs in this paper.In D-RBF-PU the value Lu is directly approximated by the PU method without any detour via local approximants s of u, say where s L are local approximants of Lu on patches .The resulting global approximation is denoted by s L which is clearly different from Ls of the standard approach (5.3).The first notable difference is that in the new method the weight functions are not required to be differentiated.However, these two methods are basically different as they fall into two different discretization categories, namely the trial space and the direct discretizations highlighted by Schaback in [42,44] in the context of scattered data approximations.See also [31,32] for developments of the second approach in the moving least squares approximation context.
Assume that the set of indices of points from X in patch is denoted by J , i.e.
In RBF based PU methods the local approximation spaces on are constructed by the shifts of RBFs to points in X (possibly) plus a polynomial space.In our case in this paper we assume where ψ j ( ; •) are Lagrange functions associated to patch .In D-RBF-PU an obvious choice for local approximants s L in (6.1) is which is used successfully in [30] for different derivative operators L defined on Euclidean domains in R d .However, as we pointed out in section 3.2 this will cause serious difficulties for handling the derivatives of normal vectors when L is, for example, the Laplace-Beltrami operator on surfaces.An alternative choice for case L = then is s L = ∇ • (I φ,X (∇ (I φ,X u))).(6.3)This repeated interpolation and differentiation approach, as described before, avoids differentiation with respect to normal vectors and highly simplifies the computation.By mimicking the relation (3.4) we have where u stands for u values at X .See [28,46] for the use of this approach in RBF-FD methods on surfaces.This approach gives the approximate Laplacian at all interpolation points in X , simultaneously, while the only approximation at the central point is used in RBF-FD methods.In contrast to this, the D-RBF-PU method benefits from the approximation of at all points in X .We note that for the surface gradient ∇ the standard interpolation s L = ∇ (I φ,X u) will be used.In both cases (L = ∇ and L = ) we may expand the local approximant in terms of function values as instead of (6.2) to obtain the global PU approximation Usually, Shepard's functions (5.2) are used as PU weights.These functions are as smooth as their generating compactly supported functions . For a test point x ∈ let I(x) be the set of indices of patches in { } which cover the point x.A PU approximation at point x is influenced by all trial points in ∪ ∈I(x) X leading to a denser (but still sparse) differentiation matrix compared with other localized methods such as RBF-FD.Nevertheless, RBF-PU is computationally less expensive than RBF-FD for setting up the differentiation matrix because the number of local linear systems to be solved is highly reduced from N, the number of test points, to N c , the number of PU patches.In additions, D-RBF-PU allows the use of discontinuous weight functions because in this method the weights are not required to be differentiated.A useful discontinuous weight can be defined by w (x) = 1, if is the nearest patch to x, 0, otherwise, ( which gives the total weight 1 to the nearest patch and null weights to the others.To identify the nearest patch we can measure and compare the distance between x and the center of all patches.This weight function results in a final matrix which is as sparse as the RBF-FD matrix because it assigns a single patch per any test point.Thus, the D-RBF-PU method with weight function (6.6) is faster than RBF-FD for setting up the differentiation matrix and it is as fast as RBF-FD for solving the sparse final linear system.We will see in numerical tests that the majority of costs for both methods is subjected to the setting up phase because fast iterative linear algebra solvers are available for solving the final sparse system.As another remark, D-RBF-PU can be viewed as a RBF-FD method in a PU setting.Comparing (4.1) and ( 6.4), one can see that for a test point is covered by patches for ∈ I(x k ), all RBF-FD approximants s L (x k ), ∈ I(x k ), are computed and then joined by to form the overall approximant s L (x k ).The important difference is that a local approximation s L is shared with all test points which are covered by patch .This is the key property that makes the D-RBF-PU much faster than the standard RBF-FD method.The situation becomes simpler and the algorithm becomes faster when the weight function (6.6) is applied.
In this case, each test point x k is subjected to a single local set X = ∩ X for an index = k ∈ {1, . . ., N c } in which k is the closest patch to x k .From (6.5) and (6.6) we have = k and 0 otherwise.Again we emphasize that the difference between this approach and the RBF-FD method is that in the new method a set (stencil) X is shared with many test points while in the RBF-FD each stencil is associated to a unique test point.This is also connected to the overlapped RBF-FD approach presented in [45].

Some points about error bounds
In our PU setting we assume that the assigned covering { } is regular which means that there exists a global constant K such that every x ∈ is covered by at most K patches out of { }, i.e. |I(x)| K for all x ∈ .On the other hand, we assume that for a given discrete set X ⊂ there exists a global constant C such that diam( ) Ch X, for 1 where diam( ) is the maximum geodesic distance between points in .In this case the covering { } is called local with respect to the discrete set X .
As it is described in [30], the analysis of D-RBF-PU method follows the theory given in [43] for nodal meshless methods.
If we are looking for the errors in nodal values u = (u(x 1 ), . . ., u(x N )) T , the consistency is analyzed by finding an upper bound for A L u − f p where • p is the p-norm in R N .According to the construction, A L (k, : Here, by A(k, :) we mean the k-th row of A. Thus, for consistency we assume there exist error bounds If the final matrix A L is full rank and admits the global C S for a constant C S independent of spacing distance between points in X , then This is a classical error analysis containing the stability and consistency terms on the bound.However, our analysis may leave some parts unproven.The first open problem concerns the invertibility of the square matrix A L for the D-RBF-PU method.One can bypass this problem by overtesting and seeking for a unique least squares solution instead of an exact solution for A L u = f (see e.g.[44,43]).If the repeated interpolation and differentiation approach is used for approximating u then the test and trial points should coincide and overtesting is not applicable.One can use the standard approach (6.2) instead but it needs manipulation with derivatives of normal vectors.However, many experiments show that the square matrix A L is invertible for different distributions of trial (and test) points X .Thus, in almost all practical situations overtesting is not required at all.The existence of the stability constant C S is the next unproven ingredient not only for the method of this paper but also for all nonsymmetric meshless methods such as RBF-FD and standard RBF-PU methods.In experiments with Matlab for the case p = ∞ we can use the fact that the condest command estimates the L 1 condition number which is the L ∞ condition number of the transpose.Thus, we can estimate C S by L .This command is very fast for sparse matrices.In case p = 2 one can compute the smallest singular value of A L via SVD to measure the 2-norm of the inverse matrix.Experiments show excellent stability bounds for the D-RBF-PU method on various surfaces.See section 7.
The consistency bounds ε L for surface operators ∇ and can be estimated by modifying the results of [21,33] for local patches ⊂ instead of itself.The theory holds for finitely smooth positive definite RBFs whose Fourier transforms decay algebraically.We assume that where K β is the modified Bessel function of second kind of order β.The native space of a positive definite function with property (6.7) on R d is the Sobolev space H τ (R d ) while the native space of the restriction of on is the Sobolev space H τ −1/2 ( ).We refer the reader to [41,56] for definitions and more details about native spaces of (conditionally) positive definite functions and to [20] for properties of restricted kernels from R d to smooth embedded submanifolds.

The above error bounds can be simply updated if we replace the compact manifold by open and bounded coverings
⊂ and X by X = X ∩ .For a fixed , since is a paracompact submanifold of R 3 there exists an atlas A = {( k , U k )}, k = 1, . . ., K of finitely slice charts for with an associated intrinsic atlas A = {( k , U k )} [27].The sets U k ⊂ R 3 are open in R 3 and cover , and each k : U k → B(0, 1) is a one-to-one smooth map.For the intrinsic atlas the sets U k ⊂ are open and k : U k → B (0, 1) = {x = (x, y, z) ∈ B(0, 1) : z = 0} are considered as copies of the open ball in R 2 .Now, let {χ k } be a partition of unity subordinate to {U k }.For a function u defined on , the projections π k (u) : R 2 → R are defined by Sobolev spaces H σ ( ) are defined by the norms , where K is the number of charts in the atlas.It is obvious that the norm depends on the special choice of atlas A and partition of unity {χ k }.But, the same space arises by different collections of these objects and the norms are equivalent.
Using this construction, the theory of [21,33] can be adapted to conclude a similar theorem with and X being replaced by and X , respectively.But since the periodicity is lost on local surfaces ⊂ , the duality trick of [33] can not be applied to double the orders for smoother functions in H β for values of β up to 2σ .The error bounds read as u H β ( ) , 3 < β σ , (6.8) for all u ∈ H β ( ).Finally we have the following theorem.Theorem 6.2.Let be a smooth 2-dimensional submanifold of R 3 and X ⊂ be a discrete and quasi-uniform set of points with sufficiently small fill-distance h X, .Let { } N c =1 be an open, bounded, regular and local (with respect to X) covering for with PU functions {w }.If s L is the direct PU approximation of Lu for L = Id, ∇ , , then (6.9) for all u ∈ H β ( ) and all x ∈ X provided that = φ( • 2 ) satisfies (6.7) with τ > k + 3/2 and σ = τ − 1/2 where k = 0, 1, 2 for L = Id, ∇ , , respectively.
Proof.Since {w } forms a partition of unity, we have for any It remains to find the local bounds ε L for = 1, . . ., N c , which are obtainable from (6.8) because we have (s h X, and u H β ( ) u H β ( ) for all = 1, . . ., N c , the final bound (6.9) results.

Cost comparison
The total computational cost of a localized RBF-based method consists of (i) the cost of setting up the final differentiation matrix, and (ii) the cost of solving the final linear system.Usually, the final system of such a method is highly sparse and can be effectively inverted using the available iterative linear algebra solvers for sparse and large systems.Thus, the majority of the computational costs arises in the first phase where on one hand a data structure should be performed to collect the indices of local sets, and on the other hand lots of local linear systems should be solved for calculating the stencil weights or constructing the local approximants.
Assume that we are given N quasi-uniform trial and N test points together with N c number of PU patches on surface .According to the construction, N c is much smaller than N. Besides, the number of trial points in each patch is independent of h X, (and thus independent of N) because the covering is assumed to be both local and regular.This means that the number of points in each patch remains unchanged, approximately, when the discretization is refined.Thus if n := |X | then there exists a global constant n P independent of the discretization distance h X, such that n n P .This is true for the size of RBF-FD stencils as well and we can assume | X k | n P for 1 k N. Assume that the set of indices of trial points in each local domain is known in advance.In D-RBF-PU N c and in RBF-FD N local linear systems should be solved to set up the final matrix.In fact, in RBF-FD the stencils are changed per test points while in D-RBF-PU a local patch is shared with many test points in a certain neighborhood.Consequently, the total cost for solving all local systems is of order n 3 P N c for the D-RBF-PU method and it is of order n 3 P N for the RBF-FD method.As in the PU approximation N c is a small fraction of N, say N c ≈ N k , the D-RBF-PU is k times faster than RBF-FD for solving all local systems.We will observe this in numerical experiments of section 7.
The cost of collecting the indices of points among all N trial points in each patch or stencil is ignored in the above cost analysis.This can be done in a preprocessing step using a boxing strategy to avoid collecting local indices for each patch or each stencil separately.The overall complexity would be of order 2N for RBF-FD and of order N + N c for D-RBF-PU.As a more relevant strategy we can use a kd-tree algorithm that takes O(N log N) time for building a tree on X and O(N 1−1/d ) time for range searchings.
Finally, in D-RBF-PU the complexity of evaluating the final solution at an evaluation set X e is dominated by the cost of building another data structure on X e to collect the indices of evaluation points in each patch .However, in RBF-FD new stencils should be formed per each evaluation point and new local systems should be solved.
Altogether, both methods possess the same complexity rate but D-RBF-PU is approximately k times faster than RBF-FD where N is approximately k times N c .

Numerical experiments
The first part of this section is devoted to some numerical experiments of the D-RBF-PU method and its comparison with the RBF-FD method for a PDE of the form on some smooth embedded submanifolds ⊂ R 3 of dimension 2. Here, b ∈ R 3 and c ∈ R are sufficiently smooth functions of spatial variable x ∈ , and f is a given source function.In numerical experiments of this part the unit sphere, a torus and a red blood cell (RBC) are used as computational domains.In Cartesian coordinates these surfaces are defined as We set c 1 = 1, r 1 = 1 3 , c 0 = 0.25, c 2 = 2, c 4 = −1.5 and r 2 = 1.For trial points and patch centers we use phyllotaxis spiral nodes based on the idea given in [38] in all experiments.These nodes are also available for download from [22].The points on RBC are obtained by mapping the points on the unit sphere using the parametric equations of this surface.Example node sets are presented in Fig. 2.
A PU covering { } is simply constructed by letting = B(ω , ρ) ∩ for N c quasi-uniform centers {ω 1 , . . ., ω N c } =: X c ⊂ .Here B(ω, ρ) is the open ball in R 3 of radius ρ and center ω.Patch centers ω on are obtained by the same procedure as described for trial points.For quasi-uniform points we set N c = N/16 to have h X c , ≈ 4h X, .The patch radii are set to be proportional to h X c , , say ρ = C ovlp h X c , , to obtain a local and regular covering.The overlapping constant C ovlp > 0 should be large enough to ensure both the inclusion ⊂ ∪ and the accuracy of local approximants.On the other side, it should be as small as permissible to get a reasonably sparse final linear system.We will illustrate some experiments to come up with a proper constant C ovlp .
As a smooth PU weight, the function 2) where ψ(r) is the C 4 compactly supported Wendland's function.However, we will apply the constantgenerated PU weight function (6.6) in most of our experiments.
In this study the restrictions of kernels Matérn (MA): + 32(εr) 3 + 25(εr) 2 + 8εr + 1 on different surfaces are used as basis functions.The first kernel is positive definite on R d for d < 2τ while the second one is positive definite on R d for d 3. The native space of the MA kernel on the all above 2 dimensional surfaces is H τ −1/2 ( ) while it is H 4.5 ( ) for the WE kernel.The parameter τ = 5 is used in all experiments.In both cases the shape parameter ε > 0 is tuned such that a reasonable accuracy is obtained without encountering any serious instability.We set ε = 3 and ε = 1/3 for MA and WE, respectively.
In all cases we use the function u(x, y, z) := cos(10z + 3) e −x 2 − sin(4x + 5 y), as a true solution of the PDE.The velocity vector b(x) ≡ (1, 2, −1) and c(x) ≡ 1 are used in all experiments.In each case, the right hand side function f is calculated, accordingly.
To adjust a good overlapping constant C ovlp we test the accuracy vs. amount of sparsity for different surfaces with both MA and WE kernels.Some results are given in Fig. 3 on sphere, torus and RBC surfaces.Relative errors (left vertical axis) and percentage of nonzero elements of the final matrix (right vertical axis) are shown in terms of the overlap constant C ovlp .The smooth PU weight is applied for the left panels and the constant-generated PU weight for the right ones.According to these results and other experiments with different numbers of points which are not presented here for the sake of brevity, we set C ovlp = 1.1 to have a balance between accuracy and sparsity.We note that in this experiment the iterative approach is applied for discretizing the Laplace-Beltrami operator.While not presented here, a similar behavior is observed with the standard approach.We further observe from Fig. 3 that the D-RBF-PU method with the constant-generated weight is about 2.5 times sparser than the case with the smooth weight function.
In Fig. 4 errors and convergence orders of D-RBF-PU with smooth and constant-generated weights and comparison with RBF-FD on three surfaces using WE and MA kernels are depicted.For results on this figure the Laplace-Beltrami operator is discretized based on the repeated interpolation and differentiation approach.The same results for standard approach are given in Fig. 5.In all plots computational convergence orders are obtained by the linear least squares fitting to error values and are written alongside the figure legends.To make the comparison fair the RBF-FD stencils are chosen to be X k = B(x k , ρ) ∩ X where ρ is the patch radius in D-RBF-PU.Both methods possess the same theoretical order although in this experiment the D-RBF-PU method with both smooth and discontinuous weights gives more accurate results and higher convergence orders.However, it seems possible to scale both methods in different ways to obtain better accuracies and convergence rates.A selection of plots for the stability constant is presented in Fig. 6 for the three surfaces with WE kernel.While not presented here, the same behavior is observed for experiments with the MA kernel.The final matrix is extremely stable as C S ≈ 1 in almost all cases.We emphasize that the stability constant C S measures the conditioning of the final sparse differentiation matrix.It is of course different from the conditioning of local RBF systems on PU patches.The problem of conditioning of local systems at small values of the shape parameter ε might be treated by some tricky techniques such as Contour-Padé or RBF-RA algorithms [57,58].Equipping the method with such stabilization techniques allows to run the algorithm for higher values of N to obtain a better accuracy.We did not pursue this issue in this paper, instead we tuned the shape parameter to obtain a reasonable accuracy without encountering any instability.
As we pointed out before, D-RBF-PU is computationally cheaper than RBF-FD.In Fig. 7 the CPU times used for setting up and for solving the final linear systems together with the total CPU times are compared.Since these systems are sparse and well-conditioned, the solving time is very small compared to the setting up time.Looking at the total time plot, although the rates are the same, D-RBF-PU is about 10 times faster than RBF-FD because in D-RBF-PU much fewer number of local systems should be solved to form the final matrix.While not presented here, cost plots for other surfaces obey the same behavior.

Application to reaction-diffusion equations
Among others, reaction-diffusion equations arise from mathematical models to describe Turing pattern formations in some biological applications and spiral waves in an excitable chemical media [3,14,34,35,51,53].The solutions of reactiondiffusion equations can exhibit pattern forming from random initial conditions.However, most results reported in the literature concern these phenomena on Euclidean domains.Growing attention has been recently paid for such models on general surfaces to investigate the effect of curvature to pattern and wave formations.Applications on simple surfaces such as spheres and cones could be found in [37,54,59] and applications on general surfaces may be found in [5,6,8,21,25,29,39,52].We consider the system of reaction-diffusion equations  where f and g are nonlinear reaction terms and δ 1 and δ 2 are diffusion parameters.Here u = u(x, t) and v = v(x, t) are chemical unknowns where x and t stand for spatial and time variables, respectively.The initial conditions are imposed by u(x, 0) = u 0 (x) and v(x, 0) = v 0 (x) for x ∈ and known functions u 0 and v 0 .
In this section we consider two examples of this time-dependent system on different surfaces.Additions to surfaces defined in the previous subsection, a Dupin's Cyclide, a tooth shape surface and a frog will be used in our experiments.The surfaces and example node sets are presented in Fig. 8.The first two surfaces are respectively defined as The parameters for Dupin's Cyclide are c 1 = 2, c 2 = 1.9, c 3 = √ 0.39 and c 4 = 1.These two surfaces are used for example in [28] to treat the same problem using a compact RBF-FD method.For the case of frog we are only given a cloud X containing 5900 points and approximate normals at each point.A subset of 530 points out of X is used as patch centers and a varying patch radius strategy is applied so that each patch contains exactly 41 interpolation points.The D-RBF-PU method is applied for discretizing the spatial variables and the semi-implicit backward differential formula of order two (SBDF2) [2] is employed for time integration.At the first time step SBDF1 is used to bootstrap SBDF2.The final linear system was efficiently solved by the GMRES algorithm with a standard incomplete LU preconditioner.It is known that as a necessary condition for time stability the spectrum of the Laplace-Beltrami differentiation matrix should fall in the left half plane.Experiments show that if the fill distance h X, is small enough then it appears that all of the eigenvalues do in fact lie in the left half plane.In few cases, small positive real parts of order 10 −16 − 10 −8 appear in the spectrum.This might be bypassed by applying a small diagonal increment on the discrete Laplacian matrix.Fig. 9 shows the eigenvalue patterns on different surfaces for some specific number of nodes which are used for numerical simulations in the following.The WE kernel is used in this and all experiments from here on.To model the Turing patterns we consider a particular form of (7.1) with for known parameters α, β, τ 1 , τ 2 and γ .The system (7.1) with the above reaction terms is known as the Brusselator equation where (u, v) = (0, 0) is a unique steady state for γ = −α [4,8,21,54].In the absence of diffusion terms (δ 1 = δ 2 = 0) and with certain parameter values the steady state (0, 0) is linearly stable while in the presence of diffusion it is unstable [51].Depending on the values of other coefficients this instability exhibits different types of patterns when the steady state is reached [4].Here we illustrate spots and stripes patterns.
In our experiments the initial condition is given by u = v = 0 except on a bound of width 0.1 around the equator of each surface where the values of u and v are chosen to be uniformly randomly distributed in [−0.5, 0.5].We use the RandStream function in Matlab to generate such values [19].We set t = 0.01 in SBDF2 and execute the simulation until a steady state solution is reached.To obtain spot patterns we set δ 2 = 0.0045, δ 1 = 0.516δ 2 , τ 1 = 0.02, τ 2 = 0.2, α = 0.899, β = −0.91,and γ = −α.For stripe patterns we use δ 2 = 0.0021, τ 1 = 3.5, τ 2 = 0 and keep the values of other parameters unchanged.These values are all motivated from the values used in [8,19].In Fig. 10 the steady state values of u, the resulting spot and stripe patterns, on RBC, torus, tooth and frog are shown.These results are in good agreement with those calculated in [19] and other sources.
Our next application of the D-RBF-PU method for solving a reaction-diffusion equation on surfaces concerns the equation (7.1) with reaction terms This is a Fitzhugh-Nagumo type system which is a simple model for dynamics of many excitable media.Here, u and v can be considered as some chemical concentrations or as membrane potential and current.The known parameters α, a and b govern the reaction kinetics.Small values of parameter α result in two separated regions on the surface each with values u = 0 and u = 1 with a thin interface separating these two regions.The solution of this system is a spiral wave evolving on the surface.Here we present some results on the sphere and on the Dupin's Cyclide.Results are presented in Fig. 11   on both surfaces.Motivating by [21], parameter values a = 0.75, b = 0.02, α = 0.02, δ 1 = 2.5(2π /50) 2 and δ 2 = 0 are assigned.We observe from Fig. 11 that the initial density u evolves into approximately symmetric spiral waves that remain intact but meander around the surfaces throughout the simulation time.

Conclusion
In this work we have presented a novel localized RBF method based on partition of unity for discretizing a large class of PDE problems on manifolds.The new D-RBF-PU method yields a numerical algorithm that is more efficient, more accurate and faster than standard RBF-PU and RBF-FD methods.It can be easily applied to any smooth orientable surface that is discretized with a set of scattered nodes and with approximations to the normal to the surface at each of the nodes.As an important property we can mention that in the new method surface PDE operators are not required to be applied on rational PU weight functions.This property simplifies the algorithm, highly, and allows the use of some discontinuous weights to develop a sparser (than RBF-PU) and faster (than RBF-FD) localized RBF technique.The consistency bounds were obtained for finitely smooth basis functions whose Fourier transforms decay only algebraically.Moreover, the stability, accuracy and efficiency properties of the scheme were tested on a general PDE defined on different surfaces.
In additions, we have shown how the new method can be easily adaptable to reaction-diffusion equations on surfaces.In particular, the method was applied on the Brusselator equation and on a Fitzhugh-Nagumo type system to obtain the formation of Turing patterns and spiral waves in an excitable media, respectively.For advancing in time, the semi-implicit BDF2 scheme was employed for time discretization.The final linear system was efficiently solved by the GMRES algorithm with a standard incomplete LU preconditioner.Experiments showed that the final matrices on different manifolds satisfy the conditions that ensure eigenvalue stability for temporal integration.
Finally, we stress that the method can be used to solve various types of PDEs on surfaces.An extension to convectiondiffusion problems is left for a future work.Moreover, the proposed method can be extended for solving a boundary value problem defined on a surface subdomain with a sufficiently smooth boundary.Imposing the boundary condition using the proposed RBF based method is straightforward because the boundary can be viewed as another submanifold with codimension 2.

Declaration of competing interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Davoud Mirzaei reports financial support was provided by Institute for Research in Fundamental Sciences.Davoud Mirzaei reports financial support was provided by Iran National Science Foundation.

Fig. 1 .
Fig. 1.Overlapping circular coverings on the unit sphere and a torus.

5 .
Partition of unity approximation Let { } N c =1 be an open and bounded covering of that means all patches are open and bounded and = N c =1 .

Fig. 3 .
Fig. 3. Accuracy (left vertical axis) and sparsity (right vertical axis) with respect to the overlap constant C ovlp in D-RBF-PU method.First row: sphere, second row: torus, third row: RBC.Left column: the smooth PU weight, right column: the constant-generated PU weight.

Fig. 4 .
Fig.4.Errors and convergence orders of D-RBF-PU with smooth and constant-generated weights and comparison with RBF-FD on sphere, torus and RBC surfaces using Mátern (first row) and Wendland (second row) kernels.The repeated interpolation and differentiation approach is applied to discretize the Laplace-Beltrami operator.Computational convergence orders are written alongside the figure legends.

Fig. 5 .
Fig.5.Errors and convergence orders of D-RBF-PU with smooth and constant-generated weights and comparison with RBF-FD on sphere, torus and RBC surfaces using Mátern (first row) and Wendland (second row) kernels.The standard approach is applied to discretize the Laplace-Beltrami operator.Computational convergence orders are written alongside the figure legends.

Fig. 6 .
Fig. 6.The behavior of the stability constant C S for D-RBF-PU and RBF-FD methods on different surfaces.Results show a perfectly stable final linear system in each case.

Fig. 7 .
Fig. 7.The CPU times used for setting up and solving the final linear system on RBC surface at different numbers of points.The rates are the same but D-RBF-PU is approximately 10 times faster.

Fig. 9 .
Fig. 9.The spectrum of the discrete Laplace-Beltrami with D-RBF-PU method and discontinuous PU weight on different surfaces.All eigenvalues fall on the left half plane.

Fig. 10 .
Fig. 10.Quasi steady Turing spots and stripe patterns resulting from solving the Brusselator equation.In all plots, red corresponds to a high concentration of u and blue to a low concentration.(For interpretation of the colors in the figures, the reader is referred to the web version of this article.)

Fig. 11 .
Fig. 11.Spiral wave patterns resulting from solving the Fitzhugh-Nagumo equation.The first row shows the excitation variable u on the sphere at times t = 0, 15, 18 and 21 (from left to right).The second row shows u on Dupin's Cyclide at times t = 0, 42, 45 and 48 (from left to right).In all plots, black corresponds to a high concentration of u and white to a low concentration.