Approximation and application of minimizing movements for surface PDE
Elliott Ginder Karel Svadlenka Takuma Muramatsu


GraphEx: A Graph-based Extraction Method for Advertiser Keyphrase Recommendation
    Kamesh Madduri
    September 9, 2024
================================================================================

empty
page-1


top=2cm
1.2


§ INTRODUCTION

We extend the applicability of minimizing movements (MM) to approximating solutions of surface partial differential equations (SPDE) and apply the technique to approximate mean curvature flow (MCF)<cit.> and hyperbolic MCF (HMCF)<cit.> on surfaces.
The MBO algorithm <cit.> and the HMBO algorithm <cit.>, both based on the level set method, are well-known approximation methods for such curvature flows.
We recall that the MBO algorithm is an approximation method for mean curvature flow and is based on solving the heat equation.
On the other hand, the HMBO algorithm is an approximation method for the hyperbolic mean curvature flow and involves solving the wave equation.
It has been shown that the MBO algorithm and HMBO algorithm can approximate curvature flow under area preservation constraints, as well as in the multiphase setting <cit.>.
Minimizing movements <cit.> are used to realize the area conservation condition, and the signed distance vector field <cit.> is used for calculations involving multiphase regions.

On the other hand, for curvature flow on curved surfaces, approximation methods for the mean curvature flow were presented in <cit.>.
The authors show that the Closest point method (CPM) <cit.> can be used to approximate mean curvature flow on surfaces.
This is done by extending the values of functions defined on a surface to the ambient space of the surface. In turn, the CPM enables the approximation of surface gradients and other differential quantities by making use of the surround Euclidean space. However, no previous studies have treated approximation of solutions to  curvature flow under area preservation involving interfaces in the multiphase and  curved surface setting.

This study develops approximation methods for mean curvature flow and hyperbolic mean curvature flow under multiphase area preservation conditions for interfaces moving on curved surfaces.

Our approach is to extend MM to the case of surface PDE and to use their framework to apply the MBO and HMBO algorithms. Similar to <cit.>, our generalizations make use the CPM which we combine with the surface-constrained signed distance vector field <cit.>.

The outline of this paper is as follows. In section <ref>, we describe the research background and our objectives. Our objectives are based on conventional approximation methods such as the MBO algorithm and HMBO algorithm. To achieve our goals, approximation methods for constrained partial differential equations on surfaces are required. We therefore demonstrate that partial differential equations on surfaces can be computed using the Closest Point Method and introduce surface-type minimizing movements to handle partial differential equations with constraints. In section <ref>, we discuss computational techniques related to partial differential equations on surfaces. In particular, we create an approximation method by combining the Closest Point Method with minimizing movements and perform numerical error analyses for heat and wave equations defined on surfaces. In section <ref>, we discuss our approximation method for mean curvature flow on surfaces and hyperbolic mean curvature flow on surfaces. This includes an explanation of the signed distance vector field, which is required to handle multiphase domains, and the area preservation condition that is achieved through the use of minimizing movements. In section <ref>, we discuss mean curvature flow and hyperbolic mean curvature flow on surfaces and describe the method for enforcing area preservation conditions in multiphase environments. We summarize the contents of this paper and discuss future challenges in section <ref>.



§ BACKGROUND

In this section, we will briefly explain our goals, and the mathematical frameworks used in our research. In section <ref>, we will touch upon our research objectives. Then, in section <ref>, we introduce the Closest Point Method (CPM), and in section <ref> we will introduce the minimizing movements (MM). Again remark that, by combining these methods, we obtain an approximation method for partial differential equations with constraints on curved surfaces.


 §.§ Objectives

As stated in Section <ref>, the goal of this research is to create an approximation method for interfacial motion on curved surfaces. Here, we will explain two representative examples of interfacial motions on curved surfaces: mean curvature flow <cit.> and hyperbolic mean curvature flow <cit.>. The approximation method for mean curvature flow that we emply is known as the MBO algorithm, which alternates between solving the heat equation and constructing level set functions <cit.>. The numerical solution method for hyperbolic mean curvature flow is known as the HMBO algorithm, which alternates between solving the wave equation and constructing level set functions <cit.>. Examples of more complex interfacial motions involving area preservation and extendtions to multiphase regions are illustrated. Here, multiphase regions refer to regions where the domain is divided into three or more regions by the interface. We use the fact that the area preserving condition can be realized by imposing a constraint on the partial differential equations used in the MBO and HMBO algorithms <cit.>. We again remark that the case of multiphase regions is possible to treat by using the signed distance vector field <cit.>.

Based on the above considerations, the purpose of this research is as follows:


  * establish an approximation method for partial differential equations with constraints on curved surfaces

  * extend the MBO and HMBO algorithma to curved surfaces and realize numerical compute mean curvature flow and hyperbolic mean curvature flow on surfaces

  * generalize our famework to treat the above problems in the multiphase setting.



 §.§ Closest Point Method

When numerically solving partial differential equations on surfaces, approximation of surface gradients (SG) on the surface are necessary. 
For a smooth surface S embedded in n-dimensional Euclidean space, the SG of a function u on the surface S is given by

    ∇_S u=∇ u-n(n·∇ u),

where n is the unit normal vector of the surface S, and ∇ is the usual gradient in the Euclidean space.

In the CPM, an approximation of the SG on the surface is obtained by smoothly extending the values of the function defined on the surface in the direction of the surface normal vector <cit.>. 
This is enabled by constructing a function that provides the closest point on the surface S to any external point of the ambient space. 
For any point x in n-dimensional space, the function C_S that gives the closest point on the surface S is defined as follows:

    C_S(x)=_y∈ S|| x-y||

For example, if S is the unit circle (n=2), then C_S is given by:

    C_S(x,y)=(x/√(x^2+y^2),y/√(x^2+y^2)),   (x,y)∈ℝ^2

If S is the unit sphere (n=3), then C_S is given by:

    C_S(x,y,z)=(x/√(x^2+y^2+z^2),y/√(x^2+y^2+z^2),z/√(x^2+y^2+z^2)),   (x,y,z)∈ℝ^3

Note that, in the case of the unit circle or sphere, the closest point to the origin is not uniquely determined. 
In this case, since the distance between the origin and the surface is constant, CPM assigns an arbitrary point on the surface as the closest point to the origin. 
Figure <ref> shows the relationship between a point x in three-dimensional space and its closest point p on the surface S. 
The closest point p is given by the closest point function C_S as p=C_S(x).

Note that, if a surface S is represented by a parameterization, it is easy to find the closest point. In such a case, we just need to solve the optimization problem (<ref>). On the other hand, if the surface S is not parameterized but is represented by a point cloud or a triangulated surface, then a bit more ingenuity is required. For example, when the surface S is discretized by a point cloud, the closest point on the surface S may involve constructing implicit surfaces defined by distance functions, or require one to accept a certain level of loss of uniquess. 
The following theorems forms the basis for approximations using the CPM <cit.>.


    Let S ⊂ℝ^3 be a smooth surface. 
    Let u: ℝ^3 →ℝ be an arbitrary smooth function that has a constant value in the normal direction of the surface S near the surface. 
    Then, on the surface S,
    
    ∇ u= ∇_S u

    holds <cit.>.
    Here, ∇_S designates the surface gradient on the surface S.


    Let v be an arbitrary smooth vector field in ℝ^3 that is tangent to the surface S or to a surface that is at a fixed distance from the surface S.
     Then, on the surface S,
    
    ∇· v=∇_S· v

    holds.

If u is a function defined on a surface S, then the function u(C_S(x)) given by the closest point function C_S is a function that has a constant value in the direction of the normal vector of the surface.
Therefore, by Theorem <ref>, we have:

    ∇ u(C_S(x)) = ∇_S u(x),  x∈ S

Furthermore, since ∇ u(C_S(x)) is always tangent to a surface that is at a short and fixed distance from the surface S, Theorem <ref> implies that:

    ∇·∇ u(C_S(x)) = ∇_S·∇_S u(x),  x∈ S

The operator ∇_S·∇_S on the right-hand side of the above equation is intrinsic to the surface S and is denoted by Δ_S (the Laplace-Beltrami operator).  
From the above theorem, it can be seen that by using the extension given by C_S, approximation methods used in Euclidean space can be applied to approximate differential operators such as ∇_S and Δ_S.


In <cit.>, examples of solving partial differential equations on surfaces using CPM and numerical methods are presented. 



As mentioned earlier, in numerical calculations of curvature flow using MBO and HMBO algorithms, an approximation method for surface partial differential equations with constraints is required to achieve the area-preservation condition.
One effective approximation method for constrained partial differential equations is the method of minimizing movements.
In the next section <ref>, we will the extention of minimizing movements to the case of surface PDE.



 §.§ Minimizing movements

Here we will explain the method of minimizing movements (MM), also known as Discrete Morse Flow <cit.>. 
Minimizing movements are a method for approximating the gradient flow of an energy functional

    ℰ(u)=∫_Ω L(∇ u(x),u(x),x) dx

by iteratively minimizing functionals of the form

    ℱ_n(u)=∫_Ω|u-u_n-1|^2/2hdx+ℰ(u)

within a suitable function space and where h>0 is a suitable time step. 
Here, Ω is a region with given boundary conditions, and u_n is an approximation of u at time t=nh. 
The Euler-Lagrange equation of each functional ℱ_n represents an approximation of the gradient flow of the energy functional. 
By changing ℱ_n(u), various approximations of solutions to partial differential equation can be obtained as the Euler-Lagrange equation of ℱ_n. 
For example, for a given α > 0, if we set,

    ℱ_n(u)=∫_Ω|u-u_n-1|^2/2h+α|∇ u|^2/2dx

we obtain an approximation of a heat equation, and if we set

    ℱ_n( u)
        = ∫_Ω| u-2 u_n-1+ u_n-2|^2/2h^2
        +α|∇  u|^2/2dx

we obtain an approximation of the wave equation.
Minimizing movements are based on energy minimizations, so it is possible to naturally handle constrained partial differential equations by adding penalty terms to the functional. 

In the next section we will create an approximation methods for constrained partial differential equations on surfaces by combining minimizing movements with the CPM.



§ NUMERICAL CALCULATION OF PDES ON SURFACES

Here we will provide an overview of an approximation method that combines CPM and MM, and explain its algorithm.
We also perform a numerical convergence analysis for the surface heat and wave equations using our method.


 §.§ Approximation of PDEs on surfaces by CPM.

We will explain the method of discretization in time when approximating solutions to the heat and wave equations on a curved surfaces using CPM.
To this end, let S be a closed smooth surface without boundary in ℝ^3.
We consider the following surface heat equation (<ref>):

    u^S_t(t,x)=αΔ_S u^S(t,x),   x∈ S,    t>0 
    
    	u^S(0,x)=f(x),   x∈ S

and surface wave equation (<ref>): 

    u^S_tt(t,x)=αΔ_S u^S(t,x),    x∈ S,   t>0 
    
            u^S_t(0,x)=V_0(x),    x∈ S           
    
            u^S(0,x)=f(x),   x∈ S

Here, α>0 is a constant, f(x) is the initial condition, V_0 is the initial velocity, and Δ_S is the Laplace-Beltrami operator on the surface S.
We remark that boundary conditions are not included in equations (<ref>) and (<ref>) because the surface S is without boundary.
For a given time step size h>0, we approximate the time derivative in (<ref>) using a forward difference, and in (<ref>) we use a centered difference approximation with respect to time.
By defining u^S_n to be an approximation of u^S at time nh where n=0,1,⋯, we obtain:

    u^S_n+1(x)=u^S_n(x)+hαΔ_S u_n^S(x), 
    
            u_0^S(x)=f(x),   x∈ S

which is a approximation scheme for equation (<ref>).
Similarly, the result for equation (<ref>) is given by:

    u^S_n+1(x)=2u^S_n(x)-u^S_n-1(x)+h^2αΔ_S u_n^S(x), 
    
            u^S_-1(x)=u^S_0(x)-hV_0(x), 
    
            u_0^S(x)=f(x),   x∈ S



Since Δ_S is included in the right-hand side of equations (<ref>) and (<ref>), they are difficult to compute in the general setting.
However, in the CPM, the following equations (<ref>) and (<ref>) are used to compute the solutions in the space Ω surrounding the surface S.
Here, u_n is a function value defined on Ω at time nh.

    u_n+1(x)=u_n(C_S(x))+hαΔ u_n(C_S(x)), 
    
            u_0(x)=f(C_S(x)),
        x∈Ω


    u_n+1(x)=2u_n(C_S(x))-u_n-1(C_S(x))+h^2αΔ u_n(C_S(x)), 
    
            u_-1(x)=u_0(C_S(x))-hV_0(C_S(x)),                                                  
    
            u_0(x)=f(C_S(x)),
        x∈Ω

Here, C_S is defined by equation (<ref>), and Δ = ∇·∇.
Since equations (<ref>) and (<ref>) do not contain Δ_S, it is possible to apply standard numerical approximation techniques in the surrounding Euclidean space calculate surface gradient quantities.
Also, since the surface is given by a point cloud, interpolation can be used to define the numerical solution restricted to the surface S or at any other point in the domain Ω.

Although explicit methods were used to discretize the time derivatives in equations (<ref>) and (<ref>), implicit methods can also be used <cit.>.
The combination of CPM and MM for the calculation of equations (<ref>) and (<ref>) are described in detail in Section <ref>.







 §.§ Combination of CPM and MM

As mentioned in Section <ref>, when performing calculations for curvature flow with an area preservation constraint via the MBO or HMBO algorithms, an approximation method for the constrained partial differential equation is necessary.
Here we explain the approximation method for the constrained partial differential equations on surfaces by combining CPM and MM.



We will explain our method for applying CPM to minimizing movement for the surface heat equation (<ref>), and the surface wave equation (<ref>).
As described in Section <ref>, applying CPM yields the approximations for the surface heat equation (<ref>) and the surface wave equation (<ref>), given by equations (<ref>) and (<ref>), respectively.
As a numerical method, utilizing the method of minimizing movements requires one to approximate functional values. In particular, for n=0,1,⋯, using a time step size h>0 and a constant α>0, the following functional values are required and can be approximated, for example, by means of the finite element method:

    ℱ_n+1(u)
        = ∫_Ω| u(x)- u_n(C_S(x))|^2/2h
        +α|∇  u(x)|^2/2dx


    ℱ_n+1(  u)
        = ∫_Ω|  u(x)-2  u_n(C_S(x))+  u_n-1(C_S(x))|^2/2h^2
        +α|∇   u(x)|^2/2dx

Here, Ω is a sufficiently large region that covers the surface S, and u_n minimizes functional ℱ_n.

In the following, we will show that the Euler-Lagrange equations for equations (<ref>) and (<ref>) lead to the implicitly discretized equations using the CPM method for partial differential equations on surfaces.

Let ϕ be an arbitrary function from C_0^∞(Ω) and ϵ be a real number.
We compute the first variation of equation (<ref>) as follows:

    .d/dϵℱ_n+1(u+ϵϕ)|_ϵ=0=0.

The first variation is:

    d/dϵℱ_n+1(u+ϵϕ)
             =d/dϵ∫_Ω| (u+ϵϕ)- u_n(C_S)|^2/2h
        +α|∇  (u+ϵϕ)|^2/2dx
        =∫_Ω(u+ϵϕ)- u_n(C_S)/hϕ
        +α∇  (u+ϵϕ)·∇ϕ dx.

Substituting ϵ=0 into equation (<ref>), we obtain:

    .d/dϵℱ_n+1(u)|_ϵ=0    =∫_Ωu- u_n(C_S)/hϕ
            +α∇ u·∇ϕ dx
        =∫_Ω( u- u_n(C_S)/h
            -αΔu)ϕ dx+α∫_∂Ω∂ u/∂νϕ dS,

where ∂Ω is the boundary of Ω, and ∂ u/∂ν is the outer normal derivative of u on ∂Ω.
Since ϕ is an arbitrary C_0^∞(Ω) function, the boundary integral in  (<ref>) is zero, and we have:

    .d/dϵℱ_n+1(u)|_ϵ=0    =∫_Ω( u- u_n(C_S)/h
        -αΔu)ϕ dx

A weak form of the Euler-Lagrange equation (<ref>) is therefore:

    ∫_Ω( u- u_n(C_S)/h
        -αΔu)ϕ dx=0

Since ϕ is arbitary, the fundamental lemma of the calculus of variations applies to obtain:

    u- u_n(C_S)/h
        -αΔu=0

which, written as an approximation scheme states:

    u=u_n(C_S)+hαΔu


Equation (<ref>) is an implicit form of the time-discrete surface heat equation (<ref>) obtained by using the CPM.

The functional (<ref>) can be treated in the same fashion. We obtain:

    d/dϵℱ_n+1(u+ϵϕ)
             =d/dϵ∫_Ω| (u+ϵϕ)- 2u_n(C_S)+u_n-1(C_S)|^2/2h^2
        +α|∇  (u+ϵϕ)|^2/2dx
        =∫_Ω(u+ϵϕ)- 2u_n(C_S)+u_n-1(C_S)/h^2ϕ
        +α∇  (u+ϵϕ)·∇ϕ dx.

Setting ϵ=0 in (<ref>), we have:

    .d/dϵℱ_n+1(u)|_ϵ=0    =∫_Ωu- 2u_n(C_S)+u_n-1(C_S)/h^2ϕ
            +α∇ u·∇ϕ dx
        =∫_Ω( u- 2u_n(C_S)+u_n-1(C_S)/h^2
            -αΔu)ϕ dx+α∫_∂Ω∂ u/∂νϕ dS.

As before, ∂ u/∂ν is the derivative of u in the direction of the outer normal vector ν on ∂Ω.
Since ϕ is an arbitrary in C_0^∞(Ω), we obtain the following equation:

    .d/dϵℱ_n+1(u)|_ϵ=0    =∫_Ω( u- 2u_n(C_S)+u_n-1(C_S)/h^2
        -αΔu)ϕ dx.

It follows that

    u- 2u_n(C_S)+u_n-1(C_S)/h^2
        -αΔu=0,

which can be written (<ref>):

    u=2u_n(C_S)-u_n-1(C_S)+h^2αΔu.

Equation (<ref>) is an implicit approximation of the time-discretized surface wave equation (<ref>) using the CPM (compare to Equation (<ref>)).

Having shown that the minimizing schemes above produce approximate solutions to the surface PDE (<ref>) and (<ref>) we now turn to discussing related numerical considerations..
Next, we introduce the computational algorithms for implementing the CPM and MM.


 §.§ Computational methods for the heat and wave equations on surfaces

Here we will explain the computational notions used in our numerical methods.
For simplicity, we will first explain in setting of the surface heat equation (<ref>) on a smooth closed surface S without boundary in three-dimensional Euclidean space. For the sake of clarity, we will also explain the detail in the setting of the surface wave equation (<ref>).

Let α>0 denote the diffusion coefficient and Δ_S denote the Laplace-Beltrami operator on S. Given a time step of h>0, the algorithm for the surface-type minimizing movement that we developed is as follows.



Surface-type minimizing movements for the surface heat equation (<ref>)

    
    
  * Create the computational domain Ω^D by preparing a sufficiently large Cartesian grid covering the surface S.
          Let x_min, x_max, y_min, y_max, z_min, and z_max be the coordinates of the grid boundaries, as shown in Figure <ref>.
          Let the grid spacing in the three spatial directions be given by Δ x, Δ y, and Δ z, respectively.
          Then Ω^D is defined as follows:
          
    Ω^D=.{x_i,j,k =(
                  x_i
    y_j
    z_k)|0≤ i ≤ N_x, 0≤ j ≤ N_y, 0≤ k ≤ N_z}

          where i, j, and k are natural numbers, and N_x, N_y, and N_z denote the number of grid points along the axes of the coordinate system.
          The grid points in the computational domain are expressed as follows:
          
    x_i=x_min+iΔ x,    y_j=y_min+jΔ y,   
                  z_k=z_min+kΔ z,

          
    N_x=x_max-x_min/Δ x,    N_y=y_max-y_min/Δ y,   N_z=z_max-z_min/Δ z

          For simplicity, we assume a uniform grid Δ x=Δ y=Δ z.
          

    
  * Using the closest point function C_S, compute and record the closest point on the surface S for each point in Ω^D.
    
  * To reduce computational cost, calculations are performed only in a vicinity near the surface S.
          This process is called banding.
          In particular, we extract a set of points from Ω^D whose Euclidean distance to the surface S is less than or equal to a constant value λ >0 and denote the region by Ω^D_λ.
          This is expressed as follows, where ||·|| represents the Euclidean norm.
          
    Ω_λ^D={x∈Ω^D |||x-C_S(x)||≤λ}

          

          Remark:
          Ω_λ^D is a point cloud; it consists of discrete points. In the continuous case, a sufficiently large region Ω⊂ℝ^3 covering the surface S is taken, and the region Ω_λ around the surface is defined as follows:
          
    Ω_λ={x∈Ω |||x-C_S(x)||≤λ}

          Remark:
          The value of λ needs to be chosen appropriately, depending on the interpolation method used in Step <ref> below.
          If a polynomial interpolation is used, λ depends on the degree of the interpolation.
          Here, we explain a method for determining λ when performing a linear interpolation in a two-dimensional space (higher dimensions can be treated analogously).
          We assume that the grid points in the computational domain Ω^D have equal spacing in both the horizontal and vertical directions (Figure <ref>(a)).
          To obtain the interpolated value at the point denoted by “⋆" in Figure <ref>(b), four points denoted by “∙" are required.
          In this case, the maximum distance between the interpolation point and the grid points is √(2(Δ x/2)^2).
          The maximum distance occurs in Figure <ref>(c), and its value is √(2(Δ x)^2).
          Therefore, λ must be larger than √(2(Δ x)^2).
          Thus, one choice is to set λ=√((Δ x)^2+2(Δ x)^2) when a linear interpolation is used in a two-dimensional space.
          This discussion can be generalized to the case of a d-dimensional pth degree polynomial interpolation, then we obtain <cit.>: 
    λ=√((d-1)(p+1/2)^2+(1+p+1/2)^2)Δ x

          Since we are considering surfaces in three-dimensional space, we select λ using the interpolation degree p as follows:
          
    λ=√(2(p+1/2)^2+(1+p+1/2)^2)Δ x

          
    
  * When approximating the gradient of a function in Ω_λ^D, information about the boundary points is necessary. 
          We define the characteristic function ϕ_i,j,k as follows:
          
    ϕ_i,j,k=
                      0,  x_i,j,k∈Ω_λ^D
    
                      1,   otherwise,

          from which we define the boundary points ∂Ω_λ^D of Ω_λ^D as follows:
          
    ∂Ω_λ^D={x_i,j,k∈Ω^D |ϕ_i,j,k|∇_Dϕ_i,j,k|≠0}

          where ∇_Dϕ_i,j,k=(ϕ_i+1,j,k-ϕ_i-1,j,k,
              ϕ_i,j+1,k-ϕ_i,j-1,k,
              ϕ_i,j,k+1-ϕ_i,j,k-1)/(2Δ x).
          We then join ∂Ω_λ^D with Ω_λ^D and define it as Ω̂_λ^D, that is,
          
    Ω̂_λ^D=Ω_λ^D∪∂Ω_λ^D.

          Figure <ref> shows the relationship between S, Ω_λ^D, and ∂Ω_λ^D.
          Figure <ref> is a schematic diagram of the section of Figure <ref>.
    
  * 
          Using the closest point function C_S, extend the initial condition given on the surface S to Ω^D as follows, where the initial condition at point x_i,j,k is denoted by u^0_i,j,k.
          
    u^0_i,j,k=
                  
                      u^S_0(C_S(x_i,j,k)),  x_i,j,k∈Ω̂_λ^D 
    
                      0,   x_i,j,k∈Ω^D ∖Ω̂_λ^D


          Remark:
          In order to simplify the calculations, the initial values of the grid points outside of Ω̂_λ^D are set to 0.
          
    
  * Obtain an approximate solution of the heat equation on Ω_λ^D using MM.
          To approximate u on Ω_λ^D in equation (<ref>), let u_i,j,k=u(x_i,j,k).
          Approximate the functional values of (<ref>) by means of an expression such as:
          
    ℱ_n(u)≈Δ x^3∑_x_i,j,k∈Ω_λ^D{|u_i,j,k-u_i,j,k^n-1|^2/2h+α(∇_D,xu_i,j,k)^2+(∇_D,yu_i,j,k)^2+(∇_D,zu_i,j,k)^2/2}

          Denote the minimizer of this functional by u_n, where u=(u_i,j,k).
          Here, Δ x^3 is the volume of the element, and ∇_D,xu_i,j,k, ∇_D,yu_i,j,k, ∇_D,zu_i,j,k are calculated by difference approximations as follows:
          
    ∇_D,xu_i,j,k=u_i+1,j,k-u_i-1,j,k/2Δ x

          
    ∇_D,yu_i,j,k=u_i,j+1,k-u_i,j-1,k/2Δ x

          
    ∇_D,zu_i,j,k=u_i,j,k+1-u_i,j,k-1/2Δ x

          Note that, as mentioned earlier, we assume Δ x=Δ y=Δ z.

          Remark:
          Various methods can be used to obtain the minimizer of (<ref>). Among them, from the viewpoint of computational cost, the L-BFGS method is often used <cit.>.

    
  * Create an interpolating function  I_n(x) defined on Ω_λ for the minimizer obtained in step <ref>. 
          Using I_n(x), define u_n^S(x)=I_n(x) for x∈ S.
          It should again be noted that I_n(x) is defined on Ω_λ. There are multiple methodologies for its construction.
          One example is to use trilinear interpolation, which is a linear interpolation in 3D <cit.>.
          The computations in this study have used polynomial interpolations.

    
  * Using the closest point function C_S, extend u_n^S onto Ω̂_λ^D as follows:
          
    u^n_i,j,k=u_n^S(C_S(x_i,j,k)),   x_i,j,k∈Ω̂_λ^D

    
  * Repeat steps <ref> to <ref> for n=1,2,⋯ until the desired final time is reached.


Next, we will explain the computational algorithm for the surface wave equation (<ref>).



Surface-type minimizing movements for the surface wave  equation (<ref>)

    
    
  * Perform the computations in Steps <ref> to <ref> of the previous algorithm.
    
  * Assign u_i,j,k^-1 using the initial velocity V_0 of equation (<ref>).
        This can be done, for example, by means of the backward difference approximation u_i,j,k^-1=u_i,j,k^0-hV_0(x_i,j,k).
          Note that u_i,j,k^-1 represents the value at the grid point x_i,j,k at time -h.
    
  *  Compute an approximate solution to the wave equation on Ω_λ^D using MM.
          Similar to the case of the heat equation, the functional values in (<ref>) can be approximated as follows, for n=1,2,⋯:
          
    ℱ_n(u)≈Δ x^3∑_x_i,j,k∈Ω_λ^D{|u_i,j,k-2u_i,j,k^n-1+u_i,j,k^n-2|^2/2h^2+α(∇_D,xu_i,j,k)^2+(∇_D,yu_i,j,k)^2+(∇_D,zu_i,j,k)^2/2}

          The minimizer of this functional is denoted by u_n, where u=(u_i,j,k).
          Here, Δ x^3 is the volume of the element, and ∇_D,xu_i,j,k, ∇_D,yu_i,j,k, ∇_D,zu_i,j,k are calculated in the same way as in Step <ref> of the previous algorithm.
    
  * Define u_n^S using the minimizer obtained in Step <ref> by employing the same procedure as in Step <ref> of the previous algorithm.
    
  * Using the closest point function C_S, extend u_n^S onto Ω̂_λ^D.
    
  * Repeat steps <ref> to <ref> for n=1,2,⋯ until the desired final time is reached.


In the next section, we perform a numerical error analysis using the above algorithms for the heat and wave equations on the surface of the unit sphere.
We will begin by treating the case of the surface heat equation.




 §.§ Numerical error analysis of MM for the heat equation on a surface

Here, we will perform an numerical error analysis of the algorithm for solving the surface heat equation, described in the previous section.
Unsing MM, we numerically solve the surface heat equation (<ref>) on the unit sphere S, and examine the error between the numerical solution and the exact solution.
We define the unit sphere S in the 3D space as follows:

    S={.(
        [ sinθcosϕ; sinθsinϕ;     cosθ ])|
        0≤θ≤π,0≤ϕ≤ 2π}


We perform two numerical experiments by changing the initial condition f in equation (<ref>) on the unit sphere S.
First, we explain the initial conditions used and their corresponding exact solutions.
The results of the numerical error analysis are presented in Section <ref>.


 §.§ MM and and surface heat equation: initial condition 1

Setting the diffusion coefficient to α=1, we take the initial condition f as

    f(θ)=cosθ

The exact solution of equation (<ref>) is then given by

    u(θ ,ϕ,t)=e^-2tcosθ,   t≥  0



 §.§ MM and and surface heat equation: initial condition 2

Setting the diffusion coefficient to α=1/42, we take the initial condition f as

    f(θ,ϕ)=Y^0_6(θ,ϕ)+√(14/11)Y^5_6(θ,ϕ)

where Y^m_l(θ,ϕ) are the eigenfunctions of the Laplacian on the unit sphere, known as spherical harmonics.
The exact solution of equation (<ref>) is then given by

    u(θ ,ϕ,t)=e^-t{Y^0_6(θ,ϕ)+√(14/11)Y^5_6(θ,ϕ)},   t≥ 0

as shown in <cit.>.


The results of the numerical error analysis using initial conditions 1 and 2 (described above) for the surface heat equation are described in the next section.


 §.§ MM numerical error analysis results (heat equation on the unit sphere)

We investigate the relationship between Δ x and the numerical error of the MM approximation to the solution of the surface heat equation. Computations follow the computational algorithm for the heat equation on surfaces (<ref>), presented in Section <ref>, where the spatial discretization Δ x is varied.
The L-BFGS method is used to minimize the discretized functional (<ref>). We implement the method using Optim.jl <cit.>, and calculate the functional gradient using automatic differentiation (ReverseDiff.jl <cit.> is used for this purpose).
The time step h is set to h = Δ x^2/6, and polynomial interpolation of order p=2 is used (see Section <ref>).
For both initial conditions 1 and 2, we calculate the maximum absolute error L_∞ on S at the closest point to each point in Ω_λ^D at time t_e=0.25. We note that, since the exact solution of the surface heat equation converges to 0, it becomes difficult to evaluate the error between the numerical solution and the exact solution. For this reason, we have chosen such a value of t_e (i.e., so that the L_∞-error of the absolute value of the exact solution at time t_e is sufficiently large for both initial conditions 1 and 2).
The L_∞-error is defined as

    L_∞`-error=sup_x∈Ω_λ^D|u(C_S(x),t_e)-û(C_S(x),t_e)|

where û is the numerical solution and u denotes the exact solution.


The results obtained for each Δ x are shown in Table <ref> and Table <ref>.
The results are plotted in Figure <ref>(heaterr) and Figure <ref>(heaterr_log) using both regular and log-log scales, respectively.
The legend in the figure denotes initial condition 1 by cond1, and initial condition 2 by cond2.
The time evolution is shown in Figure <ref> and Figure <ref>.
The results confirm that the L_∞-error decreases as Δ x decreases, except for the case where Δ x=0.0125.
Except for this case, the numerical error is roughly proportional to the square of Δ x when Δ x is reduced by a factor of 2.
The reason for the larger numerical error in initial condition 2 may be due to insufficient resolution relative to the initial condition.

These results confirm that the numerical solution obtained by MM converges to the exact solution of the surface heat equation as the spatial discretization converges to zero.
We also note that we observe the numerical error increases when Δ x becomes numerically too small (Figure <ref>(heaterr_log)).
Next, we will perform a numerical error analysis for the wave equation on a curved surface.













 §.§ Numerical error Analysis for the surface wave equation

Using the MM approximation scheme previously described, we numerically solve the wave equation on a curved surface (<ref>) and examine the error between the numerical and the exact solution.
We perform two computational experiments by changing the initial condition f of equation (<ref>) on the unit sphere S expressed by equation (<ref>).
First, we explain the initial conditions used and their exact solutions.
The results of the numerical error analysis are presented in Section <ref>.


 §.§ Surface wave equation: initial condition 1

With the constant α=1, we set the initial condition f as

    f(θ)=-cosθ

and the initial velocity V_0 as

    V_0=0

In this case, the exact solution of equation (<ref>) is given by

    u(θ,ϕ,t)=-cos(√(2)t)cosθ,   t≥  0



 §.§ Surface wave equation: initial condition 2

Let α = 1 and the initial condition f be

    f(θ ,ϕ)=Y^0_6(θ,ϕ)+√(14/11)Y^5_6(θ,ϕ)

with initial velocity V_0 = 0.
Then, the exact solution of equation (<ref>) is given by

    u(θ,ϕ,t)=cos(√(42)t){Y^0_6(θ,ϕ)+√(14/11)Y^5_6(θ,ϕ)},   t≥  0

Next, we will explain the results of the numerical error analysis using initial conditions 1 and 2 for the surface wave equation.


 §.§ Numerical error analysis results (wave equation on the unit sphere)

We performed several numerical simulations using the algorithm presented in Section <ref> for the wave equation (<ref>) and investigated the relationship between the numerical error and the spatial discretization Δ x.
We used the same optimization methods and interpolation degree p as those used for the surface heat equation.
The time step h was set to Δ x/10.
For initial condition 1, the absolute error L_∞ is calculated at the closest point on S for each point in Ω_λ^D at the time t=2π/√(2).
For initial condition 2, the absolute error L_∞ is calculated at the closest point on S for each point in Ω_λ^D at the time t=2π/√(42).
The maximum value of the absolute error is then obtained for each case.
We evaluated the L_∞ error at the time t_e such that the exact solution using initial conditions 1 and 2 have both  oscillated once over 0≤ t≤ t_e.
The L_∞ error is defined as follows:

    L_∞`-error=sup_x∈Ω_λ^D|u(C_S(x),t_e)-û(C_S(x),t_e)|

where û denotes the numerical solution, and u denotes the exact solution.

The results obtained for Δ x and L_∞`-error are presented in Table <ref> and Table <ref>.
Figure <ref>(waveerr) shows a graphical representation of the contents of Table <ref> and Table <ref>, while Figure <ref>(waveerr_log) shows a representation on log-log scale.
In the legend of the figures, cond1 corresponds to initial condition 1, and cond2 corresponds to initial condition 2.
The evolution of the solution over time is shown in Figure <ref> and Figure <ref>.
The results show that the L_∞`-error decreases with Δ x.
Except for when Δ x=0.2, the numerical error is approximately halved when Δ x is halved.
That is, it is observed that the numerical error converges proportional to Δ x.
Similar to the numerical error analysis for the heat equation, it was found that the numerical error for the initial condition 2 is greater than that for initial condition 1.
These results indicate that when the spatial grid spacing Δ x is sufficiently small, the numerical solution obtained by the developed numerical method for the wave equation converges to the exact solution.


In this section, we explained the approximation methods for partial differential equations on  curved surfaces using the CPM and MM methods, and presented the results of their numerical error analysis.
In the following sections, we will illustrate applications of the approximation methods to the simulation of interfacial motions on curved surfaces.













§ NUMERICAL SIMULATION OF INTERFACIAL MOTIONS ON SURFACES

In this section, we discuss approximation methods for realizing mean curvature flow and hyperbolic mean curvature flow on surfaces.
After explaining the interfacial motions on curved surfaces, we present the numerical results of our schemes.
We also deal with the mean curvature flow and the hyperbolic mean curvature flow on curved surfaces in the multiphase setting and under area preservation conditions.
Multiphase regions are realized by means of a signed distance vector field, which allows us to incorporate area preservation constraints into the functional minimizations of the MM approach.


 §.§ Surface-constrained interfacial motions

Here we will discuss the approximation methods for the mean curvature flow (MCF) <cit.> and the hyperbolic mean curvature flow (HMCF) <cit.> on curved surfaces.
The surface MCF and surface HMCF are described by the following nonlinear partial differential equations, respectively.


    γ^S_t(t,s)=-κ^S(t,s)ν^S(t,s), 
    γ^S(0,s)=γ^S_0(s)Surface MCF


    γ^S_tt(t,s)=-κ^S(t,s)ν^S(t,s), 
    γ^S(0,s)=γ^S_0(s),                 
    γ_t^S(0,s)=v_0ν^S(t,s)Surface HMCF

Here, S is a smooth surface, γ^S:[0,T)× [a,b]→ S is a smooth simple curve on the surface S, satisfying γ^S(t,a)=γ^S(t,b), κ^S is the curvature of the curve, γ_0 is the initial shape of the curve, v_0 is the initial velocity of the curve, and ν^S represents the outward unit normal vector of the curve on the surface S.
Here, γ^S_t=∂γ^S/∂ t,  γ^S_tt=∂γ^S/∂ tt.
Surface MCF and surface HMCF are equations that generalize the motion of surfaces following mean curvature flow or hyperbolic mean curvature flow in the Euclidean space to the setting of interfaces moving on curved surfaces.
Interfaces moving by surface MCF tend to decrease their length and smoothing their shape over time, while interfaces moving by surface HMCF tend to oscillate.

In the case that surface MCF and surface HMCF are subject to the area-preservation conditions, the interfaces should move while preserving the areas of the regions enclosed by the interfaces.
We handle such motions by extending the MBO algorithm and the HMBO algorithm to the surface PDE setting using CPM, MM, and a surface version of the signed distance vector field, to develope approximate solutions for surface MCF (surface MBO) and surface HMCF (surface HMBO).
We remark that our methods can also handle interfacial motions with area-preservation conditions  in the multiphase setting.
This is enabled by means of a signed distance vector field that is used to encod the shape of interfaces atop the surface.


 §.§ The signed distance vector field on surfaces

In this section, we discuss the signed distance vector field <cit.> and its extension to the surface setting.
The signed distance vector field (SDVF) is used to encode the shape of multiphase regions by means of vector directions.
When performing numerical calculations under area preserving conditions, the signed distance can be used in the two-phase setting.
However, for interface motions involving three or more phases and area preservation conditions, it is not possible to distinguish each phase using a single signed distance function.
On the other hand, the SDVF can be used to distinguish phase locations and shapes even in the case of three or more phases.
The SDVF is constructed by assigning a special vector to each phase of the multiphase region. Each vector is weighted by its signed distance from each interface <cit.>. The SDVF is described below.

Let K be the number of phases, ϵ > 0 be an interpolation parameter, P^i be the region of phase i, p_i be the vector from the barycenter of a (K-1)-dimensional simplex to each vertex (refer to Figure <ref> for K=3, <cit.>), d_S^i(x) be the signed distance function to phase i, and χ_E be the characteristic function of the set E.
The signed distance vector field on the surface S (Surface SDVF) z_S^ϵ is given by the following equation:

    z_S^ϵ (x)=∑_i=1^K(p_iχ_{d_S^i≥ϵ/2}+1/ϵ( ϵ/2+d_S^i)p_iχ_{-ϵ/2<d_S^i <ϵ/2}),   x∈ S

where,

    χ_E(x)=
            1  x∈ E, 
    
            0  otherwise,
        
        d_S^i(x)=inf_y∈∂ P^i||x-y||_S  x∈ P^i, 
    
            -inf_y∈∂ P^i||x-y||_S  otherwise.

Here, ||x-y||_S represents the geodesic distance on the surface between the two points x and y, and can be expressed as the value that minimizes the following length functional:

    ||x-y||_S=min_Γ⊂ SLength(Γ)

where Γ is a curve on the surface connecting x and y, and Length(Γ) is the length of Γ along the surface.

Remark:
From here on we will omit the S in z_S^ϵ and d_S^i, and denote them simply by z^ϵ and d^i, respectively.



In the next Section <ref>, we introduce an approximation method for Surface MCF, based on the MBO algorithm using the surface SDVF. We refer to appoximation method as the surface MBO algorithm.
In Section <ref>, we introduce an approximation method for the surface HMCF, which is based on the HMBO algorithm and the surface SDVF. Similarly, we refer to our approximation method as the surface HMBO.




 §.§ Surface MBO

In this section, we discuss an approximation method for surface MCF.
We begin by describing the surface MBO in the two-phase setting.
Then we explain the surface MBO with and without the area preservation condition in the multiphase setting.
The computational results of our surface MBO in the multiphase setting are presented in Section <ref>.

We use the CPM as an approximation for the surface partial differential equation on the surface.
To reduce the computational cost of the CPM, we limit the computations to a small tubular region around S, as described in Section <ref>.
This region is defined by accumulating points from the computation grid that are located within a constant distance λ from S as Ω_λ (see equation (<ref>)).
We determine the value of λ in the same way as in equation (<ref>).
We use a natural number N and the final time T>0 to define the time step τ=T/N.


  §.§.§ Surface MBO for two-phase regions

Below, we explain a method for approximating the MCF motion of an interface in a 2-phase region. Here, 
γ^S_n represents the shape of the curve at time nτ, where τ=T/N is the time step size, T>0 is the final time, and γ^S_0 is the initial curve.
Let d_n(x) be the signed distance function from the point x on the surface to the interface γ^S_n.
That is, d_n(x) is defined as follows:

    d_n(x)=inf_y∈∂ P^n||x-y||_S,     x∈ P^n, 
    -inf_y∈∂ P^n||x-y||_S,     otherwise.

where P^n is the region occupied by phase n. 

The surface MBO for a 2-phase domain is as follows:

    



    
  1. Create d_0 using Equation (<ref>) from the initial curve γ^S_0.
    
  2. Extend d_0 to Ω_λ:
        
    d_0^λ(x)=d_0(C_S(x)),   x∈Ω_λ

    
  3. Repeat the following for n=0,1,⋯,N-1:
        
            
            
  * For t∈[0,τ), solve the heat equation in Ω_λ:
                  
    u_t=αΔ u(x,t), 
    
                              u(x,0)=d^λ_n(x),
                            x∈Ω_λ,  τ>t>0

                  where α>0 is a constant representing the diffusion coefficient.
            
  * Define the new curve γ^S_n+1 as the zero level set of the solution u(x,τ) of Equation (<ref>) on the surface S:
                  
    γ^S_n+1={x∈Ω_λ|S∩{u(x,τ)=0}}

            
  * Create d_n+1 from γ^S_n+1 using Equation (<ref>).
            
  * Extend d_n+1 to Ω_λ and define d^λ_n+1 as follows:
                  
    d_n+1^λ(x)=d_n+1(C_S(x)),  x∈Ω_λ

        

The results of numerical error analysis using the proposed algorithm are shown in Section <ref>.
In the above algorithm, the sign of the signed distance function is used to extract the interface.
In the next section, we will implement the surface MBO for multiphase regions using the signed distance vector field instead of the signed distance function.


  §.§.§ Surface MBO for multiphase regions

Here, we explain a method for approximating surface MCF on interfaces consisting of K≥ 2 multiphase regions on the surface S.
Let P^i_n represent the region of phase i at time nτ, and let P_n=⋃_i=1^KP^i_n.
We provide initial regions P_0^i for each phase i as initial conditions.
Additionally, we denote the vector given to phase i as p_i, as explained in Section <ref>.
The surface SDVF, written z_n^ϵ(x), is obtain from (<ref>) using P_n.
The surface MBO for multiphase regions described as follows.

    
    
  1. Using equation (<ref>), create the surface SDVF z_0^ϵ from the initial domain P_0.
    
  2. Extend z_0^ϵ to Ω_λ and denote it as z_0^ϵ,λ as follows:
        
    z_0^ϵ,λ(x)=z_0^ϵ(C_S(x)),   x∈Ω_λ

    
  3. Repeat the following for n=0,1,⋯,N-1:
        
            
            
  * Solve the following vector-valued heat equation for t∈[0,τ):
                  
    u_t=αΔu(x,t),          
    u(x,0)=z^ϵ,λ_n(x), 
      x∈Ω_λ,  τ>t>0

                  Here, α>0 is a constant representing the diffusion coefficient.
            
  * Extract the solution u(x,τ) of equation (<ref>) on the surface S and denote it as u^S as follows:
                  
    u^S(x)=u(x,τ),  x∈ S∩Ω_λ

            
  * Obtain P_n+1 on the surface S using u^S as follows:
                  
    P_n+1=⋃_i=1^K{P^i_n+1}   P^i_n+1 ={x∈ S;u^S(x)·p_i≥u^S(x)·p_j, for all j∈{1,⋯,K}}

            
  * Update the surface SDVF z_n+1^ϵ from P_n+1 using Equation (<ref>).
            
  * Extend z_n+1^ϵ to Ω_λ:
                  
    z_n+1^ϵ,λ(x)=z_n+1^ϵ(C_S(x)),   x∈Ω_λ

        

Next, we will implement the surface MBO with a prescribed area-constraint for multiphase regions by combining MM with the above algorithm.



  §.§.§ Multiphase surface MBO with area constraints

Here, we will explain our method for approximating multiphase surface MCF under K area constraints on the surface S.
The initial vector field is constructed from the regions P_0^i for each phase i, where the vector for each phase i, denoted by p_i, is prescribed as in Section <ref>.
The sign-distance vector field z_n^ϵ(x) is then obtained via equation (<ref>) from P_n.
When approximating constrained interfacial motions, MM are often used when treating the case of area-preservation <cit.>.

In particular, MM can be used in combination with penalty for each area constraint.

We take a sufficiently large positive constant M and set h=τ/M.
Given w_m-1, we define the functional (<ref>) used in the MM as follows:

    ℱ_m(w)=∫_Ω_λ(|w-w_m-1|^2/2h+α|∇w|^2/2)dx
        +ρ∑^K_i=1|A^i-V^i_w|^2,

where we use the extension (<ref>) of the surface SDVF to Ω_λ.
In (<ref>), α>0 and ρ>0 are constants, and

    A^i=V^i_z_0^ϵ,λ,   V^i_w=∫_Ω_λH^ϵ(ϕ^i_w(x))dx,
    ϕ^i_w(x)=inf_y∈∂ Q^i_w||x-y||_S,   x∈ Q^i_w, 
    
                -inf_y∈∂ Q^i_w||x-y||_S,   otherwise,
              
            H^ϵ(u)=
                1,   u>ϵ
    1/2+u/2ϵ+1/2πsinπ u/ϵ,    -ϵ≤ u≤ϵ
    
                0,   u<-ϵ
    
    
            Q^i_w={x∈Ω_λ|w(x)·p_i≥w(x)·p_j, for all j∈{1,⋯,K}},

where z^ϵ,λ_0 is the extension of z^ϵ_0 obtained from the initial region P_0 on the surface S to Ω_λ.
The parameter ρ controls the strength of the penalty. 
Note that when ρ=0, equation (<ref>) takes the same form as the vectorized version of equation (<ref>), and applying the following algorithm results in a numerical solution of surface MCF without area preservation.

The surface MBO for realizing multiphase area-preserving surface MCF is as follows:

    

    
  1. Using equation (<ref>), create the surface SDVF z_0^ϵ from the initial domain P_0.
    
  2. Extend z_0^ϵ to Ω_λ and define:
        
    z_0^ϵ,λ(x)=z_0^ϵ(C_S(x)),   x∈Ω_λ

    
  3. Using (<ref>), obtain A^i from z_0^ϵ,λ.
    
  4. Repeat the following for n=1,2,⋯,N-1:
        
            
            
  * Set w_0 = z^ϵ,λ_n-1(x).
            
  * For m=1,⋯,M, find the minimizer w of the functional ℱ_m(w) (refer to equation (<ref>)). Denote the minimizer by w_M.
            
  * Extract w_M on the surface S and denote it by w^S:

                  
    w^S(x)=w_M(x),  x∈ S∩Ω_λ

            
  * Obtain P_n+1 on the surface S using w^S as follows:
                  
    P_n+1   =⋃_i=1^K{P^i_n+1},
    
                          P^i_n+1   ={x∈ S|w^S(x)·p_i≥w^S(x)·p_j, for all j∈{1,⋯,K}}

            
  * Create the surface SDVF z_n+1^ϵ from P_n+1 using equation (<ref>).
            
  * Extend z_n+1^ϵ to Ω_λ and denote it as z_n+1^ϵ,λ as follows:
                  
    z_n+1^ϵ,λ(x)=z_n+1^ϵ(C_S(x)),   x∈Ω_λ.

        

The numerical examples using the approximation methods presented in this section are shown in Section <ref>. 
In the next section, we will explain an approximation methods for hyperbolic mean curvature flow on surfaces.


 §.§ Surface HMBO

In this section, we discuss an approximation method for surface HMCF on a surface S. 
First, we explain the surface HMBO for two-phase regions. Then, we describe its multiphase counterpart without area preservation. Afterwards, the surface HMBO for multiphase regions with area constraints is presented. 
The numerical results using the surface HMBO for multiphase regions is presented in Section <ref>.

As before, let Ω_λ be a tubular region of S where the distance from S is within a constant value λ as in Eq. (<ref>). 
The value of λ is determined using the same method as described at  (<ref>). 
The time step used in the implemented algorithm is denoted by τ=T/N, where N is a natural number representing the number of steps and T>0 is the final time.


  §.§.§ Surface HMBO for two-phase regions

Below, we explain the method for approximating the motion of interfaces in a 2-phase region evolving by surface HMCF. 

Given an initial curve γ^S_0 and an initial velocity v_0, let  γ^S_n represent the shape of the curve at time nτ. 
Let d_n(x) be a signed distance function on the surface to the interface γ^S_n. 
That is, if we let P_n be the region enclosed by γ_n^S, d_n(x) can be expressed by (<ref>).
The surface HMBO for a 2-phase region is then described as follows:

    

    
  1. Define γ^S_1=γ^S_0+τ v_0 using the initial curve γ^S_0 and the initial velocity v_0.
    
  2. Create d_0 and d_1 from γ^S_0 and γ^S_1 using equation (<ref>).
    
  3. Extend d_0 and d_1 to Ω_λ and denote them by d^λ_0 and d^λ_1 respectively, as follows:
        
    d_l^λ(x)=d_l(C_S(x)),  x∈Ω_λ, l=0,1

    
  4. Repeat the following for n=1,2,⋯,N-1:
        
            
            
  * Solve the following wave equation: 
                  
    u_tt=αΔ u,                                                              
    
                              u(x,0)=2d^λ_n(x)-d^λ_n-1(x), 
    
                              u_t(x,0)=0
                            x∈Ω_λ,  τ>t>0

                  where α>0 is a constant.
            
  * Define the new curve γ^S_n+1 as the zero level set of the solution u(x,τ) of equation (<ref>) on the surface S:
                  
    γ^S_n+1={x∈Ω_λ|S∩{u(x,τ)=0}}

            
  * Create d_n+1 from γ^S_n+1 using equation (<ref>).
            
  * Extend d_n+1 to Ω_λ and define d^λ_n+1 as follows:
                  
    d_n+1^λ(x)=d_n+1(C_S(x)),   x∈Ω_λ.

        

The results of the numerical error analysis using the above algorithm are presented in Section <ref>. 
In the above algorithm, the signed distance function is used to detect the location of the interface. 
As before, by using the signed distance vector field instead of the signed distance function, we can implement the surface HMBO in the multiphase setting. This will be explained in the next section.


  §.§.§ Surface HMBO for multiphase regions

In the following, we describe a method for approximating surface HMCF of  interfaces separating K multiphase regions on a surface S. 
To this end, let P^i_n represent the region of phase i at time nτ, and P_n=⋃_i=1^KP^i_n. 
In addition to the initial regions P_0^i, we also provide the initial velocities v_0^i for each phase i. 
We denote the vector given to phase i as p_i, as explained in Section <ref>. 
Again, z_n^ϵ(x) denotes surface SDVF defined by equation (<ref>) and which is constructed from P_n. 

The surface HMBO algorithm for multiphase regions is as follows:

    

    
  1. Determine P_1 from the initial domain P_0 and the initial velocities v_0^i on ∂ P_0^i. For details, refer to Section <ref>.
    
  2. Using equation (<ref>), create the surface SDVFs  z_0^ϵ, and z_1^ϵ, from P_0,P_1.
    
  3. Extend z_0^ϵ, and z_1^ϵ to Ω_λ, and denote their extensions by z_0^ϵ,λ and z_1^ϵ,λ respectively, as follows:
        
    z_l^ϵ,λ(x)=z_l^ϵ(C_S(x)),   x∈Ω_λ,   l=0,1

    
  4. Repeat the following for n=1,2,⋯,N-1:
        
            
            
  * Solve the vectorial wave equation: 
                  
    u_tt=αΔu,           
                              
    u(x,0)=2z^ϵ,λ_n(x)-z^ϵ,λ_n-1(x), 
    u_t(x,0)=0,
                            x∈Ω_λ,  τ>t>0

                  where, α>0 is a constant.

            
  * Extract the solution u(x,τ) of equation (<ref>) on the surface S and denote it by  u^S:

                  
    u^S(x)=u(x,τ),  x∈ S∩Ω_λ.

            
  * Obtain P_n+1 on the surface S using u^S as follows:
                  
    P_n+1   =⋃_i=1^K{P^i_n+1},
    
                          P^i_n+1   ={x∈ S|u^S(x)·p_i≥u^S(x)·p_j, for all j∈{1,⋯,K}}

            
  * Create the surface SDVF z_n+1^ϵ from P_n+1 using equation (<ref>).
            
  * Extend z_n+1^ϵ to Ω_λ and denote it by z_n+1^ϵ,λ as follows:
                  
    z_n+1^ϵ,λ(x)=z_n+1^ϵ(C_S(x)),   x∈Ω_λ

        

In the next section, we show how MM can be combined with the above surface HMBO algorithm to realize multiphase surface HMCF of interfaces with optional area-preserving conditions.


  §.§.§ Surface HMBO for multiphase  area-preserving motions

In the following, we describe our method for approximating multiphase surface HMCF, where the area of each domain's preserved. 
We assume that the initial conditions for each phase i are given by the initial domain P_0^i and the initial velocity v_0^i. 
As in section <ref>, we denote the vector given to phase i by p_i.


Similar to section <ref>, MM are used to incorporate the area-preserving conditions. We take a sufficiently large integer M and set h=τ/M. 
Let the functions w_m-1 and w_m-2 be defined using the surface SDVF extended to Ω_λ and the MM functional be as follows:


    ℱ_m(w)=∫_Ω_λ(|w-2w_m-1+w_m-2|^2/2h^2+α|∇w|^2/2)dx
        +ρ∑^K_i=1|A^i-V^i_w|^2

Here, α>0 and ρ>0 are constants, and A^i and V_w^i are defined by (<ref>).

The definitions of w_1 and w_0 are as described before. Namely, we use z^ϵ,λ_0 (the extension of z^ϵ_0 obtained from the initial domain P_0) together with the initial velocities v^i_0  on the surface S. 
Here, ρ>0 is a sufficiently large penalty parameter. 
Note that when ρ=0, the functional (<ref>) is the same as the vectorized extension of equation (<ref>), and that applying the algorithm below would result in a numerical approximation of surface HMCF without area preservation. 

The Surface HMBO that realizes area preservation for multiphase domains is as follows.


    
     
  1. 
        Determine P_1 from the initial domain P_0 and initial velocities v_0^i on ∂ P_0^i. For details, refer to Section <ref>.
    
  2. 
        Using equation (<ref>), create the surface SDVF z_0^ϵ and z_1^ϵ from P_0 and P_1.

    
  3. Define the extensions of z_0^ϵ and z_1^ϵ to Ω_λ  as follows:
        
    z_l^ϵ,λ(x)=z_l^ϵ(C_S(x)),   x∈Ω_λ,   l=0,1

    
  4. 
    Using equation (<ref>), compute each A^i from z_0^ϵ,λ.
    
  5. Repeat the following for n=1,2,⋯,N-1:
        
            
            
  * Set w_0=w_1=2z^ϵ,λ_n(x)-z^ϵ,λ_n-1(x).
            
  * For each m=2,⋯,M, minimize ℱ_m(w) (given by (<ref>)) and denote each minimizer by w_m.
                  
            
  * Let w^S denote the restriction of  w_M to the surface S as follows:

                  
    w^S(x)=w_M(x),  x∈ S∩Ω_λ

            
  * Obtain P_n+1 on the surface S using w^S as follows:
                  
    P_n+1   =⋃_i=1^K{P^i_n+1}
     P^i_n+1   ={x∈ S|w^S(x)·p_i≥w^S(x)·p_j, for all j∈{1,⋯,K}}

            
  * Create Surface SDVF z_n+1^ϵ from P_n+1 using Equation (<ref>).
            
  * Let z_n+1^ϵ,λ be the extension of  z_n+1^ϵ to Ω_λ:
                  
    z_n+1^ϵ,λ(x)=z_n+1^ϵ(C_S(x)),   x∈Ω_λ.

        

Numerical examples using the method described in this section are presented in section <ref>.


§ NUMERICAL RESULTS AND CONSIDERATIONS


In this section, we use the approximation methods presented in sections <ref> and <ref> to numerically solve the mean curvature flow and hyperbolic mean curvature flow on surfaces under various conditions. 

The discrete approximation of Ω_λ uses a uniformly spaced orthogonal grid Ω_λ^D with a spacing of Δ x in all three directions.
Note that Ω_λ^D is obtained using the same method as in section <ref>. 
Details regarding the approximation of the MM functional values are explained in section <ref>. 
In all cases, the interpolation parameter ϵ used to construct the surface SDVF represented by Equation (<ref>) is set to ϵ=0.03.

Remark:
The interpolation parameter ϵ needs to be appropriately selected depending on the discretization of the surface S. 
In the case that the surface is discretized by a point cloud, it was found from the numerical investigations in this study that a value of 3-5 times the average distance to neighboring points within the point cloud is appropriate for the interpolation parameter ϵ.


 §.§ Regarding the initial conditions

Boundaries between regions determine the shape of the interface and hence the initial conditions used in the numerical calculations. 
The following two types of initial conditions for the numerical calculations on the unit sphere were used. In both cases, points with different colors indicate different phases.

    
  * Two-phases on the unit sphere
          
    
  * Four-phases on the unit sphere
          
          

          (Left: Figure viewed from an oblique angle. Right: Figure viewed from directly above.)




 §.§ Computational details regarding surface mean curvature flow

In the following, we introduce the parameters used in our computations of surface mean curvature flow. 
We performed computations of the surface  mean curvature flow involving interfaces in two and  four phase environments on the unit sphere, with and without area preservation.

Discussions of the corresponding computations are presented in section <ref>.

    
  Result of two-phase mean curvature flow on the unit sphere (without area preservation) 


    We use the two-phase initial condition shown in Figure <ref> and the algorithm described in section <ref> for the numerical calculation. 
    Parameters were set as follows:
    
    α=1.0,  Δ x=0.05,    h= Δ x^2/6,  τ=15h,  ρ=0

    The numerical result is shown in Figure <ref>.
    
  Result of two-phase mean curvature on the unit sphere (with area preservation) 


    We used the two-phase initial condition shown in Figure <ref> and the algorithm described in section <ref> for the numerical calculation. 
    Parameters were set as follows:
    
    α=0.05,  Δ x=0.05,    h= Δ x^2/6,  τ=100h,  ρ=10^3

    The numerical result is shown in Figure <ref>.
    
  Result of four-phase mean curvature flow on the unit sphere (without area preservation) 


    We used the four-phase initial condition shown in Figure <ref> and the algorithm described in section <ref> for the numerical calculation. 
    Parameters were set as follows:
    
    α=1.0,  Δ x=0.05,    h= Δ x^2/6,  τ=15h,  ρ=0

    The numerical result is shown in Figure <ref>.
    
  Result 1 of four-phase mean curvature flow on the unit sphere (with area preservation) 


    We used the four-phase initial condition shown in Figure <ref> and the algorithm described in Section <ref> for the numerical calculation. 
    Parameters were set as follows:
    
    α=1.0,  Δ x=0.015,    h= Δ x^2,  τ=15h,  ρ=10^5

    The numerical result is shown in Figure <ref>.
    
  Result 2 of four-phase mean curvature flow on the unit sphere (with area preservation) 

    We used the four-phase initial condition shown in Figure <ref> and the algorithm described in Section <ref> for the numerical calculation. 
    Parameter were set as follows:
    
    α=0.1,  Δ x=0.01,    h= Δ x^2,  τ=100h,  ρ=4×10^4

    The numerical result is shown in Figure <ref>.



 §.§ Computational details regarding surface hyperbolic mean curvature flow

Here, we introduce the parameters used in the numerical calculation of hyperbolic mean curvature flow on surfaces. 
Numerical calculations were carried out in the two and four-phase setting on the unit sphere, both with and without area preservation. 
The discussion of the results is presented in section <ref>.

    
  Result of two-phase hyperbolic MCF on the unit sphere (without area preservation) 

    We used the two-phase initial condition shown in Figure <ref> and the algorithm described in section <ref> for the numerical calculation. 
    Parameters were set as follows:
    
    α=0.1,   Δ x=0.05,    h= Δ x,  τ=5h,  ρ=0

    The numerical result is shown in Figure <ref>.
    
  Result of two-phase hyperbolic MCF on the unit sphere (with area preservation) 


    We used the two-phase initial condition shown in Figure <ref> and the algorithm described in section <ref> for the numerical calculation. 
    Parameters were set as follows:
    
    α=0.1,  Δ x=0.05,    h= Δ x,  τ=5h,  ρ=10^3

    The numerical result is shown in Figure <ref>.
    
  Result of four-phase hyperbolic MCF on the unit sphere (without area preservation) 

    We used the four-phase initial condition shown in Figure <ref> and the algorithm described in section <ref> for the numerical calculation. 
    Parameters were set as follows:
    
    α=1,  Δ x=0.01,    h=7.84Δ x/1000,  τ=200h,  ρ=0

    The numerical result is shown in Figure <ref>.
    
  Result of four-phase hyperbolic MCF on the unit sphere (with area preservation) 


    We used the four-phase initial condition shown in Figure <ref> and the algorithm described in section <ref> for the numerical calculation. 
    Parameters were set as follows:
    
    α=1,  Δ x=0.01,    h= 1.4Δ x/100,  τ=100h,  ρ=10^6

    The numerical result is shown in Figure <ref>.






















 §.§ Numerical error analysis of the area-preserving condition in the two-phase setting


We investigated how well the area enclosed by the interface in the two-phase setting on the unit sphere (shown in Figure <ref>), is preserved under the area-constrained mean curvature flow, and hyperbolic mean curvature flow. Here, we numerically solved the mean curvature flow and the hyperbolic mean curvature flow on the unit sphere with area preservation 
using the algorithm “Surface MBO for multiphase regions with the area preservation condition" introduced in section <ref>, 
and the algorithm “Surface HMBO for multiphase regions with the area preservation condition" introduced in section <ref>, respectively. 
Note that setting K=2 in the above algorithms yields an approximation method for two-phase regions.

Numerical errors were investigated as follows. 
The approximation of A^1 obtained 
in Step 3 of the algorithm “Surface MBO for multiphase regions with the area preservation" and in Step 4 of “Surface HMBO for multiphase regions with the area preservation" is denoted by V_0. 
The surface SDVF z_2^ϵ,λ obtained by executing the above algorithms for one step is used to calculate V^1_w using Eq. 
(<ref>). Here, z_2^ϵ,λ is calculated in Steps 4.f and 5.f of the above algorithms, respectively. 
The solution obtained by executing the algorithm for one step corresponds to the time τ. 
The approximate value of V^1_w is denoted by V_τ and the value of the error ERR is defined as follows:

    ERR=|V_0-V_τ|.

We investigated the response of ERR to changes in ρ (the penalty parameter for the area preservation) for the minimizing movements Eq. (<ref>) and Eq. (<ref>). 
The parameters were as follows:


    
    
  Two-phase MCF on the unit sphere (with area preservation) 

    
    α=0.05,  Δ x=0.05,    h= Δ x^2/6,  τ=100h,   10^-1≤ρ≤10^3.

    
  Two-phase hyperbolic MCF on the unit sphere (with area preservation) 

    
    α=0.1,  Δ x=0.05,    h= Δ x,  τ=5h,    10^-1≤ρ≤10^3.


The parameters used in the calculation of “two-phase MCF on the unit sphere (with area preservation)” were the same as those used in section <ref>, except for the values of ρ. 
The result obtained by evolving the system with ρ=10^3 is shown in Figure <ref>. 
The parameters used in the calculation of “two-phase hyperbolic MCF on the unit sphere (with area preservation)” were the same as those used in section <ref>, except for the value of ρ. 
The result obtained by evolving the system with ρ=10^3 is shown in Figure <ref>. 
Table <ref> shows the specific values of ρ used for the numerical error analysis. 
The results are presented in Table <ref> and Figure <ref>, where the x-axis uses a logarithmic scale. 
Discussions of the results is presented in section <ref>.





 §.§ Discussion

In this section, we will explain the results of the numerical calculations conducted in sections <ref> to <ref>. Our results are summarized in the following order: surface MCF, surface hyperbolic MCF, and numerical error analysis of the area preservation conditions in the two-phase setting. Following this, we discuss our observations regarding the motion of interfaces in the simulations, properties of the functionals used in our approximation methods, and issues related to area preservation.

Numerical results of surface MCF involving interfaces in the two-phase setting (Figure <ref>) without area preservation show that curves on the unit sphere disappear over time. Similar to the flat setting, their length decreases, and the interface becomes nearly circular before contracting to a single point. 
When the area preservation condition (Figure <ref>) is prescribed, the curvature of the interfaces tend to decrease over time, and the interface converges to a circular shape. 
The area remains approximately constant, and the curve approaches a stationary state at the terminal time. 
In the four-phase setting (Figure <ref>) without area preservation, the interface smooths itself while the length of the network decreases with time. Each interface moves to maintain the junctions where they intersect before contracting to a single point. 
When the area preservation condition is prescribed (see Figure <ref> and Figure <ref>), the curvature of the interfaces tend to decrease over time, while the junctions are maintained. However, even in the steady state, we note some slight irregularities near the junction (see Figure <ref>). Figure <ref> confirms that the curvature of the curve is still relatively large. The areas of each phase changed slightly compared to the initial state. Both observations can be attributed to various factors including: ambient mesh spacing, interpolation methods, and optimization stopping criteria.

Next, we explain the numerical results for the hyperbolic mean curvature flow on surfaces. In the two-phase setting without area-preservation (Figure <ref>), the length of the curve decreases while oscillating and  approaching a circular shape, before contracting to a point. Under the area-preservation condition (Figure <ref>), the shape of interface converges to a circle over time and the area remained approximately constant. After becoming approximately circular, the curve continued to move atop the surface of the sphere.

In the four-phase setting, the numerical results of the hyperbolic mean curvature flow (Figure <ref>) showed that, without the area-preservation condition, the interfaces oscillated with time while the total length of the network decreased. 
The interfaces evolved while maintaining junctions, and eventually the network contracted to a single point. 
When the area-preserving condition was applied, the interfaces oscillated while while preserving their areas and maintaining junctions over time. 
The area of each phase remained almost constant throughout the evolution, and eventually reached a nearly stationary state.


In the two-phase setting, surface MCF and hyperbolic MCF with area preservation both showed a tendency for the error ERR (equation (<ref>)) to decrease as ρ increases.
In Figure <ref>, ERR decreased as ρ increased over the range 10^-1≤ρ≤10, and after ρ>10, ERR increased and became almost constant and less than 0.001.
Overall, we observe that if we approximate the area-preserving MCF and hyperbolic mean curvature flow using the proposed numerical method, the value of ERR will decrease as ρ increases.
However, it is expected that increasing ρ beyond a certain value will not reduce ERR below a certain threshold.

In the numerical results (Figure <ref> and Figure <ref>) for the interfaces moving according to surface MCF with area preservation in the multiphase setting, a slight irregularity was observed at the point of stationary state. 
This observation indicates that there are points on the interface with relatively large curvature, which is contrary to the expected result. 
The reason for this stagnation, similar the original MBO, can be attributed to the fact that even though there are points with large curvature along the interface, the curvature may still be too small to resolve for a given threshold length.
Relatedly, setting an interface with sufficiently large curvatures along its initial curve tends to eliminate points with large curvature at the steady state. 
Alternatively, using a smaller spatial discretization tends to alleviate such constrains on the interface's motion. This, of course, leads to an increase in the computational time required by the method. 
In fact, all the methods developed in this study require a relatively fine spatial grid, and refining it further causes a significant increase in the required computation time. 
For example, changing the spatial grid width from 0.05 to 0.01 for the method “Surface MBO for multiphase regions with the area preservation” described in section <ref> increases the require computational time by a factor of 60. 
Consequently, improving the methods used in this study (especially their computation time) is an important future task.


Regarding the surface hyperbolic MCF, the oscillation of the interface tends to decrease with time (Figure <ref> and Figure <ref>). From the point of view of conservation of energy, this observation is unexpected. One of the reasons for this is thought to be the use the MM used in the algorithm created.
In MM, it is known that the energy of the obtained numerical solution decreases compared to the exact solution as time increases <cit.>.
Since the algorithm created reconstructs the interface based on the numerical solution obtained from the MM, it is understandable that the kinetic energy ofthe numerical solution of the interface decreases with time.


In this research, we have used minimizing movements to impose area constraints on surface-constrained multiphase interfacial motions.  Since minimizing movements require one to minimize a functional, we have also considered the influence of the functional used in this process and on the corresponding numerical results. Again, the functional used in our method for dealing with the area-constrained curvature flows in the multiphase setting (Equations (<ref>) and (<ref>)) is expressed as follows:

    ℱ_m(w)=∫_Ω_λ(F(w,w_m-1,w_m-2)+α|∇w|^2/2)dx
        +ρ∑^K_i=1|A^i-V^i_w|^2

where V^i_w is the V^i_w included in Equations (<ref>) and (<ref>), and F is expressed as follows.

    F(w,w_m-1,w_m-2)=
        |w-w_m-1|^2/2h,     Surface MBO
    |w-2w_m-1+w_m-2|^2/2h^2,     Surface HMBO

For simplicity, let us define

    J(w)     =α|∇w|^2/2
    
        P(w)     =ρ∑^K_i=1|A^i-V^i_w|^2

Then, the functional in equation (<ref>) can be expressed as follows:

    ℱ_m(w)=∫_Ω_λ(F(w,w_m-1,w_m-2)+J(w))dx
        +P(w)

The functional P determines the penalty of to the area constraints. If we set P(w)=0, then equation (<ref>) corresponds to the functional used in the absence of area preservation.
Emphasis of each area constraint is controlled through the value of ρ. 
However, if ρ is taken too large so that F and J are significantly smaller than P, then the minimizing scheme will tend to focus only on the penalty term. That is, during the process of minimizing the functional, the significance of F and J are diminished when compared to that of P.
As a result, the approximate solution may deviate from the expected result. One may observe a jagged interface, even after several minimizers and at the stationary state. 
On the other hand, if P is too small, the area constrain of each phase will not be satisfied at an acceptable level. 
Therefore, in order to approximate the motion of each interface following the mean curvature flow or the hyperbolic mean curvature flow while satisfying the area constraints, it may be necessary to adjust the ratio of the magnitudes of F, J, and P.  
Such an approach is would avoid large differences in the magnitudes of F, J, and P. 
However, it is not clear what ratio the magnitudes of F, J, and P should satisfy at present. 
We would like to return to this and related topics in a future study.


In imparting the area constraints a top the surface S, we numerically solved the constrained partial differential equations in the tubular region Ω_λ.
The functionals used in the surface-type MM (Eq. <ref>, Eq. <ref>) were designed to conserve volume in Ω_λ, where the width of Ω_λ is a constant λ (Eq. <ref>). 

Consider the case of an interface in the two phase setting. 
Let Q be the region enclosed by the interface on the surface S, let A be the area of Q, let R be the region obtained by extending Q in the normal direction of the surface S to Ω_λ, and let V be the volume of R. 
Figure <ref> shows a schematic diagram of the relationships between S, Ω_λ, R, and Q.
In this case, assuming that λ is sufficiently small, we note that we can approximate V as

    V≈2Aλ
 
In section <ref>, we investigated the numerical error of the area preservation for two-phase regions. 
We observed that the mean curvature flow and the hyperbolic mean curvature flow conserve area at higher precisions as ρ is increased. 
Since V is approximated by Eq. (<ref>) and λ is a constant, it is expected that for two-phase regions on a surface, increasing ρ will better conserve the area A surrounded by the interface on the surface. 


Remark:
The value of the width λ of Ω_λ used in Section <ref> is given by

    λ    =√(17)Δ x 
       ≈ 0.2

This is obtained by substituting p=3 and Δ x=0.05 into equation (<ref>).


The computational results in sections <ref> and <ref> showed that area preservation can be approximately satisfied on surfaces even in the multiphase case.
However, we have not yet performed a numerical error analysis to describe the relationship between the parameters of the computational algorithm, and the area preservation condition for cases other than the two-phase setting. We would like to treat this in a separate study.













§ SUMMARY

This study developed approximations methods for surface-constrained mean curvature flow and hyperbolic mean curvature flow of interfaces.
This was achieved by first creating approximation methods for a surface partial differential equations by combining the closest point method with minimizing movements. We then extended the methods to implement the conventional MBO and HMBO algorithms on surfaces.
In addition, we constructed the surface-signed distance vector field to distinguish multiphase geometries on surfaces.

Numerical error analyses of our methods were performed for the surface heat and wave equations, and convergence with respect to the spatial discretization was investigated. It was found that the numerical solution of the partial differential equation on the surface obtained by our  approximation methods converges to the exact solution.

By using the surface version of the signed distance vector field, we extended the MBO and HMBO algorithms to the surface-constrained setting. These were used to perform numerical calculations of mean curvature flow and hyperbolic mean curvature flow for two and four phase interfacial motions.

The numerical error of the prescribed area in the two-phase setting for mean curvature flow and hyperbolic mean curvature flow on surfaces was evaluated. Our results confirm that increasing the value of the penalty parameter ρ leads to higher precision in the area preservation.

Improvements to our approximation methods could be made by adjusting the energy functionals used in the MM method. Namely, it is known that, by using appropriate functionals, energy conservation can be realized <cit.>. Therefore, creating approximation methods that conserve energy and performing their numerical error analysis for equations such as the surface-wave equation is an important future task. This is expected to clarify questions about the energy dissipation of the interface in the HMBO algorithm. In addition, we would like to design generalized surface-type threshold dynamics which impart damping terms on target interfacial motion.



§ APPENDIX




 §.§ Surface HMBO and Initial Velocity for Multiphase Regions


In the surface HMBO for multiphase regions (see Section <ref> and Section <ref>) one needs to determine P_1 from the initial shape P_0 and initial velocities v_0^i of each phase. Here we describe the method.
Let Γ_ij be the interface between phase i and phase j. We define the following:

    Γ=⋃_i,jΓ_ij,   ℱ_ij=Γ-Γ_ij,   𝒥_ij=Γ_ij∩Σ_Γ_ij


    Σ_Γ_ij={C|{Γ_ij∩ C}≠∅ , C ∈ℱ_ij}

Γ represents the union of all the interfaces, ℱ_ij represents the interfaces other than Γij, and  𝒥_ij represents the endpoints of Γij. Also, Σ_Γ_ij represents the set of interfaces connected to Γ_ij.
In the multiphase setting, the surface HMCF is represented by the following nonlinear partial differential equation. For each interface Γ_ij, we have:

    d^2 Γ_ij/dt^2=-κ_ijn_ij,      t>0                                                                        
    Γ_ij(t=0)=Γ_ij^0                                                                                     
    dΓ_ij/dt(t=0)=V_ij^0n_ij^0                                                                                      
    Γ_ij(P)=σ(P),       t≥ 0,  σ∈Σ_Γ_ij,    P ∈𝒥_ij.

Here, κ_ij denotes the mean curvature of Γ_ij, n_ij represents the outward unit normal vector of Γij, Γ_ij^0 represents the interface between phases i and j at the initial time, V_ij^0 represents the initial velocity of Γ_ij^0, and n_ij^0 represents the outward unit normal vector of Γ_ij^0.
The fourth equation in Equation (<ref>) represents the continuity condition that the interface Γ_ij should satisfy.
Without the continuity condition, each interface may move independently over time, which could lead to the loss of junctions.

In the Surface HMBO for multiphase regions, after determining the initial shape P_0, the interface set Γ_ij^0 is defined for each interface. Then, for each interface, the set {Γ_ij^1} is defined as follows:

    Γ_ij^1=Γ_ij^0+hV_ij^0n_ij^0

Following this, the regions P_1 for each phase are determined from the set {Γ_ij^1}.


 §.§ Numerical error analysis of the surface MBO for two-phase regions

In this section, we present the results of an numerical error analysis for the mean curvature flow on surfaces using the algorithm introduced in section <ref>.
The analysis focuses on the motion of a circle a top the unit sphere.
As shown in Figure <ref>, let r denote the radius of the circle on the unit sphere.
A circular interface moving by the mean curvature flow on the surface of the unit sphere satisfies the following differential equation: <cit.>:

    dr/dt=r^2-1/r,   t>0 
     r(0)=r_0

where r denotes the radius of the circle at time t and r_0>0 denotes the radius of the circle at the initial time.
The exact solution of Eq. (<ref>) can be obtained as follows:

    r(t)=√(1-(1-r_0^2)exp(2t))


We investigated the numerical error's dependence on the grid spacing Δ x (the discretization of the ambient space surrounding the unit sphere). Calculations were done using the algorithm introduced in section <ref> and we employ minimizing movements to solve Eq. (<ref>). The computation implementation of the MM employs the technique introduced in section <ref>. We set r_0=2/3, the time step was h=Δ x^2/5, and the threshold length was τ=0.03. 
Let r̅(t) denote the average radius of the data points that approximate the interface at time t.

Remark:
The average radius of the data points approximating the interface is defined as follows. Assume that at time t, the interface is computed and represented by M points. For i=1,2,⋯,M, the coordinates of the i-th point are denoted as (x_i,y_i,z_i). Then, the average radius r̅(t) is calculated as follows: 

    r̅(t)=∑_i=1^M √(x_i^2+y_i^2)/M


The results obtained for Δ x = 0.05 and 0.025 are shown in Figure <ref>. The figure shows the exact solution r(t) and the average radius r̅(t) obtained our method. The points where the curves are interrupted in Figure <ref> correspond to the times that the interface could no longer be detected.

The exact solution at t=0.24 is approximately r(0.24) 0.31965745, while the average radius r̅(t) obtained from the numerical solution is 0.28291438 for Δ x=0.05 and 0.2955379 for Δ x=0.025.

Although the numerical errors are relatively small at the beginning of the computations, due to the coarsening of the numerical grid for small interfaces, the numerical errors tend to increase as time increases. We also note that the average radius of the numerical solution tends to be smaller than that of the exact solution. Regarding convergence, we indeed observe the tendency of numerical errors to decrease as Δ x becomes decreases.

When the time step size h is proportional to Δ x^2 and τ is fixed, reducing Δ x is expected to improve the accuracy of the  approximation.





 §.§ Numerical error analysis of surface HMBO for two-phase regions

In this section, show the results of a numerical error analysis for the hyperbolic mean curvature flow on a surface using the algorithm introduced in section <ref>.
Our analysis focuses on the motion of a circle on the surface of the unit sphere.
As shown in Figure <ref>, let r(t) denote the radius of a circle constrained to the unit sphere.

Since the hyperbolic mean curvature flow represents the motion in which the normal acceleration of the interface is proportional to the mean curvature, the test problem corresponding to Equation (<ref>) as follows:

    d^2r/dt^2=r^2-1/r,   t>0 
    .dr/dt|_t=0=V_0   
    
            r(0)=r_0

Here, V_0 is the initial speed of the interface, and r_0 is the radius of the circle at the initial time.

We assume that the exact solution represented by the numerical solution obtained by numerically solving Eq. (<ref>).
We compare this numerical solution with that obtained by our own computational algorithm.
Our computations use DifferentialEquations.jl <cit.> to numerically solve Eq.n (<ref>).
Similar to before, we investigated the error's dependence on the spatial discretization Δ x used in the surrounding space.
We use the algorithm introduced in section <ref> for our numerical calculations.
Minimizing movements were used to solve Eq. (<ref>) in the algorithm of section <ref>.
Our calculations set r_0=2/3, V_0=0 and the time step h was assigned to h=Δ x/10. The threshold time τ was set to τ=0.01.
As before, we define r̅(t) as the average radius of the data points that approximate the interface at time t using Eq. (<ref>).


We present the numerical results for Δ x=0.1, 0.05, 0.025 in Figure <ref>. The figure shows the numerical solution r(t) obtained by numerically solving equation (<ref>) and the average radius r̅(t) obtained using our own method. The points where the curves are interrupted in the lower right of Figure <ref> correspond to the time that the interface has dissappeared.

The value of r(t) at t=0.6 is approximately r(0.6) 0.50635371, while the average radius r̅(t) obtained from the numerical solution is 0.293178268 for Δ x=0.1, 0.380734809 for Δ x=0.05, and 0.474013006 for Δ x=0.025.

Except for Δ x=0.1, the average radius continues to decrease over time, and the interface can no longer be detected. However, for Δ x=0.1, the average radius starts to increase at some point. This is because the interface that initially shrinks inward eventually becomes a point and starts to expand outward. After the interface has contracted, the subsequent expansion is an interesting feather of the hyperbolic mean curvature flow. A detailed analysis of this phenomenon is a future research topic.

In all cases, the numerical error is relatively small at the beginning of the calculations, but tends to increases over time. Of course, decreasing Δ x tends to decrease the numerical error.  Setting the time step h proportional to Δ x and fixing τ should lead to improved accuracy for decreasing Δ x.




 §.§ Implementation methods

Here we explain a few important details about the implementation of the numerical algorithms introduced in section <ref> and section <ref>. The numerical algorithms described above involve the calculation of functional values and integrals, which require discretizations a computer.

Similar to section <ref> and section <ref>, we consider a smooth surface S in a three-dimensional Euclidean space and discretize a sufficiently large domain Ω_λ that covers S. As before, we refer to the discretized space as Ω_λ^D. As a result of the discretization, Ω_λ^D is divided into N_x points in the x-direction, N_y points in the y-direction, and N_z points in the z-direction, as in section <ref>. The interval between divisions is assumed to be equal in all three directions and is denoted by Δ x.

    
  Approximation of functional values

    In our method, it is necessary to compute the value of the functionals included in the first term of Eq. (<ref>) and Eq.  (<ref>). We will explain the approximation using Eq. (<ref>) as an example. The first term of Eq. (<ref>) is expressed by
    
    ℱ_m(w)=∫_Ω_λ(|w-w_m-1|^2/2h+α|∇w|^2/2)dx.

    In equation (<ref>), function w: Ω_λ→ℝ^K-1(K:number of phases) is a vector-valued and so the functional value can be computed as follows:
    
    ℱ_m(w)=∑_i=1^K-1∫_Ω_λ(|w^i-w_m-1^i|^2/2h+α|∇ w^i|^2/2)dx

    Here, w^i is the ith component of w. The integrals included in equation (<ref>) are approximated using the same technique as in equation (<ref>). For equation (<ref>), we used the same method as in the calculation of equation (<ref>).

    
  Approximation of V^i_w

 The value of V^i_w appearing in Eq. (<ref>) and Eq. (<ref>) represents the volume of the extension of phase i within Ω_λ. It can be approximated using H^ϵ and ϕ^i_w defined in Eq. (<ref>) as follows:
    
    V^i_w    =∫_Ω_λH^ϵ(ϕ^i_w)dx,                                                       
       ≈Δ x^3 ∑_x_i,j,k∈Ω_λ^DH^ϵ(ϕ^i_w(x_i,j,k))

    
  Calculation of w^S 

        The methods developed here include a step where one must extract the values of w_M at the points of S:
        
    w^S(x)=w_M(x),  x∈ S∩Ω_λ

        In the numerical calculations, Ω_λ is discretized. Therefore, an interpolation on Ω_λ^D (the discretized grid points of Ω_λ) is required in order to obtain the values on the surface S.
        In this study, we have used third-order polynomial interpolations.

    
  Creation of geodesic distance functions 

        When constructing the surface SDVF, signed distance functions on the surface are required. These calculation are not very simple, and we use a method based on the Fast Marching Method <cit.> to construct the signed distance vector field on the surface.





plain

