74
Minimizing movement approach without using
distance function for evolving spirals by the
crystalline curvature with driving force
Takeshi Ohtsuka

Division of Pure and Applied Sciences,
Graduate School of Science and Technology, Gunma University
4‐2, Aramaki‐machi, Maebashi, Gunma 371‐8510, Japan
tohtsuka@gunma‐u.ac.jp

1.

Introduction

This is a brief introduction to the numerical methods for evolving spirals having several
centers by the crystalline eikonal‐curvature flow.
An anisotropic curvature of a curve \Gamma , which is denoted by H_{\gamma} , is defined by thechanging ratio of an anisotropic perimeter functional

E_{\gamma}( \Gamma)=\int_{\Gamma}\gamma_{0}(-n)d\sigma(1.1)

with a density function \gamma_{0}:\mathbb{S}^{1}arrow(0, \infty) . Here, n\in S^{1} is a continuous normal vector fieldof \Gamma denoting the orientation of \Gamma . It is well‐known that, if \gamma(p)=|p|\gamma_{0}(p/|p|) is smoothon \mathbb{R}^{2}\backslash \{0\} and strictly convex, then the Wulff shape \mathcal{W}_{\gamma} defined by\mathcal{W}_{\gamma}=\{p\in \mathbb{R}^{2};\gamma^{\circ}(p)\leq 1\}(1.2)

has a boundary \partial \mathcal{W}_{\gamma} with constant anisotropic curvature, where \gamma^{\circ}(p)=\sup\{p\cdot q;\gamma(q)\leq

1\} is the support function of \gamma . In fact, if the normal vector for the calculation of thecurvature is oriented to the interior of \mathcal{W}_{\gamma} , then H_{\gamma}=1 on \partial \mathcal{W}_{\gamma} ; see [4] or [7] for details.
We say H_{\gamma} is the crystalline curvature if \mathcal{W}_{\gamma} is a convex polygon. Since \gamma^{\circ} is positively
homogeneous of degree 1, we here assume that \gamma^{\circ}(p) is a convex and piecewise linear
function, i.e.,\gamma^{\circ}(p)=\sup_{0\leq j\leq N_{\gamma}-1}m_{j}\cdot p
for m_{j}\in \mathbb{R}^{2} , where N_{\gamma}\geq 3 is a number of the facets of \mathcal{W}_{\gamma} . In this case, one can find that\gamma=(\gamma^{\circ})^{\circ} is also a convex and piecewise linear function by the convex analysis provided

that \gamma is convex. By the context of the above argument, we here impose the followings.

(A1) \gamma:\mathbb{R}^{2}arrow[0, \infty) is convex,(A2) \gamma is positively homogeneous of degree 1,(A3) \gamma>0 on \mathb {S}^{1},(A4) \gamma is piecewise linear.

75
See also [18] or [19] for the fundamental results from convex analysis used in the above.
Under these settiings, we consider the evolving spirals \Gamma_{t} by the crystalline eikonal‐
curvature flow with driving force of the form

V_{\gamma}=-\rho_{c}H_{\gamma}+f ,(1.3)

where \rho_{c}>0 is a constant, and f denotes the driving force of the evolution. Note thatthe normal velocity V_{\gamma} is calculated with the Finsler metricd_{\gamma}(x, y)=\gamma^{\circ}(x-y)
so that \Gamma_{t} can be evolved as polygonal spirals having parallel facets with those of \mathcal{W}_{\gamma}.
There are some pioneering works by a discrete model by the evolution of facet length

due to [3, 20], which is extended to the motion of a spiral by [8, 9]. However, the evolving
spirals having several centers may merge with each other, which causes singularities on the
curves. Then, it is quite natural to introduce an implicit formulation of evolving curves,which is, for example, level set method due to [17]. However, such a PDE approach withassumption (A4) includes L^{1} type regularization, namely, singular and nonlocal diffusion.
To overcome this difficulty, the variational approach due to [2] is one of powerful options;see [1] for the variational approach to the crystalline curvature motion. Chambolle [6]introduced an algorithm combining the variational approach and the level set method

using signed distance function. Oberman, Osher, Takei and Tsai [13] proposed an iterative
method to calculate Chambolle’s algorithm, which is based on Bregman method [5].
The aim of this paper is to extend the numerical method due to [13] to the evolutionof spirals by (1.3). The crucial difficulty of our problem lies into the fact that a spiralcurve does not divide the domain into two regions. This issue seems to be the same as

the level set method for spirals, which is overcome by [14]. However, even if we use theidea of the sheet structure function due to [11, 14], discontinuity of the signed distancefunction still remains. To overcome this issue, we introduce an idea to use a general levelset function instead of the signed distance function. Note that the idea in this paper isthe revised version of that introduced in [15].

2.

Proposed Algorithm

In this section we propose a minimizing movement approach for evolving spirals by (1.3).
We first review the algorithm by [6] briefly. Let \Omega\subset \mathbb{R}^{n} be a domain, and an interface
Let d:\Omegaarrow \mathbb{R} be a signed distance function from \Gamma , which is positive
We also impose that the normal velocity of \Gamma is positive when \Gamma evolves tothe region where d<0 . Then, consider the minimizer w^{*} of the functional\Gamma\subset\Omega be given.

inside of \Gamma .

E(w)= \int_{\Omega}\gamma(\nabla w)dx+\frac{1}{2h}\Vert w-d\Vert_{L^{2} ^{2}.
Then, w^{*} formally satisfies- div[D\gamma(\nabla w^{*})]+\frac{w^{*}-d}{h}=0,

76which implies\Gamma_{h} :=\{x\in\Omega;w^{*}(x)=0\}=\{x\in\Omega;d(x)=-hdiv[D\gamma(\nabla w^{*}(x))]\}.
Note that -div[D\gamma(Vw^{*})] denotes an anisotropic curvature of level set of w^{*} with ananisotropic energy density \gamma (see [7]). Hence, one can regard \Gamma_{h} as the evolution of \Gamma by
V=-H_{\gamma} in a short time interval [0, h].

When we apply this idea to the evolution of spirals, the crucial difficulty lies in the factthat spiral curve does not divide the domain into two regions, which causes not only theinconsistency between the inside and outside of the spiral, but also the discontinuity ofthe distance function on far away region from the spiral. To overcome these difficulty, we
here propose an approach using general level set function instead of the signed distance.

2.1.

Level set method for evolving spirals

We first review the level set method for evolving spirals due to [14] and its evolution
equation for an anisotropic mean curvature flow with driving force.
Let \Omega\subset \mathbb{R}^{2} be a bounded domain with smooth boundary. We denote the centersa_{N}\in\Omega . We choose r>0 satisfying B_{r}(a_{i})\cap B_{r}(a_{j})=\emptyset ifof spirals by a_{1}, a_{2},i\neq j , and B_{r}(a_{j})\subset\Omega so that W := \Omega\backslash \bigcup_{j=1}^{N}\overline{B_{r}(a_{j})} has smooth boundary, where
B_{r}(a)=\{x\in \mathbb{R}^{2};|x-a|<r\} for r>0 and a\in \mathbb{R}^{2} . Let m_{j}\in \mathbb{Z} be a signed number ofspirals associated with a_{j} . It means that\bullet\bullet|m_{j}| curves are attached to a_{j} as their endpoint,
if m_{j}>0 (resp. m_{j}<0 ), then these curves rotates around a_{j} with counterclockwise(resp. clockwise) rotation provided that V>0.

See [14] for the details of the definition of the signed number of spirals. In this paper,
N as “spirals” ora
we call the union of spiral curves associated with a_{\dot{j} for j=1,2,
spiral pattern” interchangeably.
To describe the spirals with level set formulation, we introduce a pre‐determined func‐tion \theta due to Kobayashi [11, 12]; define\theta(x)=\sum_{j=1}^{N}m_{j}\arg(x-a_{j}).

Let \Gamma_{t}\subset\overline{W} be a spiral pattern at time t\geq 0 , and n\in \mathbb{S}^{1} be a continuous unit normal
vector field of \Gamma_{t} denoting the orientation of the evolution of \Gamma_{t} . Then, we describe \Gamma_{t}and n by\Gamma_{t}=\{x\in\overline{W};u(t, x)-\theta(x)\equiv 0 mod 2\pi \mathbb{Z}\}, n=-\frac{\nabla(u-\theta)}{|\nabla(u-\theta)|}.
Here and hereafter, we write the above shortly as \Gamma_{t}=\{u(t)-\theta\equiv 0\} . Note that \theta should
be a multi‐valued function to describe the spirals completely. However, \nabla\theta is defined as

77a single‐valued function. Therefore, we observe that \Gamma_{t} is locally given as a usual level set
formulation by u-\theta , and then we find

H_{\gamma}=-divD\gamma(\nabla(u-\theta)).

On the normal velocity of \Gamma_{t} , we impose that the normal velocity of \Gamma_{t} should be calculatedwith the metric derived by \gamma , so that \Gamma_{t} should be evolved with the polygonal spiral curve
having parallel facets between those of \mathcal{W}_{\gamma} . Then, we have

V_{\gamma}= \frac{u_{t} {\gamma(\nabla(u-\theta) }.
Hence, we obtain the level set equation of (1.3) asu_{t}-\gamma(\nabla(u-\theta))\{\rho_{c}div[D\gamma(V(u-\theta))]+f\}=0
2.2.

in (0, T)\cross W.

(2.1)

Without distance function scheme

We propose a numerical method to calculate (2.1) due to the minimizing movement ap‐
proach with a general level set function.
Consider the situation that a spiral pattern denoted by \Gamma\subset\overline{W} is given by the levelset formulation\Gamma=\{x\in\overline{W};u(x)-\theta(x)\equiv 0 mod 2\pi \mathbb{Z}\}with u\in C(\overline{W}) . Note that u possibly is not the signed distance function. Let us consider
the minimizer w^{*} ofE(w;u)= \int_{W}\rho_{c}\gamma(\nabla(w-\theta) dx-\int_{W}fwdx+\frac{1}{2h} \Vert\frac{w-u}{\gamma(\nabla(u-\theta) }\Vert_{L^{2} ^{2}
Then, w^{*} should satisfyw^{*}=u+h\gamma(\nabla(u-\theta))\{\rho_{c}div[D\gamma(\nabla(w^{*}-\theta))]+f \}with a suitable boundary condition. We now focus on the analogy between the above andthe implicit finite difference scheme of (2.1), i.e.,u(t+h)\approx u(t)+h\gamma(\nabla(u-\theta))\{\rho_{c}div[D\gamma(\nabla(u- \theta))]+f\}.
Accordingly, we regard

S_{h}(\Gamma):=\{x\in\overline{W};w^{*}-\theta\equiv 0 mod 2\pi \mathbb{Z}\}as the result of the evolution of \Gamma by (1.3) in a short time span
every given timeh>0 .

Hence, fort>0 and the initial spiral patterm \Gamma_{0} , we obtain the solution \Gamma_{t} of theevolution of spirals by (1.3) by \Gamma_{t}=S_{h}^{[t/h]}(\Gamma_{0}) , where [a] is the Gaussian bracket of a\in \mathbb{R}.

78
3.

Split Bregman iteration

Our proposed method includes the problem finding the minimizer w^{*} of an energy func‐
tional of the form

E(w;g)= \int_{W}\rho_{c}\gamma(\nabla(w-\theta) dx-\int_{W}fwdx+\frac{1}{2h} \Vert\frac{w-g}{V}\Vert_{L}^{2},for f, g\in L^{2}(W) and \psi:\overline{W}arrow \mathbb{R} satisfying \varphi/\psi\in L^{2}(W) for \varphi\in L^{2}(W) . In this section,
we give a brief introduction to the split Bregman iteration to solve the above problem,

which is due to [13].

3.1.

Review: Bregman method

We first describe the split Bregman method due to [13] briefly.
The key idea to find w^{*} is to interpret the problem to the following constraint mini‐
mization problem: find a minimizer (w^{*}, d^{*}) of

F(w, d;g) := \int_{W}\rho_{c}\gamma(d-\nabla\theta)dx-\int_{W}fwdx+\frac{1}{2h} \Vert\frac{w-g}{\psi}\Vert_{L^{2} ^{2}subject to d=\nabla w.

For this problem, we introduce a functional with penalty term:

F_{\mu}(w, d;g) :=F(w, d;g)+ \frac{\mu}{2}\Vert d-\nabla w\Vert_{L^{2} ^{2}

= \int_{W}\rho_{c}\gamma(d-\nabla\theta)dx-\int_{W}fwdx+\frac{1}{2h} \Vert\frac{w-g}{\psi}\Vert_{L^{2} ^{2}+\frac{\mu}{2}\Vert d-\nabla w\Vert_{L^{2} }^{2}
for (w, d)\in H^{1}(W)\cross L^{2}(W;\mathbb{R}^{2}) . We now apply Bregman iteration to obtain (w^{*}, d^{*}) .
We shall give a brief review of the Bregman method with an abstract functional. See[16] for details. Let X be a reflexive Banach space. Assume that F, H:Xarrow[0, \infty] be con‐vex functionals, and H is smooth. Let us consider the following constrained minimization
problem:
J(u)=F(u)+\mu H(u) subject to H(u)=0.
We introduce the Bregman distance D_{F}^{p}(u;v) for F of the form

D_{F}^{p}(u;v)=F(u)-F(v)-\langle p, u-v\rangle,
where

of

p\in\partial F(v) and \partial F(v)=\{q\in X^{*};F(u)\geq F(v)+\langle q, u-v\rangle\} is a subdifferential

F at v\in X .

Then, one can find that the iterationu^{k+1}=\arg m\dot{{\imath}}n(D_{F}^{p^{k}}(u;u^{k})u\in X+\mu H(u)) , u^{0}=0, p^{0}=0, p^{k}\in\partial F(u^{k})yields the desired minimizer u^{*}\in X by u^{*}=\lim_{karrow\infty}u^{k} under some suitable assumptions.
Moreover, it is important to find that, if H(u)=\Vert Au-\zeta\Vert_{Y}^{2}/2 with a linear operator

79
A:Xarrow Y onto a Hilbert space Y and\zeta\in Y , then the minimizer u^{k+1} of the above isequivalent to the following iteration:u^{k+1}=\arg m\dot{ \imath} nu\in\cdot(J(u)+\frac{\mu}{2}\Vert Au-\zeta-b^{k} \Vert_{Y}^{2}),b^{k+1}=b^{k}+\zeta-Au^{k+1} , b^{0}=0.

By applying the above sketch to our problem, we consider the itaration(u^{k+1}, d^{k+1})= \arg\min_{(w,d)}[\int_{W}\rho_{c}\gamma(d-\nabla\theta)dx- \int_{W}+ \frac{1}{2h}\Vert\frac{w-g}{\psi}\Vert_{L^{2} ^{2}+\frac{\mu}{2}\Vert d- \nabla w-b^{k}\Vert_{L^{2} ^{2}]
fwdx

,

(3.1)

b^{k+1}=b^{k}+\nabla u^{k+1}-d^{k+1}, u^{0}=g, d^{0}=b^{0}=0(3.2)

with \psi=\gamma(\nabla(g-\theta)) to obtain the minimizer (u^{*}, d^{*}) of F(w, d;g) . To solve this problem,
we introduce the following alternate iteration;u^{k,j+1}= \arg m\dot{ \imath} nw[-\int_{W}fwdx+\frac{1}{2h}\Vert\frac{w-g} {\psi}\Vert_{L^{2} ^{2}+\frac{\mu}{2}\Vert d^{k,j}-\nabla w-b^{k}\Vert_{L^{2} ^{2}],d^{k,j+1}= \arg m_{d}\dot{ \imath} n[\int_{W}\rho_{c}\gamma(d-\nabla\theta)dx+ \frac{\mu}{2}\Vert d-\nabla u^{k,j+1}-b^{k}\Vert_{L^{2} ^{2}]

(3.3)
(3.4)

Then, we observe that (u^{k+1}, d^{k+1})= \lim_{jarrow\infty}(u^{k,j}, d^{k,j}) .

The first minimization (3.3) is established by just solving the following boundary value

problem of an elliptic equation;

w-h\psi\triangle w=g+h\psi(f-div(d^{k,j}-b^{k}))\frac{\partial w}{\partial\nu}=d^{k,j}-b^{k}

in W,

(3.5)

on \partial W.

(3.6)

The second minimization (3.4) is established by calculating the minimizer of integrating
function x\mapsto\rho_{c}\gamma(x-z)+(\mu/2)|x-y|^{2} , which is considered in the next subsection.

3.2.

Shrinkage function

According to [13], we calculate the polyhedral shrinkage function for the second minimiz‐ing problem of the split Bregman iteration. Moreover, our integrating function includes
the translation term \nabla\theta in \gamma . Thus, we present the calculation of polyhedral shrinkage

as in [13] with the translation briefly.
We first consider the minimizing problem of the integrating function in (3.4) withgeneral settings. Remark that \gamma(x) is represented as\gamma(x)=\sup_{p\in \mathcal{W}_{\gamma} p\cdot x

80
when \gamma is convex, where \mathcal{W}_{\gamma}=\{p\in \mathbb{R}^{2};\gamma^{\circ}(p)\leq 1\} , and \gamma^{\circ}(p)=\sup\{p\cdot q;\gamma(q)\leq 1\}.
Let us define the projection map P_{\mathcal{W}_{\gamma} :\mathb {R}^{2}ar ow \mathcal{W}_{\gamma} onto \mathcal{W}_{\gamma} , which is defined by

P_{\mathcal{W}_{\gamma} (x)= \arg\min_{y\in \mathcal{W}_{\gamma} |y-x|^{2}Let us consider the minimization problemx^{*}=x^{*}(y, z)= \arg_{x\in \mathb {R}^{2} m\dot{ \imath} n\{\frac{1}{c} \gamma(x-z)+\frac{1}{2}|x-y|^{2}\}with a constant(3.7)

c>0 with the above notations.

We first derive the necessary condition for the minimizer.

Lemma 3.1 Let x^{*}=x^{*}(y, z) be the minimizer defined by (3.7). Then,x^{*}=y- \frac{1}{c}P_{\mathcal{W}_{\gamma}}(c(y-z).

Proof. Note that, in this proof, we denote the part which has no influence to theminimizer x^{*} by C , which may change by step‐by‐step calculation.
Set \varphi(x,p)=(1/c)p\cdot(x-z)+(1/2)|x-y|^{2} and interpret the minimizing problem(3.7) to finding x^{*} which attains the following minimum:

M= \min_{x\in \mathb {R}^{2}p}\max_{\in \mathcal{W}_{\gamma} \varphi(x,p)=\min_ {x\in \mathb {R}^{2}p}\max_{\in \mathcal{W}_{\gamma} (\frac{1}{c}p\cdot(x-z)+ \frac{1}{2}|x-y|^{2}).

Note that x\mapsto\varphi(x,p) is convex, \lim_{|x|arrow\infty}\varphi(x,p)=\infty uniformly for p\in \mathcal{W}_{\gamma} , andp\mapsto\varphi(x, p) is affine. Then, by the minimax theorem due to [21] we observe that

M= \max_{p\in}\min_{\mathcal{W}_{\gamma^{x\in \mathb {R}^{2} }(\frac{1}{c} p\cdot(x-z)+\frac{1}{2}|x-y|^{2}).

By stragithforward calculation we have\varphi(x,p)=\frac{1}{2}|x|^{2}- (\begin{ar y}{l 1 y-p\overline{c} \end{ar y}) x+C= \frac{1}{2}|x- (\begin{ar y}{l 1 y-p\overline{c} \end{ar y}) |^{2}+C.
.

Hence, we obtain\min_{x\in \mathbb{R}^{2} \varphi(x,p)=\varphi(\begin{ar ay}{l 1 y-p, \overline{c} \end{ar ay})=\frac{1}{c}p (\begin{ar ay}{l 1 y-zp\overline{c} \end{ar ay}) + \frac{1}{2c^{2} |p|^{2}.

=- \frac{1}{2c^{2} |p|^{2}+\frac{1}{c}p\cdot(y-z)=-\frac{1}{2c^{2} |p-c(y-z) |^{2}+C.
Consequently, the minimum M is attained atx^{*}(y, z)=y- \frac{1}{c}p^{*}, p^{*}=P_{\mathcal{W}_{\gamma}}(c(y-z)).

\squareThe following corollary directly follows from Lemma 3.1, so that we omit its proof.

81
81
Corollary 3.2 It holds that c(y-z)\in \mathcal{W}_{\gamma} if and only if x^{*}(y, z)=z.
Let us consider the case when c(y-z)\not\in \mathcal{W}_{\gamma} . From here on, let \gamma be a convex andpiecewise linear function, i.e.,\gamma(p)=_{0\leq j}\max_{\leq N_{\gamma}-1}n_{j}\cdot p,
where n_{j}=r_{j}(\cos\theta_{j}, \sin\theta_{j}) ( r_{i}>0 is a constant) satisfies(W1) \theta_{0}<\theta_{1}<\theta_{2}<<\theta_{N_{\gamma}-1}<\theta_{0}+2\pi,(W2) \theta_{j}<\theta_{j+1}<\theta_{j}+\pi for j=0,1 ,, N_{\gamma}-1.

Note that, the index number j of n_{j} has been extended to \mathb {Z} by regarding j\in \mathbb{Z}/(N_{\gamma}\mathbb{Z}) ,in other words, n_{j+kN_{\gamma}}=n_{j} for k\in \mathbb{Z} . For the calculation of the subdifferential \partial\gamma of \gamma,
we partition \mathbb{R}^{2}\backslash \{0\} into the regions
Q_{i}= { p\neq 0;n_{i}\cdot p>n_{j}\cdot p for j\neq i },R_{i,i+1}=\{p\neq 0;n_{i}\cdot p=n_{i+1}\cdot p>0\}

for i\in \mathbb{Z} . Assume the following.

(W3) Q_{i}=\Xi_{i,i-1}\cap\Xi_{i,i+1}\neq\emptyset , where \Xi_{j,k}=\{p\neq 0;n_{j}\cdot p>n_{k}\cdot p\}.

Note that (W3) implies \partial Q_{i}=R_{i-1,i}\cup R_{i,i+1}\cup\{0\} . Moreover, \gamma is smooth on Q_{i} and has
singularities on each R_{i,i+1} . By using the above notations, we now characterize x^{*} whenc(y-z)\not\in \mathcal{W}_{\gamma}.

Lemma 3.3 Let \gamma(p)=\max_{0\leq j\leq N_{\gamma}-1}n_{j} . p with n_{j}=r_{j}(\cos\theta_{j}, \sin\theta_{j}) for r_{j}>0 and\theta_{j}\in \mathbb{R} , and assume that (W1)-(W3) holds. Let c>0 and y, z\in \mathbb{R}^{2} satisfy c(y-z)\not\in \mathcal{W}_{\gamma}.
Set\lambda_{i}=\frac{c(y-z)\cdot(n_{i}-n_{i+1})-n_{i}\cdot n_{i+1}+|n_{i+1}|^{2} {|n_{i}-n_{i+1}|^{2} ,\xi_{i}=\lambda_{i}n_{i}+(1-\lambda_{i})n_{i+1} .

For the minimizer x^{*} defined by (3.7), the followings hold.
(i) If x^{*}-z\in R_{i,i+1} , then \lambda_{i}\in[0,1] and (c(y-z)-\xi_{i})\cdot\xi_{i}\geq 0.

(ii) If x^{*}-z\in Q_{i} , then \lambda_{i}>1 and \lambda_{i-1}<0.
Proof. By (3.7) we have\frac{1}{c}\gamma(x-z)+\frac{1}{2}|x-y|^{2}\geq\frac{1}{c}\gamma(x^{*}-z)+ \frac{1}{2}|x^{*}-y|^{2}.
This implies that\gamma(x-z)\geq\gamma(x^{*}-z)-c(x^{*}-y)\cdot(x-x^{*})-\frac{c}{2}|x-x^{*} |^{2},

(3.8)
(3.9)

82
and thus -c(x^{*}-y)\in\partial\gamma(x^{*}-z) . The subdifferential \partial\gamma can be explicitly calculated as\parti l\gam a(q)=\{ begin{ar y}{l \{n_i}\ ifq\inQ_{i} forsomei\n {0,1 N_{\gam a}-1\}, \{ lambdan_{i}+(1-\lambda)n_{i+1};\lambda\in[0,1]\}ifq\inR_{i,+1} forsomei\n {0,1 N_{\gam a}-1\}. \end{ar y}

By Corollary 3.2, c(y-z)\not\in \mathcal{W}_{\gamma} implies that x^{*}-z\neq 0 , and thus x^{*}-z \in\bigcup_{i=0}^{N_{\gamma}-1}Q_{i}\cup\bigcup_{\dot{i}=0}^{N_{\gamma}-1}R_{i, +1} . By combining the above, we obtain the followings.

(i) If x^{*}-z\in R_{\dot{i}},i+1 , then x^{*}=y-(1/c)(\lambda n_{i}+(1-\lambda)n_{i+1}) for some \lambda\in[0,1].

(ii) If x^{*}-z\in Q_{i} , then x^{*}=y-(1/c)n_{i}.
Let us consider the case (i). In this case, we have-c(x^{*}-y)=\lambda n_{i}+(1-\lambda)n_{i+1}for some \lambda\in[0,1] . By taking an inner product between n_{i} or n_{i+1} and the above, we
obtain-c(x^{*}-y)\cdot n_{i}=\lambda|n_{i}|^{2}+(1-\lambda)n_{i+1}\cdot n_{i} ,-c(x^{*}-y)\cdot n_{i+1}=\lambda n_{i}\cdot n_{i+1}+(1-\lambda)|n_{i+1}|^{2}(3.10)(3.11)

Note that x^{*}-z\in R_{\dot{i}},i+1 implies n_{i}\cdot(x^{*}-z)=n_{i+1}\cdot(x^{*}-z) . Then, (3.10) and (3.11)
yield the following;c(y-z)\cdot n_{i}-\lambda|n_{i}|^{2}-(1-\lambda)n_{i+1}\cdot n_{i}=c(y-z)\cdot n_{i+1}-\lambda n_{i}\cdot n_{i+1}-(1-\lambda)|n_{i+1}|^{2}We calculate the above and obtain\lambda=\lambda_{i}=\frac{c(y-z)\cdot(n_{\dot{i} -n_{i+1})-n_{i}\cdot n_{i+1}+ |n_{i+1}|^{2} {|n_{\dot{i} -n_{i+1}|^{2} \in[0,1].
Moreover, we have x^{*}=y-(1/c)\xi_{i} , and then (x^{*}-z) . n_{i}=(x^{*}-z)\cdot n_{i+1}\geq 0 . Thisimplies(c(y-z)- \xi_{\dot{i} )\cdot\xi_{i}=(c(x^{*}-\frac{1}{c}\xi_{i}-z)-\xi_{i}).

\xi_{i}=c[\lambda_{i}(x^{*}-z)\cdot n_{i}+(1-\lambda_{i})(x^{*}-z)\cdot n_{i+1}]\geq 0.

We next consider the case (ii). In this case we have(x^{*}-z)\cdot n_{i}>(x^{*}-z)\cdot n_{j}forj\neq i.

Since x^{*}=y-(1/c)n_{i} , the formula (3.8) yields that\lambda_{i}=\frac{c(y-z)\cdot(n_{i}-n_{i+1})-n_{i}\cdot n_{i+1}+|n_{i+1}|^{2} {|n_{\dot{i} -n_{i+1}|^{2}= \frac{c(x^{*}-z+(1/c)n_{i})\cdot(n_{i}-n_{i+1})-n_{i}\cdot n_{i+1}+|n_{i+1}|^ {2} {|n_{i}-n_{i+1}|^{2}= \frac{c(x^{*}-z)\cdot(n_{i}-n_{i+1})+|n_{i}-n_{i+1}|^{2} {|n_{i}-n_{\dot{i}+ 1}|^{2} >1.

83
Similarly, we obtain\lambda_{i-1}=\frac{c(y-z)\cdot(n_{i-1}-n_{i})-n_{i-1}\cdot n_{i}+|n_{i}|^{2} {|n_{i-1}-n_{i}|^{2}= \frac{c(x^{*}-z+(1/c)n_{i})\cdot(n_{i-1}-n_{\dot{i} )-n_{\dot{i}-1}\cdot n_{\dot{i} +|n_{i}|^{2} {|n_{i-1}-n_{i}|^{2}= \frac{c(x^{*}-z)\cdot(n_{i-1}-n_{i})}{|n_{i-1}-n_{i}|^{2} <0. \square
By Lemma 3.3, either the following (A) or (B) occurs when c(y-z)\not\in \mathcal{W}_{\gamma} :(A) There exists i\in\{0,1, N_{\gamma}-1\} such that \lambda_{i}\in[0,1] and (y-(1/c)\xi_{i})\cdot\xi_{i}\geq 0.
(B) There exists i\in\{0,1, N_{\gamma}-1\} such that \lambda_{i}>1 and \lambda_{i-1}<0.
Hence, by combining Corollary 3.2 and Lemma 3.3, we obtain the following scheme tocalculate the minimizer d^{k,j+1} as (3.4). See [10] for the practical way to construct \gamma^{\circ} froma polyhedral \gamma.

\frac{Scheme(Sh)toca1cu1ated^{k,j+1}of(3.4)}{Sety=\nabla u^{k,j+1}(x)+b^{k}(x) andz=\nabla\theta(x})\bullet.

If c(y-z)\in \mathcal{W}_{\gamma} , i.e., \gamma^{\circ}(y-z)\leq 1/c , then set

d^{k,j+1}=z=\nabla\theta(x)\bullet.

If \gamma^{\circ}(y-z)>1/c, then set d^{k,j+1} as follows: Set \lambda_{i} and \xi_{i} by (3.8) and (3.9),
respectively.
‐ If there exists i\in\{0,1, , N_{\gamma}-1\} such that \lambda_{i}\in[0,1] and (y-(1/c)\xi_{i})\cdot\xi_{i}\geq 0,then setd^{k,j+1}(x)=y- \frac{1}{c}\xi_{i}=\nabla u^{k,j+1}(x)+b^{k}(x)-\frac{1}{c}\xi_ {i}.
‐ Otherwise, find i\in\{0,1, , N_{\gamma}-1\} satisfying \lambda_{i}>1 and \lambda_{i-1}<0 , and thensetd^{k,j+1}(x)=y- \frac{1}{c}n_{i}=\nabla u^{k,j+1}(x)+b^{k}(x)-\frac{1}{c}n_{i}.
4.

Numerical results

In this section, we propose some numerical results on the evolution of spirals by (1.3) byour scheme.

For a given initial curve \Gamma_{0}\subset \mathbb{R}^{n} , we give a continuous initial data u_{0}\in C(\overline{W})satisfying \Gamma_{0}=\{u_{0}-\theta\equiv 0\} . Then, by repeating the algorithm proposed in §2.2, we
obtain a sequence of minimizers u_{n} := \arg\min_{u}E(u;u_{n-1}) and a family of the spirals\Gamma_{n}=\{u_{n}-\theta\equiv 0\} . One can find two loops to obtain u_{n+1} from u_{n} ;84(Outer) Finding u_{n+1} from u_{n} by (3.1)-(3.2) ,(Inner) Finding (u_{n}^{k+1}, d^{k+1}) from (u_{n}^{k}, d^{k}) by (3.3)-(3.4) .
See Algorithm 1 for the summary of our algorithm. In the both Outer and Inner loops, we
calculate F_{\mu}(u_{n}^{k,j}, d^{k,j};u_{n}) or F_{\mu}(u_{n}^{k}, d^{k};u_{n}) to check if they are the minimizer of F_{\mu}(w, d;u_{n})with \psi=\gamma(\nabla(u_{n}-\theta)) . However, it is failed when the set \{x\in\overline{W};\gamma(\nabla(u_{n}-\theta))=0\}
has nonempty interior. To avoid this issue, we introduce an approximation of F_{\mu}(w, d;g)by cut‐off only for checking the condition breaking the loop; define

F_{\mu,\alpha}(w, d;g):= \int_{W}\rho_{c}\gamma(d-\nabla\theta)dx-\int_{W}fwdx+ \frac{1}{2h}\Vert\frac{w-g}{\max\{\gamma(\nabla(g-\theta) ,\alpha\} \Vert_{L^ {2} ^{2}+\frac{\mu}{2}\Vert d-\nabla w\Vert_{L^{2} ^{2},

with a positive constant \alpha\ll 1 . Then, to break the inner and outer loop, we check if(for Inner)(for Outer)|F_{\mu,\alpha}(u_{n}^{k,j+1}, d^{k,j+1};u_{n})-F_{\mu,\alpha}(u_{n}^{k,j}, d^{k,j};u_{n})|<\varepsilon_{in},|F_{\mu,\alpha}(u_{n}^{k+1}, d^{k+1};u_{n})-F_{\mu,\alpha}(u_{n}^{k}, d^{k}; u_{n})|<\varepsilon_{out}for some \varepsilon_{in}, \varepsilon_{out}\ll 1 . The equation (3.5)-(3.6) can be solved by SOR method. Remark
Algorithm 1: Minimizing movement without distance function
Input: \Gamma_{0}\subset\overline{W} and u_{0}\in C(\overline{W}) such that \Gamma_{0}=\{u_{0}-\theta\equiv 0\}.

Output: \Gamma(T)=\{u(T)-\theta\equiv 0\} with some function u(T):\overline{W}arrow \mathbb{R}.

(Time step) for n=0,1,[T/h]-1 do\ovalbx{tsmREJCT}In\do{iath1zeu_}^0={n,'bd\cot=.0;Seu_{n+1}^ktg=;(u_{Seb+1}^k\nalu_{ 1}-d^k+;\ovalbx{tsmREJCT}\for_{i|Fmu,alph(n}^kj+1d{,;g)-F_\mualph}({n^kj,d_\otJ;g)|<varepsiln{}hb^So,1ve(35)-6with\ps=gam(nbl,-the)oaiu_{n}^kj+1;f=.0,do
.

reak;

that we do not apply the cut‐off of \psi=\gamma(\nabla(u_{n}-\theta)) for (3.5).
Finally, we present some numerical results on the evolution of spirals by (1.3) with
f\equiv 2 and \rho_{c}=0.04 , i.e.,

V_{\gamma}=2-0.04H_{\gamma}

(4.1)

85
with Algorithm 1. Here we set \Omega=[-2.5,2.5]^{2} and the numerical lattice on \Omega as D=\{x_{i,j}=(i\triangle x, j\triangle x)\in\Omega;-100\leq i, j\leq 100\} , so that the mesh size is \triangle x=0.025 . We
impose that a_{j}\in D , and set r=(2-10^{-8})\triangle x . The time span and the number of time
steps are chosen as h=0.02\cross\triangle x and 2000 steps, and then T=1.0 . For the split
Bregman method, we set \mu=\rho_{c}, \varepsilon_{in}=10^{-2}, \varepsilon_{out}=10^{-5} , and \alpha=10^{-8}.
Figure 1 and 2 presents profiles of spirals with two centers; N=2, a_{1}=(-0.75,0) ,a_{2}=(0.75,0) . The settings for Figure 1 and Figure 2 are as follows.

(Figure 1): \gamma(p)=\gamma_{1}(p)=|p_{1}|+|p_{2}| , i.e., N_{\gamma}=4,\theta_{j}=\frac{\pi}{4}+\frac{\pi j}{2}, r_{j}=\sqrt{2},\theta(x)=\arg(x-a_{1})+\arg(x-a_{2}) , u_{0}\equiv 0 .

(Figure 2): \gamma(p)=\sqrt{2}\gamma_{\infty}(p)=\sqrt{2}\max\{|p_{1}|, |p_{2}|\} , i.e., N_{\gamma}=4,\theta(x)=\arg(x-a_{1})-\arg(x-a_{2}) , u_{0}\equiv-\pi .

(4.2)\theta_{j}=\frac{\pi j}{2}, r_{j}=\sqrt{2},(4.3)

for p=(p_{1},p_{2}) . Note that \gamma_{1}=\gamma_{\infty} and \gamma_{\infty}^{\circ}=\gamma_{1} , and then \mathcal{W}_{\gamma_{1}}=[-1,1]^{2} is a square and\mathcal{W}_{\gamma_{\infty} is a diagonal square whose vertices are (\pm 1,0) and (0, \pm 1) . As in Figure 1 and 2,
one can find the evolution of square spiral with topological change reflecting the shape of\mathcal{W}_{\gamma_{1} and \mathcal{W}_{\gamma_{\infty} , respectively. Moreover, we can choose an initial curve which is not parallelto the facets of \mathcal{W}_{\gamma} as in Figure 2. Note that such an initial curve is not admissible which

is defined in [9], although the admissibility is required in the discrete model due to [9].

Our scheme can be applied to the evolution which violates the admissibility of the curvein finite time.

t=0

t=0.2

t=0.6

t=0.8

t=0.4

t=1.0

Figure 1: Profiles of spirals evolving by (4.1) with (4.2).

86

t=0

t=0.2

t=0.6

t=0.8

t=0.4

t=1.0

Figure 2: Profiles of spirals evolving by (4.1) with (4.3).

References[1] Fred Almgren and Jean E. Taylor. Flat flow is motion by crystalline curvature forcurves with crystalline energies. J. Differential Geom., 42(1):1-22 , 1995.
[2] Fred Almgren, Jean E. Taylor, and Lihe Wang. Curvature‐driven flows: a variational
approach. SIAM J. Control Optim., 31(2):387-438 , 1993.
[3] Sigurd Angenent and Morton E. Gurtin. Multiphase thermomechanics with interfa‐
cial structure. II. Evolution of an isothermal interface. Arch. Rational Mech. Anal.,

108(4):323-391 , 1989.

[4] G. Bellettini and M. Paolini. Anisotropic motion by mean curvature in the context
of Finsler geometry. Hokkaido Math. J., 25(3):537-566 , 1996.
[5] L. M. Brègman. A relaxation method of finding a common point of convex sets andits application to the solution of problems in convex programming. Ž. Vyčisl. Mat.
i Mat. Fiz., 7:620−631, 1967.

[6] Antonin Chambolle. An algorithm for mean curvature motion.
Bound., 6(2):195-218 , 2004.

Interfaces Free

[7] Yoshikazu Giga. Surface evolution equations: A level set approach, volume 99 of
Monographs in Mathematics. Birkhäuser Verlag, Basel, 2006.

87[8] Tetsuya Ishiwata. Crystalline motion of spiral‐shaped polygonal curves with a tip
motion. Discrete Contin. Dyn. Syst. Ser. S, 7(1):53-62 , 2014.
[9] Tetsuya Ishiwata and Takeshi Ohtsuka. Evolution of spiral‐shaped polygonal curveby crystalline curvature flow with a pinned tip. Discrete Contin. Dyn. Syst. Ser. B,
to appear.

[10] Tetsuya Ishiwata and Takeshi Ohtsuka. Numerical analysis comparing ODE ap‐
proach and level set method for evolving spirals by crystalline eikonal‐curvature flow.
Hokkaido University Preprint Series in Mathematics, \# 1121:1-16 , 2019.

[11] Ryo Kobayashi.

A brief introduction to phase field method.

AIP Conf. Proc.,

1270:282−291, 2010.

[12] Hitoshi Miura and Ryo Kobayashi. Phase‐field modeling of step dynamics on growing
crystal surface: Direct integration of growth units to step front. Cryst. Growth Des.,

15(5):2165-2175 , 2015.
[13] Adam Oberman, Stanley Osher, Ryo Takei, and Richard Tsai. Numerical methods
for anisotropic mean curvature flow based on a discrete time variational formulation.
Commun. Math. Sci., 9(3):637-662 , 2011.

[14] T. Ohtsuka, Y.‐H.R. Tsai, and Y. Giga. A level set approach reflecting sheet structurewith single auxiliary function for evolving spirals on crystal surfaces.
Scientific Computing, 62(3):831-874 , 2015.

Journal of[15] Takeshi Ohtsuka. Minimizing movement approach for spirals evolving by crystalline
curvature using level set functions. Oberwolfach Reports, 14:314−317, 2017.

[16] Stanley Osher, Martin Burger, Donald Goldfarb, Jinjun Xu, and Wotao Yin. An
iterative regularization method for total variation‐based image restoration. Multiscale
Model. Simul., 4(2):460-489 , 2005.

[17] Stanley Osher and James A. Sethian. Fronts propagating with curvature‐dependent
speed:

algorithms based on Hamilton‐Jacobi formulations.

J. Comput. Phys.,

79(1):12-49 , 1988.

[18] R. Tyrrell Rockafellar. Convex analysis. Princeton Mathematical Series, No. 28.
Princeton University Press, Princeton, N.J., 1970.

[19] Rolf Schneider. Convex bodies: the Brunn‐Minkowski theory, volume 44 of Encyclo‐
pedia of Mathematics and its Applications. Cambridge University Press, Cambridge,
1993.

[20] Jean E. Taylor. Constructions and conjectures in crystalline nondifferential geometry.
In Differential geometry, volume 52 of Pitman Monogr. Surveys Pure Appl. Math.,pages 321‐336. Longman Sci. Tech., Harlow, 1991.

[21] J. v. Neumann. Zur Theorie der Gesellschaftsspiele. Math. Ann., 100(1):295-320,
1928.