A rapid numerical method for the Mullins-Sekerka flow with application to contact angle problems

The Mullins-Sekerka problem is numerically solved in $\mathbb{R}^2$ with the aid of the charge simulation method. This is an expansion of the numerical scheme by which Sakakibara and Yazaki computed the Hele-Shaw flow. We investigate a sufficient condition for the number of collocation points to ensure that the length of the generated approximate polygonal curves gradually decreases. We propose a new benchmark function for the Mullins-Sekerka flow to confirm that the scheme works well. Moreover, by changing the fundamental solutions of the charge simulation method, we are successful to establish a numerical scheme that can be used to treat the Mullins-Sekerka problem with the contact angle condition.


Introduction
The Mullins-Sekerka problem is a quasi-stationary Stefan problem with surface tension.Its solution describes the time evolution of the interface separating two-phases filled with different materials.At a time t, one material is filled in a smooth bounded region Ω t and another material is filled in the outer region R 2 \Ω t .The boundary Γ t = ∂Ω t represents the interface.The temperature of the material is denoted by u = u(t, x) for x ∈ R 2 \Γ t and t ≥ 0.Then, the Mullins-Sekerka problem asks us to find (u, Γ t ) that satisfies where n denotes the normal vector to Γ t outgoing from Ω t , and κ denotes the curvature of Γ t .The normal velocity of Γ t is denoted by V .More precisely, the motion of Γ t is governed by d dt X(t, s) = V (t, X(t, s))n(t, X(t, s)) where the interface Γ t is parameterized as X(t, s) for s ∈ [0, 2π].V describes the speed of Γ t in the direction to n.Let [ϕ] be the jump in the normal direction of a quantity ϕ across Γ t ; namely, [ϕ] (x) := lim ε→0 {ϕ(x − εn) − ϕ(x + εn)} for each x ∈ Γ t .
The Mullins-Sekerka equation is a limit of a phase-field model, which is called the Cahn-Hilliard equation where two phases are separated by a transition layer instead of a sharpe interface.The Cahn-Hilliard equation is important to understand spinodal decomposition, which explains a phenomenon that compounded solutes and solids are stable at high temperatures but become unstable at low temperatures and eventually separated by sharp interfaces.The convergence of the Cahn-Hilliard equation to the Mullins-Sekerka equation was first formally shown by Pego [22].Stoth [27] gave a rigorous proof of the convergence in the case where the domain under consideration is a ball in R 3 , and the initial data and the boundary values are all radially symmetric.In a general dimension, Alikakos et al. [2] proved that a family of smooth solutions to the Cahn-Hilliard equation tends to a smooth solution to the Mullins-Sekerka equation provided that the latter exists.
The Mullins-Sekerka problem has been well studied from an analytical point of view.Chen et al. [7] proved the existence of a classical solution to the Mullins-Sekerka problem local in time in the two dimensional case, and Escher and Simonett [9] proved it in general dimensional cases.In the literature on weak solutions, Luckhaus and Sturzenhecker [17] proposed a weak notion of the solution to the Mullins-Sekerka problem and gave its global time existence result whenever the sum of the surface area measure does not change discontinuously over time.Under the same assumption, Bronsard et al. [6] established a weak solution of a multiphase Mullins-Sekerka problem with a triple junction and showed its existence.As indicated in our experiment, this setting cannot be applied when multi particles exist and are very close to each other (see Section 4.6).Röger [25] was successful to remove this assumption in terms of geometric measure theory, although his result could be applied only to the two-phase case and excluded the multi-phase case and contact angle case.Recently, Hensel and Stinson [13] proposed a varifold solution to the Mullins-Sekerka problem based on the energy dissipation property and included a fixed contact angle condition.Julin et al. [14] revealed an asymptotic behaviour of the Mullins-Sekerka flow in the torus T 2 ⊂ R 2 .They proved that a flat flow solution to the Mullins-Sekerka problem exponentially converges to the finite union of disks whenever the perimeter of the initial data is less than 2 (see Theorem 1.3 [14]).
There are several works that treated the Hele-Shaw problem or the Mullins-Sekerka problem numerically.Though a typical Mullins-Sekerka problem is considered in a smooth bounded domain, Bates et al. [4] considered the problem (1) in R 2 and translated the original problem (1) into a corresponding boundary integral equation.Eventually, they did not have to care about the boundary of the container; they split an initial curve into several segments and regarded as a part of circles.In this way, they could calculate a discrete version of the curvatures and derive a linear system of equations whose unknown variables are normal velocities at each vertex by means of the Gibbs-Thomson law (this is the second condition of (1)).Recall that our unknown variables in the linear system are coefficients of fundamental solutions in contrast to their scheme (see (10) and ( 12)).Moreover, they adopted the semi-implicit method to stabilize their scheme with a small time step when solving the linear system.For other studies using the boundary integral method, we refer the reader to [28], [19], and [5].Barrett et al. [3] proposed a parametric finite element scheme for the Stefan problem with surface tension and proved the well-posedness and stability of their scheme.In the course of the discussion, they also revealed that the scheme applies to the Mullins-Sekerka problem.
Feng and Prohl [10] showed that their numerical scheme, the so-called the fully discrete mixed finite element scheme to construct discrete solutions to the Cahn-Hilliard equation, tends to the solution of the Mullins-Sekerka equation, provided that a global-in-time classical solution exists.Recently, Nürnberg [21] introduced a front tracking method for the Mullins-Sekerka flow which is based on the finite element method.He proved its unconditional stability and volume conservation law of the scheme.Its accuracy was confirmed in terms of two concentric circles.
The purpose of this paper is to propose a discrete scheme to solve the Mullins-Sekerka problem numerically, revealing that our scheme has some desiable properties.To this end, we follow the scheme proposed by Sakakibara and Yazaki [26].Moreover, we confirm that the proposed scheme possesses curve-shortening property (CS) and area-preserving (AP) properties.These facts are predictable because the original scheme for the Hele-Shaw problem also has such properties.However, we focus on CS and derive a discrete variant of the estimation related to the length of the curve.This outcome is reported in Corollary 1.
At this stage, we provide a brief explanation of our scheme.Suppose that a smooth curve Γ t is given for some t > 0.Then, Γ t is approximated by an N polygon Γ N t where N ≥ 3 is a positive integer.The interior domain of Γ N t is designated as Ω N t .The vertices of Γ N t are denoted by X i (1 ≤ i ≤ N ).For convenience, we adopt a periodic rule for the indices of X i , such as In addition, normal vectors N i = N i (t) and tangential vectors T i = T i (t) at X i are suitably defined.Among the standard numerical methods, we adopt the charge simulation method (CSM) that was originally developed to approximate a solution to the Laplace equation in a bounded domain with Dirichlet boundary conditions.CSM is a variant of the method of fundamental solutions (MFS), in which approximate solutions are expressed as a linear combination of fundamental solutions to partial differential equations under consideration.CSM requires choosing proper charge points y + 1 , • • • , y + N from R 2 \Ω N t and collocation points X 1 , • • • , X N on Γ t and expresses an approximate solution as follows: where E is the fundamental solution of the Laplace equation, that is This is the basic idea of CSM.For more detail about CMS, see Katsurada and Okamoto [16].The above conventional scheme was modified by Murota [20] to make the scheme possess invariance properties that original continuous problems also have.Concretely, he alternatively used the following combination.
In this scheme, we have to add an equality to a linear system because we have one more value Q + 0 to find.For instance, Murota assumed that the sum of Q + i s equals zero.In solving the Hele-Shaw problem numerically, Sakakibara and Yazaki [26] improved Murota's invariant scheme to make this additional assumption more natural.They defined dummy singular points z + i (1 ≤ i ≤ N ) and replaced the combination of fundamental solutions by It can be observed that the above function U + is invariant under translation and scaling without any additional assumptions.Hence, it is possible to impose an area-preserving requirement that seems more natural than the zero-average assumption.Taking singular points y − i and z − i from Ω N t , we also have a function U − being harmonic in R 2 \Ω N t that satisfies the Dirichlet boundary condition.We should solve an external potential problem to find such a U − .However, it is impossible by either the finite difference method or the finite element method owing to the unboundedness of the domain where the problem is considered.As imposed in the fourth equality of (1), each point x on Γ N t is required to move at the speed that is equal to the jump of the normal derivatives of U + and U − across Γ N t .Once we obtain such U + and U − , the direct differentiation of U + and U − yields the representative normal velocity Tangential velocity W i and its vector T i are required to stabilize the scheme, although they have no effect on the geometry of the curve (see Proposition 2.4 [8] for instance).Normal velocity V i and its vector N i definitely control the motion of a curve.Finally, we discretize the time variable as t = n∆t(0 ≤ n ≤ N ) and rearrange the evolution equation (5) as follows: A particular novelty of this study is treating a boundary contact case of the Mullins-Sekerka problem in the half plane R 2 + .To this end, we replace the fundamental solutions of the combination by the Green function on R 2 + .Since the curves under consideration are open, we modify the structure of the proposed scheme.Well-posedness of the Mullins-Sekerka problem with ninety contact angle condition was established by Abels et al. [1].After that, Garcke and Rauchecker [12] show stability and instability results of stationary solutions to the linearized problem that is either flat or a part of a circle.
The reminder of this paper is organized as follows.In Section 2, we rigorously state how to implement the proposed scheme.Section 3 is devoted to list our main results without the proofs.In Section 4, we give several examples of implementation of the scheme.Moreover, the accuracy of our scheme is confirmed in terms of an annulus-like domain which consists of three concentric circles and a continuous function being harmonic except on the circles.To our best knowledge, this is a new feature in the literature of the Mullins-Sekerka problem as a benckmark function that can describe the two-phase motion.We shall extend the scheme to the boundary contact cases in Section 5.In Section 7, we collect all proofs of Theorems and Propositions whose justification has been postponed.

Numerical scheme
In this section, we present a concrete procedure to construct approximate polygons.This is a natural extension of the scheme proposed in [26].

Polygonal approximation of the interface
Let N ≥ 3 be a positive number, and Ω N be an N polygonal domain in R 2 with vertices This periodic rule is always applied unless otherwise stated explicitly.The symbol [X i−1 , X i ] denotes the line segment connecting X i−1 and X i .This is called the edge in the sequel.Then, the boundary Γ N := ∂Ω N readily designs an N polygon and is expressed as For each 1 ≤ i ≤ N , we define Here we have used the notation (a, b) ⊥ := (b, −a).Moreover, the midpoint of the edge The outer angle ϕ i of Γ N at each vertex X i is expressed as follows: where the function sgn is defined by for each a, b ∈ R 2 .By using ϕ i , we set Then, we define the discrete curvature of Γ N at X * i as follows: Moreover, we define the normal vectors N i and the tangential vectors T i at X i in terms of the normal vectors n i and the tangential vectors t i of the edge [X i−1 , X i ] as follows: See the figure below to take a look at our setting.

Approximation of the normal velocity
Solve the internal problem.In this step, we solve the Laplace equation with a Dirichlet boundary condition given by κ i .To this end, we must choose singular points Then, suppose that the solution of the Dirichlet boundary problem is of the form for some clearly harmonic in Ω N due to its structure.To determine these values, we shall solve the following linear system of equations: where Solve the external problem.As in the previous step, we take singular points Then, the solution to the exterior Dirichlet boundary problem is expressed as where Derive the normal velocity.Once Q + i and Q − i are determined, we can define a representative normal velocity at X i .By setting v we define the representative normal velocities by Derive the tangential velocity.To reduce instability of the scheme, we consider the tangential vector that does not affect the shape of the polygon.To this end, we adopt a uniformly distribute method (UDM).Let us only state formulae to obtain the tangential velocities.We refer the readers to Sakakibara and Yazaki [26] for derivation of the scheme.
Tangential velocities W i (1 ≤ i ≤ N ) are defined as follows: Here we have set Ψ i := i j=1 ψ j for each 1 ≤ i ≤ N where ψ 1 := 0 and

Time evolution of the polygonal interface
We are now in the position to state our numerical scheme.Given an initial curve Γ 0 = ∂Ω 0 , we approximate it by an N polygonal curve where The domain surrounded by Γ N 0 is denoted by Ω N 0 .For t > 0, suppose that an N polygon Γ N t is an approximation of Γ t to be found.The vertices of Γ N t are denoted by are defined as above.Following the continuous version of the evolution law, we consider Finally, the time variable is discretized as t n = n∆t(n = 0, 1, • • • ) with a given time step ∆t > 0 and the forward Euler method is applied to (15).Using the notations we provide the fully discrete scheme below.
with the initial condition 16) is denoted by Γ N,n ∆t .The interior domain of Γ N,n ∆t is denoted by Ω N,n ∆t .

Properties of the scheme
As mentioned in Introduction, the proposed numerical scheme is expected to have the areapreserving and the curve-shortening properties since the original continuous problem has such properties.Here, we assume that Γ t is a simply closed curve that separates R 2 into a bounded domain Ω t and an unbounded domain R 2 \Ω t .Suppose that (u, Γ t ) is a solution to the Mullins-Sekerka problem (1).Let us distinguish u inside Ω t by writing u i , whereas u outside Ω t is presented as u e .Then, the area A t of Ω t never changes.Indeed, we compute Here, R > 0 is arbitrarily taken so large that Ω t ⊂⊂ B(0, R).Thus, d dt A t = 0.Moreover, we see that the length L t of Γ t does not increase in time.Indeed, Lemma 3).Therefore, by letting R → ∞, we see that the right-hand side of the above identity is not positive.In this section, we will observe that these fine properties are valid even in the discrete version of the solutions derived by MFS.Let us prepare several useful formulae to proceed argument.All proofs of the following results are postponed and stated in Section 7.
Let L = L(t) and A = A(t) be the length and the area of a polygon evolving in time, subject to (15), respectively.Namely, The following formulae are used to derive the tangential velocities of each vertex of the polygon.
Proposition 1.Let L and A be defined by ( 17) and (18), respectively.Then, the following formulae are valid: where Proposition 2 (Uniform boundedness of charges).Set . where Theorem 1 (Curve shortening property).Assume that L is defined by (17).Then, it holds that where C N denotes a positive constant depending on N .
Remark 1.Despite the appearance of the term O(N 2 log N ), the right-hand side of (21) becomes nonpositive for N large enough.Indeed, since Hence, a direct calculation gives: Here, we have recalled that To get a desired property, we fix N so large that Thanks to Proposition 1, we can expect that our scheme decreases the length of the polygon step-by-step.This expectation turns out to be true subsequently (see Corollary 1).
Remark 2. If the domain that we approximate by a polygon is a circle, then the term O(N 2 log N ) appearing in Theorem 1 is expected to be replaced by O( α N N ) for some 0 < α < 1.Let us write the approximate solutions of ( 9) and (11) as U + (N ) and U − (N ) , respectively.Moreover, we replace the domain in which we solve the external problem with an annulus.Then, from the results of Katsurada (Theorem 2.2 and Theorem 4.1 [15]) and Murota (Theorem 2.4 [20]), we can expect that U + (N ) and U − (N ) converge exponentially to the exact solutions U + and U − as N → ∞.We cannot apply these results directly because their selection of collocation points is slightly different from ours.Thus, let us assume that the results are available for our study.Since Γ N is a regular polygon, we see that κ i = κ for every 1 ≤ i ≤ N .Then, the AP property, which is the second constraints of (10) and (12), guarantees that Since U + (N ) and U − (N ) are radially symmetric due to their construction, • n i has the same sign for all 1 ≤ i ≤ N .Hence, we can assume that this value is nonnegative without loss of generality.Then, we can estimate as follows: Here, we note that Γ N i ⊂ Ω for each 1 ≤ i ≤ N to derive the last inequality.Moreover, there are neither singular points nor dummy singular points in the domain Ω N i sandwiched between Γ N i and ∂Ω for every 1 ≤ i ≤ N for N ∈ N large enough.To see this, recall that This obviously implies that y ± i , z ± i / ∈ Ω N i for sufficiently large N ∈ N. Hence, we obtain From this observation, we can calculate as follows: Admitting Katsurada's result(See Remark4.1 [15]), the integration of the right hand side of the above inequality tends to ∂Ω ∂U ∂n dS as N → ∞ and this quantity equals zero if U + and U − are exact solutions to (1).In this way, we can predict that it is not necessary to take N so large to ensure the CS property holds.
Theorem 2 (Area preserving property).Assume that the vertices X i are uniformly distributed, namely, r i = r i+1 for each 1 ≤ i ≤ N and these are differentiable with respect to the time variable.Then, the area A surrounded by Γ is theoretically preserved for all t; that is Ȧ = 0, provided that we adopt the UDM.
Theorem 1 indicates that our scheme has the CS property in some sense once we assume the differentiability of X i with respect to the time variable.Our next trial is to derive a timediscrete version of Theorem 1.For simplicity of notation, let us use the symbol v n i instead of v + i + v − i .Moreover, let r n i and κ n i denote the corresponding values of Γ N,n ∆t defined in ( 7) and ( 8), respectively.
Theorem 3 (Discrete version of the CS property).For each n ≥ 0, let (16).Then, it holds that for each n ≥ 0, ) We deduce from Theorem 1 and Theorem 3 that the highest degree with respect to N of negative terms on the right-hand side of ( 22) equals max {3 − α, 9  2 − 2α}.On the other hand, one of the positive terms equals 5−2α.Thus, by solving the inequality max {3 − α, 9  2 − 2α} > 5 − 2α, we see that the right-hand side of ( 22) should be negative for N ∈ N large enough whenever α > 2. Consequently, we obtain the following corollary.
Corollary 1.For any α > 2, suppose that ∆t := N −α .Then, for sufficiently large N , L n+1 ≤ L n is valid for each n ≥ 0. Remark 3. In [26], several numerical experiments were carried out with N = 100 and ∆t := 0.1N −2 .This choice seems reasonable since ∆t = N − 5 2 , and this setting is consistent with the result of Corollary 1.

Numerical examples 4.1 Pi curve
We borrowed the coordinates of the pi curve as the initial data from https://ja.wolframalpha.com/, as [26].See Figure 2 for a numerical result.To confirm the uniform boundedness of

Star
The star-shaped curve is close to a circle.As expected, this shape is deformed into a circle by our scheme.See Figure 5 for a numerical result.

Tube
Mayer [18] showed that the one-sided Mullins-Sekerka flow (Hele-Shaw flow) does not necessarily preserve the convexity of the initial curve.This is a unique feature of the Mullins-Sekerka flow.This fact was also numerically observed by Bates et al. [4] provided that an initial curve is a tube with two circular end caps, even in the 2-phase Mullins-Sekerka flow.We shall confirm that our scheme yields the same result as well.We set N = 50, ∆t = 0.1N −2 .In addition, assume that the length and thickness of the tube are 8.0 and 1.0, respectively.

Accuracy of the scheme
As mentioned in [4], two concentric circles and a function u in R 2 are available to comfirm the accuracy of the numerical scheme.Let us recall the definition of u as follows: For 0 Then, u is readily harmonic in Ω := {R 1 < |x| < R 2 } and satisfies the Gibbs-Thomson law on ∂B R 1 (0) ∪ ∂B R 2 (0).However, we cannot adopt this function as a benchmark for our scheme because this model can only describe the Hele-Shaw flow, which is one-sided.Indeed, though the normal derivative of u jumps across ∂Ω, u is constant in R 2 \Ω so that the jump coincides with −∇u • ν and this case is something special.Hence, we should prepare a new example of u whose normal derivative does not vanish in both side across the boundary.To this end, we alternatively consider the following function u and the domains Ω i and Ω e .Let 0 < R 1 < R 2 < R 3 < R and we define Though the computation will be much more complicated since Ω i is not connected similarly to Ω e , we see that the normal derivative of u from both sides of across ∂B R 2 (0) ⊂ ∂Ω i does not vanish.Observe that the function u defined by (24) rigorously satisfies the equalities (??) and the radii R 1 , R 2 and R 3 evolve, subject to the following system of ordinary differential equations: , In the sequel, let us confirm that this u can be approximated by the proposed scheme.First, to estimate the error between the results of the proposed scheme and the rigorous solution, we need to derive a numerical solution of (25).To this end, we adopt the fourth-order Runge-Kutta method.We draw three concentric circles: C 1 , C 2 , and C 3 with radii of R 1 , R 2 , and R 3 , respectively.Then, we plot N 1 , N 2 , and , respectively.Due to a technical reason, we assume that N 2 = N 3 .For singular points y ± k,j and dummy singular points z ± k,j , we set: The quantity d that indicates the distance of each charge point y ± k,j from the edge [X k,j−1 , X k,j ] is set to 1 √ N universally.However, we adjust this quantity by multiplying it with a scaling parameter.This scaling will stabilize the position of each charge point, even if the length of the edge [X k,j−1 , X k,j ] becomes quite small.Note that this modification is unnecessary if the domain approximated by a polygon is connected because the shape of the domain converges to a circle in time and the length of the edge is bounded below by a positive constant.In this experiment, the radii are chosen as R 1 := 1.0, R 2 := 5.0, and R 3 := 10.0.Solving (25) numerically for these parameters shows that the circle C 1 vanishes at approximately t = 0.56, and a topological change occurs.However, if N is too small, then the error between the approximate solution and the rigorous one is not neglectable, and we cannot proceed the scheme until t = 0.56.This difficulty may stem from observation that the approximate solution U + should be constant inside C 1 although this cannot be expected due to the structure of the approximate solution.To avoid this problem, we calculate R 1 (t) and R 3 (t) in terms of the system (25), whereas we approximate R 2 (t) by the proposed scheme.To be more accurate, we apply the fourth-order Runge-Kutta method to obtain R 1 (t), R 2 (t), and R 3 (t) for each time-step.To reduce the calculation load, we skip the step of the UDM.In other words, we always assume that the tangential velocity equals zero.Being optimistic, it is worth to consider the normal velocity because the concentric circles are symmetric with respect to the origin.In this way, we can proceed with the scheme and measure the error of R 2 until a topological change occurs.
This r < 0 is often called the experimental order of convergence (EOC).The quantity (26) was utilized in Eq. (3.1) [24] to present an experimental performance analysis of a PDE-solving system called ELLPACK.Let us define Err(N ) in our setting.We implement the time-discrete evolution of polygons until k∆t > 0.5 holds.Let M be the number of implementations and set T := M ∆t where ∆t := 0.1N −2 as usual.Then, we define Err(N ) by where R 2 (k) := |X k 2,1 | and R2 (t) denotes the rigorous solution of ( 25).This quantity is nothing but the mean value of the relative error between the radii of the approximate solution and the rigorous one.In terms of this error, let us summarize the result of the experiment.From this result, we can expect the convergence rate of the present scheme to be at least

Disappearance of particles
Bates et al. [4] treated a case where there are several particles in R 2 .According to [4], if the particles are disjoint circles, then the largest circle grows, and smaller particles shrink and eventually disappear.See Figure 5 [11] for a 3D case.The proposed scheme can also be applied to such cases.In order to observe this phenomenon, we prepare four circles in R 2 imitated by 20 regular polygons, namely, we set N = 20 * 4 = 80.While evolving, the circles will be removed from the target of the numerical calculation when their area becomes smaller than a prescribed setting.(g) step7.

Coalescence of particles
Setting the initial datum as two ovals and placing them closely, we can observe the coalescence of particles as shown in Figure 8 [4].This type of phenomenon was rigorously formulated by Röger [25] with the aid of the notion of varifolds.He was successful to establish the existence of a weak solution to the Mullins-Sekerka problem without the assumption that the loss of area never happens.However, there is no uniqueness result of the weak solution.In this experiment, we check the distance between the collocation points of the ovals step-by-step.If the distance becomes lower than a prescribed value, then we remove the collocation points and their neighboring collocation points from the ovals.Second, we artificially connect the ovals and regard them as one polygon.Each oval has 42 collocation points, that is, N = 42 * 2 = 84.For the numerical result, see Figure 11.

Dumbbell
When the initial data is a dumbbell, which is the sum set of two circles connected by a narrow bar, a pinch off phenomenon may occur.For example, let r, l and b be the circle's radius ,the bar's length, and the thickness, respectively, with r := 30, l := 5, b := 0.25.Then, the thickness of the dumbbell decreases in time and eventually pinches off.In this case, we remove the vertices of the polygon when the bar thickness becomes smaller than the prescribed value.See Figure 12.

Extension to boundary contact cases
Let us consider the Mullins-Sekerka flow with a homogeneous Neumann boundary condition in the half plane Formally, the problem we are concerned with is as follows: Here, Γ t is an open curve in R 2 + whose endpoints are bonded on ∂R 2 + for every t ≥ 0. Alternative fundamental solutions.To solve this problem numerically, we have to use another approximate solution different from ( 9) and (11).Let us explain how to construct approximate solutions to (28) step-by-step.Similar to the case where Γ t does not touch the boundary, suppose that Γ N t is an N polygon whose vertices are X i (1 ≤ i ≤ N ).The indices are numbered counterclockwise so that X 1 and X N are the endpoints of Γ N t .X 1 and X N are allowed to move in time only on ∂R 2 + , that is the x 1 -axis.Singular points i .Now, we define the approximate solutions U + and U − as follows: where we have used the notation y := (y , −y n ) when y = (y , y n ) ∈ R N −1 × R. Observe that the functions U + and U − are symmetric across the x 1 -axis and fulfill the homogeneous Neumann boundary conditions owing to their structures.The real numbers Q ± i (0 ≤ i ≤ N ) are determined with the aid of the Dirichlet boundary condition and the AP property.We should note that κ 2 and κ N are affected by the contact angle between Γ N t and ∂R 2 + which is not originally geometry of Γ N t .Thus, we replace the definition of κ 2 and κ N with Modified normal velocity at the end points.We allowed the endpoints X 1 and X N to move in time but bonded them to the x 1 -axis.Thus, it is necessary to modify the normal velocity at the points what we have defined in (13).Suppose that the contact angle between Γ N t and the x 1 -axis at X 1 equals θ.When the edge [X 1 , X 2 ] moves at the speed V in the direction of n 2 , the movement of the endpoint X 1 should equal V sin θ .Taking this observation into account, let us define the new normal velocity vectors at the endpoints as follows: Here, we have used the notation that e 1 := (1, 0) ∈ R 2 .
Modified UDM.We change the way to implement the method of uniform distribution for the vertices X i because X 1 and X N never move along the tangential vector of Γ t , and the edge [X N , X 1 ] is not included in Γ N t .Tangential velocity W i (1 ≤ i ≤ N ) are required to satisfy the following linear system: By these changes from the case where Γ t is a simply closed curve, we can determine the velocity vectors Ẋi for each 1 ≤ i ≤ N .We exhibit the numerical result of a simulation for the Mullins-Sekerka flow by means of the way proposed above.An initial open curve is an L-shaped and contacts the x 1 -axis at 90°angle at two endpoints.The collocation points are equidistantly placed and at least at the corners of the character L. Though the contact angles become different from 90°as soon as we start the program, they eventually converge to 90°, and the shape of the curve tends to a semicircle in finite time.Remark 4. As mentioned in Section 1, local well-posedness of the problem (28) has been shown in [1] whenever Γ t is perpendicular to the boundary ∂R 2 + , that is to say, ∠(Γ t , ∂R 2 + ) = 90°.We should note that our scheme yields contact angles that are not equal to 90°, and this result is different from that of [1].We predict that the L-shaped curve cannot be presented by a height function over a fixed reference surface Σ, which also intersects the boundary of the container at 90°.Thus, we could not detect 90°numerically because of this extremely short time existence.

Acknowledgements
The author has learned how to express curves on the plane and background of the method of fundamental solutions from Yazaki's monograph.The author is grateful to Professor Norikazu Saito for his suggestion to apply the charge simulation method to this problem.He reviewed the manuscript and gave me a lot of helpful comments to improve the paper.Professor Yoshikazu Giga encouraged the author to consider the case where the heat diffusion coefficients are different and examine it.

Sequence of the proofs
In this section, we list the proofs of the results presented in Section 3 that have been postponed.
Lemma 1.Let the initial curve Γ 0 be a Jordan curve in R 2 with at least C 2 regularity.Assume that dummy singular points z i are taken as z i = X * i + d 2 n i and collocation points X i are equidistantly placed, that is to say r i = L N .Then, for each 1 ≤ i, j ≤ N , More strongly, it holds that where S(N ) denotes the set of all bijections between {1, • • • , N }.
Proof.Since Γ does not cross itself, we can choose an N so large that |X * i −y j | ≥ |X * i −y i | and |X * i − z j | ≥ |X * i −z i hold for every 1 ≤ i, j ≤ N .Indeed, viewing the triangle X * i y i+1 X * i+1 , we deduce from the cosine theorem that Moreover, we also take a large N ∈ N to obtain the last equality if necessary.By a similar argument, we obtain Fix 1 ≤ i, j ≤ N with i < j.Applying the Taylor expansion to g(t for some s ∈ (0, 1).It readily follows that For the second term on the right-hand side of (31), we compute: Let us estimate the third term.To this end, note that holds due to an elementary geometric observation.Suppose that Then, the third term is bounded by 1 2 On the other hand, if |X * i − y j + (y j − z j )s| ≥ |X * i − z j |, then the third term on the right-hand side of (31) is bounded by 1 2 By these observation, we see that Next, we consider the function x → x |x| 2 and again utilize the Taylor expansion to obtain Thus, we have Finally, let us estimate H j .By the definition, Before proving the main theorems, we prepare some lemmas.
Lemma 2. For each 1 ≤ i ≤ N , let S i be defined by holds where Proof.Direct differentiation shows that Here, the symbol ∇ 2 denotes the Hessian.The first two terms of (33) have been estimated in Lemma 1 and the proof of Theorem 2 at [26].Since the structure of U − is the same as U + , the same argument works well to estimate the third and fourth terms.We omit the proof.Lemma 3. Let f : R 2 → R be partially differentiable and satisfy ∇f Proof.From the assumption for ∇f , there exists R > 0 such that for some C > 0, it holds Fix any x ∈ R 2 \B R and set x n := (R + n) x |x| for each n ∈ N ∪ {0}.We construct a sequence { xn } n by apply Taylor's expansion around x n−1 to find xn ∈ (x n−1 , x n ) for which valid from the construction of {x n } n and {x n } n .Moreover, we again use the Taylor expansion around x = x n 0 to obtain x ∈ (x n 0 , x) such that f (x) = f (x n 0 ) + ∇f (x) • (x − x n 0 ).Then, we compute Note that the right-hand side of the above inequality is finite and independent of the choice of x ∈ R 2 \B R .Therefore, we can conclude that f is bounded in R 2 .Proof.The first assertion is straightforward from the following calculation: The second assertion follows immediately by setting α := 2 in Lemma 3.
Remark 5.An argument similar to (35) can be found in [4] in which the boundary integral method was applied to construct an approximate solution to (1).To show decay of the derivative of solutions as |x| → ∞, they used the area-preserving restriction, namely the mean value of the jump of normal derivatives across the phase interface equals zero.Due to the definition (11), we do not have to use the second equality of (12).where α > 0 is a positive constant which is independent of N due to Lemma 1.We should note that S(N ) = N ! to get the last order.By contrast, the cofactors of A N can be estimated as follows: For 2 ≤ i, j ≤ N + 1, These estimations tell us that while the absolute sum of the first columns has the order O(N !N ), that of second columns and beyond have the order O(N !N 2 ).This implies nothing but A N 1 = O(N !N 3 2 ).Therefore, we complete the proof of the first value by Recalling the relationship A N Q N = K N , we have Q N = A −1 K N .Since Γ is smooth, the curvature of Γ has a global maximum so that κ i = O(1) as N → ∞ for each 1 ≤ i ≤ N .Then, we compute for each 0 ≤ i ≤ N , = O(1).
Proof of Theorem 2. As shown in Proposition 1, the derivative of A with respect to the time variable is expressed as Owing to the UDM, it immediately follows that errA = 0. On the other hand, the second equations of ( 10) and ( 12) require that N i=1 v + i r i = N i=1 v − i r i = 0. Therefore, we have Ȧ = 0. Lemma 5.For each 1 ≤ i ≤ N , T i = cos i t i − sin i n i = cos i t i+1 + sin i n i+1 , N i = cos i n i + sin i t i = cos i n i+1 − sin i t i+1 .
Proof.This is easily seen if one notes that T i is derived by either rotating t i by the angle ϕ i 2 counterclockwise or rotating t i+1 with the angle ϕ i 2 clockwise.The second formula immediately follows by rotating the first formula by the angle π 2 .
We combine the formula (5) with Lemma 5 to get the following formulae:

Figure 1 :
Figure 1: An N polygon mimicking a curve.

Figure 3 :
Figure 3: Time evolution of length, area and curvatures for pi-curve.

and 1 N
N j=1 |Q − j | for a large N , we extract the first step of the scheme for each N and draw their values as a graph.We infer from Figure4that all of them may be dominated by constants and depend only on the initial curve.

Figure 7 :
Figure 7: 3 concentric circles and the position of charge points.

Figure 8 :
Figure 8: Relative error of R 2 (t) for each N .

Figure 14 :
Figure 14: Evolution of L-shaped open curve.

Proof of Theorem 1 .(− 1 )
Take an N ∈ N so large that A N is regular.To confirm that the first value is finite, let us estimate det A N and A N separately.The cofactor expansion of A N givesdet A N = tσ H σ(i) G 1,σ(1) • • • G i−1,σ(i−1) G i+1,σ(i+1) • • • G N,σ(N )

Table 1 :
EOC with respect to N .