A localized meshless method for diffusion on folded surfaces

A localized meshless method for diffusion on folded surfaces Ka Chun Cheunga,∗, Leevan Linga , Steve Ruuthb a Department of Mathematics, Hong Kong Baptist University, Kowloon Tong, Hong Kong. of Mathematics, Simon Fraser University, BC, Canada. b Department Abstract Partial differential equations (PDEs) on surfaces arise in a variety of application areas including biological systems, medical imaging, fluid dynamics, mathematical physics, image processing and computer graphics. In this paper, we propose a radial basis function (RBF) discretization of the closest point method. The corresponding localized meshless method may be used to approximate diffusion on smooth or folded surfaces. Our method has the benefit of having an a priori error bound in terms of percentage of the norm of the solution. A stable solver is used to avoid the ill-conditioning that arises when the radial basis functions (RBFs) become flat. Keywords: Power function, radial basis function (RBF), closest point method, diffusion, surface. 1. Introduction A variety of numerical methods have been proposed to approximate the solution of partial differential equations on surfaces. These include discretizations on parameterized surfaces, finite difference and finite element methods on triangulated surfaces, and embedding methods. Each class of methods is characterized by certain advantages and disadvantages. Parametrization methods can be very efficient, but these methods require a parametrization of the surface. Typically, the construction of a parameterization leads to distortions in the surface as well as singularities [1, 2, 3, 4]. We may avoid this complication by discretizing on triangulated surfaces. Unfortunately, finite difference discretizations on triangulated surfaces are more difficult to construct on curved surfaces than on the plane and it can be difficult to compute geometric quantities such as the curvature and the normal to the surface ∗ Corresponding author Email address: kachun.charles@gmail.com (Ka Chun Cheung) 1 This project was supported by a CERG Grant of the Hong Kong Research Grant Council, an FRG Grant of HKBU, and a grant from NSERC Canada. 1 A localized meshless method for diffusion on folded surfaces 2 [5, 6, 7]. On the other hand, finite element methods on triangulated surfaces can be effective when applied to problems of elliptic or parabolic type [8, 9]. In the embedding class of methods, the surface and the corresponding surface PDE are represented in the underlying d-dimensional embedding space (e.g., a problem posed on a two-dimensional curved surface might be represented and solved in a three-dimensional neighbourhood of the surface) thereby enabling the use of standard Cartesian grid methods in Rd . The most popular finite difference [5, 6] and RBF discretizations [10, 11] of this type project the surface Laplacian to the tangent space of the surface. Often, a level set representation of the surface is used, in which case there is no direct generalization to open surfaces or surfaces without orientation. Moreover, embedding methods typically introduce artificial boundary conditions at the boundary of the computational domain, leading to an addition error and first-order accuracy for diffusion problems. These complications can be overcome by using the closest point method [12], which is a high-order accurate method that embeds the surface within the Euclidean embedding space by means of a closest point mapping. The standard explicit closest point method is comprised of an evolution step and an interpolation step: the method must interpolate function values on the surface at each time step in order to replace derivatives intrinsic to the surface with standard Cartesian operators. In [13], Piret proposed to discretize the closest point method using radial basis functions. As a consequence of using radial basis functions, the method is not restricted to Cartesian grids. This allows for two fundamental advantages over the standard closest point method. First, it is straightforward to use irregular or refined grids according to the needs of the problem. Second, the freedom of selecting nodal points naturally allows one to eliminate the interpolation step of the closest point method by making an appropriate node selection. By avoiding the interpolation step, the computational error and cost associated with interpolation is eliminated. The method proposed in this paper also combines the closest point method and the radial basis function projection method to solve PDEs on surfaces. As in Piret’s approach, our approach does not involve an interpolation step, thereby eliminating a source of error and opening the possibility of significant savings in computational cost. Apart from this, we will show that our proposed method is error optimal, i.e., the spatial approximation has a minimal error bound among all meshless methods. Moreover, the error bound is of a-priori type: it depends only on the distribution of nodes and the type of kernel or basis function. On the other hand, the discretization of the spatial variables in our approach leads to a differentiation matrix W such that ut = W u. While the differentiation matrix W is error optimal, it is also very ill-conditioned when the radial basis functions become flat or when the nodal points are very dense. Consequently, as we refine our discretization, the differentiation matrix becomes ill-conditioned [14, 15, 16]. This trade-off is inherent in meshless methods. See [17, 18] for a strategy to handle the ill-conditioning by changing the basis function to an orthogonal basis, thereby stabilizing the computation for flat radial basis functions. Throughout this paper, the differentiation matrix W arising in numerical experiments is solved by this RBF-QR A localized meshless method for diffusion on folded surfaces 3 strategy. In addition to the development of an error-optimal RBF discretization, this paper introduces a generalization of the closest point method appropriate for diffusion on folded surfaces. The standard closest point method is built on the assumption that the closest point function is sufficiently smooth within the computational band. As a consequence, the standard method cannot be used to compute the solutions of PDEs on curves with corners or surfaces with folds. In [19], a transformation was introduced to give a smooth, transformed problem which is amenable to computation via the standard closest point method. Unfortunately, the transformed problem has increased dimensionality, and consequently is computationally expensive. In this paper, we develop an extension of the closest point method for solving diffusion on curves with corners and on folded surfaces. The method computes within the original embedding space and uses node placement along the folds. This makes a radial basis function discretization particularly appropriate. Numerical experiments demonstrate that the proposed method retains high-order accuracy for smooth curves and surfaces, while giving first order convergence for computations on geometries with folds. 2. Surface Embedding The closest point method, originally proposed by Ruuth and Merriman [12], is an embedding method for solving PDEs on surfaces. The method embeds the surface operator to the Euclidean space Rn using a closest point mapping. Such an embedding is called a closest point extension. By alternating the closest point extension with a traditional PDE solver on Rn , we can obtain an explicit method for solving a surface PDE in the embedding space. Let S ⊂ Rn be a smooth surface and ∇ be the standard Euclidean gradient in Rn . Then, the surface gradient, denoted by ∇S , is defined as: ∇S := ∇ − n̂(n̂ · ∇) , (1) where n̂ is the unit normal to the surface S. Consider a C 1 vector field φ : Rn → Rn that is constant along the normal direction to S. The directional derivative satisfies d (u ◦ φ) = n̂ · ∇(u ◦ φ), 0= dn̂ which we combine with (1) to obtain ∇S (u ◦ φ) = ∇(u ◦ φ) − n̂ · ∇(u ◦ φ) n̂ = ∇(u ◦ φ), on S. (2) Thus, for any C 1 function that is constant along the normal direction to S, we observe that, on the surface, standard Euclidean gradients (defined on the embedding space) coincide with the desired surface gradients. A particularly convenient vector field φ : Rn → S is the closest point mapping φS (x) = arg min |x − s|, for all x ∈ S, s∈S (3) 4 A localized meshless method for diffusion on folded surfaces where S ⊆ Rn is a neighborhood of S over which the closest point mapping φS is C 1 . Generally, the maximum possible S is determined by the curvature of S. For any C 1 surface, we can properly embed S ⊂ S. For example, if S ∈ R2 is circle, then the maximum possible S will be the whole plane R2 except the center of S. Similar conditions can be derived for the surface divergence for any differentiable vector field ϕ : Rn → Rn . Similarly, if ϕ is tangent to the surface and all the surfaces displaced by a fixed distance from S, then the surface divergence of ϕ coincides with the usual divergence at the surface ∇S · ϕ = ∇ · ϕ, on S. (4) A combination of these two principles demonstrates that we may evaluate the traditional Laplacian ∆(u ◦ φS ) = ∇ · ∇(u ◦ φS ) in S instead of the surface Laplacian ∆S (u ◦ φS ) = ∇S · ∇S (u ◦ φS ) on S for u ∈ C 2 . In fact, any nonlinear diffusion operator can be embedded to a traditional one in S; see [12]. 3. S–Orthogonal Data Distribution For the sake of clarity, let uS be the restriction to the surface S for any function u defined on Rn . Consider the heat equation intrinsic to the surface S, i.e., ∂t uS = ∆S uS on S. (5) After embedding the surface PDE to S, we get an embedded PDE: ∂t u = ∆(u ◦ φS ) in S. (6) n We can use any solver developed for the embedded PDE posed on R to solve the surface PDE (5). Meshless methods are methods that do not require any mesh structure. When applied to the closest point method with an appropriate mesh, meshless methods give a trivial discretization of the closest point extension step. In this work, we adopt the localized meshless approach in [20, 21, 13, 22, 23, 24]. First, let X = {xi }N i=1 ⊂ S be the set of data points used to discretize S; see Figure 1 for a schematic demonstration. Using the set of normal vectors n(X) = {n̂(xi ) : xi ∈ X}, the set X can be extended to a discretization for S. We define Gh (X) : S → S such that Gh (X) = X + hn(X) = {xi + h n̂(xi ) : xi ∈ X, h ∈ R}. It is easy to see that Gh (X) is a set of points extended from X by a distance h along its normal direction n(X). For simplicity, assume h > 0 is sufficiently small (w.r.t. k) such that, G : S → S and X = G(X) := X ∪ G±h (X) ∪ · · · ∪ G±kh (X) ⊂ S. (7) With N = |X| points on S, we now have a set of M = (2k + 1)N data points in S. By construction, trivially, we also have a discretized closest point mapping (3), that maps u|X → u|X . The tradeoff is that the data points in X are not regular. 5 A localized meshless method for diffusion on folded surfaces S X X\X xi S Xi \xi X i \Xi Figure 1: A schematic demonstration for a data point xi on surfaces S ⊂ R2 along with its symmetric 3-point local neighbourhood Xi and S-orthogonal extension X i . 3.1. Point distribution for corners and folded surfaces S Suppose for now that our surface S is nonsmooth and S = Sj , where each Sj is smooth. Define the singular curve Γ to be the intersection of these subsurfaces, i.e., the set on which S is nonsmooth. It is desirable to have data on Γ to enhance the resolution of the numerical approximation. We now describe mesh selection on singular surfaces. For each data point xi ∈ X ⊂ S, we either have xi inside a smooth subsurface or on the singular curve. For xi ∈ Sj \ ∂Sj (for some j), its local neighbourhood should be restricted to the same subsurface; i.e., Xi ⊂ Sj . This local neighbourhood can then be extended in the S-orthogonal direction to yield X i . Now, forTxi ∈TΓ, xi must belong to two or more smooth subsurfaces. That is, xi ∈ Si1 · · · Sik ⊆ Γ where k ≥ 2. Notice that our discretization requires smoothness of the closest point function in a band around each subsurface, a requirement that can be satisfied by taking h sufficiently small. This contrasts with the original closest point method that places a more strict, global restriction over S. Figure 2 gives a simple case of two smooth curves meeting at a corner point xi . At the corner point xi , asymmetric local neighbourhoods are used to satisfy Xi1 ⊂ S1 and Xi2 ⊂ S2 . In our method, these local neighbourhoods will lead to multiple sets of weighings for the corner xi . An arithmetic average will be taken as the final weighting. Note that the S-orthogonal extensions X i1 can extend past subsurface S2 , and vice versa for an acute angled corner with our method. These overlapping grids will not cause any trouble for our method. Figure 3(a) shows an example of a dented sphere S ⊂ R3 composed of S1 , a part of the unit sphere, and S2 , the dented top. Figure 3(b) is a bird’s-eye 6 A localized meshless method for diffusion on folded surfaces S1 S2 X xi Xi2 \xi X i2 \Xi1 S1 S2 X xi Xi1 \xi X i1 \Xi1 Figure 2: A schematic demonstration of the local neighbourhoods (Xi1 and Xi2 ) and extenT sions (X i1 and X i2 ) for a corner point xi ∈ S1 S2 . S2 S1 (a) 1 Γ X xi 0.8 Xi2 Xi2 ∩ Γ 0.6 1 0.4 0.2 0.2 0.2 0.4 0.6 0.8 (b) Xi1 Xi1 ∩ Γ 0.6 0.4 0 Γ X xi 0.8 0 0.2 0.4 0.6 0.8 (c) Figure 3: A nonsmooth manifold in R3 and a schematic data point distribution. 7 A localized meshless method for diffusion on folded surfaces view of some data points X near a part of Γ. In this projected annular region belongs to S1 and the inner region belongs to meet at Γ (represented by the solid line). For a data point xi and X i2 are obtained via neighborhood searches. See Figures for illustrations of these sets. view, the outer S2 . ST1 and S2 ∈ S1 S2 , X i1 3(b) and 3(c) 3.2. Reduced semi-discrete form For each xi ∈ X ⊂ S and some ε > 0, let the ith local neighbourhood be denoted by Xi = {x ∈ X ⊂ S : |x − xi | < ε}. Then, this neighbourhood on S can be extended to S by X i = G(Xi ) ⊂ X. Assuming we have a set of generalized finite-difference weightings wij ∈ R, such that for each xi ∈ X, X wij u(ξ j ), ξ j ∈ X i , Lu(xi ) ≈ (8) j where L is a linear functional. For data xi ∈ Γ on theSsurface S singularity, the linear combination in (8) involves all points ξ j ∈ X i1 · · · X ik . A detailed discussion on the generalized finite-difference method is presented in the next section. For now, suppose (8) is available for all data points in X. Collecting all these weights in an M × M sparse matrix W , as in the standard finite difference method, we obtain a semi-discrete formulation, ∂t u(X) = W (u ◦ φS )(X). (9) Advancing (9) via any time discretization scheme will yield a set of new nodal values at all points in X. For simplicity, we consider the forward Euler scheme which yields u(X, t + ∆t) = (u ◦ φS )(X, t) + ∆tW (u ◦ φS )(X, t). (10) It remains to discretize the closest point mapping u ◦ φS . This will yield a standard difference equation to update the approximate solution of the embedded PDE (6). Recall that our target is the surface PDE (5) rather than the embedded (6). In other words, we are only interested in the function values of u ◦ φS or, equivalently, u values at points X belonging to S. Those at points X \ X are specified by the closest point mapping, i.e., u ◦ φS (X) = u(G −1 (X)). Notice, however, that our data points give the function values outside S as copies of those on the surface. Thus, in (10), the values of u(X \X) and the corresponding weightings in W are unnecessary. Applying this observation, the closest point mapping can be written in matrix form as (u ◦ φS )(X, t) = Eu(X, t), for E = 1 ⊗ IN ∈ {0, 1}M×N , where ⊗ is the kronecker product and E is an upsampling with the (2k + 1) vector 1 of all ones, assuming points in X are listed first in X as in (7). We also define the corresponding downsampling E † = eT1 ⊗ IN ∈ {0, 1}N ×M that A localized meshless method for diffusion on folded surfaces 8 extracts values at X. Applying operators E and E † , the M × M system (10) can be reduced to an N × N system: u(X, t + ∆t) = = = = E † u(X, t + ∆t),
E † (u ◦ φS )(X, t) + ∆tW (u ◦ φS )(X, t) ,
E † Eu(X, t) + ∆tW Eu(X, t) ,
IN + ∆t E † W E u(X, t), (11) where the N × N matrix E † W E only requires N sets of finite difference weights wij , given by (8), for the points X instead of all M sets for X. Equation (11) gives the forward Euler approximation of the time discretization of the surface PDE ∂t uS = LS uS , (5), on data points X ⊂ S. Other time-stepping methods may also be considered. Similar to the forward Euler case, such methods should be applied to the reduced semi-discrete form for the surface PDE, ∂t uS (X) = (E † W E) uS (X). (12) 4. Norm-Optimal Approximation We now address how the generalized finite difference weightings (8) can be obtained for all xi ∈ X ⊂ S to evaluate the matrix E † W E in the reduced semidiscrete surface PDE (12). Recall that xi ∈ X is any given data point on the surface S. The subset of neighboring points Xi ∈ S of xi is then obtained by some neighborhood search technique over the set X, such as a kd-tree search. Then, the subset X i ⊂ X is obtained by extending points in the normal direction, as in (7), and we denote ni = |X i |. Throughout this section, the index i will be used to denote the current point under consideration and will be omitted when no ambiguity arises. Generally speaking, we need a method to look for the weighting {wij } in (8) such that the linear functional δxi L is well-approximated by a linear combination of the function values of u at X i , X wij u(ξj ), ξ j ∈ X i . δxi L[u] = Lu(xi ) ≈ (13) j Equivalently, in matrix form, it suffices to find the matrix W in the following linear system Lu(X) ≈ E † W u(X), where E † W is of size N ×M with entries [E † W ]ij = wij in (13). For the singular points xi ∈ Γ, the linear functional δxi L has multiple approximations, one for each subsurface to which xi belongs. To select the weighting so that the discretization on a set of given nodes is error optimal, we introduce a reproducing kernel Hilbert space. As will be shown, a priori error estimation is possible with our approach. A localized meshless method for diffusion on folded surfaces 9 Let K : Ω × Ω → R be a symmetric, strictly positive definite kernel with Ω ⊃ X i and let HK (Ω) = span{K(·, y) : y ∈ Ω}. We define an associated bilinear form h ·, ·iK for HK by + * X X XX cj K(·, xj ), dk K(·, yk ) = cj dk K(xj , yk ). j k j K k Thus, HK (Ω) is a pre-Hilbert space with h·, ·iK as an inner product. Completion of HK (Ω) yields a Hilbert space which is called the native space, denoted by H(Ω) for simplicity. The space H(Ω) satisfies hK(x, ·), K(y, ·)iH , hf, K(x, ·)iH , K(x, y) = f (x) = hδxi L, δxj LiH∗ δxi Lx δxj Ly K(x, y) = Lx Ly K(xi , xj ), = for all x, y ∈ Ω, f ∈ H and δx L ∈ H∗ , where H∗ is the dual space of H and the superscript in δxi Lx is used to emphasize its action on the first variable of K(x, y). Suppose we want to find the weights in (8) to minimize the absolute error [25]. The approximation error is given by Lu(xi ) − X wj u(ξ j ) = ξj ∈X i ≤ δx i L − X wj δξj H∗  δx i L − u H X  wj δξj [u] , , (14) where wj := wij is the weight in (8). We define the power function εδx L,X i ,w := i P δxi L − wj δξj for a linear functional δxi L, node point X i , and weighting w. Note that the square of the norm of the power function is a bilinear form such that: X 2 2 P 2 (w) := εδx L,X i ,w H∗ = δxi L − wj δξj H∗ , i ξj ∈X i XX wj wk hδxk , δxj iH∗ , wj hδxi L, δxj iH∗ + X XX = Lx Ly K(xi , xi ) − 2 wj Lx K(xi , xj ) + wj wk K(xk , xj ).(15) = hδxi L, δxi LiH∗ − 2 X Finding an error optimal formula is equivalent to minimizing P 2 . It is easily shown that the solution w∗ of the following system minimizes P 2 : X wj∗ K(xk , xj ) = Lx K(xi , xk ), xk ∈ X i . (16) xj ∈X i Equivalently, in matrix form a minimizer is obtained if w∗ satisfies K(X i , X i )w∗ = Lx K(xi , X i ). (17) A localized meshless method for diffusion on folded surfaces 10 Applying this result to (15), we observe that the minimum of P 2 (i.e., the square of absolute error of (8) at the point xi ) is given by X wj∗ Lx K(xi , xj ). min P 2 (w) = P 2 (w∗ ) = Lx Ly K(xi , xi ) − w xj ∈X i Consequently, the solution of (16) is the weighting such that the approximation (13) has the smallest error in the sense of the native space norm k · kH . Note that the native space H is a reproducing kernel Hilbert space. For each type of kernel K, there is a corresponding native space NK . In our prob2 1 lem, we take K(x, y) = e−( c kx−yk2 ) . For the Matérn kernel K(x, y) :=
m−d/2 kx−yk2 Km−d/2 1c kx − yk2 , x, y ∈ Rd , where Kv (r) is the modified Bessel function of the second kind, the corresponding native space is the Sobolev space W2m (Rd ) for m > d/2. For the native space of other kernels, see [26, 27]. In both cases, the shape parameter c controls the shape of the kernel. The tradeoff property for the radial basis function method [28, 29] is that the accuracy improves as the shape parameter c grows large or as the nodal point distribution becomes dense. Computationally, either case comes at the cost of more ill-conditioned matrices. Once the kernel K is fixed, one can directly compute the power function and error bound without any knowledge of the exact solution. Thus, we can obtain an error bound for a given node before computing any numerical solution. Note that the error bound (14) can also be interpreted as the relative error of the discretization in terms of the H-norm of the exact solution u. 4.1. The meshless finite difference method The meshless finite difference method, a.k.a. the localized meshless method [20, 22, 23, 24], has been developed for the same purpose to find the weighting w in (13). Here, we shall show that the weighting found by the meshless finite difference method coincides in exact arithmetic with ours in (16). Suppose we are given a set of data sites X i , but with unknown function values u(X i ). Meshless methods, with a suitable kernel function K, give an interpolant of ũ in the form X ũ(x) = wj K(x, xj ). xj ∈X i Imposing interpolation conditions ũ = u at all x ∈ X i , we express the weighting as w = K(X i , X i )−1 u(X i ), where u(X i ) is a vector of (unknown) function values and the matrix entries are defined by [K(X i , X i )]jk = K(xj , xk ) for xj , xk ∈ X i . Note that the matrix K(X i , X i ) is nonsingular provided the kernel K is positive definite. 11 A localized meshless method for diffusion on folded surfaces The meshless finite difference method uses information from the interpolant to approximate the unknown function. For our discussion we have Lu(xi ) ≈ = Lũ(xi ), X wj Lx K(xi , xj ), xj ∈X i = = = (Lx K(xi , X i ))T · w, (Lx , K(xi , X i ))T K(X i , X i )−1 u(X i ), X w · u(X i ) = wj u(xj ). xj ∈X i Therefore, the meshless finite difference method returns the solution of K(X i , X i )w = Lx K(xi , X i ), as the weighting w, which is identical to our previous results (16) and (17). Thus, the meshless finite difference method is error optimal. Note, however, that the reproducing kernel Hilbert space approach [26, 25, 30] has the advantage of allowing us to measure error by using the corresponding power function. 4.2. Finding optimal-weights numerically To numerically obtain the generalized finite difference weighting for point xi ∈ X ⊂ S, we solve the ni × ni matrix system (17). This gives an ni = |X i | point interpolant over a neighbourhood of xi . For any commonly used global kernel, such as the Gaussian and Matérn kernels, the kernel matrix K := K(X i , X i ) is full and, generally speaking, ill-conditioned. To solve (17) directly, it is necessary to invert the kernel matrix, which gives rise to numerical instability. In our 2 1 numerical computation, we adopt the Gaussian kernel K(x, y) = e−( c kx−yk2 ) and apply the stable RBF-QR method by Fornberg et al. [17]. Note that the method can also be applied to the Matérn kernel [31, 17] in which case one will be working with error bounds in some standard Sobolev spaces, which are the native spaces of the Matérn kernel [26]. It is not difficult to see why the matrix K(X i , X i ) becomes ill-conditioned for large c-values. In the limit c → ∞, the Gaussian kernels become flat and form a numerically unstable basis. To avoid this problem, Gaussian radial basis functions may be changed to another set of stable basis functions sharing the same 2 1 span. It is shown in [17] that the Gaussian kernel K(x, xk ) = e−( c kx−xk k2 ) has an expansion in 2D, j−p K(x, xk ) = ∞ X 2 X j−p c dj,m cj,m (xk )Tj,m (x) j=0 m=0 for 0 ≤ m ≤ ⌊ 2j ⌋ = c Tj,m (x) = e−r s dj,m sj,m (xk )Tj,m (x), j=0 m=1−p j−p 2 2 + ∞ 2 X X (18) with p = mod2 j. In this expression, the functions /c2 2m r
Tj−2m (r) cos (2m + p)θ , for 2m + p = 0 (19) 12 A localized meshless method for diffusion on folded surfaces and s Tj,m (x) = e−r 2 /c2 2m r
Tj−2m (r) sin (2m + p)θ , for 2m + p 6= 0, (20) form a linearly independent set and are defined using the Chebyshev polynomials Tm (r) where x = (r, θ) is the polar coordinate of x. The coefficients in (18) depend on the point distribution of X i only and are given by  j + 2m + p
j − 2m − p
−1 dj,m = c−2j 2j−2m−1 ! ! , 2 2
2 2 cj,m (xk ) = b2m+p tj−2m e−rk /c rkj cos (2m + p)θk 1 F2 (αj,m , βj,m , c−4 rk2 ),
2 2 sj,m (xk ) = b2m+p tj−2m e−rk /c rkj sin (2m + p)θk 1 F2 (αj,m , βj,m , c−4 rk2 ), where b0 = 1 and bm = 2, m > 0, t0 = 12 and tj = 1, j > 0. The funcand βj,m = tion 1 F2 is the hypergeometric function with αj,m = j−2m+p+1 2   . j − 2m + 1 j+2m+p+2 2 Despite the complexity, if we truncate the infinite sum in (18) after J > ni terms, the kernel matrix can be analytically decomposed as K(X i , X i )T = CDV T = (QR)DV T , (21) where C = QR is an ni × J matrix consisting of coefficients cj,m (X i ) and sj,m (X i ) and its full QR-factorization, D, is a J × J diagonal matrix consisting c of entries dj,m . The n × J matrix V = V (x) consists of the functions Tj,m (x) s and Tj,m (x). The QR-factorization of the matrix C gives C = QR = Q [R1 R2 ] where R1 contains the left ni × ni block of R and R2 contains the remaining columns. Let D1 be the ni × ni principal submatrix of the diagonal matrix D and denote the remaining nonzero (J − ni ) × (J − ni ) submatrix by D2 . We observe that T QT K(X i , X i ) = RDV T , = [R1 R2 ] = [R1 D1 T D1−1 R1−1 QT K(X i , X i ) = [Ini  D1 D2 R2 D2 ]V T ,  V T, D1−1 R1−1 R2 D2 ]V T . Once the above expansion is known, we define the new basis functions as H(X i , X i ) = K (QR1−T D1−T ) = V [Ini {z } | | =:G D1−1 R1−1 R2 D2 ]T . {z } (22) Note that the diagonal elements of D are computed analytically, whereas R1−1 R2 is computed by backward substitution. The relation between the new kernel matrix and the old kernel matrix is given by H(X i , X i ) = K(X i , X i ) G. A localized meshless method for diffusion on folded surfaces 13 Instead of solving the original system, K(X i , X i )w = Lx K(xi , X i ), we find w in the new (discrete) basis via H(X i , X i )w = Lx H(xi , X i ). (23) Note that the right-hand vector can be obtained analytically from (19) and (20), and that the interpolation matrix with respect to the new kernel, H(X i , X i ), is no longer symmetric. In practice, some polynomial terms should be added to the expansion (13) in order to guarantee polynomial accuracy, or when the kernel K is conditionally positive definite. See [26, 17] for further details. 5. Numerical Demonstrations To begin, we consider diffusion on various surfaces S ⊂ R2 . Three point local neighbourhoods Xi will be used to extend one layer outwards and one inwards, as shown in Figures 1 and 2. Throughout this section, we will use the (scaled) 2 1 Gaussian kernel K(x, y) = e−( c kx−yk2 ) , with various shape parameters c, to evaluate the generalized finite difference weights in (13) and the power function P in (15). The reduced semi-discrete system (12) is evolved in time by the forward Euler scheme. As our interest is the spatial discretization of the surface embedding, small time steps are used: ∆t = 0.001~2 where ~ is the fill distance [26]. Our first example evolves diffusion on the unit circle. The convergence of the relative errors ε and the power function P for three different shape parameters, c = 0.5, 1, and 2 are given in Figure 4. Convergence rates are approximately second order for all three tested shape parameters (rates are estimated by the last 5 data points). Notably, the proposed scheme is convergent for all tested shape parameters. Moreover, the power function P also exhibits second order convergence. To test the performance on nonsmooth curves we apply the same method to a diffusion problem on the square [−1, 1]2 . As shown in Figure 5, the proposed method is still convergent but the rate drops from second order to first. As in the case of the unit circle, we observe a close agreement between the ε– and P–convergence rates. Provided the temporal discretization error is negligible, the numerically accessible convergence rate of P provides a good estimate for the convergence rate of the error. The power function can be viewed as the local truncation error in the finite difference method. We emphasize that the evaluation of P is pointwise; that is, at the generalized finite difference stencils (8) and (13), the power function can be evaluated by the stencils and assigned weights for any SPD kernel, independent of the solution. The value of P yields an upper bound for the spatial discretization error: when the solution is unknown, which is usually the case, the discretization error is bounded by P times an unknown factor kukH. In our convergence plots, we observe that the ratio ε/P changes with the shape parameter c. Since the native 14 A localized meshless method for diffusion on folded surfaces −1 −2 10 ε P 10 (slope = 2.002) (slope = 1.983) 0 ε P 10 (slope = 2.002) (slope = 2.000) ε P (slope = 2.005) (slope = 2.000) −2 10 −3 −1 10 10 −3 10 −4 −2 10 10 −4 10 −5 −5 10 −2 10 −1 10 0 10 −3 10 −2 10 c = 2.0 −1 10 0 10 10 −2 10 c = 1.0 −1 10 0 10 c = 0.5 Figure 4: Convergence profiles for diffusion on the unit circle S. The relative error ε and power function P are plotted for various grid spacings and shape parameters c. 1 0 10 ε P 10 (slope = 1.001) (slope = 1.000) P −1 10 (slope = 0.971) (slope = 1.000) −2 10 −1 0 10 −3 10 −2 10 −1 10 −4 10 −3 −1 10 c = 2.0 0 10 10 −2 10 (slope = 0.861) (slope = 0.999) 1 10 10 ε P 0 10 10 −2 10 2 ε −2 −1 10 c = 1.0 0 10 10 −2 10 −1 10 0 10 c = 0.5 Figure 5: Convergence profiles for diffusion on the square S = [−1, 1]2 . The relative error ε and power function P are plotted for various grid spacings and shape parameters c. 15 A localized meshless method for diffusion on folded surfaces 4 4 3 3 2 2 1 1 0 −2 0 1 0 0 −5 5 0 ε P 10 (slope = 0.995) (slope = 1.000) −1 10 (slope = 0.994) (slope = 1.000) −1 −2 −2 10 −3 −1 0 10 10 −2 10 (slope = 0.995) (slope = 1.000) −1 −2 10 5 10 10 −3 ε P 10 10 0 0 ε P 10 10 −2 10 2 0 −5 2 0 10 3 −3 −1 10 0 10 10 −2 10 −1 10 0 10 Figure 6: Convergence profiles for diffusion on various triangles. The relative error ε and power function P are plotted for various grid spacings with the shape parameter c = 2.0. space H is generated by the c-dependent Gaussian kernel in our demonstrations, different shape parameters yield different spaces H and hence different kukH . Determining the optimal shape parameter is challenging and is not the purpose of this work. Nonetheless, the introduction of the power function allows an easy way to search for the optimal c by brute force. The values of P shown in Figures 4 and 5 are the maximum over all stencils used in a particular run. Intuitively, the first order spatial convergence observed in Figure 5 is due to the error at the corners (this conjecture can be confirmed numerically). We may determine the optimal shape parameters for a problem with corners by minimizing the power function at the corners. Experimentally, we find c = 2 gives good results, so we take it for all computations going forward. Figure 6 gives results for diffusion over 3 different folded, triangulated curves. All three cases give similar errors and first order convergence. The observed quadratic convergence on smooth surfaces and linear convergence on folded surfaces is also found in 3D. In Figure 7, the convergence rates (as measured by the max-norm error) for ε and P are shown for the unit sphere and a punch-in sphere, which is constructed by reflecting all points above z = 0.5 in the place z = 0.5 (thereby allowing an analytic solution to be derived). For the unit sphere, data points X are quasi-uniformly placed on the manifold [32]. For the punch-in sphere we place points along the fold. The local neighborhood 16 A localized meshless method for diffusion on folded surfaces 0 10 ε P ⇐ (slope = 1.951) (slope = 2.069) 0 10 ε P (slope = 0.971) (slope = 0.936) −1 10 −2 10 −1 10 ⇒ −3 10 −2 10 −1 10 0 10 −2 10 −1 10 0 10 Figure 7: Convergence profiles for diffusion on the unit sphere S and a punch-in sphere. The relative error ε and power function P are plotted for various grid spacings and the shape parameter c = 2.0. Middle: snapshot at t = 1 for the same initial condition. of xi on each smooth (sub-) surface is chosen to contain 13 points: 5 points are chosen on the fold and 8 on each smooth subsurface. Similar to the treatment in 2D, we average the finite difference weighings over the two subsurfaces to yield the final spatial discretization for points on the fold. The unit sphere and the punch-in sphere share identical solutions for any given initial condition despite the existence of a fold (we can map between the unit sphere and the punch-in sphere via a reflection in z = 0.5). The numerical schemes used for these two manifolds are clearly different and we observe quadratic spatial convergence for the unit sphere and linear convergence for the punch-in sphere. Figure 8 shows how the difference between the two numerical schemes evolves over time for one fixed spatial discretization. The top two subfigures are the respective initial conditions: they are identical if we unfold the punch-in sphere. From the snapshots in Figure 8, we can see that a disagreement between the two numerical solutions emerges and diffuses away from the fold, which is not surprising since the two numerical recipes are different at the fold. We can see that the relative difference grows gradually with time due to the decay of the solution, whereas the absolute difference is rather steady in time. Clearly, the difference is solely due to the treatment of the fold. Note that our initial conditions vary along the meridian lines. This choice was made to avoid over-simplifying the example via symmetric data at the fold. Not all manifolds are easily unfolded. As an example, consider a flattop sphere constructed by replacing the z-coordinates of the unit sphere by min(z, 0.5). Figures 9 and 10 present two simulated results with different initial conditions. The first simulation begins with two hot spots on the large smooth surface symmetric about the x-axis. As time evolves, we observe the 17 A localized meshless method for diffusion on folded surfaces Initial condition for both unit sphere and punch-in sphere Max 1 1 0 0 Mean −1 −1 0.5 0 −0.5 0.5 1 0 0 −0.5 −1 1 0 Min −1 Relative difference of the approximation on unit sphere and punch-in sphere 0 −0.5 −1 0.5 t = 0.400 t = 0.200 t = 0.000 Max 0 −0.5 −1 0 −0.5 0.5 1 −1 0 Abs. diff. ≤ 00 Rel. diff. ≤ 0.000% 0 −0.5 −1 0 −0.5 0.5 0 0 −0.5 Abs. diff. ≤ 2.78e−03 Rel. diff. ≤ 0.415% 1 −1 −1 1 0 Abs. diff. ≤ 3.50e−03 Rel. diff. ≤ 0.783% 0 −0.5 −1 0.5 t = 1.000 t = 0.800 t = 0.600 Mean 0 −0.5 −1 0 −0.5 1 −1 −1 0.5 0 0 −0.5 Abs. diff. ≤ 3.66e−03 Rel. diff. ≤ 1.222% 0 −0.5 1 −1 0.5 0 0 −0.5 Abs. diff. ≤ 3.63e−03 Rel. diff. ≤ 1.808% −1 1 0 Abs. diff. ≤ 3.53e−03 Rel. diff. ≤ 2.615% Min Figure 8: Top: initial conditions. Middle, Bottom: relative difference of the numerical approximation on the unit sphere and the punch-in sphere at various times t with “identical” initial conditions. A localized meshless method for diffusion on folded surfaces 0.5 0 −0.5 −1 0.5 0.5 t = 0.400 t = 0.200 t = 0.000 0.5 0 −0.5 1 −1 0.5 0 max = +1.0000 min = −1.0000 Max 0 −0.5 −1 −1 0 −0.5 18 0 −0.5 0.5 1 −1 0 max = −0.2096 min = −0.7323 0 −0.5 −1 1 0 max = −0.4107 min = −0.5751 Mean 0 −0.5 −1 0.5 0.5 t = 1.000 0.5 t = 0.800 t = 0.600 0.5 0 −0.5 −1 0 −0.5 1 −1 0 max = −0.4669 min = −0.5248 0.5 0 −0.5 −1 0 −0.5 1 −1 0 max = −0.4836 min = −0.5087 0.5 0 −0.5 −1 1 0 max = −0.4888 min = −0.5026 Min Figure 9: Simulation of diffusion on a flat-top sphere with two hot spots in the x-symmetric initial profile. diffusion progress across the fold to the flat-top and across the large spherical surface. For each snapshot, the maximum and minimum values of the numerical solution are also presented to confirm that the maximum principle is satisfied at consecutive snapshots. The second simulation is initialized with a hot spot located near the edge of the flat-top in order to highlight the method’s treatment of diffusion across a fold. Snapshots are presented at nonuniform times to better visualize the rapid initial flow. We observe that heat flows across the fold at a rate similar to the heat flow within the flat-top. At the final time, the location of hottest temperature coincides with the hot spot for the initial condition. Finally, we note that it will take longer for this simulation to reach steady state due to its more isolated initial profile. 6. Conclusions We propose an embedding method to solve diffusion on manifolds with folds. Our discretization is based on meshless methods. This approach has the advantage of eliminating the interpolation step in the underlying embedding method, thereby eliminating a source of error and computational cost. Moreover, the flexibility of data point placement is central to our approach for representing folds and the associated differential operators. In future work we will develop extensions of the method to more general operators and more general surfaces. Higher order accurate methods will also be developed. A localized meshless method for diffusion on folded surfaces 0.5 0 −0.5 −1 0.5 0.5 t = 0.200 t = 0.100 t = 0.000 0.5 0 −0.5 1 −1 0.5 0 max = +1.0000 min = −1.0000 Max 0 −0.5 −1 −1 0 −0.5 19 0 −0.5 0.5 1 −1 0 max = −0.1344 min = −1.0000 0 −0.5 −1 1 0 max = −0.4461 min = −0.9985 Mean 0 −0.5 −1 0.5 0.5 t = 1.000 0.5 t = 0.500 t = 0.300 0.5 0 −0.5 −1 0 −0.5 1 −1 0 max = −0.5842 min = −0.9906 0.5 0 −0.5 −1 0 −0.5 1 −1 0 max = −0.7096 min = −0.9596 0.5 0 −0.5 −1 1 0 max = −0.8095 min = −0.8936 Min Figure 10: Simulation of diffusion on a flat-top sphere with an off-centered initial hot spot at the top. References [1] M. S. Floater, K. Hormann, Surface parameterization: a tutorial and survey, in: In Advances in Multiresolution for Geometric Modelling, Springer, 2005, pp. 157–186. [2] L. Lui, Y. Wang, T. Chan, Solving PDEs on manifolds with global conformal parameterization, in: N. Paragios, O. Faugeras, T. Chan, C. Schnrr (Eds.), Variational, Geometric, and Level Set Methods in Computer Vision, Vol. 3752 of Lecture Notes in Computer Science, Springer Berlin Heidelberg, 2005, pp. 307–319. [3] J. Stam, Flows on surfaces of arbitrary topology, ACM Trans. Graph. 22 (3) (2003) 724–731. [4] A. Witkin, M. Kass, Reaction-diffusion textures, SIGGRAPH Comput. Graph. 25 (4) (1991) 299–308. [5] M. Bertalmo, 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. [6] 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. A localized meshless method for diffusion on folded surfaces 20 [7] G. Turk, Generating textures on arbitrary surfaces using reaction-diffusion, SIGGRAPH Comput. Graph. 25 (4) (1991) 289–298. [8] G. Dziuk, Finite elements for the Beltrami operator on arbitrary surfaces, in: S. Hildebrandt, R. Leis (Eds.), Partial Differential Equations and Calculus of Variations, Vol. 1357 of Lecture Notes in Mathematics, Springer Berlin Heidelberg, 1988, pp. 142–155. [9] G. Dziuk, C. Elliott, Surface finite elements for parabolic equations, Journal of Computational Mathematics 25 (4) (2007) 385 – 407. [10] N. Flyer, G. B. Wright, A radial basis function method for the shallow water equations on a sphere, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 465 (2106) (2009) 1949–1976. [11] E. Fuselier, G. Wright, Scattered data interpolation on embedded submanifolds with restricted positive definite kernels: Sobolev error estimates, SIAM Journal on Numerical Analysis 50 (3) (2012) 1753–1776. [12] S. J. Ruuth, B. Merriman, A simple embedding method for solving partial differential equations on surfaces, J. Comput. Phys. 227 (3) (2008) 1943– 1961. [13] 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. [14] E. Larsson, B. Fornberg, Theoretical and computational aspects of multivariate interpolation with increasingly flat radial basis functions (2003). [15] B. Fornberg, G. Wright, Stable computation of multiquadric interpolants for all values of the shape parameter, Computers and Mathematics with Applications 48 (56) (2004) 853 – 867. [16] E. Larsson, B. Fornberg, A numerical study of some radial basis function based solution methods for elliptic PDEs, Comput. Math. Appl 46 (2003) 891–902. [17] B. Fornberg, E. Larsson, N. Flyer, Stable computations with Gaussian radial basis functions, SIAM J. Sci. Comput. 33 (2) (2011) 869–892. [18] E. Larsson, E. Lehto, A. Heryudono, B. Fornberg, Stable computation of differentiation matrices and scattered node stencils based on Gaussian radial basis functions, SIAM Journal on Scientific Computing 35 (4) (2013) A2096–A2119. [19] P. Rockstroh, T. Märx, S. Ruuth, An embedding technique for the solution of reaction-diffusion equations on algebraic surfaces with isolated singularities, ArXiv preprint arXiv:1211.5298. A localized meshless method for diffusion on folded surfaces 21 [20] T. Cecil, J. Qian, S. Osher, Numerical methods for high dimensional Hamilton-Jacobi equations using radial basis functions, Journal of Computational Physics 196 (2004) 327–347. [21] M. Li, W. Chen, C. S. Chen, The localized RBFs collocation methods for solving high dimensional PDEs, Eng. Anal. Bound. Elem. 37 (10) (2013) 1300–1304. [22] C. Shu, H. Ding, K. Yeo, Local radial basis function-based differential quadrature method and its application to solve two-dimensional incompressible Navier-Stokes equations, Computer Methods in Applied Mechanics and Engineering 192 (7-8) (2003) 941 – 954. [23] A. Tolstykh, On using RBF-based differencing formulas for unstructured and mixed structured-unstructured grid calculations, Proceedings of the 16 IMACS World Congress, Lausanne, 2000, pp. 1–6. [24] G. Wright, Radial basis function interpolation: Numerical and analytical developments, Ph.D. thesis, University of Colorado (2003). [25] R. Schaback, Direct discretizations with applications to meshless methods for PDEs, Dolomites Research Notes on Approximation 6 (2013) 27–50, proceedings of DWCAA12. [26] G. F. Fasshauer, Meshfree Approximation Methods with MATLAB, World Scientific Publishing Co., Inc., River Edge, NJ, USA, 2007. [27] H. Wendland, Scattered Data Approximation, Cambridge University Press, 2004. [28] R. Schaback, Error estimates and condition numbers for radial basis function interpolation, Advances in Computational Mathematics 3 (3) (1995) 251–264. [29] R. Schaback, Multivariate interpolation and approximation by translates of a basis function, World Scientific Publishing Co., Inc, 1995, pp. 491–514. [30] R. Schaback, A computational tool for comparing all linear PDE solvers, Advances in Computational Mathematics (2014) 1–23. [31] B. Fornberg, E. Larsson, N. Flyer, Stable computations with Gaussian radial basis functions in 2-D, Tech. Rep. 2009-020, Department of Information Technology, Uppsala University (Aug. 2009). [32] L. Ling, Pointonsphere, MATLAB Central File Exchange (2005). URL http://www.mathworks.com/matlabcentral/fileexchange/6977