An RBFFD closest point method for solving PDEs on surfaces
Abstract
Partial differential equations (PDEs) on surfaces appear in many applications throughout the natural and applied sciences. The classical closest point method (Ruuth and Merriman, J. Comput. Phys. 227(3):19431961, [2008]) is an embedding method for solving PDEs on surfaces using standard finite difference schemes. In this paper, we formulate an explicit closest point method using finite difference schemes derived from radial basis functions (RBFFD). Unlike the orthogonal gradients method (Piret, J. Comput. Phys. 231(14):46624675, [2012]), our proposed method uses RBF centers on regular grid nodes. This formulation not only reduces the computational cost but also avoids the illconditioning from point clustering on the surface and is more natural to couple with a grid based manifold evolution algorithm (Leung and Zhao, J. Comput. Phys. 228(8):29933024, [2009]). When compared to the standard finite difference discretization of the closest point method, the proposed method requires a smaller computational domain surrounding the surface, resulting in a decrease in the number of sampling points on the surface. In addition, higherorder schemes can easily be constructed by increasing the number of points in the RBFFD stencil. Applications to a variety of examples are provided to illustrate the numerical convergence of the method.
keywords:
closest point method, embedding method, radial basis functions, finite differencescomma,sort
1 Introduction
Many applications in the natural and applied sciences require the solution of partial differential equations (PDEs) on surfaces. Image processing applications include the placement of an image on a surface Turk (1991), the restoration of a damaged pattern on a surface Bertalmio et al. (2001) and the segmentation and the denoising of images on surfaces Tian et al. (2009); Biddle et al. (2013). In biology, applications include the formation of patterns on animal coats Murray (2001) and the wound healing process Olsen et al. (1998). In computer graphics, applications are found in the topic of real time fluid visualization on surfaces Auer et al. (2012).
Various numerical methods have been developed to approximate the solution of PDEs on surfaces. These include methods applied on parametrized surfaces, on triangulated surfaces and on surfaces embedded in a higher dimensional space. Solution of PDEs on parametrized surfaces can be efficient for surfaces where a parametrization is possible Lui et al. (2005); Floater and Hormann (2005), however a parametrization of a surface often leads to distortions of the surface and singularities Floater and Hormann (2005). Triangulated surfaces avoid these issues. Finite difference methods can be applied to solve PDEs on triangulated surfaces Turk (1991), but there are difficulties in the calculation of geometric quantities, including the normal vector and the curvature of a surface Bertalmıo et al. (2001). On the other hand, finite element methods on triangulated surfaces are effective in solving parabolic or elliptic PDEs Dziuk and Elliott (2007). Methods using surfaces embedded into a dimensional space extend the surface PDE in the embedding space and solve the extended PDE using standard Cartesian methods.
A popular method employs a level set representation of surfaces and a projection operator to solve surface PDEs Greer (2006); Flyer and Wright (2009). Typically, the computational domain consists of points in a neighborhood of the surface. This may lead to the introduction of artificial boundary conditions at the boundary of the computational domain which can degrade the accuracy of the method Cheung et al. (2015). Meshfree approximations using radial basis functions (RBFs) are also becoming popular within the embedded surfaces class Fuselier and Wright (2013); Shankar et al. (2015).
The closest point method Ruuth and Merriman (2008) is an embedding method which uses a closest point representation of the surface to solve PDEs on surfaces. In the classical formulation Ruuth and Merriman (2008); Macdonald and Ruuth (2008), the discretization is carried out in a neighborhood of the surface using standard finite difference schemes and barycentric Lagrangian interpolation. The implicit closest point method was introduced in Macdonald and Ruuth (2009) to provide a stable approximation of the LaplaceBeltrami and other higherorder surface differential operators. Application of the implicit closest point method to the solution of eigenvalue problems appears in Macdonald et al. (2011). See also Marz and Macdonald (2012) for a study of the theoretical foundation of the closest point method.
The extension via the closest point mapping is also used as part of the development for other methods for solving surface PDEs. In Piret (2012), the author uses the closest point mapping to derive an RBF method for solving surface PDEs. The method gives a highorder approximation to the solutions of surface PDEs in a variety of examples. See also Cheung et al. (2015) for a related RBF method that carries out a local approximation of surface differential operators to solve PDEs on folded surfaces. In addition, the extension via the closest point mapping is used in the computation of integrals over curves and surfaces Kublik and Tsai (2016) and in the solution of PDEs on closed, smooth surfaces using volumetric variational principles Chu and Tsai (2017).
In this paper, we propose a new method using a closest point representation of a surface and finite difference stencils derived from radial basis functions (RBFFD). Notably, the use of RBFFD leads to a method (RBFCPM) that evaluates derivatives on the surface, rather than in the embedding space. This eliminates the interpolation step in the evaluation of derivatives, thereby eliminating a potential source of error and computational cost. In addition, standard RBFFD and global RBF methods for surface PDEs may suffer illconditioning due to small separating distances in the surface points Wendland (2004). Our method uses RBFFD stencils on regular Cartesian grid nodes, thus allowing irregular collocation points on the surface and avoiding the illconditioning that may arise due to point clustering on the surface. Due to the regularity of the RBFFD stencil, the collocation matrix associated with the calculation of the RBFFD weights is independent of the surface, and its inverse can be accurately calculated locally. Due to repeated patterns on the RBF centers, only a small number of collocation matrix inverses need to be calculated, thus reducing the computational cost over existing RBF methods.
In our method, secondorder accuracy in can be achieved with smaller computational domains and fewer points on the surface relative to the classical closest point method. Furthermore, higherorder schemes are obtained simply by increasing the number of points in the finite difference stencil.
The paper unfolds as follows. In Section 2, we review the classical closest point method and RBF approximation. Section 3 gives our new method and studies the selection of parameters and the computational domain. Section 4 considers the performance of the method using a variety of convergence studies and numerical examples in two and three dimensions. Finally Section 5 concludes and explores potential future work.
2 Numerical methods review
2.1 The closest point method
The classical closest point method Ruuth and Merriman (2008) is a simple numerical method for approximating the solution of PDEs on surfaces. In this section, we review the method and its components.
2.1.1 Surface representation
To begin, the closest point to the surface is defined:
Definition 1
Let be a surface and be some point in the embedding space . Then,
is the closest point of to the surface .
In a neighborhood of the surface, will be smooth for a smooth surface Marz and Macdonald (2012).
To discretize, a Cartesian grid is introduced in the embedding space, over a neighborhood of the surface. Typically, this neighborhood includes all grid nodes whose Euclidean distance to the surface is less than or equal to some constant . Following a common convention (e.g., Leung and Zhao (2009)), we refer to this localized computational domain as the computational tube, and the corresponding radius as the computational tube radius. To avoid introducing discontinuities into , the computational tube radius should satisfy , where is an upper bound on the curvatures of Chu and Tsai (2017). The grid points and their closest points together form a closest point representation of the surface .
The method used to determine the closest point function depends on the type of surface under consideration. For simple surfaces, such as the sphere and the torus, an analytical formula for the closest point function is available, and the preferred approach is to simply evaluate the formula. On the other hand, for ellipsoids, the Möbius strip, and many other interesting shapes, the surface may be given in parameterized form. Here, standard numerical optimization techniques can be applied to find the closest point on the surface (cf. Merriman and Ruuth (2007)). Finally, we consider surfaces in triangulated form. In this case, we follow Macdonald and Ruuth (2008) and loop over the list of triangles. For grid nodes near a triangle (specifically, nodes that are within a Euclidean distance of ), we compute and store the closest point on . For each grid node, the closest point over all stored possibilities is the closest point on the surface. See Macdonald and Ruuth (2008) for further details on this procedure. Another approach which follows the causality of the eikonal equation can also be applied Tsai (2002).
2.1.2 The equivalence principles
In the closest point method, we do not solve the surface PDE problem on the surface directly. Instead, we solve a suitable differential equation defined over the computational domain. The key property of this differential equation (the embedding equation) is that its solution on the surface must agree with the solution to the original surface PDE. Values off the surface do not directly give the solution to the PDEonsurface problem.
To form the embedding equation, surface derivatives are replaced with closest point operators and standard Cartesian derivatives. The foundation of this procedure consists of two principles Ruuth and Merriman (2008); Cheung and Ling (2017): the equivalence of gradients and the equivalence of divergence.
Principle 1
Let be any function on that is constant along normal directions of . Then, at the surface, intrinsic gradients are equivalent to standard gradients, .
Principle 2
Let be any vector field on that is tangent to and tangent to all surfaces displaced by a fixed distance from . Then, at the surface, .
General surface differential operators can be replaced with the corresponding Cartesian differential operators by combining the two principles. Of particular interest to us is the composition of the divergence and gradient operators: Let be any function on that is constant along normal directions of . Then, on the surface, the LaplaceBeltrami operator is equivalent to the standard Laplacian operator, for all . This property is referred to as the equivalence of the Laplacian. As a consequence, heat flow intrinsic to a surface can be approximated by alternating constantalongnormal extension with standard heat flow in the underlying embedding space. Other, similarly straightforward combinations of the gradient and divergence principles lead to replacements for nonlinear diffusive flows such as curvature motion intrinsic to a surface Merriman and Ruuth (2007); Ruuth and Merriman (2008). Higherorder operators can also be approximated in the embedding space, although additional extension operators may be needed. See, for example, Macdonald and Ruuth (2009) where two extension operators are used as part of the approximation of fourthorder operators.
2.1.3 The closest point method
Evolution of the embedding equation may be carried out in a similar fashion, to yield the explicit closest point method Ruuth and Merriman (2008). Specifically, given a closest point representation of a surface , the explicit closest point method alternates between the following two steps:

Closest point extension. Carry out a constantalongnormal extension of to yield by for each in the tubular computational domain .

Evolution. The PDE is solved on the tubular computational domain in the embedding space for one time step (or one stage of a RungeKutta method).
Note that the closest point extension step is an interpolation step since is not necessarily a grid point. In Ruuth and Merriman (2008), barycentric Lagrange interpolation is used with polynomial degree , where is the order of finite differences schemes and is the order of the derivatives. Localization of the computation is accomplished by computing over a computational tube surrounding the surface Ruuth and Merriman (2008); Macdonald and Ruuth (2009). For secondorder finite differences and secondorder derivatives, it is sufficient to choose a computational tube radius of where
(1) 
in a dimensional embedding space Ruuth and Merriman (2008).
2.2 Global Radial Basis Function (RBF) approximation
RBF approximation is a powerful tool for approximating smooth functions on a variety of geometries. Following Fasshauer (2007); Fornberg and Flyer (2015a), given an RBF (see Table 1 for some RBF options) and a set of scattered points called RBF centers, the RBF interpolant has the form
(2) 
with coefficients . For the RBF interpolation of any smooth function , the coefficients can be found by interpolation conditions at , i.e. by solving the symmetric linear system
(3) 
for , where and for in an orderly sense.
Name of RBF  Abbreviation  Definition 

Smooth RBFs  
Gaussian  GA  
Multiquadratic  MQ  
Inverse multiquadratic  IMQ  
Inverse quadratic  IQ  
Piecewise smooth RBFs  
Cubic  CU  
Thin plate spline  TPS 
Then, we can use the interpolant to approximate the derivatives of . If is a differential operator, then the quantity at some point can be approximated as
(4) 
in which acts upon the variable in the basis function . In matrix form, we can express (4) as
(5) 
where . This expression provides a radial basis function pseudospectral method, which can be easily localized to give a RBF finite difference discretization (RBFFD) for the differential operator evaluated at the data site . Some compact RBFFD stencils can be found in Fornberg and Flyer (2015b) for applications in geosciences. Details on RBFFD stencils that use a given number of nearest neighbors are found in Flyer et al. (2016).
3 A closest point method for solving PDEs on surfaces using RBFFD
In this section, we introduce an explicit RBF closest point method for solving PDEs on surfaces. Our method uses RBFFD stencils that consist of closest neighboring grid points. As part of our method, we provide an approach for the calculation of a computational tube around the surface.
For illustration purposes, consider the heat equation
(6) 
intrinsic to a surface , where is a function defined solely on the surface. Using a set of specially constructed surface data points , our aim is to design an efficient, accurate, and robust finite difference scheme in order to spatially discretize equation (6), i.e.
where is a differentiation matrix that takes the vector of function values to the approximated values of at .
In the literature, this can be done by the orthogonal gradients method Piret (2012) and projection methods Fuselier and Wright (2013) using RBFs, in which the geometry of plays a key role in the numerical stability. In the orthogonal gradients method, the RBF centers are constructed by extending the surface points in the embedding space in the normal direction. The geometry of the surface and the points determine that of , which should ideally be quasiuniform. Projection methods work solely on with the corresponding analysis carried out in Sobolev spaces on . Like other typical kernel approximation theories, convergence comes when gets dense, i.e. as the fill distance . Both approaches involve solving interpolation problems, whose conditioning depends on the minimum separating distance of Wendland (2004). Thus, it is common to require that the mesh ratio of the surface points is bounded. In short, quasiuniform points need to be used on surfaces, which may not be an easy task.
Finite difference schemes, on the other hand, work on regular grids with mesh ratio exactly equal to 1. Yet, extra work is required to apply finite differences to surfaces in general (e.g., the original closest point method Ruuth and Merriman (2008)). We propose a new RBF kernel based formulation to get the best of both worlds.
3.1 Description of the method
We start with a collection of Cartesian grid points on a small tubular domain containing the surface . Then, we define surface data points via to form a set ; see Figure 1 for a schematic demonstration. Applying the equivalence of the Laplacian property yields the relation
which holds for any , and where we denote the constantalongnormal extension of by .
We now deploy the methodology locally to obtain an RBFFD approximation. For each surface point for some and , let denote the nearest neighborhood of . Locally, we take the –local interpolant of using basis function at centers , denoted by below, as the function in Section 2.2. Then, (5) gives an approximation scheme
In other words, the nonzero RBFFD weight is given as the product of a row vector and the inverse matrix of for . Using for , we can assemble the RBFFD matrix such that
(7) 
In addition, we also construct an RBFFD projection matrix such that
(8) 
simply by replacing with the identity map in the computation. The matrix will appear later as part of our timestepping scheme.
Note that all approximations are done in the embedding space, which is independent of the geometry of . By using a regular , we obtain collocation matrices that depend on the geometry of , but are independent of the surface . Thus, we can precompute all , and the evaluation of each RBFFD weight requires only matrixvector multiplications. As a consequence, we can employ accurate and expensive solvers to compute the matrices without harming the overall performance of the proposed method. This greatly improves computational efficiency over the existing methods in which no matrix inverse can be reused due to the lack of repeated pattern in data point geometry.
The RBF finite difference stencil used in this paper consists of the closest grid points , to each surface point . Methods that use RBFFD stencils of a given number of nearest neighbors on scattered nodes require the use of quasiuniform nodes on the surface Shankar et al. (2015). By using RBFFD stencils on the grid nodes, we avoid the illconditioning of the RBFFD matrices that arises from the clustering of nodes on the surface. Figure 1 shows an example of an RBFFD stencil that consists of the closest grid points to a particular surface point (displayed using a squared red dot).
We are now ready to temporally discretize the heat equation
intrinsic to the surface . Using the forward Euler scheme with spatial discretization at as above, we have
(9) 
for . Recall that in our proposed setup is regular, whereas could be highly nonuniform. It is desirable to work on the nodes.
Let be the approximated values of on and at time . This timestepping scheme is initialized using the initial condition . It is equivalent to consider the discrete equation of the constantalongnormal extended function, and (9) becomes
(10) 
We use (7) and (8) to approximate the right hand side from the stored approximated solution values . Note that we do not use pointwise projection, i.e., , because contains discretization and approximation errors. Using (8) to approximate introduces some averaging into the approximation and increases numerical stability. Indeed, Figure 2 shows the eigenvalues of the discretization of equation (10) on the unit circle for points, a grid size of and a time stepsize of . The scheme that uses pointwise projection leads to an unstable system (eigenvalues larger than ) whereas a stable approximation can be achieved using the projection operator.
By design, we have , whose approximated values will be used to define . With all of the above considered, the approximate solution can be updated from time to by
(11) 
In (11), the exterior nodes are also used implicitly in the construction of RBFFD matrices and defined in (7) and (8) as RBF centers. In other words, the proposed RBFCPM runs solely based on RBF interpolations with the regularly placed nodes. To evaluate the numerical approximation, say at for simplicity, one can evaluate . Other timestepping schemes are also possible. In our experiments, we use the thirdorder, threestage SSP RungeKutta scheme Shu and Osher (1988) in advectiondominant problems due to its good linear stability along the imaginary axis. For simplicity, and to differentiate from other methods, we shall refer to our RBF discretization as the RBFCPM.
Given a collection of Cartesian grid points in a small tubular domain that contains the surface , the algorithm of the RBFCPM for a time dependent PDE consists of the following steps:

Compute the set of surface points via the closest point representation of the surface : , for , .

Solve the surface PDE using an explicit timestepping scheme, e.g. (11), until the final time.
Our method allows precomputation for solving local linear systems; provided that two local neighborhoods and share the same geometrical arrangement, the corresponding interpolation matrices are identical, i.e, . Therefore, an expensive but accurate method can be used.
3.2 Parameters
The convergence order of the RBFFD schemes is limited by the smoothness of the employed kernel and the geometry of the RBF centers Davydov and Schaback (2016). In this paper, we employ Gaussian RBFs so that the smoothness of the kernel will not be a limiting factor. To avoid any potential problem of illconditioning, we use the stable RBFGA method which provides an accurate and stable algorithm and is a cheaper stabilization method over RBFQR Fornberg et al. (2013). These stabilization techniques are independent of the condition number of the matrix (see Section 3.1) using a direct calculation using Gaussian RBFs Fornberg et al. (2011, 2013).
There are two parameters appearing in the RBFCPM. These are the shape parameter of the Gaussian RBFs and the number of points in the stencil used locally for the RBF interpolation. While this section considers the dependence of the numerical method on both parameters, our emphasis will be on the number of points . The parameter was found to have little effect on our results.
We consider two test problems to measure the error of the discrete LaplaceBeltrami operator in comparison to the exact. The first test problem (P1) approximates the LaplaceBeltrami operator applied to the function on the unit circle. The relative error in this problem can computed using the known, exact solution . In our second problem (P2), the LaplaceBeltrami operator is applied to the function on the unit sphere
Here, the relative error can be computed using the known, exact solution .
We apply the RBFFD method in Equation (7) to test problems P1 and P2 for various mesh spacings and selected values. See Figure 3 for the corresponding max norm relative errors. We find that the observed order of the method increases as the number of points in the finite difference stencil increases. Due to the regularity of the grid spacing, and the smoothness of the kernel, our choice of parameter has little effect on the results. In particular, the observed orders of convergence for , and are the same for problems P1 and P2. For this reason, we simply choose in the numerical experiments presented in Section 4.
3.3 Computational tube
The RBFFD stencils used in the RBFCPM, i.e., for , are formed using the nearest neighboring regular grid points with a predetermined spacing . For such problems, the matrixbased formulation of the closest point method Macdonald and Ruuth (2009) can be used to obtain the minimalsized computational tube. In this paper, we can take a simpler approach to identify the computational tube by identifying a sufficiently large tube radius .
Consider the Gauss circle problem Weisstein (2004). In its standard form, it is posed as follows: Find the number of integer lattice points inside a circle with radius centered at the origin. We shall consider a related formulation: Find the number of ordered pairs , with integers , such that
where the radius of the circle is chosen as , for a fixed integer . Generalizations to higher dimensions are also available. In three dimensions, the problem uses a sphere centered at the origin. In this case, we find the number of ordered triplets , with integers , such that
where the radius of the sphere is , for an integer .
The integer solutions of the Gauss circle problem and its generalization to three dimensions are given by sequences A057655 (two dimensions) and sequences A117609 (three dimensions) in the OnLine Encyclopedia of Integer Sequences oei (2011). Table 2 shows some of the integer solutions for the Gauss circle problem in two and three dimensions.
(2D)  (3D)  

0  1  1 
1  5  7 
2  9  19 
3  9  27 
4  13  33 
5  21  57 
For a stencil that uses the closest grid nodes to a surface point, a circle can be constructed centered at the surface point that contains these grid nodes. In order to construct a computational tube around the surface using the Gauss circle problem, we need to find a sufficiently large circle independent of the position of the surface point relative to the surrounding grid nodes. In the optimal case, the surface point lies on a grid node; see Figure 4 (left). In such an occurrence, the Gauss circle problem can be applied directly, scaled properly with the grid size . Otherwise, the surface point does not lie on a grid node. In order to use the Gauss circle problem, the distance between the surface point and its closest grid node needs to be added to the radius of the Gauss circle problem. The worst case appears in Figure 4, where the surface point lies midway between all four surrounding grid nodes.
To guarantee a sufficiently large circle/sphere radius, we must allow for the worst case. This leads to a relation between the computational tube radius (corresponding to the circle/sphere radius ) and the number of points in the RBFFD stencil; see Table 3.
(2D)  (2D)  (3D)  (3D) 

(  9  (  27 
(2+  13  (2+  33 
21  57  
25  81  
93 
4 Numerical experiments
In this section, we test our method, the RBFCPM, on a number of examples in two and three dimensions. For problems involving the LaplaceBeltrami operator, we choose in two dimensions and in three dimensions, which provides a second order approximation for smooth solutions (cf. Figure 3). The radius of the computational tube is set according to the values specified in Table 3. Finding the optimal RBF shape parameter is out of the scope of this paper, and thus we set (unscaled Gaussian RBFs are used). Unless stated otherwise, forward Euler with a time stepsize is used for the discretization of the time derivatives.
4.1 Examples in two dimensions
First, we perform numerical experiments for applications in two dimensions to test the convergence of the proposed method.
4.1.1 Heat equation on a circle
In our first experiment, we consider the heat equation
intrinsic to the unit circle . Following Ruuth and Merriman (2008), for an initial profile , the exact solution is
Using an analytic closest point representation of the unit circle, the surface heat equation is discretized and solved using Equation (11). Table 4 shows the relative errors as well as the convergence rates for different grid sizes and number of points on the computational domain. The convergence rate here appears to be at least secondorder. Using the original closest point method, the number of points in the computational tube that are required for a second order approximation with is , whereas our method uses points. This corresponds to a reduction of .
Rel. error ()  Conv. rates  

0.2  172  7.15   
0.1  336  1.22  2.55 
0.05  688  2.23  2.46 
0.025  1376  5.15  2.11 
0.0125  2708  1.35  1.93 
0.00625  5464  3.15  2.10 
4.1.2 Heat equation on a semicircle
Next, we consider the surface heat equation on the unit semicircle with homogeneous Dirichlet boundary conditions. Given an initial profile , the exact solution is
Following Macdonald et al. (2011), we introduce a modified closest point mapping which equals for points that map to the interior of the semicircle. Grid nodes that satisfy are called ghost points . At such points, the function is extended by . Figure 5 shows the computational tube used in this example.
Table 5 presents the relative errors as well as the convergence rates for different grid sizes and number of points on the computational domain. Secondorder convergence is observed. Using the original closest point method, the number of points in the computational tube that are required for a second order approximation at is , whereas the RBFCPM uses points.
Rel. error ()  Conv. rates  

0.2  110  7.38   
0.1  192  1.14  2.69 
0.05  368  2.12  2.43 
0.025  712  5.02  2.08 
0.0125  1378  1.34  1.91 
0.00625  2756  3.13  2.09 
4.1.3 Advection equation on an ellipse
The next example is the advection equation on an ellipse. Following Ruuth and Merriman (2008), the equation
with being the arclength, is imposed on an ellipse with major axis along the yaxis and minor axis along the xaxis. By Ruuth and Merriman (2008), application of the closest point principles to the surface PDE leads to the embedding PDE
with
For an initial profile , the exact solution for subsequent times is
where is the length of the perimeter of the ellipse. Using a parametrization for the ellipse, the closest point representation is calculated via optimization techniques. Due to the generic centered nature of the RBFFD stencils used for approximating the firstorder derivatives, the TVDRK3 scheme Gottlieb and Shu (1998) is chosen for the time discretization with a time stepsize . In this example, the computational tube radius is with points in the finite difference stencil. Table 6 shows the error at the final time and the estimated order of convergence of the method for various grid spacings and number of points in the computational domain. Secondorder convergence is observed. The number of points in the computational tube that are required for a second order approximation at using the original closest point method is . Our method uses fewer points.
Rel. error ()  Conv. rates  

0.2  136  8.99   
0.1  272  9.80  3.20 
0.05  552  2.25  2.12 
0.025  1080  5.59  2.01 
0.0125  2168  1.40  2.00 
0.00625  4332  3.51  1.99 
0.003125  8652  8.75  2.00 
4.1.4 Advectiondiffusion equation on an ellipse
Our next example considers advectiondiffusion on an ellipse. The equation
with being the arclength, is imposed on an ellipse with major axis (along the yaxis) and minor axis (along the xaxis). Similar to the previous example, application of the closest point principles to the surface PDE gives
with
For an initial profile , the exact solution has the form
where is the length of the perimeter of the ellipse. Table 7 shows the relative error of the approximate solution compared to the exact as well as the estimated order of convergence. Secondorder convergence is observed.
Rel. error ()  Conv. rates  

0.2  172  9.66   
0.1  348  1.34  2.85 
0.05  692  4.88  1.46 
0.025  1384  1.25  1.97 
0.0125  2792  2.72  2.20 
0.00625  5552  6.86  1.99 
0.003125  11100  1.68  2.03 
4.2 Examples in three dimensions
Next, we apply our proposed method to examples in three dimensions.
4.2.1 Heat equation on a sphere
For our first three dimensional example, consider the heat equation
on the unit sphere . For an initial profile , the exact solution for all times is
An analytic closest point representation of the unit sphere is used and the surface heat equation is discretized and solved using Equation (11). Table 8 shows the relative errors as well as the convergence rates for different grid sizes and number of points in the computational tube. Secondorder convergence is observed. Using and a second order finite difference discretization in the original closest point method leads to points in the computational domain, while the RBFCPM uses points. This corresponds to a reduction of .
Rel. error ()  Conv. rates  

0.2  2240  8.18   
0.1  8072  2.21  1.89 
0.05  31416  5.42  2.03 
0.025  125216  1.36  1.99 
0.0125  498392  3.40  2.00 
4.2.2 Advection on a torus
In this example, the solution of the advection equation on a torus is approximated. Following Greer (2006), for a torus defined as
the advection equation is given by
For an initial profile
where
the exact solution at time is
Rel. error ()  Conv. rates  

0.1  11392  1.76   
0.05  45464  2.99  2.56 
0.025  181480  4.88  2.62 
0.0125  725200  9.52  2.36 
0.00625  2901248  1.80  2.40 
In this example, an analytic closest point representation of the torus is used. The computational tube radius is and a stencil with points is chosen. The stable TVDRK3 scheme is chosen for the discretization in time with stepsize . Table 9 shows the error and the convergence rate of the method for various grid sizes . Convergence is at least secondorder. The number of points in the computational tube that are required for a second order approximation at using the original closest point method is . Our method uses fewer points.
4.2.3 Image denoising on a sphere
Our next example concerns image denoising for a textured image on the unit sphere. Following Biddle et al. (2013), we apply the PeronaMalik model to denoise surface images. The equation has the form
where is the diffusion coefficient given by
and is a coefficient. The parameter and the final time of the computation control the denoising of an image.
To construct the initial image, we add Gaussian noise (zero mean and 0.2 standard deviation) to the image of two birds. The noisy image is scaled to the interval . Figure 6 shows the results after denoising the image using and time steps with and . We find this choice for the parameter is sufficiently small to preserve edges in the denoised image. In this example, the computational tube radius is chosen to be with points in the stencil. The number of points in the image is .
4.2.4 Reactiondiffusion systems
Our final example evolves the GrayScott reactiondiffusion model Gray and Scott (1984) on a triangulated surface. The Gray Scott model describes the chemical reaction
where , and are chemicals. The corresponding surface model has the form
where , are the concentrations of the chemicals, , are the diffusion rates, is the conversion rate from to and is the feed rate of .
For parameter choices of and , a variety of patterns are observed as and are varied. Figure 7 shows two of these patterns on the surface of the Stanford Bunny Turk and Levoy (1994). The closest point representation to the triangulated surface is calculated according to the method described in Section 2.1.1. In this example, the final time is 15000 and the time stepsize is for a spatial grid size of . See Figure 7 for patterns arising for two different choices of the parameters and McGough and Riley (2004).
5 Summary
In this paper, an explicit closest point method is introduced that uses finite differences derived from radial basis functions (RBFFD). In our method, an RBFFD approximation of surface derivatives is formed using the grid points closest to a surface point. Localization of the computation is accomplished by computing over a tube whose radius is obtained from the solution to the Gauss circle problem. An advantage of our algorithm relative to the standard finite difference CPM is a reduction of the computational tube radius, leading to the reduction of the grid points in the computational domain and their corresponding closest points on the surface. Also, higherorder schemes are easily constructed by increasing the number of points in the finite difference stencil. When compared to RBF methods, our algorithm does not require quasiuniform distribution of points on the surface. In addition, the repeated patterns in our computational geometry allows us to use an algorithm to invert (a small number of) collocation matrices, thereby reducing computational cost over other existing methods. Numerical experiments are provided to validate the method for different types of PDEs on surfaces.
The RBFFD discretization introduced in this paper solves surface PDEs using explicit time stepping methods. Implicit RBFFD schemes which allow for large time steps for stiff problems are also needed, and are a focus of our current work. Related to this, the approximation of the eigenvalues of surface operators using the RBFCPM method is particularly interesting (cf. Macdonald et al. (2011)). Another focus of our work is the solution of PDEs on moving surfaces. In moving closest point representations, grid node deactivation may occur Leung and Zhao (2009). Methods based on RBFFD discretizations accommodate irregular stencils and are therefore particularly attractive for such problems. For a discussion on the issue of grid node deactivation, and an initial method using the original closest point method, see Petras and Ruuth (2016).
Acknowledgements
The first and third authors were partially supported by an NSERC Canada grant (RGPIN 227823). The first author was partially supported by the Basque Government through the BERC 20142017 program and by Spanish Ministry of Economy and Competitiveness MINECO through BCAM Severo Ochoa excellence accreditation SEV20130323 and through project MTM201569992R BELEMET. The second author was partially supported by a Hong Kong Research Grant Council GRF Grant (HKBU 11528205) and a Hong Kong Baptist University FRG Grant.
References
 Turk (1991) G. Turk, Generating textures on arbitrary surfaces using reactiondiffusion, vol. 25, ACM, 1991.
 Bertalmio et al. (2001) M. Bertalmio, A. L. Bertozzi, G. Sapiro, Navierstokes, fluid dynamics, and image and video inpainting, in: Computer Vision and Pattern Recognition, 2001. CVPR 2001. Proceedings of the 2001 IEEE Computer Society Conference on, vol. 1, IEEE, I–355, 2001.
 Tian et al. (2009) L. L. Tian, C. B. Macdonald, S. J. Ruuth, Segmentation on surfaces with the closest point method, in: Image Processing (ICIP), 2009 16th IEEE International Conference on, IEEE, 3009–3012, 2009.
 Biddle et al. (2013) H. Biddle, I. von Glehn, C. B. Macdonald, T. Marz, A volumebased method for denoising on curved surfaces, in: Image Processing (ICIP), 2013 20th IEEE International Conference on, IEEE, 529–533, 2013.
 Murray (2001) J. D. Murray, Mathematical Biology. II Spatial Models and Biomedical Applications Interdisciplinary Applied Mathematics V. 18, SpringerVerlag New York Incorporated, 2001.
 Olsen et al. (1998) L. Olsen, P. K. Maini, J. A. Sherratt, Spatially varying equilibria of mechanical models: Application to dermal wound contraction, Mathematical biosciences 147 (1) (1998) 113–129.
 Auer et al. (2012) S. Auer, C. B. Macdonald, M. Treib, J. Schneider, R. Westermann, Realtime fluid effects on surfaces using the closest point method, in: Computer Graphics Forum, vol. 31, Wiley Online Library, 1909–1923, 2012.
 Lui et al. (2005) L. M. Lui, Y. Wang, T. F. Chan, Solving PDEs on manifolds with global conformal parametriazation, in: Variational, Geometric, and Level Set Methods in Computer Vision, Springer, 307–319, 2005.
 Floater and Hormann (2005) M. S. Floater, K. Hormann, Surface parameterization: a tutorial and survey, Advances in multiresolution for geometric modelling 1 (1).
 Bertalmıo et al. (2001) M. Bertalmıo, L.T. Cheng, S. Osher, G. Sapiro, Variational problems and partial differential equations on implicit surfaces, Journal of Computational Physics 174 (2) (2001) 759–780.
 Dziuk and Elliott (2007) G. Dziuk, C. M. Elliott, Surface finite elements for parabolic equations, Journal of Computational MathematicsInternational Edition 25 (4) (2007) 385.
 Greer (2006) J. B. Greer, An improvement of a recent Eulerian method for solving PDEs on general geometries, Journal of Scientific Computing 29 (3) (2006) 321–352.
 Flyer and Wright (2009) N. Flyer, G. B. Wright, A radial basis function method for the shallow water equations on a sphere, in: Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, The Royal Society, rspa–2009, 2009.
 Cheung et al. (2015) K. C. Cheung, L. Ling, S. Ruuth, A localized meshless method for diffusion on folded surfaces, Journal of Computational Physics 297 (2015) 194–206.
 Fuselier and Wright (2013) E. J. Fuselier, G. B. Wright, A highorder kernel method for diffusion and reactiondiffusion equations on surfaces, Journal of Scientific Computing 56 (3) (2013) 535–565.
 Shankar et al. (2015) V. Shankar, G. B. Wright, R. M. Kirby, A. L. Fogelson, A radial basis function (RBF)finite difference (FD) method for diffusion and reaction–diffusion equations on surfaces, Journal of scientific computing 63 (3) (2015) 745–768.
 Ruuth and Merriman (2008) S. J. Ruuth, B. Merriman, A simple embedding method for solving partial differential equations on surfaces, Journal of Computational Physics 227 (3) (2008) 1943–1961.
 Macdonald and Ruuth (2008) C. B. Macdonald, S. J. Ruuth, Level set equations on surfaces via the closest point method, Journal of Scientific Computing 35 (23) (2008) 219–240.
 Macdonald and Ruuth (2009) C. B. Macdonald, S. J. Ruuth, The implicit closest point method for the numerical solution of partial differential equations on surfaces, SIAM Journal on Scientific Computing 31 (6) (2009) 4330–4350.
 Macdonald et al. (2011) C. B. Macdonald, J. Brandman, S. J. Ruuth, Solving eigenvalue problems on curved surfaces using the closest point method, Journal of Computational Physics 230 (22) (2011) 7944–7956.
 Marz and Macdonald (2012) T. Marz, C. B. Macdonald, Calculus on surfaces with general closest point functions, SIAM Journal on Numerical Analysis 50 (6) (2012) 3303–3328.
 Piret (2012) C. Piret, The orthogonal gradients method: A radial basis functions method for solving partial differential equations on arbitrary surfaces, Journal of Computational Physics 231 (14) (2012) 4662–4675.
 Kublik and Tsai (2016) C. Kublik, R. Tsai, Integration over curves and surfaces defined by the closest point mapping, Research in the mathematical sciences 3 (1) (2016) 3.
 Chu and Tsai (2017) J. Chu, R. Tsai, Volumetric variational principles for a class of partial differential equations defined on surfaces and curves, arXiv preprint arXiv:1706.02903 .
 Wendland (2004) H. Wendland, Scattered data approximation, vol. 17, Cambridge university press, 2004.
 Leung and Zhao (2009) S. Leung, H. Zhao, A grid based particle method for moving interface problems, Journal of Computational Physics 228 (8) (2009) 2993–3024.
 Merriman and Ruuth (2007) B. Merriman, S. J. Ruuth, Diffusion generated motion of curves on surfaces, Journal of Computational Physics 225 (2) (2007) 2267–2282.
 Tsai (2002) Y.h. R. Tsai, Rapid and accurate computation of the distance function using grids, Journal of Computational Physics 178 (1) (2002) 175–195.
 Cheung and Ling (2017) K. C. Cheung, L. Ling, A kernelbased embedding method and convergence analysis for surfaces PDEs, SIAM J. Sci. Compt., Accepted 2017.
 Fasshauer (2007) G. E. Fasshauer, Meshfree approximation methods with MATLAB, vol. 6, World Scientific, 2007.
 Fornberg and Flyer (2015a) B. Fornberg, N. Flyer, Solving PDEs with radial basis functions, Acta Numerica 24 (2015a) 215–258.
 Fornberg and Flyer (2015b) B. Fornberg, N. Flyer, A primer on radial basis functions with applications to the geosciences, vol. 87, SIAM, 2015b.
 Flyer et al. (2016) N. Flyer, B. Fornberg, V. Bayona, G. A. Barnett, On the role of polynomials in RBFFD approximations: I. interpolation and accuracy, Journal of Computational Physics .
 Shu and Osher (1988) C.W. Shu, S. Osher, Efficient implementation of essentially nonoscillatory shockcapturing schemes, J. Comput. Phys. 77 (2) (1988) 439–471.
 Davydov and Schaback (2016) O. Davydov, R. Schaback, Optimal stencils in Sobolev spaces, arXiv preprint arXiv:1611.04750 .
 Fornberg et al. (2013) B. Fornberg, E. Lehto, C. Powell, Stable calculation of Gaussianbased RBFFD stencils, Computers & Mathematics with Applications 65 (4) (2013) 627–637.
 Fornberg et al. (2011) B. Fornberg, E. Larsson, N. Flyer, Stable computations with Gaussian radial basis functions, SIAM Journal on Scientific Computing 33 (2) (2011) 869–892.
 Weisstein (2004) E. W. Weisstein, Gauss’s circle problem, From MathWorld – A Wolfram Web Resource, URL http://mathworld.wolfram.com/GausssCircleProblem.html, 2004.
 oei (2011) OEIS Foundation Inc., The online encyclopedia of integer sequences, URL http://oeis.org, accessed: 20161025, 2011.
 Gottlieb and Shu (1998) S. Gottlieb, C.W. Shu, Total variation diminishing RungeKutta schemes, Mathematics of computation of the American Mathematical Society 67 (221) (1998) 73–85.
 Gray and Scott (1984) P. Gray, S. Scott, Autocatalytic reactions in the isothermal, continuous stirred tank reactor, Chemical Engineering Science 39 (6) (1984) 1087 – 1097.
 Turk and Levoy (1994) G. Turk, M. Levoy, Zippered polygon meshes from range images, in: Proceedings of the 21st annual conference on Computer graphics and interactive techniques, ACM, 311–318, 1994.
 McGough and Riley (2004) J. S. McGough, K. Riley, Pattern formation in the Gray–Scott model, Nonlinear analysis: real world applications 5 (1) (2004) 105–121.
 Petras and Ruuth (2016) A. Petras, S. Ruuth, PDEs on moving surfaces via the closest point method and a modified grid based particle method, Journal of Computational Physics 312 (2016) 139–156.