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