Numerical analysis comparing ODE approach and level set method for evolving spirals by crystalline eikonal-curvature flow

TL;DR

This paper compares the numerical results of ODE and level set methods for simulating the evolution of crystalline spiral curves, finding minimal differences in their area calculations despite different evolution laws.

Contribution

It provides a numerical comparison between ODE and level set models for crystalline spiral evolution, highlighting their close agreement in results.

Findings

Small area differences between models

Differences in evolution laws are minor

Both methods produce consistent results

Abstract

In this paper, the evolution of a polygonal spiral curve by the crystalline curvature flow with a pinned center is considered with two view points, discrete model consist of an ODE system of facet lengths and a level set method. We investigate the difference of these models numerically by calculating the area of the region enclosed by these spiral curves. The area difference is calculated by the normalized L1 norm of the difference of step-like functions which are branches of arg(x) whose discontinuities are only on the spirals. We find the differences of the numerical results considered in this paper are very small even though the evolution laws of these models around the center and the farthest facet are slightly different.

Equations

$$E_{\gamma }(\Gamma )={\int }_{\Gamma }\gamma (\mathbf{n})d\sigma$$

$${\mathcal{W}}_{\gamma }=\{p\in {\mathbb{R}}^2;p\cdot q\le \gamma (q)\text{mbox}forq\in {\mathbb{S}}^1\},$$

$$\beta V_{\gamma }=U-{\rho }_c{H}_{\gamma }\text{mbox}on{\Gamma }_t,$$

$${\Gamma }_L(t)=\{x;u(t,x)-\theta (x)\equiv 0\mathrm{mod}2\pi \mathbb{Z}\},\mathbf{n}=-\frac{\nabla (u-\theta )}{|\nabla (u-\theta )|}$$

$$V_{\gamma }=\frac{u_t}{\gamma (-\nabla (u-\theta ))},H_{\gamma }=-\mathrm{div}\{\xi (-\nabla (u-\theta ))\},$$

$$\tilde{\beta }(\nabla (u-\theta ))u_t-\tilde{\gamma }(\nabla (u-\theta ))\left[\mathrm{div}\{\tilde{\xi }(\nabla (u-\theta ))\}+U\right]=0,$$

$$\mathcal{D}(t)=\frac{1}{|W|}{\int }_W\frac{{\theta }_D(t,x)-{\theta }_L(t,x)}{2\pi }dx,$$

$${\mathbf{N}}_j=(\cos {\phi }_j,\sin {\phi }_j),{\mathbf{T}}_j=(\sin {\phi }_j,-\cos {\phi }_j).$$

$$L_{j}(t)=\left\{\begin{array}[]{ll}\{\lambda y_{j}(t)+(1-\lambda)y_{j-1}(t);\ \lambda\in[0,1]\}&\mbox{for}\ j=k,k-1,\ldots,1\\ \{y_{0}(t)+\lambda{\bf T}_{0};\ \lambda>0\}&\mbox{if}\ j=0\end{array}\right.$$

$${\mathbf{T}}_j=\frac{y_{j-1}(t)-y_j(t)}{|y_{j-1}(t)-y_j(t)|}.$$

$${\dot{d}}_k=c_k^{-}\left(U-\frac{{\rho }_c{\ell }_{k-1}}{d_{k-1}}\right),$$

$$\{\begin{array}{rl}{\dot{d}}_{k-1}& =-b_{k-1}\left(U-\frac{{\rho }_c{\ell }_{k-1}}{d_{k-1}}\right)+c_{k-1}^{-}\left(U-\frac{{\rho }_c{\ell }_{k-2}}{d_{k-2}}\right),\\ {\dot{d}}_j& =-b_j\left(U-\frac{{\rho }_c{\ell }_j}{d_j}\right)+c_j^{+}\left(U-\frac{{\rho }_c{\ell }_{j+1}}{d_{j+1}}\right)+c_j^{-}\left(U-\frac{{\rho }_c{\ell }_{j-1}}{d_{j-1}}\right)\\ & \text{mbox}forj=2,3,\dots ,k-2,\\ {\dot{d}}_1& =-b_1\left(U-\frac{{\rho }_c{\ell }_1}{d_1}\right)+c_1^{+}\left(U-\frac{{\rho }_c{\ell }_2}{d_2}\right)+c_1^{-}U,\end{array}$$

$$b_j=\frac{1}{{\beta }_j}\left(\frac{1}{\tan ({\phi }_{j+1}-{\phi }_j)}+\frac{1}{\tan ({\phi }_j-{\phi }_{j-1})}\right),c_j^{\pm }=\pm \frac{1}{{\beta }_{j\pm 1}\sin ({\phi }_{j\pm 1}-{\phi }_j)}$$

$$y_k(t)=O,y_{j-1}(t)=y_j(t)+d_j(t){\mathbf{T}}_j\text{mbox}forj=k,k-1,k-2,\dots ,1.$$

$${\Gamma }_D(0)=L_1(0)\cup L_0(0)=\{\lambda {\mathbf{T}}_0;\lambda \ge 0\}.$$

$$T_{k+1}=\sup \{T>T_k;d_k(t)<{\rho }_c{\ell }_k/U\text{mbox}fort\in [T_k,T]\}.$$

$${\Gamma }_L(t)=\{x\in \overline{W};u(t,x)-\theta (x)\equiv 0\mathrm{mod}2\pi \mathbb{Z}\},\mathbf{n}=-\frac{\nabla (u-\theta )}{|\nabla (u-\theta )|},$$

$$H_{\gamma }=-\mathrm{div}\xi (-\nabla (u-\theta ))$$

$${\mathcal{W}}_{\gamma }=\{p\in {\mathbb{R}}^2;{\gamma }^{\circ }(p)\le 1\}$$

$$d_{\gamma }(x,y)={\gamma }^{\circ }(x-y)$$

$$\text{mbox}with{\gamma }^{\circ }({\mathbf{N}}_j)=1\text{mbox}forj=0,1,2,\dots ,N_{\gamma }-1.$$

$$V_{\gamma }=\frac{u_t}{\gamma (-\nabla (u-\theta ))}$$

$$\tilde{\beta }(\nabla (u-\theta ))u_t-\tilde{\gamma }(\nabla (u-\theta ))\left\{{\rho }_c\mathrm{div}\tilde{\xi }(\nabla (u-\theta ))+U\right\}=0$$

$$\vec{\nu }\cdot \nabla (u-\theta )=0$$

$$\gamma (p)=\underset{0\le j\le N_{\gamma }-1}{\max }{n}_j\cdot p=\sum _{j=0}^{N_{\gamma }-1}(n_j\cdot p){\chi }_{Q_j}(p)$$

$$\chi_{Q}(x)=\left\{\begin{array}[]{ll}1&\mbox{if}\ x\in Q,\\ 0&\mbox{otherwise}\end{array}\right.$$

$$\xi (p)=\sum _{j=0}^{N_{\gamma }-1}{\chi }_{Q_j}(p)n_j,$$

$$\sigma(z;p_{1},p_{2}):=\left\{\begin{array}[]{ll}{\displaystyle\frac{z}{\sqrt{z^{2}+\varepsilon^{2}(|p_{1}|+|p_{2}|)^{2}}}}&\mbox{if}\ z\neq 0,\\[5.0pt] 0&\mbox{otherwise}\end{array}\right.$$

$${\chi }_Q(x)\approx \zeta (f(x);p_1,p_2),\text{mbox}with\zeta (z;p_1,p_2):=\frac{\sigma (z;p_1,p_2)+1}{2}$$

$$\xi (p)\approx \sum _{j=0}^{N_{\gamma }-1}\zeta (f_j(x);p_1,p_2)n_j$$

$${\mathcal{W}}_{\gamma }=\{p\in {\mathbb{R}}^2;{\gamma }^{\circ }(p)\le 1\}$$

$${\gamma }^{\circ }(p):=\underset{0\le j\le N_{\gamma }-1}{\max }{m}_j\cdot p,m_j={\eta }_j(\cos {\psi }_j,\sin {\psi }_j)$$

Full text

The first author is partly supported by JSPS KAKENHI Grant Number 15H03632 and 16H03953.

Tetsuya Ishiwata

Department of Mathematical Sciences, Shibaura Institute of Technology

Fukasakuk 309, Minuma-ku, Saitama 337-8570, Japan

Takeshi Ohtsuka

Division of Pure and Applied Science, Faculty of Science and Technology, Gunma University

Aramaki-machi 4-2, Maebashi, 371-8510 Gunma, Japan

1. Introduction

The crystalline curvature of a curve $\Gamma$, which is denoted by $H_{\gamma}$, is defined by the changing ratio of an anisotropic surface energy functional

[TABLE]

for a singular density function $\gamma\colon\mathbb{R}^{2}\to[0,\infty)$ with respect to the volume of a region enclosed by $\Gamma$, where ${\bf n}$ is a continuous unit normal vector field of $\Gamma$ and $d\sigma$ is the line element. Here, singular means that the Wulff shape

[TABLE]

which satisfies $H_{\gamma}=1$ on $\partial\mathcal{W}_{\gamma}$, is a convex polygon. See [7] for details of the crystalline curvature. Such a singular energy expresses the surface energy of the polygonal structure of interfaces like as crystal surface. The typical example of $\gamma$ is $\ell^{1}$ norm. For describing general settings, we here assume that

1. (A1)

$\gamma$ is convex, 2. (A2)

$\gamma$ is positively homogeneous of degree 1, i.e., $\gamma(\lambda p)=\lambda\gamma(p)$ for $p\in\mathbb{R}^{2}$ and $\lambda>0$, 3. (A3)

$\gamma>0$ on $\mathbb{S}^{1}$ 4. (A4)

$\gamma$ is piecewise linear.

Note that (A2) is for the level set formulation of curves mentioned later. Moreover, (A4) is a sufficient condition to the singularity of $\mathcal{W}_{\gamma}$ for the crystalline curvature, since $\mathcal{W}_{\gamma}=\{p\in\mathbb{R}^{2};\ \gamma^{\circ}(p)\leq 1\}$ and $(\gamma^{\circ})^{\circ}=\gamma$ if $\gamma$ is convex, where $\gamma^{\circ}(p):=\sup\{p\cdot q;\ \gamma(q)\leq 1\}$ is a support function of $\gamma$. See [13] for details of the properties of $\gamma$ and $\gamma^{\circ}$.

In this paper we consider the evolution of a convex polygonal spiral by

[TABLE]

where $V_{\gamma}$ is an anisotropic normal velocity under the Finsler metric defined by $\mathrm{dist}_{\gamma}(x,y)=\gamma^{\circ}(x-y)$, and $U>0$ and $\rho_{c}>0$ are assumed to be constants. (Note that we do not assume the symmetricity of this metric.) For this evolution of a pinned spiral, the authors of this paper introduce a discrete model by an ODE system of the facet lengths in [9], due to the idea of [1, 14, 8], see also [7] for details.

On the other hand, Tsai, Giga and the second author [11, 10] introduced a level set formulation for evolving spirals with fixed centers. According to their formulation, an evolving spiral curve with a fixed center at the origin is given as

[TABLE]

with an auxiliary function $u(t,x)$ and a pre-defined multivalued function $\theta(x)=\arg x$. Then, $V_{\gamma}$ and $H_{\gamma}$ are interpreted as

[TABLE]

where $\xi=D\gamma$. Hence, we obtain the level set equation for (1) of the form

[TABLE]

where $\tilde{\beta}(p)=\beta(-p)$, $\tilde{\gamma}(p)=\gamma(-p)$, and $\tilde{\xi}(p)=\xi(-p)$.

The aim of this paper is to show the numerical difference between the spirals calculated by the discrete model due to [9] and the level set method due to [10]. To measure the difference between these spirals, we calculate the area of the region enclosed by their spirals. It is established by calculating

[TABLE]

where $\theta_{D}$ and $\theta_{L}$ are branches of $\theta$ whose discontinuities are only on the spiral curves $\Gamma_{D}(t)=\sum_{j=0}^{k}L_{j}(t)$ obtained by the discrete algorithm and $\Gamma_{L}(t)$ by the level set method, respectively. A practical way to construct $\theta_{L}$ from solution $u$ of the level set equation is provided in [10]. Thus, we shall give a way to construct $\theta_{D}$ in §3.2. Note that the discrete model in [9] is constructed from $\mathcal{W}_{\gamma}$, we shall give a way to construct $\gamma$ from $\gamma^{\circ}$ in §3.1 to obtain the level set equation corresponding to the discrete model.

2. Models

In this section, we recall the discrete model due to [9] and the level set method due to [10]. To compare the evolving spiral curves from these models, we have to give a Wulff shape $\mathcal{W}_{\gamma}$ for the discrete model and corresponding surface energy density $\gamma$ for the level set method. In this section, we consider the situation $\mathcal{W}_{\gamma}$ and corresponding $\gamma$ are already given. A practical way to obtain $\gamma$ from $\mathcal{W}_{\gamma}$ will be discussed in §3.1. We briefly review mathematical results on these models.

2.1. Discrete model

We recall the ODE model by [9].

We first prepare some notations for $\mathcal{W}_{\gamma}$. Let $\mathcal{W}_{\gamma}$ be a $N_{\gamma}$ sided convex polygon. The $j$-th facet of $\mathcal{W}_{\gamma}$ has an outer unit normal vector $\mathcal{N}_{j}$ with angle $\varphi_{j}$ for $j=0,1,2,\ldots,N_{\gamma}-1$. Set the unit tangential vector ${\bf T}_{j}$ of the $j$-th facet as well as the definition of the Frenet frame, i.e.,

[TABLE]

We assume the followings for expressing the convexity of $\mathcal{W}_{\gamma}$.

1. (W1)

$\varphi_{0}<\varphi_{1}<\varphi_{2}<\cdots<\varphi_{N_{\gamma}-1}<\varphi_{0}+2\pi$. 2. (W2)

$\varphi_{j}<\varphi_{j+1}<\varphi_{j}+\pi$ for $j=0,1,2,\ldots,N_{\gamma}-1$.

Note that $\varphi_{N_{\gamma}}=\varphi_{0}$. We denote the length of the $j$-th facet of $\mathcal{W}_{\gamma}$ by $\ell_{j}>0$.

We next prepare the notation of an evolving polygonal spiral. We denote an evolving polygonal spiral curve by (1) by $\Gamma_{D}(t)=\bigcup_{j=0}^{k}L_{j}(t)$. According to [9], we here consider the evolution of a positive convex polygonal spiral. Assume that the $j$-th facet $L_{j}(t)$ is given as

[TABLE]

with vertices $y_{j}(t)$ ($j=0,1,2,\ldots,k-1$) and the center $y_{k}(t)=O$. Assume that

[TABLE]

We have extended the number $j$ of ${\bf T}_{j}$ from $j=0.1,2,\ldots,N_{\gamma}-1$ to $\mathbb{Z}$; let ${\bf T}_{j+nN_{\gamma}}={\bf T}_{j}$ for $j=0,1,2,\ldots,N_{\gamma}-1$ and $n\in\mathbb{Z}$. Then, the evolution of $\Gamma_{D}(t)$ by (1) with fixed center $y_{k}(t)=O$ is expressed by an ODE system for $d_{j}(t)=|y_{j}(t)-y_{j-1}(t)|$ of the form

[TABLE]

where $b_{j}\in\mathbb{R}$ and $c^{\pm}_{j}>0$ are numerical constants defined by

[TABLE]

and $\beta_{j}=\beta({\bf N}_{j})$. Tracking the evolution of $\Gamma_{D}(t)$ is established by drawing $\Gamma_{D}(t)$ with setting

[TABLE]

See Figure 1 for details of $\Gamma_{D}(t)$ described with the above notations.

In this paper, we give an initial curve as $k=1$ with $d_{1}(0)=0$, i.e., $y_{1}(0)=y_{0}(0)=O$ and

[TABLE]

For evolution of a “spiral”, a new facet should be generated as the resultant of the evolution of present facets. Let $T_{1}=0$ and inductively set the generation time of facet $L_{k+1}(t)$ as

[TABLE]

When $t=T_{k+1}$, we add a new facet $L_{k+1}(T_{k+1})$ with $y_{k+1}(T_{k+1})=O$ and $d_{k+1}(T_{k+1})=0$. Then, change the spiral center to $y_{k+1}(t)$ from $y_{k}(t)$.

In summary, the algorithm of our discrete model for evolving polygonal spiral by (1) is as follows:

1. (I)

The generation time $T_{k}$ and curve $\Gamma_{D}(T_{k})=\bigcup_{j=0}^{k}L_{j}(T_{k})$ (with $d_{k}(T_{k})=0$) are given. 2. (II)

Solve (2)–(3) on $[T_{k},T_{k+1}]$ to obtain the evolution of $\Gamma_{D}(t)$. 3. (III)

When $t=T_{k+1}$, add a new facet $L_{k+1}(T_{k+1})$ with $y_{k+1}(T_{k+1})=O$ (then $d_{k+1}(T_{k+1})=0$) as the fixed center of $\Gamma_{D}(t)$. Then, return to (I).

The existence and uniqueness of solution to (2)–(3), the existence of the sequence $\{T_{k}\}_{k=1}^{\infty}$ of the generation times, $\lim_{k\to\infty}T_{k}=\infty$, and the intersection-free result of $\Gamma_{D}(t)$ are obtained by [9]; see it for details of the mathematical results.

2.2. Level set method

We recall the level set method [10] for a evolving spiral corresponding to the discrete model explained in the previous section.

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded domain with a smooth boundary. Consider the evolution of a single spiral by (1), and have set the center of a spiral at the origin. We give such a spiral curve and its direction of the evolution, which is denoted by ${\bf n}\in S^{1}$, with the level set method due to [10] as

[TABLE]

where $W=\{x\in\Omega;\ |x|>\rho\}$ for a constant $\rho>0$, and $\theta=\arg x$. According to [5], we obtain the anisotropic curvature $H_{\gamma}$ of $\Gamma_{L}(t)$ as

[TABLE]

with $\xi=D\gamma$ and $\gamma\in C^{2}(\mathbb{R}^{2}\setminus\{0\})$ satisfying (A1)–(A3). It is well-known that

[TABLE]

with $\gamma^{\circ}(p)=\sup\{p\cdot q;\ \gamma(q)\leq 1\}$, and $H_{\gamma}=1$ on $\mathcal{W}_{\gamma}$; see [2] for details. Moreover, from the context of derivation of (2)–(3) as in [9], one can find a self-similar solution with extension of $\mathcal{W}_{\gamma}$ for the motion of closed curve by $V=1$ ((1) with $U=1$ and $\rho_{c}=0$), which means that we measure the normal velocity with the Finsler metric

[TABLE]

Then, the normal velocity in this case should be given by

[TABLE]

since $\gamma(D\gamma^{\circ}(p))=1$ for $p\in\mathbb{R}^{2}\setminus\{0\}$ under some additional regularity and convexity assumptions on $\gamma$ and $\gamma^{\circ}$; see [2] for details.

As a boundary condition of the evolution with (1), we impose the right angle condition between $\Gamma_{L}(t)$ and $\partial W$. Then, the level set equation of the motion of spirals by (1) is of the form

[TABLE]

where $\vec{\nu}\in\mathbb{S}^{1}$ is the outer unit normal vector field of $\partial W$, and $\tilde{\beta}(p)=\beta(-p)$, $\tilde{\gamma}(p)=\gamma(-p)$ and $\tilde{\xi}(p)=\xi(-p)$. See [5] for details of the level set method.

Mathematical analysis for (6)–(7) with $\gamma\in C^{2}(\mathbb{R}^{2}\setminus\{0\})$ and $\beta\in C(\mathbb{R}^{2}\setminus\{0\})$ is established in [11]. For given initial data $u_{0}\in C(\overline{W})$, there exists a unique global viscosity solution $u\in C([0,\infty)\times\overline{W})$ to (6)–(7) with $u(0,\cdot)=u_{0}$. Moreover, the uniqueness of evolution of $\Gamma_{L}(t)$ is established in [6]; if there are continuous viscosity solutions $u$ and $v$ to (6)–(7) satisfying $\Gamma_{L}^{u}(0)=\Gamma_{L}^{v}(0)$ with the same orientations, then $\Gamma_{L}^{u}(t)=\Gamma_{L}^{v}(t)$ for $t>0$, where $\Gamma_{L}^{u}(t)=\{x\in\overline{W};\ u(t,x)-\theta(x)\equiv 0\mod 2\pi\mathbb{Z}\}$. Hence, we may give an arbitrary $u_{0}\in C(\overline{W})$ to obtain the motion of $\Gamma_{L}(t)$. In this paper, we give $u_{0}$ for (4) as $u_{0}\equiv\varphi_{0}$ due to [10].

Recall that we consider the situation such that $\mathcal{W}_{\gamma}$ is a convex polygon. The assumption (A4) is imposed for such a situation. Then, $\gamma$ is now given as

[TABLE]

with some $Q_{j}\subset\mathbb{R}^{2}$ for $j=0,1,2,\ldots,N_{\gamma}-1$, where

[TABLE]

for $Q\subset\mathbb{R}^{2}$. The crucial problem for solving (6)–(7) is how to treat $\mathrm{div}\tilde{\xi}(\nabla(u-\theta))$. For this problem, approximation of $\xi$ by the analogy of the stability result as in [4] is a simple option. From (8), we formally obtain

[TABLE]

so that we approximate $\chi_{Q_{j}}$ with the method as in [3] to remove the singularities. More precisely, we use the function

[TABLE]

with $\varepsilon\ll 1$ to approximate the sign function $z/|z|$. This function is also used in [9] when we approximate $\xi=D\gamma$ of $\gamma(p)=\|p\|_{1}=|p_{1}|+|p_{2}|$ or $\gamma(p)=\|p\|_{\infty}=\max\{|p_{1}|,|p_{2}|\}$ for $p=(p_{1},p_{2})$. In general, consider the case when $Q$ is given as a level set of a continuous function $f$, i.e., $Q=\{x\in\mathbb{R}^{2};\ f(x)>0\}$ and $\mathbb{R}^{2}\setminus\overline{Q}=\{x\in\mathbb{R}^{2};\ f(x)<0\}$. Then, we approximate $\chi_{Q}$ by

[TABLE]

for a suitable parameter $(p_{1},p_{2})$. (We often choose $(p_{1},p_{2})=\nabla f$ like as in [3], or $(p_{1},p_{2})=(1,0)$ for simplicity.) Hence, we obtain the approximation

[TABLE]

by a level set functions $f_{j}\in C(\mathbb{R}^{2})$ for $Q_{j}$.

3. Measuring difference

3.1. Crystalline energy density

Let us consider the situation such that the Wulff shape $\mathcal{W}_{\gamma}$ and a support function $\gamma^{\circ}\colon\mathbb{R}^{2}\to[0,\infty)$ satisfying

[TABLE]

are given. Note that $\gamma^{\circ}=\sup\{p\cdot q;\ \gamma(q)\leq 1\}$ is a convex and positively homogeneous of degree 1. According to these facts and that $\mathcal{W}_{\gamma}$ is a convex polygon, we assume that $\gamma^{\circ}$ is given as

[TABLE]

with $\eta_{j}>0$ and $\psi_{j}\in\mathbb{R}$. Assume that

1. ($\gamma$1)

$\psi_{0}<\psi_{1}<\psi_{2}\cdots<\psi_{N_{\gamma}-1}<\psi_{0}+2\pi$, 2. ($\gamma$2)

$\psi_{j}<\psi_{j+1}<\psi_{j}+\pi$ for $j=0,1,2,\ldots,N_{\gamma}-1$. 3. ($\gamma$3)

$P_{j}=\{p\in\mathbb{R}^{2};\ m_{j}\cdot p\geq m_{k}\cdot p\ \mbox{for}\ k=0,1,2,\ldots,N_{\gamma}-1\}=\Xi_{j,j-1}\cap\Xi_{j,j+1}\neq\emptyset$ for $j=0,1,2,\ldots,N_{\gamma}-1$, where $\Xi_{j,k}=\{p\in\mathbb{R}^{2};\ m_{j}\cdot p\geq m_{k}\cdot p\}$.

(Note that $\psi_{j+nN_{\gamma}}=\psi_{j}$ for $n\in\mathbb{Z}$.) We now propose a practical way to reconstruct a convex and piecewise linear $\gamma\colon[0,\infty)\to[0,\infty)$ from the above settings. Note that we do not impose the normalizing assumption (5) in this section.

We first remark that

[TABLE]

when $\gamma$ is convex, since $\gamma=(\gamma^{\circ})^{\circ}$. Then, by ($\gamma$3), we find $\theta_{0},\theta_{1},\theta_{2},\ldots,\theta_{N_{\gamma}-1}$ such that $\theta_{j}<\theta_{j+1}<\theta_{j}+2\pi$ and

[TABLE]

for $j=0,1,2,\ldots,N_{\gamma}-1$. It should be calculated by

[TABLE]

In fact, by (9) we have

[TABLE]

Let $a_{j}$, $b_{j}$ be constants defined by

[TABLE]

and $c_{j}$ be a constant satisfying

[TABLE]

Then, (10) yields that

[TABLE]

Let us consider the formula $p\cdot q\leq\gamma^{\circ}(q)$ with $q=(\cos\theta,\sin\theta)$ and $p=(x,y)$. If $\theta\in[\theta_{j},\theta_{j+1}]$, then we observe that

[TABLE]

Then, one can find that

[TABLE]

where

[TABLE]

Moreover, one can find $\Pi_{j,j+1}=\Pi_{j+1,j+1}$. In fact, by definition of $\Pi_{j,j+1}$ and (11), we observe that

[TABLE]

which implies $\Pi_{j,j+1}=\Pi_{j+1,j+1}$. Hence, we obtain

[TABLE]

Set $r_{j}=\left[\eta_{j}(\cos\theta_{j}\cos\psi_{j}+\sin\theta_{j}\sin\psi_{j})\right]^{-1}=\left[\eta_{j}\cos(\theta_{j}-\psi_{j})\right]^{-1}$, and

[TABLE]

Then, we observe that

[TABLE]

which implies that

[TABLE]

Hence, we observe that $\gamma(p)=\max_{0\leq j\leq N_{\gamma}-1}n_{j}\cdot p$.

Summary. Assume that ($\gamma$1)–($\gamma$3) hold. Let $\gamma^{\circ}\colon\mathbb{R}^{2}\to[0,\infty)$ be given as

[TABLE]

Set

• •

$\theta_{j}=c_{j}+\pi/2$ with $c_{j}\in\mathbb{R}$ such that

[TABLE]

• •

$r_{j}=\left[\eta_{j}\cos(\theta_{j}-\psi_{j})\right]^{-1}$

for $j=0,1,2,\ldots,N_{\gamma}-1$. Then,

[TABLE]

Remark 1**.**

1. (i)

When we give only the parameters of $\ell_{j}$ and $\psi_{j}$ for $\mathcal{W}_{\gamma}$, then we have to set the location of the origin $O\in\mathcal{W}_{\gamma}$ to determine $\gamma^{\circ}$. Note that $\eta_{j}$ depends on the location of the origin in $\mathcal{W}_{\gamma}$. 2. (ii)

There is a case that $P_{j}=\emptyset$ and thus $P_{k}\neq\Xi_{k,k-1}\cap\Xi_{k,k+1}$ for some $j,k\in\{0,1,2,\ldots,N_{\gamma}-1\}$ when $N_{\gamma}\geq 4$, even if $\psi_{j}$ satisfies ($\gamma$1). In fact, $\gamma^{\circ}(p)=\max_{0\leq j\leq 3}n_{j}\cdot p$ with

[TABLE]

implies that $P_{1}=\emptyset$. 3. (iii)

Notice that the above way also can be applied to construct $\gamma^{\circ}$ from a given $\gamma$.

3.2. Difference function

Once we obtain $\Gamma_{D}(t)=\bigcup_{j=0}^{k}L_{j}(t)$ or $\Gamma_{L}(t)$, then we compare with $\Gamma_{D}(t)$ and $\Gamma_{L}(t)$ by calculating the measure of the region enclosed by $\Gamma_{D}(t)$ and $\Gamma_{L}(t)$. It is established as follows; we construct the height functions

[TABLE]

of the stepwise surface at $\Gamma_{D}(t)$ or $\Gamma_{L}(t)$ with step height $h_{0}=1$, respectively. Note that $\theta_{D}(t,x)$ or $\theta_{L}(t,x)$ is a branch of $\theta(x)$ whose discontinuity is only on $\Gamma_{D}(t)$ or $\Gamma_{L}(t)$, respectively. According to [10], the practical way to construct $\theta_{L}(t,x)$ is given in [10]. Hence, we here give a practical way to construct $\theta_{D}(t,x)$.

1. (i)

We first pick up the rotation number $n\in\mathbb{N}$ for the facet number $k\in\mathbb{N}$ of $\Gamma(t)=\bigcup_{j=0}^{k}L_{j}(t)$, i.e., $k=\bar{k}+nN_{\gamma}$ with $\bar{k}\in\{0,1,2,\ldots,N_{\gamma}-1\}$. 2. (ii)

Then, we now set

[TABLE]

a branch of $\arg x$ whose discontinuity is only on

[TABLE]

(See Figure 2(2).) 3. (iii)

Let us set

[TABLE]

(gray regions in Figure 2(3)). To remove a discontinuity on a dash line in $\partial R_{k-1}(t)$, we set

[TABLE] 4. (iv)

We inductively set

[TABLE]

to remove illegal discontinuities of $\Theta_{\ell-1}$ from $\ell=1$ to $\ell=k$, where

[TABLE]

for $j=0,1,\ldots,k-1$ (see Figure 2(4) for $R_{k-2}(t)$).

Consequently, we set

[TABLE]

Hence, we can define the difference $\mathcal{D}(t)$ between $\Gamma_{D}(t)$ and $\Gamma_{L}(t)$ as

[TABLE]

4. Numerical results

In this section, we present some numerical simulations measuring the difference between $\Gamma_{D}(t)$ and $\Gamma_{L}(t)$ evolving by

[TABLE]

i.e., (1) with $\beta\equiv 1$ and $U=1$ for some kinds of $\gamma$. The initial curve is chosen as

[TABLE]

and then (4) for the discrete model, and $u(0,x)=\varphi_{0}$ for the level set method. Throughout this section, we set

[TABLE]

for some $s\in\mathbb{N}$, and then $\Delta x=0.02/s$. In the following subsections, we will choose time intervals of the numerical simulations so that the curves $\Gamma_{D}(t)$ does not touch to the outer boundary $\partial\Omega$. In other words, we avoid the situation that the boundary condition on $\partial\Omega$ makes difference between $\Gamma_{D}(t)$ and $\Gamma_{L}(t)$. Note that, however, the difference of the boundary condition at the center and the evolution law of the first facet $L_{0}(t)$ are still remains.

We calculate the ODE system (2)–(3) by 4-th order Runge-Kutta method with the time span $\Delta t=10^{-6}$. From these numerical results, we construct $h_{D}(t,x)$ on each numerical mesh $D_{s}$ to compare the results with those from the level set method. On the other hand, the level set equation (6)–(7) is calculated by the explicit finite difference scheme as in [10] with the time span $\Delta t=0.1\times\Delta x^{2}$. See also [10] for the way to construct $h_{L}(t,x)$ with the step height $h_{0}=1$. To draw a graph of $\mathcal{D}(t)$, we pick up the data $\mathcal{D}(t_{k})=\mathcal{D}(kT/20)$ for $0\leq k\leq 20$ on the calculating time interval $[0,T]$.

We now recall the difference between the discrete model in §2.1 and the level set method in §2.2.

1. (i)

The domain of the level set method has a “center” $B_{\rho}=\{x\in\mathbb{R}^{2};\ |x|\leq\rho\}$ with a finite radius $\rho>0$. However, the discrete has the center at the origin as a point (null set). 2. (ii)

The boundary conditions are different:

• •

[Discrete model] $L_{0}(t)$ evolves by $V=1$ since $d_{0}(t)=\infty$. On the other hand, the behavior of the facets associated with center is imposed with fixing and the generation rule of new facets.

• •

[Level set method] The right angle conditions, in particular, $\Gamma_{L}(t)\perp\partial B_{\rho}$ and $\Gamma_{L}(t)\perp\partial\Omega$ are imposed by (7).

Because of the above differences, we have no conjectures of convergence between $\Gamma_{L}(t)$ and $\Gamma_{D}(t)$ now. Moreover, from the numerical results of the isotropic case in [10, 12], not only tending the approximation parameters to zero but also letting $\rho\to 0$ is required for numerical accuracy. Thus, we shall check the numerical results with fixed radius $0<\rho\ll 1$ and reducing radius $\rho=O(\Delta x)$.

4.1. Square spiral

The first examination is the square spiral case, i.e.,

[TABLE]

Thus, we define the parameters of $\mathcal{W}_{\gamma}$ for the discrete model as

[TABLE]

For the level set equation, since $\gamma^{\circ}(p)=\max\{|p_{1}|,|p_{2}|\}$ for $p=(p_{1},p_{2})$, we observe that

[TABLE]

We calculate the ODE system (2)–(3) and the level set equation (6)–(7) for

[TABLE]

on the time interval $[0,1]$. See [9, §4] for details to approximate $\xi$ of the above $\gamma$. Figure 3 are profiles of the diagonal spiral at $t=1$ with the above setting. Note that, in this and following sections, the profile of spirals by the level set method is calculated with $\rho=0.02-10^{-8}$ and $\Delta x=0.0050$ ($s=4$).

The left figure of Figure 4 presents the graph of $\mathcal{D}(t)$ for $s=2,3,4,5,6$ with a fixed center radius $\rho=0.02-10^{-8}$. One can find that the differences are less than $4\%$ of the area $|W|$ for all cases, although the value of $\mathcal{D}(t)$ becomes worse when we choose smaller $\Delta x$. The best one is the case with $\Delta x=0.010$ ($s=2$).

On the other hand, we obtain better results when $\rho=O(\Delta x)$. The right figure of Figure 4 presents the graph of $\mathcal{D}(t)$ for $s=2,3,4,5,6$ with the center size $\rho=(2-10^{-8})\Delta x$, i.e., the setting $\rho\to 0$ as $\Delta x\to 0$. Note that the cases of $\Delta x=0.010$ ($s=2$) in both figures of Figure 4 is the same. One can find that the differences are less than 2.5% of the area $|W|$ for all cases, and $\mathcal{D}(t)$ of the cases with $s\geq 3$ are smaller than that of $s=2$, although the smallest $\mathcal{D}(t)$ is the case $\Delta x=0.0067$ ($s=3$).

4.2. Diagonal spiral

The second examination is the diagonal spiral case, i.e., the $\pi/4$ rotation of the first case; and then

[TABLE]

for the discrete model. In this case, one can find that

[TABLE]

and thus

[TABLE]

For the level set equation, we set

[TABLE]

According to [9], it is represented as

[TABLE]

and thus

[TABLE]

See [9] for the approximation of the above $\xi$. We calculate the ODE system (2)–(3) and the level set equation (6)–(7) for (13) on the time interval $[0,1]$. Figure 5 are profiles of the diagonal spiral at $t=1$ with the above setting.

The left figure of Figure 6 is a graphs of $\mathcal{D}(t)$ for $s=2,3,4,5,6$ with a fixed center radius $\rho=0.02-10^{-8}$. One can find that $\mathcal{D}(t)$ is reduced by choosing smaller $\Delta x$, and the smallest $\mathcal{D}(t)$ is the case $\Delta x=0.0033$ ($s=6$). Our numerical simulations show that the differences are less than 4% of $|W|$ if $s\geq 4$. Note that $\rho\approx 4\Delta x$ when $s=4$.

Because of the above results, we choose $\rho\approx 4\Delta x$ for accurate simulations with a reduced center radius $\rho=O(\Delta x)$. The right figure of Figure 6 presents graphs of $\mathcal{D}(t)$ for $s=2,3,4,5,6$ with $\rho=(4-10^{-8})\Delta x$. One can find that the differences are less than $5\%$ for all cases, although the worst one is that with $\Delta x=0.0033$ ($s=6$). Note that the cases of $\Delta x=0.0050$ ($s=4$) in both figures of Figure 6 are the same.

4.3. Triangle spiral

Finally, we examine a triangle spiral as an asymmetric case of $\gamma$ or $\gamma^{\circ}$. To give its settings, we first give $\gamma^{\circ}$. Because of the normalizing assumption (5), we set

[TABLE]

Then, $\mathcal{W}_{\gamma}=\{p\in\mathbb{R}^{2};\ \gamma^{\circ}(p)\leq 1\}$ implies that

[TABLE]

since $\mathcal{W}_{\gamma}$ is an equilateral triangle whose vertices are at $(1,\pm\sqrt{3})$ and $(-2,0)$. On the other hand, from the computation as in §3.1 we obtain

[TABLE]

Note that $Q_{j}$ in (8) is given as

[TABLE]

Then, we obtain

[TABLE]

where $f_{j}(p)=\min_{k\neq j}(g_{j}(p)-g_{k}(p))$. We calculate (2)–(3) or (6)–(7) for the evolution equation

[TABLE]

on the time interval $[0,0.8]$ with the above anisotropic setting. This time interval is chosen so that $\Gamma_{D}(t)$ does not touch to $\partial\Omega$ for $t\in[0,0.8]$. Figure 7 are profiles of the diagonal spiral at $t=0.8$ with the above setting.

The left figure of Figure 8 presents graphs of $\mathcal{D}(t)$ with $s=2,3,4,5,6$ with a fixed center radius $\rho=0.02-10^{-8}$. One can find that the differences are less than $4\%$ for the all cases except $\Delta x=0.0067$ ($s=3$), and the best one is that with $\Delta x=0.0040$ ($s=5$).

From the analogy of the diagonal spiral case, we choose $\rho\approx 4\Delta x$ as a reducing center radius $\rho=O(\Delta x)$. The right figure of Figure 8 presents graphs of $\mathcal{D}(t)$ with $s=2,3,4,5,6$ with a fixed center radius $\rho=(4-10^{-8})\Delta x$. One can find that the differences are less than $5\%$ when $\Delta x\leq 0.0050$ ($s\geq 4$). Note that the cases of $\Delta x=0.0050$($s=4$) in both figure of Figure 8 are the same.

5. Conclusion

In this paper, we compared the discrete model as in [9] and the level set method as in [10] for evolving spirals by the crystalline eikonal-curvature flow (1). Note that the level set equation includes the derivative of a piecewise linear energy density function. For this problem, we introduced an approximation of level set equation for the crystalline curvature flow, which is established with the approximation of the characteristic function as in [3]. To measure the difference between the two curves obtained by the discrete model and the level set method, we introduced an area difference function defined by (12). It is consist of the $L^{1}$ difference of the height function as in [10] with the step-height $h_{0}=1$. Note that the discrete and level set models are slightly different on the boundary condition at the center and the outer boundary of the domain. However, we found that the area differences of these models are less than $5\%$ of the area of the domain for square, diagonal and triangle spirals as in §4 when the resolution of the numerical lattice is enough high and the radius of the center for the level set method is suitably small.