An RBF-FD closest point method for solving PDEs on surfaces

# An RBF-FD closest point method for solving PDEs on surfaces

A. Petras L. Ling S.J. Ruuth Department of Mathematics, Simon Fraser University, Burnaby, British Columbia, Canada V5A1S6 BCAM-Basque Center for Applied Mathematics, Bilbao, Basque Country, Spain 48015 Department of Mathematics, Hong Kong Baptist University, Kowloon Tong, Hong Kong
###### 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):1943-1961, [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 (RBF-FD). Unlike the orthogonal gradients method (Piret, J. Comput. Phys. 231(14):4662-4675, [2012]), our proposed method uses RBF centers on regular grid nodes. This formulation not only reduces the computational cost but also avoids the ill-conditioning 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):2993-3024, [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, higher-order schemes can easily be constructed by increasing the number of points in the RBF-FD 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 differences
journal: Journal of Computational Physics\biboptions

comma,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 Laplace-Beltrami and other higher-order 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 high-order 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 (RBF-FD). Notably, the use of RBF-FD leads to a method (RBF-CPM) 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 RBF-FD and global RBF methods for surface PDEs may suffer ill-conditioning due to small separating distances in the surface points Wendland (2004). Our method uses RBF-FD stencils on regular Cartesian grid nodes, thus allowing irregular collocation points on the surface and avoiding the ill-conditioning that may arise due to point clustering on the surface. Due to the regularity of the RBF-FD stencil, the collocation matrix associated with the calculation of the RBF-FD 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, second-order accuracy in can be achieved with smaller computational domains and fewer points on the surface relative to the classical closest point method. Furthermore, higher-order 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,

 cpΓ(z)=argminx∈Γ∥x−z∥2

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 PDE-on-surface 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 Laplace-Beltrami 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 constant-along-normal 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). Higher-order 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 fourth-order 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:

1. Closest point extension. Carry out a constant-along-normal extension of to yield by for each in the tubular computational domain .

2. Evolution. The PDE is solved on the tubular computational domain in the embedding space for one time step (or one stage of a Runge-Kutta 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 second-order finite differences and second-order derivatives, it is sufficient to choose a computational tube radius of where

 γCPM=√(d−1)(p+12)2+(1+p+12)2Δx (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

 s(x)=nZ∑j=1λjϕ(∥x−zj∥) (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

 A(Z,Z)λ=v(Z), (3)

for , where and for in an orderly sense.

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

 Lv(x)≈Ls(x)=nZ∑j=1λjLϕ(∥x−zj∥), (4)

in which acts upon the variable in the basis function . In matrix form, we can express (4) as

 Lv(x)≈B(x,Z)A(Z,Z)−1v(Z), (5)

where . This expression provides a radial basis function pseudo-spectral method, which can be easily localized to give a RBF finite difference discretization (RBF-FD) for the differential operator evaluated at the data site . Some compact RBF-FD stencils can be found in Fornberg and Flyer (2015b) for applications in geosciences. Details on RBF-FD 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 RBF-FD

In this section, we introduce an explicit RBF closest point method for solving PDEs on surfaces. Our method uses RBF-FD 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

 ut=ΔΓu (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.

 [ΔΓu](X)≈Wu(X),

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 quasi-uniform. 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, quasi-uniform 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

 ΔΓu(X)=Δ~u(X),

which holds for any , and where we denote the constant-along-normal extension of by .

We now deploy the methodology locally to obtain an RBF-FD 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

 Δ~u(xj) ≈ Δsj(xj) = B(xj,Zj)A(Zj,Zj)−1~u(Zj) =: wj~u(Zj).

In other words, the nonzero RBF-FD weight is given as the product of a row vector and the inverse matrix of for . Using for , we can assemble the RBF-FD matrix such that

 ΔΓu(X)=Δ~u(X)≈W~u(Z). (7)

In addition, we also construct an RBF-FD projection matrix such that

 u(X)=~u(X)≈P~u(Z), (8)

simply by replacing with the identity map in the computation. The matrix will appear later as part of our time-stepping 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 RBF-FD weight requires only matrix-vector 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 RBF-FD stencils of a given number of nearest neighbors on scattered nodes require the use of quasi-uniform nodes on the surface Shankar et al. (2015). By using RBF-FD stencils on the grid nodes, we avoid the ill-conditioning of the RBF-FD matrices that arises from the clustering of nodes on the surface. Figure 1 shows an example of an RBF-FD 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

 ut=ΔΓu

intrinsic to the surface . Using the forward Euler scheme with spatial discretization at as above, we have

 u(X,tn+1)=u(X,tn)+ΔtΔΓu(X,tn)+O(Δt2), (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 time-stepping scheme is initialized using the initial condition . It is equivalent to consider the discrete equation of the constant-along-normal extended function, and (9) becomes

 u(X,tn+1) = ~u(X,tn)+ΔtΔ~u(X,tn)+O(Δt2),n∈N. (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 step-size 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

 ˜Un+1Z:=(P+ΔtW)˜UnZ,n∈N. (11)

In (11), the exterior nodes are also used implicitly in the construction of RBF-FD matrices and defined in (7) and (8) as RBF centers. In other words, the proposed RBF-CPM runs solely based on RBF interpolations with the regularly placed nodes. To evaluate the numerical approximation, say at for simplicity, one can evaluate . Other time-stepping schemes are also possible. In our experiments, we use the third-order, three-stage SSP Runge-Kutta scheme Shu and Osher (1988) in advection-dominant 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 RBF-CPM.

Given a collection of Cartesian grid points in a small tubular domain that contains the surface , the algorithm of the RBF-CPM for a time dependent PDE consists of the following steps:

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

2. Compute the RBF-FD matrices, i.e. the matrices and in (7)-(8). For each surface point , :

1. Find the closest grid points to .

2. Compute the RBF-FD weight at the surface node .

3. Solve the surface PDE using an explicit time-stepping scheme, e.g. (11), until the final time.

Our method allows pre-computation 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 RBF-FD 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 ill-conditioning, we use the stable RBF-GA method which provides an accurate and stable algorithm and is a cheaper stabilization method over RBF-QR 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 RBF-CPM. 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 Laplace-Beltrami operator in comparison to the exact. The first test problem (P1) approximates the Laplace-Beltrami 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 Laplace-Beltrami operator is applied to the function on the unit sphere

 Γ={x:x(θ,ϕ)=(cos(θ)cos(ϕ),sin(θ)cos(ϕ),sin(ϕ)),−π≤θ<π,−π2≤ϕ≤π2}.

Here, the relative error can be computed using the known, exact solution .

We apply the RBF-FD 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 RBF-FD stencils used in the RBF-CPM, i.e., for , are formed using the -nearest neighboring regular grid points with a predetermined spacing . For such problems, the matrix-based formulation of the closest point method Macdonald and Ruuth (2009) can be used to obtain the minimal-sized 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

 x2+y2≤q

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

 x2+y2+z2≤q

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 On-Line Encyclopedia of Integer Sequences oei (2011). Table 2 shows some of the integer solutions for the Gauss circle problem in two and three dimensions.

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 RBF-FD stencil; see Table 3.

## 4 Numerical experiments

In this section, we test our method, the RBF-CPM, on a number of examples in two and three dimensions. For problems involving the Laplace-Beltrami 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 step-size 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

 ut=ΔΓu

intrinsic to the unit circle . Following Ruuth and Merriman (2008), for an initial profile , the exact solution is

 u(θ,t)=e−tsinθ.

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 second-order. 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 .

#### 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

 u(θ,t)=e−tsinθ.

Following Macdonald et al. (2011), we introduce a modified closest point mapping which equals for points that map to the interior of the semi-circle. 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. Second-order 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 RBF-CPM uses points.

#### 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

 ut+us=0,

with being the arclength, is imposed on an ellipse with major axis along the y-axis and minor axis along the x-axis. By Ruuth and Merriman (2008), application of the closest point principles to the surface PDE leads to the embedding PDE

 ut+T(x,y)⋅∇u=0

with

 T(x,y)=(−y/b2,x/a2)√y2/b4+x2/a4.

For an initial profile , the exact solution for subsequent times is

 u(s,t)=sin3(2π(s−t)/L),

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 RBF-FD stencils used for approximating the first-order derivatives, the TVD-RK3 scheme Gottlieb and Shu (1998) is chosen for the time discretization with a time step-size . 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. Second-order 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.

#### 4.1.4 Advection-diffusion equation on an ellipse

Our next example considers advection-diffusion on an ellipse. The equation

 ut+us=uss,

with being the arclength, is imposed on an ellipse with major axis (along the y-axis) and minor axis (along the x-axis). Similar to the previous example, application of the closest point principles to the surface PDE gives

 ut+T(x,y)⋅∇u=Δu

with

 T(x,y)=(−y/b2,x/a2)√y2/b4+x2/a4.

For an initial profile , the exact solution has the form

 u(s,t)=e−2πt/Lsin(2π(s−t)/L),

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. Second-order convergence is observed.

### 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

 ut=ΔΓu

on the unit sphere . For an initial profile , the exact solution for all times is

 u(θ,ϕ,t)=e−2tsinϕ.

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. Second-order 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 RBF-CPM uses points. This corresponds to a reduction of .

#### 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

 Γ={x:x(ϕ,θ)=((12cos(ϕ)+1)cos(θ),(12cos(ϕ)+1)sin(θ),12sin(ϕ)),−π≤θ,ϕ≤π},

the advection equation is given by

 ut+uϕ=0.

For an initial profile

 u(ϕ,θ,0)=f(ϕ)=⎧⎨⎩g(ϕ+ππ),−π≤ϕ≤0,g(π−ϕπ),0<ϕ<π,

where

 g(s)=e1/(s−1)−e−1/se−1/s+e1/(s−1),

the exact solution at time is

 u(ϕ,θ,t)=f(ϕ−t).

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 TVD-RK3 scheme is chosen for the discretization in time with step-size . Table 9 shows the error and the convergence rate of the method for various grid sizes . Convergence is at least second-order. 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 Perona-Malik model to denoise surface images. The equation has the form

 ut=∇Γ⋅(g(|∇Γu|)∇Γu)

where is the diffusion coefficient given by

 g(s)=11+(s/λ)2

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 Reaction-diffusion systems

Our final example evolves the Gray-Scott reaction-diffusion model Gray and Scott (1984) on a triangulated surface. The Gray Scott model describes the chemical reaction

 U+2V⟶3V,
 V⟶P,

where , and are chemicals. The corresponding surface model has the form

 ut=F(1−u)−uv2+DuΔSu
 vt=−(F+k)v+uv2+DvΔSv

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 step-size 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 (RBF-FD). In our method, an RBF-FD 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, higher-order 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 quasi-uniform 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 RBF-FD discretization introduced in this paper solves surface PDEs using explicit time stepping methods. Implicit RBF-FD 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 RBF-CPM 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 RBF-FD 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 2014-2017 program and by Spanish Ministry of Economy and Competitiveness MINECO through BCAM Severo Ochoa excellence accreditation SEV-2013-0323 and through project MTM2015-69992-R 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 reaction-diffusion, vol. 25, ACM, 1991.
• Bertalmio et al. (2001) M. Bertalmio, A. L. Bertozzi, G. Sapiro, Navier-stokes, 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 volume-based 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, Springer-Verlag 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, Real-time 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 Mathematics-International 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 high-order kernel method for diffusion and reaction-diffusion 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 (2-3) (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 kernel-based 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 RBF-FD approximations: I. interpolation and accuracy, Journal of Computational Physics .
• Shu and Osher (1988) C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing 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 Gaussian-based RBF-FD 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 on-line encyclopedia of integer sequences, URL http://oeis.org, accessed: 2016-10-25, 2011.
• Gottlieb and Shu (1998) S. Gottlieb, C.-W. Shu, Total variation diminishing Runge-Kutta 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.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters