A KERNEL-BASED LEAST-SQUARES COLLOCATION METHOD FOR SURFACE DIFFUSION arXiv:2109.03409v1 [math.NA] 8 Sep 2021 MENG CHEN∗, KA CHUN CHEUNG†, AND LEEVAN LING‡ Abstract. There are plenty of applications and analysis for time-independent elliptic partial differential equations in the literature hinting at the benefits of overtesting by using more collocation conditions than the number of basis functions. Overtesting not only reduces the problem size, but is also known to be necessary for stability and convergence of widely used unsymmetric Kansa-type strong-form collocation methods. We consider kernelbased meshfree methods, which is a method of lines with collocation and overtesting spatially, for solving parabolic partial differential equations on surfaces without parametrization. In this paper, we extend the time-independent convergence theories for overtesting techniques to the parabolic equations on smooth and closed surfaces. Key words. meshfree method, Kansa method, radial basis function, method of lines, parabolic PDEs, convergence analysis. AMS subject classifications. 65D15, 65N35, 65N40, 41A63. 1. Introduction. For closed manifold S ⊂ Rd and some given square integrable functions f : S×[0, t] → R and g : S → R, we consider parabolic PDEs for some time-dependent surface scalar function u : S × [0, T ] → R in the form of (1.1a) (1.1b) u̇(y, t) + LS u(y, t) = f (y, t) for (y, t) ∈ S × [0, T ] u(y, 0) = g(y) for y ∈ S, with a second-order uniformly elliptic operator in divergence form (1.2)
LS u(y, t) := −∇S · A(y, t)∇S u(y, t) for (y, t) ∈ S × [0, T ]. Assumptions and definitions required to make sense of (1.1) will be provide in Sect. 2. We focus on strong-form collocation kernel-based spatial discretization. In [1, 2], PDE (1.1) were solved by the method of Rothe, in which the PDE is discretized first in time, and then by kernel-based collocation method in space. In this paper, we focus on the method of lines that is theoretically more completed for parabolic PDEs on spheres [3] using Galerkin formulations and also in bulk domains, see [4, 5]. An advantage of using collocation method is that we can completely remove surface integrations from our algorithms (but not theories). Readers can find our theoretical study on the convergence of semi-discretized solution ∗ Department of Mathematics, Nanchang University, Nanchang, China AI Technology Center (NVAITC), NVIDIA, USA ‡ Department of Mathematics, Hong Kong Baptist University, Kowloon Tong, Hong Kong † NVIDIA 1 in Sect. 3 (see Thm. 3.6). Algorithms for computing the fully discretized solution and the corresponding error analysis (see Thm 4.1) were provided in Sect. 4. Numerical examples in Sect. 5 verify the convergence behaviour of the proposed methods and robustness in simulating solutions to Allen-Cahn equations. 2. Notations and preliminaries. Let S ⊂ Rd be a closed, connected, orientable, and complete Riemannian manifold of dimension dim S := dS = d − 1. We further suppose that S is of class C m+1 for some integer m with bounded geometry and boundary regularity as required in the theories in [6, 7]. Under these assumptions, there exists δ > 0 such that Euclidean closest point retraction map (2.1) cp(x) := arg inf ky − xkL2 (Rd ) : Ω → S y∈S is well-defined and C m -smooth in the narrow band domain (2.2) Ω := {x ∈ Rd : kx − ykL2 (Rd ) < δ, y ∈ S}. 2.1. Differential operators without parametrization. We adopt a set of parametrizationfree definitions in [8] for differential operators on S. Let n : S → Rd be a smooth vector field that spans the normal space. We define the orthogonal projection matrix P (y) := Id − n(y)n(y)T , (2.3) that projects onto the tangent space Ty S of S at y ∈ S. For any continuously differentiable scalar surface function v : S → R, the surface gradient operator ∇S is defined by (2.4)

∇S v(y) := P (y)∇ v ◦ cp (y) = ∇ v ◦ cp (y), y ∈ S. For any continuously differentiable surface vector field g : S → Rd , the surface divergence operator ∇S · is defined by (2.5)


∇S · g(y) := P (y)∇ · g ◦ cp (y) = ∇ · g ◦ cp (y), y ∈ S. Note that these definitions are equivalent [6] to their intrinsic counterparts defined by local parametrization and Riemannian metric tensor. 2.2. Diffusion tensor. The diffusion tensor A(y, t) in (1.2) needs to satisfy some symmetric positive definiteness assumption when restricted to the tangent space Ty S of S. A SSUMPTION 1. The symmetric surface intrinsic diffusion tensor A : S ×[0, t] → Rd×d 2 in (1.2) satisfies, for all y ∈ S and t ∈ [0, T ], that (1) Ty S is an invariant subspace of A, i.e., ξ ∈ Ty S ⇒ A(y, t)ξ ∈ Ty S, (2.6a) ∀ ξ ∈ Ty S, (2) A(y, t) is uniformly positive definite on Ty S, i.e., there exists ν > 0 such that νkξk22 ≤ ξ T A(y, t)ξ ≤ ν −1 kξk22 (2.6b) ∀ ξ ∈ Ty S, and, (3) all components aij (y, t) := [A(y, t)]ij , 1 ≤ i, j ≤ d, of A are sufficient smooth with bounded derivatives. Equivalently, (2.6a)–(2.6b) ensure the existence of a set of orthonormal eigenvectors S {τi (y, t)}di=1 ⊂ Rd of the diffusion tensor A(y, t) that spans Ty S and all the associated eigenvalues are strictly positive. Moreover, the surface vector field A(y, t)∇S u(y, t) is tangent to the surface and lies in Ty S. Also, by (2.6b) and [6, Thm. 4.2–4.3], we can recast the diffusion operator in (1.2) by the Cartesian gradient and divergence operators as (2.7)
LS u(y, t) = −∇ · A(cp(y), t)∇(u(cp(y), t)) for (y, t) ∈ S × [0, T ], in terms of the cp-mapping (2.1). 2.3. Sobolev spaces and norms. Hilbert spaces on complete manifolds with bounded geometry S is defined to by H k (S) := (I − ∆S )−k/2 L2 (S), see [9], which is norm equivalent [10] to the Sobolev spaces characterized by localization by an atlas and a subordinate partition of unity [11, Def. 4.4]. Using some suitable partition of unity {Ui }i , we decompose P v = i vi with supp(vi ) ⊂ Ui and express S ∩ Ui by, say, xd = ai (x1 , . . . , xdS ) without loss of generality. Let the multi-indexed derivative with respect to the Cartesian coordinate x1 , . . . , xd be (2.8) Ddα := 1 ∂xα 1 ∂ |α| , d · · · ∂xα d |α| = α1 + · · · + αd and, we also define DdαS as in (2.8), but, with respect to x1 , . . . , xd−1 = xdS . Let the surface element with respect to variables x1 , . . . , xdS be dσi =  1 1+  dS  X ∂ai 2 2 j=1 ∂xj dx1 · · · dxdS . k We define the non-standard Hcp (S)-norm via cp-operator by (2.9) kvk2Hcp k (S) = X Z |α|≤k 2 Dα v dσ := S X i 3 kvi k2Hcp k (U ) i with (2.10) kvi k2Hcp k (U ) := i X Z |α|≤k Ui   Ddα (vi ◦ cp) 2 xd =ai (x1 ,...,xdS ) dσi , k ∈ N. We explicitly include the cp-operator in definition to go along with (2.7). For functions k (S)), we adopt the following space-time norms notation f ∈ L2 (0, T ; Hcp kf kL2 (0,T ;Hcp kf kH k (S) k (S)) = L2 (0,T ) = T Z 0 kf k2Hcp k (S) dt 1/2 . L EMMA 2.1. For any C m+1 surface S satisfying the assumptions in Sect. 2 and all k k ≤ m, the Hcp (S)-norm in (2.9) is equivalent to the W2k (S)-norm defined via atlas and subordinate partition of unity. Proof. For v ∈ W2k (S), [11, Sect. 4.2] shows that the W2k (S)-norm (i.e., defined by an atlas and a subordinate partition of unity) is equivalent to a surface norm defined by means of coordinate invariant surface integrals as in (2.9), but with localized counterparts (2.11) kvi k2W k (Ui ) = 2 X Z |α|≤k  DdαS vi Ui xd =ai (x1 ,...,xdS )  2 dσi . Note that this W2k (Ui )-norm in (2.11) is exactly the one that can be deduced by using metric tensors of Riemannian manifold with the parameterized equation
T r(x1 , x2 , . . . , xdS ) = x1 , x2 , . . . , xdS , ai (x1 , . . . , xdS ) : RdS → S ⊂ Rd , which defines a set of basis of the tangent space {τj }1≤j≤dS :=  ∂r ∂xj  = 1≤j≤dS (" ej ∂ai ∂xj #) , 1≤j≤dS where ej ∈ RdS is the j-th standard unit vector. Without loss of generality, we assume ∂ai /∂xj 6= 0 for some j, or else xd = const. In this case, the norms in (2.10) and (2.11) are trivially equivalent. Since vi = vi ◦ cp for all y = (x1 , . . . , xd ) ∈ S, we now work on the function in the
integrand of (2.11) and define V (x1 , . . . , xdS ) := vi ◦ cp : RdS → R. By xd =ai (x1 ,...,xdS ) implicit differentiation, we know that the first derivatives were connected by ∂ ∂ ∂ai ∂ V = (vi ◦ cp) + (vi ◦ cp), ∂xj ∂xj ∂xj ∂xd 4 j = 1, . . . , dS . Rewriting in matrix form yields  T ∇dS V (x1 , . . . , xdS ) = τ1 , . . . , τdS ∇d (vi ◦ cp)(y). (2.12) From the fact that ∇d (vi ◦ cp)(y) ∈ Ty S, we can uniquely express dS X ∇d (vi ◦ cp)(y) = (2.13) γj τj j=1 with some coefficients γ = [γ1 , . . . , γdS ]T ∈ RdS . Putting (2.13) into (2.12) yields ∇dS V (x1 , . . . , xdS ) = dS X γj τi · τj =: GS (y)γ. i,j=1 By assumption, the eigenvalues of the Riemannian metric tensor GS (y) of S at y ∈ S are positive and bounded for all y ∈ S. For some constants 0 ≤ c1 ≤ c2 , we have c1 kγk22 ≤ k∇dS V k22 ≤ c2 kγk22 , and, from here and on, we write ∼ to denote such norm equivalency. Recognize that γ ∈ RdS is the coefficient of the vector ∇d (v ◦ cp)(y) ∈ Ty S S of Ty S. Thus, we have k∇dS V k22 ∼ kγk22 ∼ k∇d (v◦cp)k22 , with respect to the basis {τj }dj=1 or equivalently, we write X (2.14) |α|≤1  DdαS vi xd =ai (x1 ,...,xdS )  2 X  ∼ |α|≤1  Ddα (vi ◦ cp) 2 xd =ai (x1 ,...,xdS ) . So, (2.10) and (2.11) are equivalent for |α| ≤ 1. We finish the proof by induction. Suppose (2.14) holds for |α| ≤ 1, . . . , k − 1. Then, the statement also holds true for |α| ≤ k as demonstrated as follows: X |α|≤k  DdαS vi xd =ai (x1 ,...,xdS ) ∼  X + =  X +  X + ∼ X  |α|≤k d X X  d X X  =  Ddα (vi ◦ cp) +  DdαS ∂j vi DdαS d X X  j=1 |α|=k−1  X |α|≤k−1 j=1 |α|=k−1 |α|≤k−1 |α|≤k−1 2 j=1 |α|=k−1 |α|≤k−1 ∼  |α|=k  DdαS vi xd =ai (x1 ,...,xdS ) h ∂v Ddα X  j ∂xj nh ∂v j ∂xj xd =ai (x1 ,...,xdS )  xd =ai (x1 ,...,xdS ) 2 2 + ∂vj ∂ai i ∂xd ∂xj + o ∂vj ∂ai i ◦ cp ∂xd ∂xj 2  xd =ai (x1 ,...,xdS )  2 2 xd =ai (x1 ,...,xdS ) , provided all functions (including cp and ai ) are smooth enough for the partial derivatives to 5 exit.  3. Semi-discretized trial solution and its convergence. We can apply the method of lines to discretize the surface diffusion equation (1.1) spatially by some kernel-based trial space. In this paper, we focus on manifold kernels that can be obtained by restricting [12, 13] some global, symmetric positive definite, and Sobolev space H m+1/2 (Rd ) reproducing [14] kernels Φm+1/2 : Rd × Rd → R to S. Fourier transforms of such kernels Φ̂τ decay like (3.1) c1 (1 + kωk22 )−(m+1/2) ≤ Φ̂m+1/2 (ω) ≤ c2 (1 + kωk22 )−(m+1/2) for all ω ∈ Rd , for some constants 0 < c1 ≤ c2 . Define the manifold kernels Ψm : S × S → R by Ψm (·, ·) := Φm+1/2 (·, ·)|S×S . (3.2) Then, for any m > dS /2 = (d − 1)/2, this manifold kernel Ψm reproduces H m (S), see [12]. In practice, one can use the standard Whittle-Matérn-Sobolev kernels [14] or the Wendland compactly supported kernels [15] with smoothness order m + 1/2, as is, for implementation. 3.1. Trial spaces. Let Z = {z1 , . . . , znZ } ⊂ S be the set of trial centers. The fill distance hZ and the separation distance qZ are defined respectively as (3.3) hZ := sup inf dist(ζ, η) and qZ := ζ∈S η∈Z 1 inf dist(zi , zj ). 2 zi 6=zj ∈Z The mesh ratio is ρZ = hZ /qZ ≥ 1. We say that surface points in Z are quasi-uniform if, as nZ := |Z| increases, Z satisfies qZ ≤ hZ = ρZ qZ ≤ ρqZ , for some constant ρ > 0 independent of nZ . Because of equivalence between them on some smooth and compact surfaces [12, Thm. 6], we can simply use Euclidean distant as dist(·, ·) instead of geodesic distant. We work on finite-dimensional trial spaces UZ , that is a span of the translation-invariant Ψm in (3.2) to Z, defined by (3.4) UZ,S,Ψm := span{Ψm (· , zj ) | zj ∈ Z}. Any (time-dependent) trial function u ∈ UZ,S,Ψm can be expressed by a linear combination with a set of unknown coefficients (function of time) {λj } (3.5) uZ = nZ X λj Ψm (·, zj ) =: Ψm (·, Z)λZ . zj ∈Z The associated reproducing kernel Hilbert space (a.k.a native space) NΨm (S)-norm is given 6 by nZ X kuZ k2NΨm (S) = (3.6) Ψm (zi , zj )λi λj =: λTZ Ψm (Z, Z)λZ . zi ,zj ∈Y m (S) as By [16, 17] and Lem. 2.1, we know that kuZ kNΨm (S) ∼ kuZ kW2m (S) ∼ kuZ kHcp m (S)-norm in place of native space norm defined in (2.9). Throughout the paper, we use Hcp for simplicity and norm equivalency will be taken care by generic constants. 3.2. Strong form collocation and overtesting. Let Y = {y1 , . . . , ynY } ⊂ S, nY > nZ , be a sufficiently dense set of quasi-uniform collocation points. We use strong from collocation at Y as test conditions in order to identify a numerical solution form the trial space UZ,S,Ψm . Let the discrete ℓ2 (Y )-norm be kf k2ℓ2 (Y ) := (3.7) X |f (y)|2 . y∈Y Now, we formally define the time-dependent semi-discretized trial solution, for some regularization parameter α, by (3.8) uZ,α (·, t) := h2Y arg inf u(·,t)∈UZ,S,Ψm Z T ku̇ + LS u − f k2ℓ2 (Y ) dt 0  +ku(·, 0) − gk2ℓ2 (Y ) + α2 ku(·, 0)k2H m (S) . cp Here, we focus on regularized least-squares initial condition. For sufficiently smooth g, one can use the interpolant of g from UZ,S,Ψm in (3.4) as initial condition and define the numerical solution only without initial values in (3.8). This modification will not affect the convergence analysis in this paper. Putting the trial function in the form of (3.5) into the PDE residual, u̇ + LS u − f , and evaluating at yi ∈ Y yields the following equations nZ X Ψm (yi , zj ) zj ∈Z nZ X d [LS Ψm ](yi , zj )λj (t) − f (yi , t) for t ∈ (0, T ], λj (t) + dt zj ∈Z or in an overdetermined matrix as Ψm (Y, Z)λ̇Z (t) + [LS Ψm ](Y, Z)λZ (t) − f (Y, t) for t ∈ (0, T ]. Thus, we can recast (3.8) in terms of the unknown coefficient functions as (3.9a) λZ,α (t) := arg inf λ∈((L∞ ∩H 1 )[0,T ])|Z| Z T Ψm (Y, Z)λ̇ + [LS Ψm ](Y, Z)λ − f (Y, t) 0 7 2 dt ℓ2 (R|Y | ) subject to initial condition λZ,α (0) := arg inf Ψm (Y, Z)λ − g(Y ) (3.9b) λ∈R|Z| 2 ℓ2 (Y ) + α2 λT Ψm (Z, Z)λ. Note that the initial condition for λ̇Z,α (0) can be determined by differentiating (3.9b) or by
other means, say u̇Z,α (·, 0) = Π f (·, 0) − LS g with some appropriate projection Π to the trial space. Then, the unique solvability of (3.8) is guaranteed by theorems in calculus of variations. In the remaining of this section, we will prove convergence estimates by proving a regularity estimate for surface diffusion (1.1a) in Sect. 3.3, an error estimate of regularized initial condition (3.9b) in Sect. 3.4, and finally in Sect. 3.5, the convergence of the trial function uα (·, t) := Ψm (·, Z)λZ,α (t), which is defined via the solution of (3.9a), to the exact solution u∗ (·, t). We summarize the main result in Thm. 3.6. 3.3. Regularity estimates for surface diffusion. We will walk through the details in proof to identify all necessary assumptions so that 1 (S) + kukL2 (0,T ;H µ (S)) + ku̇kL2 (0,T ;H 0 (S)) (3.10) ess sup kukHcp cp cp 0≤t≤T   0 (S)) + kgkH 1 (S) , ≤ C kf kL2 (0,T ;Hcp cp for µ = 1, 2, holds for the solution u to the surface diffusion in (1.1) with Sobolev norms defined in Sect. 2.3. After we prepare the manifold versions of all required (in)equalities, regularity estimate (3.10) can be shown by applying some standard arguments. We begin with the Green’s first identity. L EMMA 3.1. Let S ⊂ Rd be a C 3 -smooth closed manifold and v, w : S → R some C 2 -smooth scalar value functions. Then we have Z v∆S (Aw) dσ = − S Z ∇S v · (A∇S w) dσ. S for any diffusion tensor A satisfying Assumption 1. Proof. Without loss of generality, we can ignore the time-dependency in the diffusion tensor A in this proof. We cp-extend all functions to the narrow band domain Ω ⊃ S in (2.2) and denote vcp := v ◦ cp, wcp := w ◦ cp, and Acp := A ◦ cp. We start with a vector calculus inequality1 : (3.11) 1∇ ∇ · (vcp Acp ∇wcp ) = vcp ∇ · (Acp ∇wcp ) + ∇vcp · (Acp ∇wcp ). · (f F ) = f ∇ · F + (∇f ) · F 8 All smoothness requirements together ensure (3.11) is well-defined and vcp Acp ∇wcp is a smooth vector field. Applying divergence theorem on a closed manifold [18, Thm. 1.1], i.e., R X dσ = 0, to (3.11) yield S 0= Z ∇ · (vcp Acp ∇wcp ) dσ = Z vcp ∇ · (Acp ∇wcp ) + ∇vcp · (Acp ∇wcp ) dσ. S S The proof is completed after simplification by definitions (2.4) and (2.5).  Using Lem. 3.1, we test (1.1a) by u to obtain 1 d 2 dt Z Z 2 u dσ + k∇S uk2A S S dσ = Z uf dσ, S where k∇S uk2A := ∇S u · (A∇S u). Following by some standard arguments, i.e., integrate over time Z  Z tZ Z Z tZ 1 uf dσdτ, u2 dσ − g 2 dσ + k∇S uk2A dσdτ = 2 S S 0 S 0 S and apply Cauchy inequality ab ≤ 2ǫ a2 + Z tZ (3.12) 0 uf dσdτ ≤ T S T Z 1 2 2ǫ b Z with ε = 2T , yields f 2 dσdτ + S 0 1 ess sup 4 t∈[0,T ] Z u2 dσ. S We obtain an energy estimate (3.13) ess sup 0≤t≤T Z u2 dσ + S T Z Z S 0 k∇S uk2A dσdτ ≤ C Z g 2 dσ + S Z T 0 Z S  f 2 dσdτ , for some constant C > 0 independent of the solution u to (1.1). Moreover, due to Assump0 (S) are equivalent with constants depending on tion 1, we know that k∇S ukA and k∇S ukHcp the eigenvalues of A restricted to the tangent space Ty S. To get an improve regularity, we follow [19, Thm.5, Ch.7] and suppose that the diffusion tensor A does not depend on t. Testing the PDE by u̇ yields d u̇ dσ + dt S Z 2 Z S 1 k∇S uk2A dσ = 2 Z u̇f dσ, S which, after an integration in time, will lead us to the estimate (3.14) ess sup k∇S uk2A 0≤t≤T + Z 0 T Z 2 u̇ dσdτ ≤ C Z S S k∇S gk2A dσ + Z 0 T Z S  f 2 dσdτ . For strong solution, the last required estimate comes form the fact that LS u(·, t) ∈ L2 (S) for 9 a.e. 0 ≤ t ≤ T and, by elliptic regularity [2, Lem. 2.1 & 3.1] of the surface operator LS and (3.14), we have Z 0 T kuk2Hcp 2 (S) dτ ≤ C Z T kLS uk2Hcp 0 (S) dτ 0 ≤C Z T 0 ≤C (3.15) Z 0 T   2 kf k2Hcp dτ 0 (S) + ku̇kH 0 (S) cp kf k2Hcp 0 (S) dτ + Using the trivial bound k · k2L2 (0,T ;H 0 cp (S)) Z S k∇S gk2A dσ. ≤ CS ess sup0≤t≤T k · k2H 0 cp (S)) to combine (3.13)–(3.15) completes the proof of the regularity estimate in (3.10). L EMMA 3.2. Assume the data functions in the surface diffusion equation (1.1) with time-independent diffusion tensor are sufficiently smooth, f ∈ L2 (0, T ; L2 (S)), and g ∈ H 1 (S), then the regularity estimate (3.10) holds on the unique strong solution u∗ ∈ (L2 ∩ L∞ )(0, T ; H 2 (S)) and u̇∗ ∈ L2 (0, T ; L2 (S)) to (1.1). 3.4. Regularized discrete least-squares initial condition. In this section, we show that the solution to the following time independent regularized approximation problem gZ,α := arg inf ku − gk2ℓ2 (Y ) + α2 kuk2Hcp m (S) (3.16) u∈UZ,S,Ψm 1 (S)-convergent so that, later in our convergence analyin (3.8) with some appropriate α is Hcp 1 (S) in the right handed side of sis, we can show consistency by some bound on kgZ,α − gkHcp (3.10). To do so, we extend the convergence estimate in [20] for flat geometry to manifolds. L EMMA 3.3. For any m ≥ ⌊1 + dS /2⌋ and α ≥ 0, suppose that the sets Y, Z ⊂ S of sufficiently dense discrete data points satisfy (3.21). Let gZ,α ∈ UZ,S,Ψm be the regularized m (S). Then, approximant in (3.16) for g ∈ Hcp   −2 2m−dS 2m−3 kgZ,α − gk2Hcp + hdYS −2 α2 + hZ kgk2Hcp 1 (S) ≤ C hY hZ m (S) holds for some constant C independent of g. m−dS /2 When α ≥ α∗ := hY , we have   2m−3 −2 2m−dS kgZ,α∗ − gk2Hcp + hY2m−3 + hZ kgk2Hcp 1 (S) ≤ C hY hZ m (S) , without any extra denseness restriction imposed on Z and Y . m Proof. We begin with a sampling inequality for manifold functions g ∈ Hcp (S), see [21, Lem. 3.1] (with m = m + 1/2 and s = 3/2); for any discrete set Y ⊂ S with sufficiently 10 small mesh norm hY , the following holds (3.17)   dS −2 2 S kgk2Hcp kgk kgk2ℓ2 (Y ) + h2m−d m 1 (S) ≤ ChY Hcp (S) Y Let IZ g be the interpolant of g from the trial space UZ,S,Ψm . We have the following convergence estimate in [22, Cor.13] and [2, Sect.3]: (3.18) m−k−dS (1/2−1/q)+ kIZ g − gkWqk (S) ≤ ChZ kgkNΨm (S) for 0 ≤ k ≤ ⌈m − dS (1/2 − 1/q)+ ⌉ − 1. Using (3.18) with k = 1, it is sufficient to demonstrate the convergence gZ,α → IZ g in the trial space UZ,S,Ψm , within which we know a Bernstein’s inverse inequality [7, Thm. 10] −k kukHcp k (S) ≤ Ch Z kukL2 (S) , (3.19) 0≤k≤m holds for all u ∈ UZ,S,Ψm . m (S) to For any α ≥ 0, we use (3.17) and (3.19) on e := gZ,α − IZ g ∈ UZ,S,Ψm ⊂ Hcp establish a chain of upper bounds:   dS −2 S kek2Hcp kek2Hcp kek2ℓ2 (Y ) + h2m−d m (S) 1 (S) ≤ ChY Y   2m−dS = ChdYS −2 kek2ℓ2 (Y ) + α2 kek2Hcp − α2 )+ kek2Hcp m (S) m (S) + (hY   2m−dS − α2 )+ h−2m kek2L2 (S) ≤ ChdYS −2 kek2ℓ2 (Y ) + α2 kek2Hcp m (S) + (hY Z   (3.20) ≤ ChdYS −2 kek2ℓ2 (Y ) + α2 kek2Hcp m (S) , m−dS /2 either if α ≥ α∗ := hY (3.21) , or under the deseness and regularization constraint −2m S ChYdS −2 (h2m−d − α2 )+ hZ < Y 1 2 Add-in subtract-out the function g in (3.20); a direct consequence of the optimality in (3.16) is that (3.22)   dS −2 kek2Hcp kIZ g − gk2ℓ2 (Y ) + α2 kIZ g − gk2Hcp 1 (S) ≤ 2ChY m (S) . A scattered zero lemma [22, Lem.10] suggests that S kIZ g − gk2ℓ2 (Y ) ≤ nY kIZ g − gk2L∞ (S) ≤ CnY h2m−d kIZ g − gk2Hcp m (S) , Z provided m ≥ ⌊1 + dS /2⌋. Because IZ g and IZ g − g are mutually orthogonal with respect m m (S) and kIZ g − to the native space norm, i.e., Hcp (S)-norm, we know that both kIZ gkHcp m (S) ≤ kgkH m (S) . Under the assumption that Y is quasi-uniform and therefore nY ≤ gkHcp cp 11 S cqY−dS ≤ Ch−d , we arrive at Y S 2m−dS hZ kIZ g − gk2ℓ2 (Y ) ≤ Ch−d kgk2Hcp m (S) Y (3.23) and complete the proof.  3.5. Convergence estimate for semi-discretized solution. m Define a residual functional Eµ : L2 (0, T ; Hcp (S)) → R using the left handed side of the regularity estimate (3.10) by (3.24) 2 2 Eµ (u) := ess sup kuk2Hcp 1 (S) + kukL2 (0,T ;H µ (S)) + ku̇kL2 (0,T ;H 0 (S)) , cp cp for µ = 1, 2, 0≤t≤T and restate (3.10) as   2 Eµ (u) ≤ C ku̇ + LS uk2L2 (0,T ;Hcp 0 (S)) + ku(·, 0)kH 1 (S) . cp (3.25) We begin by showing that there is a good comparison function in the trial space. For any t ∈ [0, T ], let s(·, t) := IZ u∗ (·, t) be the interpolant of the solution u∗ to the PDE (1.1) from the trial space UZ,S,Ψm . Because (3.26) s(zj , t) = IZ u∗ (zj , t) = Ψm (zj , Z)[Ψm (Z, Z)]−1 u∗ (Z, t) for all zj ∈ Z and t ∈ [0, T ], differentiating with respect to t shows that ṡ(zj , t) = Ψm (zj , Z)[Ψm (Z, Z)]−1 u̇∗ (Z, t) for all zj ∈ Z and t ∈ [0, T ], and therefore ṡ(·, t) ∈ UZ,S,Ψm is the unique interpolant of the first order time derivative u̇∗ (·, t). Using (3.25) and standard interpolation theories [22, Cor. 13], we get the following error estimate: Eµ (s − u∗ )   ∗ 2 ∗ 2 ≤ C kṡ − u̇∗ k2L2 (0,T ;Hcp 0 (S)) + kLS (s − u )kL2 (0,T ;H 0 (S)) + k(s − u )(·, 0)kH 1 (S) cp cp   2m−4 ∗ 2 ∗ 2 2 ∗ 2 ≤ ChZ ku̇ kL2 (0,T ;H m−2 (S)) + ku kL2 (0,T ;Hcp m (S)) + hZ ku (·, 0)kH m (S) . cp cp Our next task is to show that the numerical solution defined by (3.8) converges to the comparison function s, i.e., Eµ (uZ,α − s) → 0 as hZ → 0, in the trial space. Below are two essential stability estimates before we can study consistency. The first comes from (3.20). m−dS /2 C OROLLARY 3.4. Let m ≥ ⌊1 + dS /2⌋ and α ≥ 0. If α < α∗ := hY 12 , further suppose that the sets Y, Z ⊂ S of sufficiently dense discrete data points satisfy (3.21). Then,   dS −2 kuk2Hcp kuk2ℓ2 (Y ) + α2 kuk2Hcp 1 (S) ≤ ChY m (S) , holds for all trial function u ∈ UZ,S,Ψm . L EMMA 3.5. Let m ≥ ⌊3 + dS /2⌋. Suppose that the sets Y, Z ⊂ S of sufficiently dense discrete data points satisfy (3.28). Let Eµ (u) be the residual functional defined in (3.24). Then, the estimate Z  Eµ (u) ≤ C hdYS 0 T ku̇ + LS uk2ℓ2 (Y ) dt + ku(·, 0)k2Hcp 1 (S)  holds for all trial function u ∈ L2 (0, T ; UZ,S,Ψm ) and µ = 1, 2. Proof. For any t ∈ [0, T ], applying the sampling inequality in [21, Lem. 3.1] (with m = m−2 m + 1/2 and s = 1/2) to u̇(·, t) + LS u(·, t) ∈ Hcp (S) for any time t ∈ [0, T ] in (3.25) at m (S) due m−2 the set Y , then inverse inequality (3.19) and the estimate kLS ukHcp (S) ≤ kukHcp to boundedness of PDE coefficients, yields   dS 2 S −4 ku̇ + LS uk2ℓ2 (Y ) + h2m−d k u̇ + L uk ku̇ + LS uk2Hcp 0 (S) ≤ ChY m−2 S Y Hcp (S) 
 −2m −2m+4 S −4 ≤ ChdYS ku̇ + LS uk2ℓ2 (Y ) + h2m−d . kuk2Hcp ku̇k2Hcp hZ 0 (S) 0 (S) + hZ Y (3.27) Under the denseness constraints ( (3.28) −2m ChY2m−4 hZ < −2m+4 ChY2m−4 hZ < 1 4 1 4 if hZ ≤ 1, otherwise, we integrate (3.27) from t = 0 to T to finish the proof.  Finally, we arrive the main result for the convergence of the proposed semi-discretized solution to the method of lines and kernel-based strong form collocation method. m T HEOREM 3.6. For some m ≥ ⌊3+dS /2⌋, let u∗ ∈ (L2 ∩L∞ )(0, T ; Hcp (S)) and u̇∗ ∈ 2 m−2 L (0, T ; Hcp (S)) be the solution to (1.1). For some surface S satisfying the assumptions in Sect. 2, kernel Ψm satisfying (3.1), and two sets of quasi-uniform points Y, Z ⊂ S satisfy m−dS /2 assumptions mentioned before and (3.28). If α < α∗ := hY , further suppose that Y and Z satisfy (3.21). Let uZ,α ∈ UZ,S,Ψm be the numerical solution defined in (3.8) and 13 Eµ (u) the residual functional defined in (3.24). Then, the following estimate holds 
2m−4−dS Eµ (uZ,α − u∗ ) ≤ C hZ ku̇∗ k2H m−2 (S) + ku∗ k2Hcp m (S) cp + 2m−2 hZ 
2m−dS + h−2 + hdYS −2 α2 ku∗ (·, 0)k2Hcp m (S) Y hZ for µ = 1, 2 and some constant C independent of u∗ . Proof. Consider Eµ (uZ,α − u∗ ) ≤ Eµ (sZ,α − s) + Eµ (s − u∗ ); the latter were analyzed in (3.27). It remains to show consistency in the first term via interpolation theories. Applying Lem. 3.5, Cor. 3.4, followed by a triangle inequality to the trial function (uZ,α − s) yields Z  Eµ (uZ,α − s) ≤ C hdYS T k(u̇Z,α − ṡ) + LS (uZ,α − s)k2ℓ2 (Y ) dt + k(uZ,α − s)(·, 0)k2Hcp 1 (S) 0 ≤C  hdYS Z 0 T k(u̇Z,α − ṡ) + LS (uZ,α − s)k2ℓ2 (Y ) dt 2 +hYdS −2 kuZ,α − sk2ℓ2 (Y ) + α2 kuZ,α k2Hcp m (S) + kskH m (S) cp
 . By Add-in and subtract-out f = u̇∗ + LS u∗ and u∗ (·, 0) into the first and second norms on the right handed side, we obtain an estimate that only depends on the interpolant s in (3.26) to u∗ Z  Eµ (uZ,α − s)≤ 2C hdYS 0 (3.29) T kṡ + LS s − f k2ℓ2 (Y ) +hdYS −2 ks(·, 0) − gk2ℓ2 (Y ) + α2 ks(·, 0)k2Hcp m (S)
 by the optimality (3.8) in the numerical solution. Using the convergence estimate in (3.18) with q = ∞, we have (3.30) kṡ + LS s − f k2ℓ2 (Y ) ≤ kṡ − u̇∗ k2ℓ2 (Y ) + kLS s − LS u∗ k2ℓ2 (Y )   ∗ 2 ∗ 2 S + kL s − L u k ≤ Ch−d k ṡ − u̇ k ∞ ∞ S S L (S) L (S) Y   −dS 2m−4−dS ≤ ChY hZ ku̇∗ k2H m−2 (S) + ku∗ k2Hcp m (S) . cp Next, (3.23) suggests (3.31) S 2m−dS ks(·, 0) − gk2ℓ2 (Y ) ≤ Ch−d hZ kgk2Hcp m (S) , Y and the orthogonality of interpolant in native space [23, Ch.10] suggests (3.32) 2 ks(·, 0)k2Hcp m (S) ≤ kgkH m (S) . cp Putting (3.30)–(3.32) into (3.29) completes the proof. 14   4. Fully discretized solution. Recall the semi-discretized problem in terms of coefficients in (3.9a) is in the from of (4.1) λZ,α (t) := arg inf λ∈((L∞ ∩H 1 )[0,T ])|Z| Z T Aλ̇ + Bλ − C 0 2 dt ℓ2 (R|Y | ) for 0 ≤ t ≤ T subject to some predetermined initial condition λZ,α (0). Suppose we numerical solve (4.1) at some partition T := {tj } of [0, T ] for λj that approximate λZ,α (tj ) by some order-p scheme. Denote the corresponding fully discretized solution at any tj ∈ T by (4.2) UT ,Z,α (·, tj ) := Ψ(·, Z)λj . Let ζj := λj − λZ,α (tj ) ∈ R|Z| and we have kζk22 = O(h2p T ). We can measure difference between fully and semi-discretized solutions by kUT ,Z,α (·, tj ) − uZ,α (·, tj )k2L2 (S) = ζjT hZ S i Ψ(·, Z)T Ψ(·, Z) dσ ζj = O(h2p T ). We now derive an error estimate by a sequence of comparison kUT ,Z,α (·, tj ) − u∗ (·, tj )k2ℓ∞ (T ;L2 (S)) ≤ kUT ,Z,α (·, tj ) − uZ,α (·, tj )k2ℓ∞ (T ;L2 (S)) +kuZ,α (·, tj ) − u∗ (·, tj )k2ℓ∞ (T ;L2 (S)) ≤ kUT ,Z,α (·, tj ) − uZ,α (·, tj )k2ℓ∞ (T ;L2 (S)) + ess sup kuZ,α (·, t) − u∗ (·, t)k2H 1 (S) 0≤t≤T (4.3) ≤ O(h2p T ) + Eµ (uZ,α − u∗ ), allowing Theorem 3.6 to be applied to the fully discretized solution. 4.1. Temporal discretization. Despite being linear, the second order Euler-Lagrange ODE for (4.1) involves multiple products of kernel matrices, which is ill-conditioned. We propose two variants of algorithms for solving (4.1): (i) difference equation and (ii) ODE approach. We focus on uniform time discretization at (4.4) |T | |T | T = {tj }j=0 := {jhT }j=0 . Further suppose that the integrand of (4.1), i.e., squares of norms of PDE residuals, is of class C p [0, T ]. This assumption imposes temporal smoothness requirements on f in (1.1) and hence to u∗ . Then, there exists a order-p backward finite difference scheme DhT , see [24, Tab. 15 3] for coefficients, such that −1  X 1  γk λk+j + O(hpT ), λ˙ (tj ) = DhT λj + O(hpT ) = γ0 λj + hT k=−p and an order-p numerical quadrature rule {wj } on T , say composite Newton-Cotes rules, so that, given λZ,α (0), we can approximate the set of solutions (4.1) at T by λZ,α (T ) := ( arg inf 1 )[0,T ])|Z| λ∈((L∞ ∩Hcp = (4.5)     arg inf  |Z|   λj ∈R 1≤j≤|T | X Z T Aλ̇ + Bλ − 0 wj ADhT λj + Bλj − C j tj ∈T 2 C 2 |Y | dt ℓ (R ) )     2 dt ℓ2 (R|Y | )  T + O(hpT ).   Finally, we arrive at a sum-of-squares problem in (4.5), in which we found standard leastsquares problems with respect to λj . Explicitly, we write the solution to (4.5) by the following recursive formula: (4.6a) λZ,α (T ) = ( arg inf ADhT λj + Bλj − λj ∈R|Z| (4.6b) =    arg inf  λj ∈R|Z| (4.6c) =   γ 0  hT γ 0 hT A+B A + B λj + † C j )|T | + O(hpT ) j=1   2 C j 2 |Y | ℓ (R ) 1 A hT −1  X k=−p 2  γj λk+j − C j ℓ2 (R|Y | ) |T | −1   1  X − A γj λk+j  + O(hpT ),  hT k=−p |T |   + O(hpT )   j=1 j=1 subject to initial condition λ0 = λZ,α (0) ∈ R|Z| . We summarize the above findings. T HEOREM 4.1. Suppose that the assumptions in Thm. 3.6 hold. For some integer p ≥ 1, m we further suppose that u∗ ∈ C p+1 (0, T ; Hcp (S)). For some p-th order backward finite difference scheme on uniform time grid T in (4.4) for the first derivative with finite difference coefficients γ−p , . . . , γ0 ∈ R. Let λ0 = λZ,α (0) ∈ R|Z| as in (3.9b) and (4.7)   −1 †  X γ  1 0 λj = Ψm (Y, Z) + [LS Ψm ](Y, Z) f (Y, tj ) − Ψm (Y, Z) γj λk+j  , hT hT k=−p for 1 ≤ j ≤ |T |, be the coefficients of the fully discretized solution UT ,Z,α (·, tj ) in (4.2). Then, the error estimate in Thm. 3.6 for semi-discretized solution extends to the fully dis16 cretized solution with an additional p-th order temporal error term, i.e., ∗ kUT ,Z,α (·, tj ) − u∗ (·, tj )k2ℓ∞ (T ;L2 (S)) ≤ O(h2p T ) + Eµ (uZ,α − u ) for holds. Note that the difference equations in (4.6c) and (4.7) were derived without any ODE in the process. To benefit from the vast library of ODE solvers with adaptive time stepping, we conclude our theoretical work with an equivalent ODE. C OROLLARY 4.2. Thm. 4.1 remains to hold if the difference solution (4.7) were replaced by some p-th order approximation to the solution of the following |Z| × |Z| systems of ODE   λ˙ (t) = Ψm (Y, Z)† f (Y, t) − [LS Ψm ](Y, Z)λ (t) (4.8)
|Z| . for λ ∈ C p+1 [0, T ] Proof. Using notations in (4.6c), the least-squares solution (4.7) solves the normal equation γ 0 hT A+B T   −1  γ   X  0 A + B λj + 1 A γj λk+j − C j  = 0 hT hT k=−p for any give partition T . By reverse-manipulation from (4.6b) back to (4.6a), we can rewrite the normal equation as below and take limit of hT = T /|T | → 0: lim hT →0  γ0 A + hT B T lim hT →0   ADhT λj + Bλj − C j = 0, 1 ≤ j ≤ |T |, to obtain   AT Aλ˙ (t) + Bλ (t) − C (t) = 0, which is the normal equation of (4.8). 0 < t ≤ T,  5. Numerical examples. This section contains examples with self-explanatory title designed for various numerical verifications and aims. We use y = (x1 , x2 , x3 ) ∈ S to define functions throughout the section. Example 1: Comparing with a meshless Galerkin method. Our first aim is to compare the proposed method with a weak formulation in [3]. We consider an example there: u̇ − 0.1∆S + 3u = f on the unit sphere with u∗ (y, t) = exp(x1 + 1/(1 + t)) for t ∈ [0, 1]. We solve the problem by m = 4 ≥ ⌊3 + dS /2⌋ Sobolev kernel and the difference equations in (4.7) in Thm. 3.6 with p = 2 for a fair comparison with the Crank-Nicolson (CN) scheme used in the Galerkin results. Using the same number of data points |Z|, time stepping 17 TABLE 1 Exmp. 1: Relative errors of discrete L2 norm at T = 1 and the corresponding estimated order of convergence (eoc) of (4.7) with p = 2 and m = 4. hT 0.06 0.04 0.02 0.01 Proposed method |Z| = 961 |Z| = 3721 |Y | = 1, 153 |Y | = 4, 465 1.198428E-4 1.198429E-4 5.602314E-5 5.602329E-5 1.250487E-5 1.250503E-5 2.927779E-6 2.927935E-6 eoc 2.09 2.16 1.88 Galerkin [3, Tab. 2] |Z| = 961 |Z| = 3721 |Y | = 23, 042 |Y | = 256, 002 1.251055E-4 1.246096E-4 6.055402E-4 5.953572E-4 1.869301E-5 1.507104E-5 1.185031E-5 4.258655E-6 hT , and same ways to calculate the relative error, we report in Tab. 1 the |Y | ≈ 120%|Z| oversampling results, which is selected to match the best accuracy in [3, Tab. 2]. It is worth noting that we use orders of magnitude fewer collocation points than quadrature points in the Galerkin formulation. The estimated orders of convergence (eoc) support the theoretical temporal convergence rate in Thm. 4.1. Example 2: Nonconstant diffusive tensor and kernels’ smoothness. Using the projection matrix P in (2.3), let A(y) = P (y) Diag[x21 + 1, 1, 1] ∈ R3×3 be the diffusion tensor. We solve the diffusion equation (1.1) on the unit sphere. We use u∗ in Exmp. 1 as exact solution. This time around, we solve the PDE via ODE (4.8) in Cor. 4.2. We perform a shortterm temporal integration by CN with a small time-step size hT = 1E − 12 in order to study the spatial convergence behaviour of the proposed method with respect to various kernel’s smoothness 2 ≤ m ≤ 7. Relative L2 (S) error by using |Z| = 100 to 1000 trial centers and 120% oversampling were shown in Fig. 1. Generally speaking, larger m yields faster (inital) eoc before error stagnation. Although m = 2 and 3 are not covered by Thm. 3.6, convergence is still observed. Also, it is not uncommon to observe faster than theoretically predicted eoc when u∗ and S are both of high order of smoothness, see numerical experiment in [1, Exmp. 2]. Example 3: Simulating Allen-Cahn equations. In this numerical experiment, we consider the Allen-Cahn equation u̇ = ∆S u+ ε12 u(1−u2 ), which is a reaction diffusion equation that models phase separation of two fluids [25] with ε > 0 being the width of the diffusion interface between two fluids. We set up the test problem on the unit sphere as in [3, 26]. For some R0 > 0, we define initial condition u(0) to be +1 within radius R0 to the northpole p and −1 otherwise. Then, R(t) = 1 − (1 − R02 )e2t is the radius of the shrinking spherical cap of the solution for 0 ≤ t ≤ T := − 21 log(1 − R02 ). Using R0 = 0.717, ε = 0.05, and |Z| = 3721 as in [3] for another round of comparison with the Galerkin method, we solve the semi-discrete equation (4.8) by the MATLAB built-in ODE45 solver. Firstly, we compute the econ-QR factorization of Ψm (Y, Z) = QR where Q and R are of size |Y | × |Z| and |Z| × |Z|. Replacing Ψm (Y, Z)† in (4.8) by QT , we call 18 Radius, R(t) Rel. L2 (S) error #timestep ❍ ❥ ❍ t hX F IG . 1. Exmp. 2: Spatial convergence profiles F IG . 2. Exmp. 3: Radius the spherical cap and with respect to kernel smoothness. number of time stepping used by ODE45. t = 5E − 4 0.01 0.20 0.30 0.35 F IG . 3. Exmp. 3: Snapshots of the shrinking cap towards the north pole. ODE45 with mass matrix R. We test oversampling ratios 120% and 200% with |Y | = 4465 and 7442; Fig. 2 (Left-y) shows the analytical and numerical radius of the spherical cap, which should be compared with [3, Fig. 3]. First of all, both kernel-based methods of lines show good accuracy for a long time. If we want to be really picky, • the proposed method shows small oscillations in radius at small t due to fewer numbers of collocation points, but • the Galerkin solutions with |Y | = 40962 and 92162 have visually different (without zoom-in) final cap-vanishing times. Fig. 2 (Right-y) shows the numbers of time steps required and Fig. 3 shows a few snapshots of the shrinking spherical cap. We end this example with an Allen-Cahn phase separation simulation on a torus2 . Using the generator in [27] with hZ1 = 0.1333 and hZ2 = 0.1, and hX = 0.066, we obtain point sets of sizes |Z1 | = 864, |Z2 | = 1696, and |X| = 3856 to repeat the same calculations. Random initial values between [−0.5, 0.5] were assigned to X and we present in Fig. 4 a selected initial condition that leads to steady state solutions with non-vanishing diffuse in2 Torus: x2 + y 2 + z 2 + 12 − (1/3)2 )2 − 4(x2 + y 2 ) = 0 19 0.02 0.5 1.0 4.0 |Z2 | = 1696 |Z1 | = 864 t = 0.006 F IG . 4. Exmp. 3: Snapshots of Allen-Cahn solutions based on the same random initial condition and collocation points with size |X| = 3856, but different sets of trial centers. terface. Comparing results in the eyeball norm, two trial solutions look obviously different at t = 0.5. We also note that the solution of Z2 with more trial centers arrives steady state earlier (c.f. t = 1 and 4). 6. Conclusions. We theoretically deduce a method of lines based on a discrete leastsquares approximation to the initial condition and strong-form collocation with sampling for diffusion equations on smooth surfaces. Convergence estimate of the (spatially) semidiscretized solution is given in Thm. 3.6. The fully discretized problem is a difference equation (see Thm. 4.1) that were connected to a system of ODE (see Cor. 4.2). Numerically, provided examples mainly compare the proposed method with a meshless Galerkin method and show that our method can reproduce similar accuracy with fewer collocation points (quadrature points in the Galerkin case). We remind readers that our strong-form theories require solutions with higher smoothness, i.e., m ≥ ⌊3 + dS /2⌋. For PDEs with lower regularity, meshless Galerkin methods should be the theoretically sound method of choice. Acknowledgements. This work was supported by the Hong Kong Research Grant Council GRF Grants and the National Science Funds of China (Grant No. 12001261). REFERENCES [1] K. C. Cheung, L. Ling, A kernel-based embedding method and convergence analysis for surfaces PDEs, SIAM J. Sci. Comput. 40 (1) (2018) A266–A287. [2] M. Chen, L. Ling, Extrinsic meshless collocation methods for PDEs on manifolds, SIAM J. Num. Anal. 58 (2) (2020) 988–1007. [3] J. Kúnemund, F. J. Narcowich, J. D. Ward, H. Wendland, A high-order meshless Galerkin method for semilinear parabolic equations on spheres, Numer. Math. 142 (2019) 383–419. [4] Y. C. Hon, R. Schaback, M. Zhong, The meshless kernel-based method of lines for parabolic equations, Comput. Math. Appl. 68 (12A) (2014) 2057–2067. [5] Y. Hon, R. Schaback, Direct meshless kernel techniques for time-dependent equations, Appl. Math. Comp. 258 (2015) 220–226. [6] T. Maerz, C. B. Macdonald, Calculus on surfaces with general closest point functions, SIAM J. Numer. Anal. 50 (6) (2012) 3303–3328. [7] T. Hangelbroek, F. J. Narcowich, C. Rieger, J. D. Ward, Direct and inverse results on bounded domains 20 [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] for meshless methods via localized bases on manifolds, in: J. Dick, F. Kuo, H. Wozniakowski (Eds.), Contemporary Computational Mathematics - A Celebration of the 80th Birthday of Ian Sloan, Springer, 2018, pp. 517–543. S. J. Ruuth, B. Merriman, A simple embedding method for solving partial differential equations on surfaces, J. Comput. Phys. 227 (3) (2008) 1943–1961. R. S. Strichartz, Analysis of the Laplacian on the complete Riemannian manifold, J. Funct. Anal. 52 (1983) 48–79. N. Yoshida, Sobolev spaces on a Riemannian manifold and their equivalence, J. Math. Kyoto Univ. 32 (3) (1992) 621–654. J. Wloka, Partial Differential Equations, Cambridge University Press, 1987. E. Fuselier, G. B. Wright, Scattered data interpolant on embedded submanifolds with restricted positive definite kernels: Sobolev error estimates, SIAM J. Numer. Anal. 50 (3) (2012) 1753–1776. F. J. Narcowich, X. Sun, J. D. Ward, Approximation power of RBFs and their associated SBFs: A connection, Adv. Comput. Math. 27 (1) (2007) 107–124. B. Matérn, Spatial variation, Vol. 36, Springer Science & Business Media, 2013. H. Wendland, Error estimates for interpolation by compactly supported radial basis functions of minimal degree, J. Approx. Theory 93 (2) (1998) 258–272. M. D. Buhmann, Radial basis functions: Theory and implementations, Vol. 12 of Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, Cambridge, 2003. H. Wendland, Error estimates for interpolation by compactly supported radial basis functions of minimal degree, J. Approx. Theory 93 (2) (1998) 258–272. S. Pigola, A. G. Setti, Global divergence theorems in nonlinear PDEs and geometry, Ensaios Matemáticos 26 (2014) 1–77. L. C. Evans, Partial differential equations, Graduate studies in mathematics, American Mathematical Society, Providence (R.I.), 1998, reimpr. avec corrections : 1999, 2002. S. Li, L. Ling, K. C. Cheung, Discrete least-squares radial basis functions approximations, Appl. Math. Comput. 355 (2019) 542–552. K. C. Cheung, L. Ling, R. Schaback, H 2 –convergence of least-squares kernel collocation methods, SIAM J. Numer. Anal. 56 (1) (2018) 614–633. E. J. Fuselier, G. B. Wright, Scattered data interpolation on embedded submanifolds with restricted positive definite kernels: Sobolev error estimates, SIAM J. Numer. Anal. 50 (3) (2012) 1753–1776. H. Wendland, Scattered data approximation, Vol. 17 of Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, Cambridge, 2005. B. Fornberg, Generation of finite difference formulas on arbitrarily spaced grids, Math. Comput. 51 (184) (1988) 699–706. S. M. Allen, J. W. Cahn, A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening, Acta Metallurgica 27 (6) (1979) 1085–1095. Y. Choi, D. Jeong, S. Lee, M. Yoo, J. Kim, Motion by mean curvature of curves on surfaces using the AllenCahn equation, Int. J. Eng. Sci. 97 (2015) 126–132. P.-O. Persson, G. Strang, A simple mesh generator in MATLAB, SIAM Rev. 46 (2) (2004) 329–345. 21