Kernel-based collocation methods for heat transport on evolving surfaces

TL;DR

This paper introduces kernel-based meshless collocation algorithms for solving heat and surfactant transport PDEs on evolving surfaces, capable of handling high curvature and discontinuities.

Contribution

It presents two novel kernel-based meshless collocation methods tailored for PDEs on evolving surfaces, including direct PDE collocation and pseudospectral surface operator approximation.

Findings

Successfully solves PDEs on high-curvature surfaces

Handles discontinuous initial conditions accurately

Demonstrates convergence and robustness of methods

Abstract

We propose algorithms for solving convective-diffusion partial differential equations (PDEs), which model surfactant concentration and heat transport on evolving surfaces, based on intrinsic kernel-based meshless collocation methods. The algorithms can be classified into two categories: one collocates PDEs directly and analytically, and the other approximates surface differential operators by meshless pseudospectral approaches. The former is specifically designed to handle PDEs on evolving surfaces defined by parametric equations, and the latter works on surface evolutions based on point clouds. After some convergence studies and comparisons, we demonstrate that the proposed method can solve challenging PDEs posed on surfaces with high curvatures with discontinuous initial conditions with correct physics.

Tables

Table 1. Table 1: Example  4.1 . The L ∞ ​ ( [ 0 , T ] ; L 2 ​ ( 𝒮 ​ ( t ) ) ) superscript 𝐿 0 𝑇 superscript 𝐿 2 𝒮 𝑡 L^{\infty}\big{(}[0,T];L^{2}(\mathcal{S}(t))\big{)} errors and the corresponding estimated order of convergence (eoc) by ( 17 ) with respect to h Z subscript ℎ 𝑍 h_{Z} , obtained by Algorithm  2 with overdetermined formulas ( n X > n Z subscript 𝑛 𝑋 subscript 𝑛 𝑍 n_{X}>n_{Z} ), under the same other settings as in Figure  3 .

$h_Z$	$n_Z$	$n_X$	$L^{\infty }([0,T];L^2\left(𝒮\left(t\right)\right))$	eoc
$\sqrt{2}$	$5$	$7$	$2.31823$e-$02$	$-$
$2^{-1}\sqrt{2}$	$9$	$14$	$2.29316$e-$03$	$3.33762$
$2^{-2}\sqrt{2}$	$18$	$27$	$5.44931$e-$04$	$2.07319$
$2^{-3}\sqrt{2}$	$26$	$54$	$1.34549$e-$04$	$2.01794$

Table 2. Table 2: Example  4.1 . The L ∞ ​ ( [ 0 , T ] ; ℋ 2 ​ ( 𝒮 ​ ( t ) ) ) superscript 𝐿 0 𝑇 superscript ℋ 2 𝒮 𝑡 L^{\infty}\big{(}[0,T];\mathcal{H}^{2}(\mathcal{S}(t))\big{)} and L ∞ ​ ( [ 0 , T ] ; L 2 ​ ( 𝒮 ​ ( t ) ) ) superscript 𝐿 0 𝑇 superscript 𝐿 2 𝒮 𝑡 L^{\infty}\big{(}[0,T];L^{2}(\mathcal{S}(t))\big{)} errors and the corresponding convergence rates (eoc) with respect to h ℎ h , obtained by Algorithm  2 with exactly determined formulas and the unfitted FEM [ 15 ] , under the same other settings as in Figure  3 .

$h$	$n$	Algorithm 2 with $X=Z$				Unfitted FEMs [15]
		($△t=h/4$)				($△t=2h^2$)
		$L^{\infty }([0,T];$	eoc	$L^{\infty }([0,T];$	eoc	$L^{\infty }([0,T];$	eoc
		${\mathscr{H}}^2\left(𝒮\left(t\right)\right))$		$L^2\left(𝒮\left(t\right)\right))$		$L^2\left(𝒮\left(t\right)\right))$
$\sqrt{2}$	$5$	$1.14643$e-$1$	$-$	$3.85031$e-$2$	$-$	$-$	$-$
$2^{-1}\sqrt{2}$	$9$	$8.38177$e-$3$	$3.77376$	$2.60321$e-$3$	$3.88661$	$1.15457$e-$1$	$-$
$2^{-2}\sqrt{2}$	$18$	$2.00492$e-$3$	$2.06371$	$5.47333$e-$4$	$2.24980$	$3.25344$e-$2$	$1.82732$
$2^{-3}\sqrt{2}$	$26$	$5.23047$e-$4$	$1.93853$	$1.34576$e-$4$	$2.02400$	$8.64172$e-$3$	$1.91258$
$2^{-4}\sqrt{2}$	$72$	$1.30694$e-$4$	$2.00075$	$3.36858$e-$5$	$1.99820$	$2.13241$e-$3$	$2.01883$
$2^{-5}\sqrt{2}$	$143$	$3.31576$e-$5$	$1.97878$	$8.44284$e-$6$	$1.99634$	$5.42960$e-$4$	$1.97357$

Table 3. Table 3: Example  4.2 . The corresponding L ∞ ​ ( [ 0 , T ] ; ℋ 2 ​ ( 𝒮 ​ ( t ) ) ) superscript 𝐿 0 𝑇 superscript ℋ 2 𝒮 𝑡 L^{\infty}\big{(}[0,T];\mathcal{H}^{2}(\mathcal{S}(t))\big{)} and L ∞ ​ ( [ 0 , T ] ; ℋ 1 ​ ( 𝒮 ​ ( t ) ) ) superscript 𝐿 0 𝑇 superscript ℋ 1 𝒮 𝑡 L^{\infty}\big{(}[0,T];\mathcal{H}^{1}(\mathcal{S}(t))\big{)} errors and the estimated order of convergence (eoc) by ( 17 ) using h = △ ​ t ℎ △ 𝑡 h=\triangle t , compared with the L ∞ ​ ( [ 0 , T ] ; ℋ 1 ​ ( 𝒮 ​ ( t ) ) ) superscript 𝐿 0 𝑇 superscript ℋ 1 𝒮 𝑡 L^{\infty}\big{(}[0,T];\mathcal{H}^{1}(\mathcal{S}(t))\big{)} results in [ 8 , Table 4] , under the same settings as in Figure  4 .

$h$	$n$	Algorithm 2				cGdG		cGcG
		$L^{\infty }([0,T];$	eoc	$L^{\infty }([0,T];$	eoc	$L^{\infty }([0,T];$	eoc	$L^{\infty }([0,T];$	eoc
		${\mathscr{H}}^2\left(𝒮\left(t\right)\right))$		${\mathscr{H}}^1\left(𝒮\left(t\right)\right))$		${\mathscr{H}}^1\left(𝒮\left(t\right)\right))$		${\mathscr{H}}^1\left(𝒮\left(t\right)\right))$
$3/4$	$36$	$5.5926$e-$3$	$-$	$3.7442$e-$3$	$-$	$-$	$-$	$-$	$-$
$1/2$	$78$	$1.9828$e-$3$	$2.557$	$9.2059$e-$4$	$3.460$	$0.574872$	$-$	$1.12954$	$-$
$1/4$	$312$	$4.3551$e-$4$	$2.187$	$1.8276$e-$4$	$2.333$	$0.303749$	$0.920$	$0.550169$	$1.04$
$1/8$	$1206$	$9.6060$e-$5$	$2.181$	$4.0135$e-$5$	$2.187$	$0.157731$	$0.945$	$0.249556$	$1.14$
$1/16$	$4836$	$2.3407$e-$5$	$2.037$	$9.7715$e-$6$	$2.038$	$0.0803239$	$0.974$	$0.112138$	$1.15$

Equations

$$\begin{array}[]{ll}\partial_{t}u+\bm{v}\bm{\,\cdot\,}\nabla u+(\nabla_{\mathcal{S}}\bm{\,\cdot\,}\bm{v})u-\varepsilon\Delta_{\mathcal{S}}u=f&\mbox{\quad in }\mathcal{S}(t)\times[0,T]\end{array}$$

$$u(\cdot ,0)=u_0\text{mbox}on\mathcal{S}(0),$$

$${\nabla }_{\mathcal{S}}:=(I_d-\mathbf{n}{\mathbf{n}}^T)\nabla$$

$${\Delta }_{\mathcal{S}}:={\nabla }_{\mathcal{S}}\cdot {\nabla }_{\mathcal{S}}$$

$$\frac{\text{d}}{\text{d}t}{\int }_{\mathcal{S}(t)}u\text{d}\mathcal{S}=0.$$

$$\mathcal{S}(t):=\Big{\{}\bm{x}(\bm{\varphi},t)\in\mathbb{R}^{d},\ \text{for}\;t\in[0,T],\ \bm{\varphi}\in\mathcal{B}\subset\mathbb{R}^{d-1}\Big{\}},$$

$$\mathbf{x}(\phi ,t)=[x_1(\phi ,t),x_2(\phi ,t)]=[r(t)\cos \phi ,r(t)\sin \phi ]\in \mathcal{S}(t)\subset {\mathbb{R}}^2.$$

$$\mathbf{v}(\mathbf{\phi },t)=\frac{\partial }{\partial t}\mathbf{x}(\mathbf{\phi },t)={\left[\begin{array}{c}{\partial }_t{x}_i(\mathbf{\phi },t)\end{array}\right]}_{i=1}^d.$$

$${\mathbf{x}}^m=\mathbf{x}(\mathbf{\phi },t^m)\subset {\mathcal{S}}^m,\mathbf{\phi }\in \mathcal{B}.$$

$$u^{m}:=u(\bm{x}^{m},t^{m})=u\Big{(}\bm{x}(\bm{\varphi},t^{m}),t^{m}\Big{)}\mbox{\quad and \quad}f^{m}:=f(t_{m}).$$

$$\frac{\partial }{\partial t}u+\mathbf{v}\cdot \nabla u$$

$$\frac{u^m-u^{m-1}}{△t}+(1-\theta ){\mathcal{A}}^m{u}^m+\theta {\mathcal{A}}^{m-1}{u}^{m-1}=(1-\theta )f^m+\theta f^{m-1},$$

$${\mathcal{A}}^m:=\mathbf{v}\cdot {\nabla }_{{\mathcal{S}}^m}-\epsilon {\Delta }_{{\mathcal{S}}^m}$$

$${\mathcal{L}}^{m}u^{m}:=\Big{(}1+{(1-\theta)}\triangle t{\mathcal{A}}^{m}\Big{)}u^{m}=g^{m},$$

$$g^{m}:=\Big{(}1-{\theta}\triangle t{\mathcal{A}}^{m-1}\Big{)}u^{m-1}+\triangle t\Big{(}(1-\theta)f^{m}+\theta f^{m-1}\Big{)}.$$

$$\Psi (\cdot ,\cdot ):{\mathcal{S}}^m\times {\mathcal{S}}^m\to \mathbb{R}$$

$$u^m=\Psi (\cdot ,Z^m){\mathbf{\lambda }}_Z^m:=\sum _{z_i\in Z^m}{\lambda }_i^m\Psi (\cdot ,z_i),$$

$$[{\mathcal{L}}^m\Psi (X^m,Z^m)]{\mathbf{\lambda }}_Z^m=g_{|X}^m.$$

$${\mathbf{x}}^{m+1}\approx {\mathbf{x}}^m+\mathbf{v}△t,$$

$$\mathbf{n}=\frac{\nabla \Psi (\cdot ,Z)\Psi (Z,Z{)}^{-1}{1}_{n_Z}}{\Vert \nabla \Psi (\cdot ,Z)\Psi (Z,Z{)}^{-1}{1}_{n_Z}\Vert }:{\mathcal{S}}^m\to {\mathbb{R}}^d,Z\subset {\mathcal{S}}^m$$

$$u^m=\Psi (\cdot ,Z)\Psi (Z,Z{)}^{-1}{u}_{|Z}^m.$$

$$({\tilde{\nabla }}_{\mathcal{S}}{u}^m{)}_{|Z}={\tilde{\nabla }}_{\mathcal{S}}\Psi (Z,Z)\Psi (Z,Z{)}^{-1}{u}_{|Z}^m.$$

$${\tilde{\nabla }}_{\mathcal{S}}{u}^m=\Psi (\cdot ,Z)\Psi (Z,Z{)}^{-1}({\tilde{\nabla }}_{\mathcal{S}}{u}^m{)}_{|Z},$$

$$\widetilde{\Delta}_{\mathcal{S}}u^{m}=\big{(}\widetilde{\nabla}_{\mathcal{S}}\bm{\,\cdot\,}\Psi\big{)}(X,Z)\Psi(Z,Z)^{-1}(\widetilde{\nabla}_{\mathcal{S}}u^{m})_{|Z}.$$

$$P(\mathbf{x})=I_d-\mathbf{n}(\mathbf{x})\mathbf{n}(\mathbf{x}{)}^T=:[p_1(\mathbf{x}),\dots ,p_d(\mathbf{x})]\in {\mathbb{R}}^{d\times d}$$

$$G_k(\mathbf{x},Z):=p_k(\mathbf{x}{)}^T[\nabla \Psi (\mathbf{x},Z)][\Psi (Z,Z){]}^{-1}.$$

$$[{\tilde{\nabla }}_{\mathcal{S}}{u}^m{]}_k$$

$${\tilde{\Delta }}_{\mathcal{S}}{u}^m$$

$$\underset{[0,T]}{\sup }||u^{\ast }-u|{|}_{L^2(\mathcal{S}(t))}.$$

$$||\nabla_{\mathcal{S}}u^{*}-\nabla_{\mathcal{S}}u||_{L^{\infty}\big{(}[0,T];L^{2}(\mathcal{S}(t))\big{)}},$$

Full text

1 Introduction

Heat transport on moving surfaces are modelled by convection-diffusion partial differential equations (PDEs) that frequently appear in physical and biological fields. Applications include dealloying by surface dissolution [1], pattern formation on evolving biological surfaces [2], modelling geometric biomembranes [3], and cell motility and chemotaxis [4]. The formal problem statement is as follows. Let $\mathcal{S}(t)\subset\mathbb{R}^{d}$ be a continuously evolving surface with codimension one; for each $t\in[0,T]$, we also assume that the surface is smooth, closed, connected, and complete. Furthermore, let $\bm{{n}}=\bm{{n}}(\bm{\,\cdot\,},t)$ and $\bm{\tau}=\bm{\tau}(\bm{\,\cdot\,},t)$ represent the unit normal and the tangent vector of $\mathcal{S}(t)$.

We consider convective-diffusion equations defined on $\mathcal{S}(t)$ in the form of

[TABLE]

subject to homogeneous Dirichlet initial conditions

[TABLE]

where $\varepsilon\geq 0$ is the diffusive coefficient, and $\bm{v}:=\bm{v}_{n}+\bm{v}_{\tau}$ is the velocity of surface motion whose normal and tangent components are denoted by $\bm{v}_{n}$ and $\bm{v}_{\tau}$ respectively. The surface differential operators in (1) are defined as

[TABLE]

and

[TABLE]

where $I_{d}$ is the $d$-dimensional identity matrix and $\nabla$ is the standard gradient operator in Euclidian spaces.

The problem we considered includes several physical models in the literature. In the absence of diffusive flux and source, i.e., $\varepsilon=0$ and $f=0$, the viscosity solution $u$ to (1) as $\varepsilon\rightarrow 0^{+}$ satisfies a mass conservation law [5], namely

[TABLE]

Readers can find more details in [2, 6, 7, 8].

To solve the target PDE (1) on the evolving surface, a surface finite element method (SFEM) were proposed with the use of weak and variational formulas based on triangular meshes intrinsically defined on surfaces [6, 9]. Embedding techniques for static problems can also be extended to deal with surface PDEs posed in higher-dimensional [10, 11, 12, 13, 14]. Some domain-type time-dependent FEMs have been implemented in narrow bands containing moving surfaces [7, 8, 15, 16]. Implicit time-stepping schemes have also been used for temporal discretization in the FEM literature mentioned above, except in [8] where weak forms were utilized in both space and time. In this paper, we apply kernel-based meshless collocation techniques to solve (1). The proposed method is intrinsic, which means that no embedding is required. One variant of our method naturally works with point clouds, which makes it capable to handle solution-driven surface evolutions.

The rest of this paper is organized as follows. Firstly, we deal with the case where the evolutions of surfaces are prescribed by some parameterizations that are known a prior. Section 2 contains a novel algorithm for this case. Next, we focus on cases where the surfaces are defined by some point clouds. In Section 3, we propose a second algorithm that does not require analytic information about normals and velocities of surface motions. In Section 4, we first test accuracy and convergence of our first algorithm. Next, we test the proposed algorithms against mass balance in the case of (2). Lastly, as a pilot study, we examine the robustness of our method in dealing with merging surfaces, which involves with high curvatures. We conclude the results and observations in Section 5.

2 A collocation method for PDEs posed on parameterized surfaces with prescribed evolution

To begin, we consider the case that the surface $\mathcal{S}(t)$ and its motion are both given by some time-dependent parameterization. Assume the evolving surface is defined as

[TABLE]

where $\bm{\varphi}\in\mathcal{B}$ are parameters in certain parameter space that define surface points \bm{x}(\bm{\varphi},t)=\big{[}x_{i}(\bm{\varphi},t)\big{]}_{i=1}^{d} for any time $t\in[0,T]$. The following example aims to clarify our notations.

Example:

Consider an evolving circle with changing radius. Using polar coordinate, we can pick $\varphi=\bm{\varphi}\in\mathcal{B}=[0,2\pi]\subset\mathbb{R}$ and for some function $r(t)>0$ to define surface points on $\mathcal{S}(t)$ at any time $t\geq 0$ by

[TABLE]

Obviously, the velocity for this surface motion is given by the partial derivative of $\bm{x}$ with respect to $t$.      $\Box$

In our notations, we have velocity vector in (1) given by

[TABLE]

With all required analytic information about the evolving surface ready, we can start discretizing (1) first in time and then in space.

Firstly, we semi-discretize (1) based on some partition $\{t_{m}\}_{m=0}^{N}$ of the interval $[0,T]$. For simplicity, we assume the partition is equispaced with time stepping size $\triangle t$ and we will employ the standard $\theta$-scheme. Note that the problem (1) in hand is given in the form of ${\partial_{t}u}=F(u)$, which is not yet ready for applying the $\theta$-scheme. Instead, we implicitly discretize the material derivatives [7, 15] of the solution to (1) as follows.

We abbreviate $\mathcal{S}(t^{m})=:\mathcal{S}^{m}$ the snapshot of our evolving surfaces at time $t_{m}$ for $m=1,\ldots,T/\triangle t$, on which we have surface points

[TABLE]

Let us also denote the solution $u$ and right-hand function $f$ at $t_{m}$ as

[TABLE]

We approximate the material derivative of $u$ and obtain

[TABLE]

Thus, the time discretization of (1) by $\theta$-method is given as

[TABLE]

where the surface elliptic operator ${\mathcal{A}}^{m}$ is defined as

[TABLE]

on the surface $\mathcal{S}^{m}$ and the velocity $\bm{v}=\bm{v}(\bm{\,\cdot\,},t^{m})$ is to be evaluated exactly. We can further simplify (5) to isolate the unknown solution $u^{m}$ and get a time-independent second-order surface PDE:

[TABLE]

where the right-hand function $g^{m}$, which depends on the known functions $u^{m-1}$ and $f$, is known and given by

[TABLE]

Aiming towards fully discretized problems, we employ our previously proposed overdetermined Kansa method [17] to solve surface PDEs. We set up an overdetermined kernel-based collocation method [18] to solve surface PDE ${\mathcal{L}}^{m}u^{m}=g^{m}$ on surfaces $\mathcal{S}^{m}$ and $\mathcal{S}^{m-1}$ (due to ${\mathcal{L}}^{m}$ and $g^{m}$) for solution $u^{m}$.

Based on the theories in [19], we can take some commonly used symmetric positive definite (SPD) kernel and restrict it to any one of the surfaces $\mathcal{S}^{m}$ to obtain an SPD surface kernel

[TABLE]

with a certain smoothness order $\mu>(d-1)/2$. Let $\varphi_{Z}\subset\mathcal{B}$ be a set of $n_{Z}$ points in the parameter space; similarly, let $\varphi_{X}\subset\mathcal{B}$ be another set of $n_{X}>n_{Z}$ points. For each $m$, let $Z^{m}=\bm{x}(\varphi_{Z},t^{m})$ and $X^{m}=\bm{x}(\varphi_{X},t^{m})$ be the set of trial centers and collocation points. We shall use denser collocation points because of available convergence theories for collocation methods in domains [18] and on surfaces [12]. In both works, they require $X$ to be sufficiently denser with respect to $Z$ resulting in overdetermined formulas, in order to establish stability estimates222In this context, stability means that errors can be bounded in continuous norms with discrete residuals..

Trial functions at time $t^{m}$ are expressed by a linear combination

[TABLE]

where $\bm{\lambda}_{Z}^{m}:=\{\lambda_{i}^{m}\}_{i=1}^{n_{Z}}$ is the unknown coefficient vector to be determined. As in the original Kansa method, we collocate (7) and (8) at set $X^{m}$ to yield an $n_{X}\times n_{Z}$ matrix system

[TABLE]

In the Euclidean counterpart, overdetermined formulations guarantee stability as mentioned before , which ultimately leads to error estimates [12, 18].

Given the parametric equation (3), the surface differential operator in (6) can be rewritten without any implicit dependency on the surface. Thus, $\mathcal{L}^{m}\Psi:\mathcal{S}^{m}\times\mathcal{S}^{m}\to\mathbb{R}$ can also be analytically evaluated for each $m=1,\ldots,T/\triangle t$. The right hand vector $g^{m}_{|X}=g^{m}(X^{m})$ is the nodal value of (8) evaluated at $X^{m}$. The overdetermined system (10) should be solved in the least-squares sense. We summarize with Algorithm 2.

Algorithm 1 : For solving PDE (1) posed on a parametrized surface (3)

${\bf Initialization:}$

$\qquad set~{}t=t^{0},\ select~{}0\leq\theta\leq 1;$

$\qquad define~{}\varphi_{Z},\varphi_{X}\subset\mathcal{B},\ compute~{}Z^{0},X^{0}\subset\mathcal{S}^{0};$

$\qquad interpolate\ u_{0}\ at\ Z^{0}\ to\ obtain\ \bm{\lambda}_{Z}^{0};$

$\qquad symbolically~{}derive~{}analytic~{}formulas~{}of~{}{\mathcal{A}}^{m}~{}in~{}{\rm(\ref{eq_lst})};$

$\qquad m=1;$

${\bf While}~{}(t\leq T)$

$\qquad compute~{}\{X^{m},Z^{m}\}\subset\mathcal{S}^{m};$

$\qquad assemble~{}and\ solve~{}{\rm(\ref{kansa})}\ in\ the\ least$-$squares\ sense;$

$\qquad m=m+1;$

$\qquad t=t^{0}+(m+1)\triangle t;$

${\bf End}$

3 Another collocation method for point clouds

In the case of surfaces defined by point clouds, some approximated method is required to discretize space without analytic formulas for the normal vectors. In this section, we assume that the initial surface $\mathcal{S}^{0}$ is given by a set of points in $\mathbb{R}^{d}$. We also assume the surface evolution velocity $\bm{v}=\bm{v}(\bm{x},t)$ is independent to the PDE solution. As surface points $\bm{x}(\bm{\varphi},t)$ evolving from $\mathcal{S}^{m}$ to $\mathcal{S}^{m+1}$ ($m\geq 0$) with velocity $\bm{v}$, we can approximate new points ${\bm{x}}^{m+1}$ on $\mathcal{S}^{m+1}$ by

[TABLE]

as shown in Figure 1. This process allows us to track all surfaces and place data points $X$ and $Z$ via local interpolation on them as in Algorithm 2; see Figure 2 for a schematic demonstration. The problem here is that ${\mathcal{A}}^{m}$ in (6) can no longer be obtained analytically. Our solution is to take the RBF pseudospectral approach [20] to approximate the surface Laplacian. The ultimate goal is to obtain an approximated version of overdetermined linear system as in (10).

Recall that $\nabla_{\mathcal{S}}:=(I_{d}-\bm{{n}}\bm{{n}}^{T})\nabla$. We can approximate $\bm{{n}}=\bm{{n}}(\bm{x})$ for $\bm{x}\in\mathcal{S}^{m}$ from the point cloud $\{\bm{x}^{m}_{k}\}_{k=1}^{K}$ specifying $\mathcal{S}^{m}$ by

[TABLE]

where $\mathbf{1}_{n_{Z}}$ is the all-ones vector of size $n_{Z}$, see [21, 22]. When implementing (12), we use subsets of $Z$ surrounding the points of interest for the sake of efficiency. Once we have all normal vectors at $X$ approximated, the computation of entries in the approximated gradient matrix $\widetilde{\nabla}_{\mathcal{S}}(X,Z)$ is straightforward. The approximation error here is due to error in $\bm{{n}}$ and the exactness of data point locations on $\mathcal{S}^{m}$.

The same technique is undesirable for the approximation of surface Laplacian $\Delta_{\mathcal{S}}:=\nabla_{\mathcal{S}}\bm{\,\cdot\,}\nabla_{\mathcal{S}}$, as the surface divergence operator will act upon $\bm{{n}}$ and we simply want to avoid approximating second derivatives from point cloud. In the RBF pseudospectral approach, we rewrite (9) in terms of the unknown nodal data $u_{|Z}^{{m}}:={u^{m}}(Z)$ of $u^{m}$ at trial centers $Z\subset\mathcal{S}^{m}$. This yields

[TABLE]

We can apply the approximation and evaluate the surface gradient at $Z$ while keeping $u_{|Z}^{m}$ unknown to yield the following $n_{Z}\times d$ matrix

[TABLE]

At this point, we have approximated nodal values of each component of $\nabla_{\mathcal{S}}$ at $Z$, which still depends on the unknown $u_{|Z}^{m}$. Without explicitly writing out the index of these components for simplicity, we can express them similar to (13) by its interpolant

[TABLE]

to which we can apply an approximate surface divergence at $X\subset\mathcal{S}^{m}$ to yield

[TABLE]

This apparently complicated procedure did come with some neat simplifications [23]. Let $\bm{x}\in\mathcal{S}^{m}$ be a surface point on which the (approximated) normal vector is denoted by $\bm{{n}}(\bm{x})$. Also let

[TABLE]

be the (approximated) projection operator to $\bm{x}\in\mathcal{S}^{m}$ with columns $p_{k}$, $k=1,\ldots,d$. For each $k$, we define row-vector functions

[TABLE]

Then, we have

[TABLE]

or, alternatively, we can use $(\widetilde{\nabla}_{\mathcal{S}}u^{m})_{|X}$ as data to approximate surface Laplacian by

[TABLE]

These equations, along with (13), are sufficient for obtaining an approximated version of the overdetermined linear system in (10) for solving (7) in terms of unknown $u_{|Z}$. Using the right most inverse of the interpolation matrix of $\Psi$ at $Z$ and the relationship $\bm{\lambda}_{Z}^{m}=[\Psi(Z,Z)]^{-1}u_{|Z}^{m}$, we can recast the system with unknown $\bm{\lambda}_{Z}^{m}$. Below is a pseudocode that summarizes the algorithm above: Algorithm 3, specifically denoted as Algorithm 3a by (15) and Algorithm 3b by (16).

Algorithm 2 a/b: For solving PDE (1) posed on point cloud

${\bf Initialization:}$

$\qquad set~{}t=t^{0},\ select~{}0\leq\theta\leq 1;$

$\qquad define~{}~{}Z^{0},X^{0}\subset\mathcal{S}^{0}\ based\ on\ some\ point\ cloud$

$\qquad compute~{}\bm{{n}}(X^{0})\ using\ {\rm(\ref{find n})};$

$\qquad interpolate\ u_{0}\ at\ Z^{0}\ to\ obtain\ \bm{\lambda}_{Z}^{0};$

$\qquad m=1;$

${\bf While}~{}(t\leq T)$

$\qquad compute~{}\bm{v}~{}and~{}update~{}\{X^{m},Z^{m}\}\subset\mathcal{S}^{m};$

$\qquad approximate~{}entries~{}in~{}{\rm(\ref{kansa})}~{}using\ {\rm(\ref{eq_intS2})}-{\rm(\ref{LL1})},~{}and~{}{\rm(\ref{LL2a})}/{\rm(\ref{LL2b})};$

$\qquad solve~{}{\rm(\ref{kansa})}\ in\ the\ least$-$squares\ sense;$

$\qquad m=m+1;$

$\qquad t=t^{0}+(m+1)\triangle t;$

${\bf End}$

4 Numerical examples

We provide five examples to numerically study the proposed algorithms for solving convective-diffusion equations (1) on evolving surfaces. In all examples, we use the restricted Whittle-Matérn-Sobolev kernels and $\theta=0.5$ in (7) for the Crank-Nicolson method (CN).

To quantify the accuracy and convergence of our numerical approximations, we employ two error measures as follows. Let $u^{*}$ denote the exact solution. The L^{\infty}\big{(}[0,T];L^{2}(\mathcal{S}(t))\big{)} error is defined by

[TABLE]

Similarly, the L^{\infty}\big{(}[0,T];\mathcal{H}^{1}(\mathcal{S}(t))\big{)} semi-norm error is given by

[TABLE]

and lastly, the L^{\infty}\big{(}[0,T];\mathcal{H}^{2}(\mathcal{S}(t))\big{)} semi-norm error is defined by

[TABLE]

At any given time $t^{m}\in[0,T]$, $m\geq 1$, by using one of the above errors $e^{m}$ and $e^{m-1}$ at two consecutive time steps corresponding to fill distances $h_{X^{m-1}}$ and $h_{X^{m}}$ for collocation point sets $X^{m-1}\subset\mathcal{S}^{m-1}$ and $X^{m}\subset\mathcal{S}^{m}$ respectively, we estimate the order of convergence by

[TABLE]

In cases of smooth initial conditions (ICs) on surfaces without (extremely) high curvatures, we can alternatively employ overdetermined expressions with $n_{X}>n_{Z}$ or exactly determined formulas with the same point sets ($X=Z$) for both Algorithms 2 and 3. Having that said, oversampling in Algorithm 2 can still result in better accuracy; see Example 4.1. The situations are quite different when imposing discontinuous ICs on surfaces with high curvatures, in which oversampling is essential. In our concluding simulations in Examples 4.4 and 4.5, we will demonstrate the robustness of our proposed methods.

Example 4.1

Evolving curve in $\mathbb{R}^{2}$.

In the first example, we compare Algorithm 2 with two time-dependent FEMs [8, 15] for solving (1) on two moving surfaces. We use the same diffusivity $\varepsilon=1$ and the source terms $f(\bm{x},t)$ by the symbolic mean by providing the same exact solutions $u^{*}(\bm{x},t)$ in the FEM papers. The smoothness order is set to be $\mu=4$.

In [15], weak formulations of (1) were discretized in some narrow bands containing the evolving surface, and an implicit temporal scheme was used to solve (1), which yields an unfitted finite element method. We aim to examine the same experiment as in [15, Section 6.3] with the use of Algorithm 2, imposed on the following evolving curve,

[TABLE]

First, we test the Algorithm 2 with oversampling. In Figure 3, the numerical solutions and their error functions are plotted by color on curves at some time points from $t=0$ to $0.5$, when using $h=2^{-5}\sqrt{2}$.

Table 1 shows the corresponding L^{\infty}\big{(}[0,T];L^{2}(\mathcal{S}(t))\big{)} errors under several different oversampling settings of $n_{X}>n_{Z}$ with the corresponding time step $\triangle t=h_{Z}/4$. With the increasing numbers of both centers and collocation points, the L^{\infty}\big{(}[0,T];L^{2}(\mathcal{S}(t))\big{)} errors decrease. We then compare results in Table 1 with those listed in Table 2 obtained by the corresponding exactly determined settings $X=Z$. We find that using $h_{X}<h_{Z}$ can improve accuracy slightly. Under such smoothness order for the kernel, oversampling or not both yield the expected second-order L^{\infty}\big{(}[0,T];L^{2}(\mathcal{S}(t))\big{)} convergence. Notably, even for settings with very few data points, say $n_{Z}=5$, we still observe accuracy of $10^{-2}$ in L^{\infty}\big{(}[0,T];L^{2}(\mathcal{S}(t))\big{)} norms that outperforms the FEM mehtod.

For comparison with the results of the unfitted FEM as found in [15], Table 2 also lists the L^{\infty}\big{(}[0,T];L^{2}(\mathcal{S}(t))\big{)} errors with $\triangle t=2h^{2}$, as well as the estimates in L^{\infty}\big{(}[0,T];\mathcal{H}^{2}(\mathcal{S}(t))\big{)} error of our Algorithm 2 with $X=Z$. It also presents convergence orders as computed by (17). Algorithm 2 achieves higher-order accuracy in L^{\infty}\big{(}[0,T];L^{2}(\mathcal{S}(t))\big{)} and L^{\infty}\big{(}[0,T];\mathcal{H}^{2}(\mathcal{S}(t))\big{)} norms, using a relatively large $\triangle t=h/4$. The magnitudes of our corresponding L^{\infty}\big{(}[0,T];L^{2}(\mathcal{S}(t))\big{)} errors are lower by a factor of about 100 relative to the unfitted FEM with many more iterations ($\triangle t=2h^{2}$).

Example 4.2

Evolving surface in $\mathbb{R}^{3}$.

This is a 3D version of the previous example, but we focus on the exactly determined case with $X^{m}=Z^{m}$ at all time steps. We use kernels with the smoothness order $4$. In [8], two Eulerian FEMs (the cGdG and cGcG methods) employed Galerkin methods both in time and space with narrow bands to solve (1). The test problem in [8, Section 11.3] is posed on a moving surface given by

[TABLE]

We solve the same problem by Algorithm 2 with $\triangle t=1/8$. Based on the initial points with $h:=h_{Z^{0}}=\triangle t$, we move each point by point-to-point motion based on the analytical formula in polar coordinate, that is, the fill distance $h_{Z^{m}}$ slightly varies in time periodically indeed. In Figure 4, we present a few snapshots of the evolved surfaces and numerical solutions by colormap.

To be compared with [8, Figures 13-14], Figure 5 illustrates the convergence profile of Algorithm 2 in L^{\infty}\big{(}[0,T];L^{2}(\mathcal{S}(t))\big{)} errors. First, Figure 5(a) is obtained by three fixed $h$-values as $\triangle t$ varies. It can be seen that errors for $h=1/4$ are almost as accurate as those for $h=1/8$, under the same $\triangle t$.

Our comparison with the cGdG and cGcG approaches in [8] is as follows. With respect to time, our convergence rates and those of the two FEMs are all around order 2 due to CN. In terms of accuracy, we obtain errors less than $10^{-5}$ by using $h=1/4$, while their best FEM result shown is approximately $10^{-4}$ using a much smaller $h=1/64$. For spatial convergence, we fix a several $\triangle t$ and let $h$ is allowed to vary to obtain Figure 5(b). Generally speaking, our numerical solutions are more accurate based on [8, Figures 13-14], which used smaller $h$-values for each $\triangle t$. Our method shows spatial convergence of $6.4$, whereas that in [8] is around $2$.

Table 3 further lists the L^{\infty}\big{(}[0,T];\mathcal{H}^{2}(\mathcal{S}(t))\big{)} and L^{\infty}\big{(}[0,T];\mathcal{H}^{1}(\mathcal{S}(t))\big{)} estimates obtained by Algorithm 2, and the L^{\infty}\big{(}[0,T];\mathcal{H}^{1}(\mathcal{S}(t))\big{)} results of two Eulerian FEMs in [8] using various $h=\triangle t$. Note that Algorithm 2 remains second order while the FEM methods drop to first order.

Example 4.3

Mass conservation on an expanding sphere.

Aiming to verify the mass conservation law in (2), we solve PDE (1) without the source term, i.e., $f=0$. We choose a small diffusion coefficient $\varepsilon=10^{-3}$ to approximate the viscosity solutions to (1). The tested surface is defined by

[TABLE]

We set the initial solution to be $u^{*}(\bm{x},0)=0.5+x_{1}x_{2}x_{3}$.

We use $X=Z$ in both Algorithms 2 and 3 (Algorithm 3a and b are practically identical here) to solve this problem with the kernels of order $4$ and $h=\triangle t=0.02$. Solutions of two algorithm agree up to some order of $10^{-6}$; their difference on $\mathcal{S}(t)$ is shown in Figure 6 for a particular parameter setup.

Now we can see Figure 7 for the masses of the solutions $\int_{\mathcal{S}(t)}u\ \text{d}\mathcal{S}$ over time, which are approximated by $\sum_{i}^{n_{X}}u(\bm{x}_{i})\triangle\mathcal{S}_{i}$ at $X$ based on some triangle mesh $\mathcal{T}$. It is easy to see that $u(\bm{x}_{i})\triangle\mathcal{S}_{i}=u(\bm{x}_{i})\mbox{Area}(\mathcal{T}_{i})/3$ where $\mathcal{T}_{i}$ is the triangle containing $\bm{x}_{i}$. For Figure 7(a), it is clear that both algorithms conserve mass as desired. When we subtract the exact mass computed using $u^{*}(\bm{x},0)$ to yield errors in Figure 7(b), we can see the different mass-conservation behaviour of the proposed algorithms. The analytically exact Algorithm 2 shows gradually increasing error in time, which is typical for accumulative error. On the other hand, the nearly-constant error curve of Algorithm 3 in Figure 7(b) suggests that the error in numerical mass comes from the approximation steps in formulation.

Example 4.4

Continuous initial conditions on point cloud with high curvatures.

The last two examples aim to prepare for problems on merging surfaces. In both cases, we consider the mean-curvature motion of a kissing-spheres, as in Figure 2, which has high initial curvature at the contact point. Here, we use the FDM proposed in [24] to obtain point clouds in $\mathbb{R}^{3}$ at each discrete time. Without analytic formula for $\mathcal{S}(t)$, we can only use Algorithm 3.

In order to test the effect of geometry, we first impose a continuous IC

[TABLE]

on the surface. The arrows in Figure 8 show mean-curvature velocities and hence the directions of motion. The diffusion coefficient and time step size are set as $\varepsilon=1$ and $\triangle t=0.001$ respectively. First, we solve (1) with $f=0$ by using the exactly determined formulas with ($n_{X}=n_{Z}=1196$). Figure 9 presents two numerical results at $t=0.05$ and $0.2$. From Figure 9(a), some oscillations near the region with high curvature can be observed. Once diffusion starts, the oscillatory solution smooth out to the whole surface, see Figure 9(b). We postpone the comparison with oversampling to the next example.

Example 4.5

A discontinuous initial condition imposed on point clouds evolving by mean-curvature motion.

We continue the study on the initial high-curvature surface in Example 4.4 but use a discontinuous IC defined by

[TABLE]

as shown in Figure 10. In the following, we will examine the robustness of Algorithms 3a and 3b.

The main results were shown in Figure 11. Because of the discontinuity, it makes sense to include some common regularization techniques into our comparison. We apply a regularized RBF interpolant [25] to interpolate the IC; then, we update it in time with $X=Z$. The employed regularization parameter $p=0.2488$ is chosen by L-curves in $L^{2}(\mathcal{S}(t))$ errors [26]. From Figure 11(a)-(b), this common-sense setup yields undesirable solutions and we clear see the oscillatory solutions propagating towards the whole surface over time.

Now we focus on the results from Algorithms 3a and b shown in Figures 11(c)-(d) and (e)-(h) respectively. As Algorithm 3b is a computationally more expensive method, we use it with fewer $Z$ nodes to compare fairly. To sum up, Algorithm 3a still show the spurious oscillations initially that will be corrected and the solutions get smoother with increasing time. Solutions of Algorithm 3b, in contrast, are physically correct without oscillations.

Before example ends, we want to make sure that Algorithm 3 really is honestly more suitable than regularization. We now apply regularization at each time, instead of to the initial time only, hoping to avoid oscillatory solutions. Various regularization parameters $p$ were used, all of which yield non-oscillatory smooth wrong solutions. We make this conclusion based on mass conservation in Figure 12.

5 Conclusions

We propose some numerical algorithms based on intrinsic kernel-based meshless techniques for solving parabolic PDEs on evolving surfaces. The proposed algorithms can be implemented either by an analytic projection or a pseudospectral approximation for differential operators on surfaces. Extension into embedding spaces is unnecessary.

The analytic approach (Algorithm 2) can achieve high-order accuracy and convergence in space when normal vectors of the surface can be differentiated analytically. In case of otherwise, the approximated meshless method (Algorithm 3) becomes handy and it can handle PDEs defined on point clouds. Another issue of study is the oversampling strategy, viz., using denser collocation points than trial centers, to yield overdetermined matrix systems to be solved by least-squares. For smooth problems (including initial conditions, surfaces, and solutions), the need of oversampling is not obvious. We demonstrate that, when the initial condition is discontinuous and when the (domain) surface has high curvature, oversampling is essential for getting physically correct solutions.