# Numerical methods for motion of level sets by affine curvature

Numerical methods for motion of level sets by affine curvature Abstract We study numerical methods for the nonlinear partial differential equation that governs the motion of level sets by affine curvature. We show that standard finite difference schemes are nonlinearly unstable. We build convergent finite difference schemes using the theory of viscosity solutions. We demonstrate that our approximate solutions capture the affine invariance and morphological properties of the evolution. Numerical experiments demonstrate the accuracy and stability of the discretization. 1. Introduction 1.1 The affine curvature partial differential equation The affine curvature evolution is one of the most fundamental geometric evolution equations after the mean curvature evolution. It was introduced by Sapiro & Tannenbaum (1994) and Angenent et al. (1998), and has applications in mathematical morphology, edge detection, image smoothing and image enhancement (see Sapiro, 2006). We are motivated by the recent work of Calder & Smart (2016), which provides an application of the affine curvature partial differential equation (PDE) to the statistics of large data sets. The convex hull peeling algorithm (Chazelle, 1985) provides an affine invariant notion of the median and the quantiles of multidimensional probability distributions. Although there is more than one way to measure data depth (Barnett, 1976), affine invariance is an important property for such measures (Liu et al., 1999). The level sets of the solution of the affine curvature PDE, with the right-hand side given by the probability density $$\rho$$, also gives an affine invariant notion of the depth of $$\rho$$. According to Calder & Smart (2016), these two notions of depth are equivalent: the rescaled data depth layers of $$N$$ data points sampled from the density, $$\rho$$, given by convex hull peeling algorithm converge, in the limit $$N \to \infty$$, to the levels given by the solution of the PDE. Compared with convex hull peeling, the PDE characterization is efficient in terms of the number of data points $$N$$: an efficient density estimation method can be used to approximate $$\rho$$, and afterward the PDE solve does not depend on $$N$$. This kind of limiting PDE approach has already been shown to be effective for nondominated sorting (Calder et al., 2014). The planar motion of level sets by affine curvature is governed by the nonlinear PDE $$\label{ACPDE} u_t = F[u] := \left\vert{\nabla u}\right\vert \left (k[u]\right)^{1/3}.$$ (AC) Here $$u = u(x,y):\mathbb{R}^2 \to \mathbb{R}$$, $$\nabla u = (u_x,u_y)$$ denotes the gradient of $$u$$, and $$k[u]$$ denotes the curvature of the level set of $$u$$ $$\label{curvature} k[u] = \text{div}\left(\frac{\nabla u}{\left\vert{\nabla u}\right\vert}\right) = \frac{u_{xx}u_y^2-2u_x u_y u_{xy}+u_{yy}u_x^2}{(u_x^2+u_y^2)^{3/2}}.$$ (1.1) The affine curvature PDE is closely related to the well-known PDE for motion of level sets by mean curvature $$u_t = {\it {\Delta}}_1 u := \left\vert{\nabla u}\right\vert k[u]$$ (MC) studied in the seminal article Osher & Sethian (1988). However, the PDE (AC) exhibits instabilities not found in the mean curvature PDE (MC), as demonstrated below. To resolve these instabilities, we propose a Lipschitz regularization of the operator. The regularized operator is also a geometric PDE, and viscosity solutions converge to solutions of the affine curvature PDE in the limit as the regularization parameter goes to zero [Giga, 2006, Theorem 4.6.1]. The advantage of the regularized operator is that it allows us to build stable, convergent explicit solvers that are not otherwise available. Moreover, the resulting discretization can be combined into a filtered scheme that achieves the higher accuracy of the otherwise unstable standard finite difference scheme. In addition, the numerical solutions exhibit the affine invariance and morphology properties of the evolution. Our approach to the convergent discretization In this paragraph, we present an overview of our approach to building an elliptic discretization for clarity. The details and supporting theory can be found in the sections that follow. Elliptic discretizations are available for $${\it {\Delta}}_1 u$$, rather than for $$k[u]$$, so we rewrite \begin{align}\label{AffPQ} F[u] & = A(\left\vert{\nabla u}\right\vert,{\it {\Delta}}_1 u) = (\left\vert{\nabla u}\right\vert^2 {\it {\Delta}}_1 u)^{1/3}, && \text{where} A(p,q) = \left(p^2 q \right)^{1/3}. \end{align} (1.2) The goal is to make use of available elliptic discretizations of $$\pm \left\vert{\nabla u}\right\vert$$ and $${\it {\Delta}}_1 u$$ to build an elliptic discretization of the full operator $$F[u]$$. However, simply inserting these operators into the function $$A(p,q)$$ is not sufficient: elliptic schemes are built by composing nondecreasing maps with elliptic operators, and $$A(p,q)$$ fails to be nondecreasing. Furthermore, explicit time discretizations require Lipschitz continuous operators, and $$A(p,q)$$ also fails to be Lipschitz continuous. Because the properties of the nonlinear and singular function $$A(p,q)$$ are so important, we first study a model equation in one dimension. Define $$\label{AffPQ1D} F^{1D}[u] := A(u_x, u_{xx}) = \left(\left\vert{u_x}\right\vert^2 \:u_{xx}\right)^{1/3},$$ (1.3) so that the $$F^{1D}[u]$$ operator has the same homogeneity in the first and second derivatives as the $$F[u]$$ operator. Similar to the higher dimensional PDE, the model equation exhibits the instability of standard finite differences. To build an elliptic scheme for the one-dimensional equation, what is needed is a nondecreasing representation of the function $$A(p,q)$$, which is consistent with $$F^{1D}[u]$$. Furthermore, we need a Lipschitz continuous approximation of $$A(p,q)$$, with Lipchitz constant $$K^h$$, to build a monotone discretization of the time-dependent PDE, using a time step $${\rm d}t \le 1/K^h$$. Once this modified function is available, we proceed by inserting the discretization of the two-dimensional operators into the modified function, which results in a convergent scheme. Related numerical work There are by now a large number of numerical methods for the level set mean curvature PDE (MC) (see review papers Deckelnick et al., 2005; Ciomaga et al., 2011). In particular, Catté and Dibos proposed a morphological scheme that satisfies the comparison principle (Cattè et al., 1995) and the first author presented a convergent wide stencil finite difference scheme (Oberman, 2004) using a median formula. For the affine curvature evolution, the recent article Esedoglu et al. (2010) gives a Bence–Merriman–Osher (Merriman et al., 1994) thresholding scheme. Similar to the scheme proposed here, it is monotone, but introduced a different regularization of the cube root, which was needed for theoretical purposes, but not in practice. Thresholding methods are effective for moving a given curve by the evolution and allow for large time steps to be taken. Alvarez and Guichard proposed a local scheme that lacks the affine invariance property (Guichard, 1994). A morphological scheme that generalized Cattè et al. (1995) was proposed by Guichard and Morel for affine curvature (Guichard & Morel, 1997). This inf–sup scheme, although morphologically invariant, has some limitations on the speed at which the level set curves move. In Moisan (1998), a nonlocal geometric morphological scheme is presented for the evolution of a single curve. These are approximated by polygons, and the scheme is shown to be monotone and affine invariant for polygons. Novelty of the work The scheme presented here uses the level set representation (AC). Thus, it moves every level set by the evolution, unlike thresholding methods, but requires a much smaller time step. It is monotone, a desirable property for the convergence, but also for its applications. Moreover, it is shown numerically that it captures the affine invariance and morphological properties. 1.2 Background mathematics Euclidean and Affine Curvature We begin with a parametric description of the affine curvature evolution and make a comparison with the more familiar mean curvature evolution. We refer to Sapiro (2006, Chapter 2) for more details. Consider a curve described parametrically $$\mathcal{C}(s):[a,b] \in \mathbb{R} \to \mathbb{R}^2$$, where $$s$$ parameterizes the curve. If the curve is parameterized by the Euclidean arc length, $$\left\vert{\frac{\text{d} \mathcal{C}} {\text{d} s}}\right\vert = 1,$$ then the curvature, $$k$$, is defined, up to a sign, by the equation $$\mathcal{C}_{ss} = k\mathbf{n}$$, where $$\mathbf{n}$$ denotes the Euclidean normal of the curve. Under mean curvature, the curve evolves with velocity $$k\mathbf{n}$$. The affine curvature arises from a different parameterization of the curve. Define the parameter, $${r}$$ by the condition that the vectors $$C_{r}$$ and $$C_{{r} {r}}$$ form a parallelogram of area $$1$$, $$[\mathcal{C}_{r},\mathcal{C}_{{r} {r}}] = 1$$. (Here the brackets denote the determinant of the matrix whose columns are given by those vectors.) Then the affine curvature of the curve is defined by $$\mu = [\mathcal{C}_{{r} {r}},\mathcal{C}_{{r} {r} {r}}]$$. Affine differential geometry is not defined for nonconvex curves. However, we can still define the affine curvature evolution using the Euclidean curvature by taking the velocity to be $$k^{1/3}\mathbf{n}$$. The affine curvature is the simplest nontrivial affine invariant of the curve (Su, 1983). Under the affine curvature evolution, any convex curve remains convex; any smooth curve remains smooth; any convex smooth curve evolves to an ellipse until it collapses to a single point and any smooth curve becomes convex after a certain time. Moreover, the affine curvature evolution is invariant under the class of special affine transformations, which are defined by matrices with determinant $$1$$ (see Theorem 2.10 below). When compared with mean curvature evolution, this shrinks curves to circles (Gage, 1984) and is invariant under the smaller class of orthogonal transformations. Level Set PDE formulation The level set method (Sethian, 1999; Osher & Fedkiw, 2003) for the affine curvature evolution results in the PDE (AC). The level set method has the following advantages compared with parameterized curve evolution: (i) it provides a natural generalization of the flows when the curve becomes singular and notions such as normals are not well defined; (ii) there is no need to track topological changes, because they are discovered when the corresponding level set is computed; and (iii) it can be discretized on a uniform grid, which is convenient for many applications. In the level set method, a curve is represented implicitly as the level set of the auxiliary function $$u(x,y,t):\mathbb{R}^2\times \mathbb{R} \to\mathbb{R}$$, that is, $\mathcal{C}(t) = \left\{(x,y)\in \mathbb{R}^2 \mid u(x,y,t) = c\right\}$ for some arbitrary constant $$c\in \mathbb{R}$$. If $$u$$ satisfies $$u_t = \left\vert{\nabla u}\right\vert \beta[u(\cdot,t)]$$, for some function $$\beta$$, which depends on the level set of $$u$$, then all its level sets move in the normal direction with speed $$\beta$$. For example, choosing $$\beta = 1$$, we obtain the time-dependent eikonal equation, $$u_t = \left\vert{\nabla u}\right\vert$$. Taking $$\beta = k[u]$$ or $$\left(k[u]\right)^{1/3}$$, we obtain (MC) and (AC), respectively. 2. The affine curvature PDE 2.1 Definition of viscosity solutions Viscosity solutions (Crandall et al., 1992) provide the correct notion of weak solution to a class of degenerate elliptic PDEs, which includes (AC) and (MC). The article Evans & Spruck (1999) focuses on the mean curvature equation. The book Giga (2006) established existence and uniqueness of viscosity solutions in a bounded domain, with Neumann boundary conditions, see Section 3.6 for the Comparison Principle, and Theorem 4.6.1 for the convergence of approximations. It is also shown that the Lipschitz constant of the solution is preserved by the evolution (Section 3.5). See also the book Cao (2003) for viscosity solutions for geometric evolution equations and numerical methods. Definition 2.1 The function $$F: \mathcal{S}^{n} \times \mathbb{R}^n \times \mathbb{R} \times {\it {\Omega}} : \to \mathbb{R}$$ is proper and degenerate elliptic (in the sense of Crandall et al., 1992) if $F(M,p,r,x) \le F(N,p,s,x), \quad \text {for all} M \succeq N, r \le s, \text{and all} x \in {\it {\Omega}}, p \in \mathbb{R}^n,$ where $$Y \preceq X$$ if $$d^\intercal Yd \le d^\intercal X d$$ for all $$d\in \mathbb{R}^n$$. Remark 2.2 By this convention, the Laplacian operator is elliptic when written $$F(M) = -\text{tr}(M)$$. Some authors use the other convention, without the minus sign. Lemma 2.3 The operator $$-F[u]$$ is degenerate elliptic. Proof. We start with the observation that using (1.1), we can write ${\it {\Delta}}_1 u = u_{\mathbf{t}\mathbf{t}}, \qquad \mathbf{t} = \frac{(-u_y,u_x)}{(u_x^2+u_y^2)^{1/2}},$ where $$\mathbf{t}$$ is the (Euclidean) unit tangent. Then, it is clear that $$-{\it {\Delta}}_1$$ is degenerate elliptic. Using the representation (1.2) the degenerate ellipticity of $$-F$$ follows. □ We now give the definition of viscosity solutions on a domain $${\it {\Omega}} \subset \mathbb{R}^2$$. We follow Giga (2006). Let $$T > 0$$. We are interested in the Cauchy problem $$\label{cauchyproblem} \begin{cases} u_t = F[u] & \text{in} (x,t) \in {\it {\Omega}} \times (0,T),\\ B(x,p) := \nu(x)^\intercal p = 0 & \text{on} (x,t) \in \partial {\it {\Omega}} \times (0,T),\\ u(x,0) = u_0(x) & \text{in} x \in \overline{{\it {\Omega}}}, \end{cases}$$ (2.1) where $$\nu$$ is the outward unit normal of $$\partial {\it {\Omega}}$$. For completeness, we recall the definition of upper and lower semicontinuous functions. Definition 2.4 A function $$u:U\to\mathbb{R}$$ is upper (resp. lower) semicontinuous at $$x\in U$$ provided $\limsup_{y \in U\to x} u(y) \leq u(x) \quad \left(\text{resp.} \liminf_{y \in U\to x} u(y) \geq u(x)\right)\!.$ We denote by $$USC(U)$$ (resp. $$LSC(U)$$) the collection of functions that are upper (resp. lower) semicontinuous at all points of $$U$$. Definition 2.5 A function $$u \in USC(\overline{{\it {\Omega}}} \times [0,T])$$ is called a viscosity subsolution of (2.1) if $$u(x,0) \leq u_0(x)$$ for $$x \in \overline{{\it {\Omega}}}$$ and for any $$\varphi \in C^2(\overline{{\it {\Omega}}}\times[0,T])$$ such that $$u-\varphi$$ has a local maximum at $$(x,t) \in \overline{{\it {\Omega}}} \times (0,T)$$ then $\begin{cases} \varphi_t(x,t) - F[\varphi](x,t) \leq 0 & \text{if}~ x \in {\it {\Omega}},\\ \min\{\varphi_t(x,t) - F[\varphi](x,t), B(x,\nabla \phi(x,t))\} \leq 0 & \text{if}~ \in \partial {\it {\Omega}}. \end{cases}$ Similarly, $$u \in LSC(\overline{{\it {\Omega}}} \times [0,T])$$ is called a viscosity supersolution of (2.1) if $$u(x,0) \geq u_0(x)$$ for $$x \in \overline{{\it {\Omega}}}$$ and for any $$\varphi \in C^2(\overline{{\it {\Omega}}}\times[0,T])$$ such that $$u-\varphi$$ has a local minimum at $$(x,t) \in \overline{{\it {\Omega}}} \times (0,T)$$ then $\begin{cases} \varphi_t(x,t) - F[\varphi](x,t) \geq 0 & \text{if}~ x \in {\it {\Omega}},\\ \max\{\varphi_t(x,t) - F[\varphi](x,t), B(x,\nabla \phi(x,t))\} \geq 0 & \text{if}~ \in \partial {\it {\Omega}}. \end{cases}$ Finally, we call $$u$$ a viscosity solution of (2.1) if $$u$$ is both a viscosity subsolution and a viscosity supersolution. We have the following uniqueness result, from Giga (2006). Theorem 2.6 (Comparison Principle) Suppose that $${\it {\Omega}}$$ is convex with $$C^2$$ boundary $$\partial {\it {\Omega}}$$. Let $$u$$ and $$v$$ be a viscosity subsolution and supersolution of (2.1). Then $$u \leq v$$ in $${\it {\Omega}} \times (0,T)$$. Remark 2.7 We will also be interested in the static equation $$F[u] = f$$. The definition of viscosity solution and comparison principle presented here for the parabolic equation generalize to the static equation. 2.2 Invariance properties of the PDE We present some properties of (AC), starting with some properties of the affine curvature operator $$F[u]$$. The proofs can be found in Oberman & Salvador (2016). Lemma 2.8 Let $$u \in C^2(\mathbb{R}^2)$$. Then the operator $$F$$ has the following properties: (1) Rescaling: for $$h > 0$$, define $$v(x,y) := u(x/h,y/h)$$ $F[v] = h^{-4/3}F[u].$ (2) Morphology: for $$g \in C^1(\mathbb{R})$$, $F[g \circ u] = g^\prime(u) F[u].$ (3) Affine invariance: for any affine map $$\phi(\mathbf{x}) = A\mathbf{x}+\textbf{b}$$, $F[u \circ \phi] = (\det A)^{2/3} F[u] \circ \phi.$ Remark 2.9 For $${\it {\Delta}}_1 u$$, rescaling gives $${\it {\Delta}}_1 v = h^{-2} {\it {\Delta}}_1 u$$ and invariance only holds for orthogonal transformations. Using the Lemma, the properties can be generalized for the affine curvature evolution (AC). Theorem 2.10 Consider $${\it {\Phi}}_t:u_0\mapsto u(\cdot,t)$$ the solution map of (AC). Then $${\it {\Phi}}_t$$ satisfies the following properties: (1) Monotonicity$$:$$$$u \leq v \Rightarrow {\it {\Phi}}_t(u) \leq {\it {\Phi}}_t(v)$$. (2) Morphology/Relabelling: for any monotone scalar function $$g$$, ${\it {\Phi}}_t(g \circ u) = g \circ {\it {\Phi}}_t(u).$ (3) Affine invariance: for any affine map $$\phi(\mathbf{x}) = A\mathbf{x}+\textbf{b}$$, ${\it {\Phi}}_t(u\circ \phi) = \left({\it {\Phi}}_{t(\det A)^{2/3}}(u)\right) \circ \phi.$ Remark 2.11 Note that our rescaling factor property $$(iii)$$ differs from the one in Moisan (1998); however, both formulas agree (and give unity) for special affine transformations that have determinant $$1$$. 2.3 An exact solution We now give an example of an exact solution for (AC). A proof can be found in Oberman & Salvador (2016). See also Giga (2006, Section 1.7.4). Lemma 2.12 Given $$a,b > 0$$, let $$u:\mathbb{R}^2\times [0,\infty)$$ be given by the shifted ellipse $u(x,y,t) = t + \frac{3}{4}(ab)^{2/3}\left(\left(\frac{x}{a}\right)^2+\left(\frac{y}{b}\right)^2\right)^{2/3} = t + \frac{3}{4}\left(\frac{b}{a}x^2+\frac{a}{b}y^2\right)^{2/3}.$ Then $$u$$ is a classical solution of (AC). Moreover, we conclude that ellipses remain ellipses of fixed eccentricity under the motion by affine curvature. 3. Finite difference equations and convergence 3.1 Elliptic finite difference schemes The framework developed in Barles & Souganidis (1991) provides conditions for when approximation schemes converge to the unique viscosity solution of a PDE. The article Oberman (2006) defines elliptic finite difference schemes, which are monotone and stable. Theorem 3.1 Consider an elliptic PDE that satisfies a comparison principle for viscosity solutions. A consistent, stable and monotone approximation scheme converges locally uniformly to the (unique) viscosity solution. Consider the domain $${\it {\Omega}} \subset D$$ where $$D$$ is an $$n$$-cube. We use a uniform grid of spacing, $$h$$, in $$D$$, and define the finite difference grid, $$G^h = \{x \in h \mathbb{Z}^n ~\mid ~ x\in D \}$$. Definition 3.2 (Elliptic finite difference equation) Let $$C(G^h)$$ denote the set of grid functions, $$u: G^h \to \mathbb{R}$$. A finite difference operator is a map $$F^h: C(G^h) \to C(G^h)$$, which has the following form, $F^h[u](x) = F^h(x, u(x), u(x) - u(\cdot)),$ where $$u(\cdot)$$ indicates the values of the grid function $$u$$. It has stencil width$$W$$ if $$F^h(x, u(x), u(x) - u(\cdot))$$ depends only on values $$u(y)$$ for $$\Vert{y-x}\Vert_\infty/h \le W$$. A solution of the finite difference scheme is a grid function that satisfies the equation $$F^h[u](x) = 0$$ for all $$x \in G^h$$. The finite difference operator is elliptic if $r \le s, ~ v(\cdot) \le w(\cdot) \implies {F^h}\left(x,r,v(\cdot)\right) \leq {F^h}(x,s,w(\cdot)).$ Definition 3.3 (Consistent) The scheme $${F^h}$$ is consistent with the equation $$F$$ if for any smooth function $$\phi$$ and $$x\in{{\it {\Omega}}}$$, $\lim_{h \to 0,y\to x} {F^h}[\phi](y) = F(D^2 \phi(x), \nabla \phi(x),\phi(x),x).$ The scheme is accurate to order $$k$$ if for any smooth function $$\phi$$ in a neighborhood of $$G^h$$ and $$x \in G^h \cap {\it {\Omega}}$$ ${F^h}[\phi](x) = F(D^2 \phi(x), \nabla \phi(x),\phi(x),x) + \mathcal{O}(h^k).$ Remark 3.4 Consistency and accuracy are usually verified by a Taylor series argument. Definition 3.5 (Discrete Comparison Principle) Given the finite difference operator $${F^h}: C(G^h) \to C(G^h)$$, the comparison principle holds for $${F^h}$$ if $$\label{comparison} {F^h}[u](x) \le {F^h}[v](x)~ \text{for all} ~x \implies u(x) \le v(x) ~\text{for all}~ x.$$ (Comp) Remark 3.6 In the discrete comparison principle, the boundary conditions are encoded in $${F^h}$$, the assumption $${F^h}[u] \leq {F^h}[v]$$ means $$u \le v$$ at Dirichlet boundary points. Uniqueness of solutions follows from the discrete comparison principle, because if $$u,v$$ are solutions, then $${F^h}(u) = {F^h}(v) = 0$$, so $$u \le v$$ and $$u\ge v$$, and thus $$u=v$$. Definition 3.7 (Lipschitz constant of the scheme) The finite difference operator $${F^h}: C(G^h) \to C(G^h)$$ is Lipschitz continuous with constant $$K^h$$ if $$K^h$$ is the smallest constant such that $\left\vert{{F^h}\left(x,r,v(\cdot)\right) - {F^h}\left(x,s,w(\cdot)\right)}\right\vert \le K^h \max(\vert{r-s}\vert, \Vert{v - w}\Vert_\infty), ~\text {for all}~ x \in G^h.$ Remark 3.8 In the definition above, the maximum on the right-hand side can be replaced with a maximum over the neighboring grid values without changing the Lipschitz constant. In Oberman (2006), it is shown that the Euler map $$u \to u - \rho F[u]$$ with constant $$\rho \le 1/K^h$$ is monotone and a (nonstrict) contraction. By adding an arbitrarily small multiple of $$u$$ to the scheme, the comparison principle holds, and the Euler map is a strict contraction. The resulting scheme provides a convergent discretization of the solution of the parabolic PDE $$u_t + F[u]=0$$, because the discretization is monotone and stable. Filtered schemes were first proposed in Froese & Oberman (2013) in the context of the Monge–Ampère equation to overcome the limited accuracy of the wide-stencil elliptic schemes. The idea is to blend a stable elliptic convergent scheme with an accurate scheme and retain the advantages of both: stability and convergence of the former and higher accuracy of the latter. The only requirement on the accurate scheme is consistency. Filtered schemes have also been used for Hamilton–Jacobi equations in Oberman & Salvador (2015). We require that the difference between the filtered scheme and the elliptic scheme is uniformly bounded. Under this assumption, the proof of the Barles–Souganidis theorem can be modified to obtain convergence results for filtered schemes (Froese & Oberman, 2013). 3.2 Example finite difference schemes Define the standard finite difference operators for $$u_x, u_y$$ and $$u_{xy}$$. Definition 3.9 (Standard finite differences for $$u_x$$, $$u_y$$ and $$u_{xy}$$) \label{standard_fd} \begin{aligned} & u_x^h := \frac{u(x+h,y)-u(x-h,y)}{2h} = u_x(x,y) + \mathcal{O}(h^2),\\ & u_{xx}^h := \frac{u(x+h,y)-2u(x,y)+u(x-h,y)}{h^2} = u_{xx}(x,y) + \mathcal{O}(h^2),\\ & u_{xy}^h := \frac{u(x+h,y+h)+u(x-h,y-h)-u(x+h,y-h)-u(x-h,y+h)}{4h^2} = u_{xy}(x,y) + \mathcal{O}(h^2), \end{aligned} (3.1) and, similarly for $$u_y^h$$, and $$u_{yy}^h$$ in the $$y$$ coordinate. The operator $$-u_{xx}^h$$ is elliptic, the others are not. Definition 3.10 Define the backward and forward first-order derivatives \begin{align*} D^-_x[u](x,y) &:= \frac{u(x,y)-u(x-h,y)}{h} = u_x(x,y) + \mathcal{O}(h), \\ -D^+_x[u](x,y) &:= \frac{u(x,y)-u(x+h,y)}{h} = -u_x(x,y) + \mathcal{O}(h), \end{align*} and, similarly for $$D^-_y[u]$$ and $$D^-_y[u]$$. The operators $$D^-_x$$ and $$-D^+_x$$ are elliptic. 3.3 Nondecreasing functions and elliptic discretizations To build elliptic difference schemes, we study the properties of nondecreasing maps. This is required, because, in general, elliptic schemes are built composing nondecreasing maps with elliptic terms. Moreover, for some nonlinear elliptic PDEs, the domain of ellipticity is restricted, and thus we need to build nondecreasing extensions of the underlying functions. In this section, we define elliptic schemes for $$\left\vert{\nabla u}\right\vert$$ and $$-\left\vert{\nabla u}\right\vert$$. This is done both as an illustration of the general principle and because these schemes will be needed for our discretization of the two-dimensional affine curvature operator. Definition 3.11 (Nondecreasing functions) For $$x,y \in \mathbb{R}^N$$ we say $$x \leq y$$ if $$x_i \leq y_i$$ for all $$i =1,\dots, N$$. The function $$F:\mathbb{R}^N \to \mathbb{R}$$ is nondecreasing, $$F \in ND(\mathbb{R}^N)$$, if $x \leq y \implies F(x) \leq F(y).$ Write $$\mathbb{R}^+ = \{x \in \mathbb{R} \mid x \geq 0 \}$$ and $$\mathbb{R}^- = \{x \in \mathbb{R} \mid x \leq 0 \}$$. Furthermore, if $$F \in ND(\mathbb{R}^N)$$ and $$F: \mathbb{R}^N \to \mathbb{R}^+$$ (resp. $$F: \mathbb{R}^N \to \mathbb{R}^-)$$ we write $$F \in ND^+(\mathbb{R}^N)$$ (resp. $$F \in ND^-(\mathbb{R}^N)$$). Remark 3.12 When $$f \in C^1(\mathbb{R}^N)$$, if $$f$$ is nondecreasing in each variable, i.e., $$f_{x_i} \geq 0$$ for all $$i =1,\dots, N$$, then $$f \in ND(\mathbb{R}^N)$$. Example 3.13 The function $$x^+ = \max(x,0) \in ND^+(\mathbb{R})$$ and $$x^- = \min(x,0)$$ is in $$ND^-(\mathbb{R})$$. Write $$N(x,y) = \sqrt{x^2+y^2}$$, and define $$N^+(x,y) := N(x^+,y^+)$$ and $$N^-(x,y) := -N(x^-,y^-)$$. Then $$N^+ \in ND^+(\mathbb{R}^2)$$, and $$N^- \in ND^-(\mathbb{R}^2)$$. Furthermore, $$N^+ = N$$ on $$\{x, y \geq 0 \}$$, $$N^- = -N$$ on $$\{x, y \leq 0 \}$$. Example 3.14 (Upwinding) More generally, if we consider the operator $$a(x)u_x$$, then the upwind discretization $a^+ D^-_x[u] + a^-D^+_x[u]$ is first-order accurate and elliptic. The first example of using elliptic discretizations as building blocks is given by the following. Example 3.15 Define $\left\vert{u_x^h}\right\vert^+ = \max\left\{-D^+_xu,D^-_xu,0\right\}, \quad -\left\vert{u_x^h}\right\vert^- = \min\left\{-D^+_xu,D^-_xu,0\right\}$ that approximate $$\vert{u_x}\vert$$ and $$-\vert{u_x}\vert$$to first order. The operators $$\left\vert{u_x^h}\right\vert^+$$ and $$-\left\vert{u_x^h}\right\vert^-$$ are elliptic, the former is non-negative and the latter is nonpositive. Example 3.16 Define $\vert{\nabla u^h}\vert^+ = N \left(|u_x^h|^+,|u_y^h|^+\right), \qquad -\vert{\nabla u^h}\vert^- = -N\left(-|u_x^h|^-,-|u_y^h|^-\right)$ which are elliptic, consistent with $$\left\vert{\nabla u}\right\vert$$, and $$-\left\vert{\nabla u}\right\vert$$, respectively, and first-order accurate. 4. Numerical methods for the model equation In this section, we build convergent discretizations to the one-dimensional model of our equation defined by $$\label{1Dmodel} F^{1D}[u] := A(u_x,u_{xx}) = \left(u_x^2 \:u_{xx}\right)^{1/3} = f, \quad x \in (-1,1),$$ (AC-1D) as well as the parabolic PDE, $u_t = F^{1D}[u]-f, \quad \text{for} (x,t) \in (-1,1) \times (0,\infty),$ along with initial and boundary conditions. To do so, we need to build elliptic discretizations for $$F^{1D}$$, which is achieved by studying the nonlinear function $$A(p,q)$$. A naive approach would suggest to simply substitute the elliptic discretizations for $$\left\vert{u_x}\right\vert$$ and $$u_{xx}$$ into $$A(p,q)$$. However, this is not possible since $$A(p,q) \not\in ND(\mathbb{R}^2)$$: $$dA/dp = 2/3(q/p)^{1/3}$$ and so $$A$$ fails to be nondecreasing when $$pq < 0$$. Our approach is the following. First, in Section 4.1, we write $$A(p,q)$$ as a sum of two nondecreasing functions in terms of $$\left\vert{p}\right\vert$$, $$-\left\vert{p}\right\vert$$, $$q$$. These terms are replaced by $$\vert{u^h_x}\vert$$, $$-\vert{u^h_x}\vert$$, $$u_{xx}^h$$, which are elliptic and consistent, and thus a consistent and elliptic discretization for $$F^{1D}$$ is build. Then, in Section 4.2, we present a Lipschitz regularization of the function $$A(p,q)$$ and proceed similarly. The Lipschitz regularization is needed to build a provably convergent explicit scheme for the parabolic PDE as discussed in Section 3.1. Finally, in Section 4.3, we present the convergence results. To save space, we omit some of the proofs. They can be found in Oberman & Salvador (2016). 4.1 An elliptic discretization of the one-dimensional operator In the next lemma, we decompose $$A(p,q)$$ into the sum of two nondecreasing functions. Lemma 4.1 Define $$A^+(p,q) \,{=}\, A(p^+,q^+)$$ and $$A^-(p,q) \,{=}\, A(p^-,q^-)$$. Then $$A^+ \,{\in}\, ND^+(\mathbb{R}^2)$$, $$A^- \,{\in}\, ND^-(\mathbb{R}^2)$$ and $-A(p,q) = A(p,-q) = A^+(\vert{p}\vert,-q) + A^-(-\vert{p}\vert,-q), \quad \text{for all}~ p,q.$ We can now build an elliptic scheme for $$-F^{1D}[u]$$. Lemma 4.2 Define the finite difference operator \begin{align} -F^{1D,e}[u] = A^+\left(\left\vert{u_x^h}\right\vert^+, -u_{xx}^h \right) + A^-\left(-\left\vert{u_x^h}\right\vert^-,-u_{xx}^h \right), \end{align} (AC)1D,e where the finite difference operators involved are defined above. Then $$-F^{1D,e}$$ is elliptic and consistent with the $$-F^{1D}$$. Proof. The operators $$|u_x^h|^+$$, $$-|u_x^h|^-$$, $$-u_{xx}^h$$ are elliptic. Then, due to Lemma 4.1, the operator $$-F^{1D,e}[u]$$ consists of the composition of the nondecreasing functions $$A^+$$ and $$A^-$$ with the three preceding operators, with the first two taking the place of $$\vert{p}\vert$$ and $$-\vert{p}\vert$$ and the last one taking the place of $$-q$$. Because the operators are elliptic and the functions are nondecreasing, $$-F^{1D,e}[u]$$ is elliptic. Consistency follows from the consistency of each of the schemes used. □ Remark 4.3 (Accuracy) One challenge in forming a discretization of (AC-1D) is that the accuracy breaks down near $$u_{xx} = 0$$. For example, if we use second-order accurate approximations of $$u_{xx}$$, the accuracy of finite differences for the operator is only $$\mathcal{O}(h^{2/3})$$. Consider the expression $A(p+h^k,q+h^2) = ((p + h^k)^2(q+ h^2))^{1/3},$ where $$k = 1$$ or $$2$$. If $$p\not=0, q=0$$ we get $A(p+h^k,q+h^2) = h^{2/3} (p + h^k)^{2/3} = \mathcal{O}(h^{2/3}).$ So regardless of the accuracy of the discretization of the $$u_x$$ term, the overall accuracy of the scheme cannot be better than $$\mathcal{O}(h^{2/3})$$. In fact, a similar argument shows that the order of accuracy is $$2k/3$$ near $$p=0, q\not=0$$. Our elliptic discretization, which uses $$k=1$$, will have an order of accuracy $$2/3$$. 4.2 Lipschitz regularization of the one-dimensional operator The function $$A(p,q)$$ fails to be Lipschitz continuous near the axis $$p=0$$, $$q=0$$. Hence, the elliptic scheme presented in the previous section is not Lipschitz, a property that is required to build a monotone convergent scheme for the time-dependent problem. Thus, using the notation for the sign function $$\text{sgn}(q) = q/\vert{q}\vert$$ for $$q\not=0$$ and $$\text{sgn}(0) = 0$$ otherwise, we regularize $$A(p,q)$$ as follows. Definition 4.4 Define for $$K= K(\delta), L = L(\delta) > 0$$, the regularized function $$\label{Adelta} A^\delta(p,q) = \text{sgn}(q)\min(\vert{A(p,q)}\vert,K\vert{p}\vert,L\vert{q}\vert).$$ (4.1) Naturally, we can then define the regularized PDE operator as $F^{1D,\delta}[u] = A^\delta(u_x,u_{xx}).$ Defining an elliptic and consistent scheme for $$F^{1D,\delta}[u]$$ is accomplished by noticing that the discretization (AC)1D,e generalizes when we replace $$A$$ with the regularized version $$A^\delta$$ in each term: just like for $$A(p,q)$$ in Lemma 4.1, we can decompose $$A^\delta(p,q)$$ into the sum of two nondecreasing functions. This is achieved in the following lemma. Lemma 4.5 Define $$A^{\delta,+}(p,q) = A^\delta(p^+,q^+)$$ and $$A^{\delta,-}(p,q) = A^\delta(p^-,q^-)$$. Then $-A^\delta(p,q) = A^{\delta,+}(\vert{p}\vert,-q)+A^{\delta,-}(-\vert{p}\vert,-q),$ where $$A^{\delta,+} \in ND^+(\mathbb{R}^2)$$ and $$A^{\delta,-} \in ND^-(\mathbb{R}^2)$$. Ellipticity and consistency of the regularized scheme with respect to $$-F^{1D,e,\delta}[u]$$ follows now easily, just like in Lemma 4.2. For this reason, we omit the proof. Lemma 4.6 For $$K = K(\delta)$$, $$L=L(\delta) > 0$$, define the finite difference scheme $$\label{elliptic_scheme_regular1D} -F^{1D,e,\delta}[u] = A^{\delta,+}\left(\left\vert{u_x^h}\right\vert^+, -u_{xx}^h \right) + A^{\delta,-}\left(-\left\vert{u_x^h}\right\vert^-, -u_{xx}^h \right).$$ (AC)1,D,e,δ Then $$-F^{1D,e,\delta}[u]$$ is elliptic and consistent with $$-F^{1D,\delta}[u]$$. Proving consistency of the regularized scheme with respect to $$-F^{1D}[u]$$ requires extra work as the parameters $$K$$, $$L$$ need to be chosen carefully. The next theorem accomplishes precisely that. Theorem 4.7 Asssume $$K = h^{-1/3}$$ and $$L = h^{-4/3}$$. Let $$x \in {\it {\Omega}}$$ be a reference point on the grid and $$\phi$$ be a smooth function that is defined in a neighborhood of the grid. Then the scheme $$F^{1D,e,\delta}$$ defined by (AC)1D,e,δ is consistent with $$F^{1D}$$ and has accuracy $F^{1D,e,\delta}[\phi](x) - F^{1D}[\phi](x) = \mathcal{O}(h^{2/3}).$ Moreover, $$F^{1D,e,\delta}$$ is Lipschitz continuous with constant $$C^h$$ given by $$\label{Lipconst1D} C^h = \frac{K}{h}+\frac{2L}{h^2} = h^{-4/3} + 2h^{-10/3}.$$ (4.2) Remark 4.8 As a result of the proof, we show as well that $$F^{1D,e,\delta}$$ is consistent with $$F^{1D,\delta}$$ with accuracy $F^{1D,e,\delta}[\phi](x) - F^{1D,\delta}[\phi](x) = \mathcal{O}\left(h^{1-\alpha} + h^{2-\beta}\right),$ if $$K = \mathcal{O}(h^{-\alpha})$$ and $$L = \mathcal{O}(h^{-\beta})$$. The optimal choice is $$\alpha = 1/3$$ and $$\beta = 4/3$$, in which case the overall accuracy of $$F^{1D,e,\delta}$$ and $$F^{1D,e}$$ with respect to $$F^{1D}$$ is the same (see Remark 4.3), meaning that no accuracy is lost due to the regularization. 4.3 Convergence theorems for the one-dimensional model Having proved the ellipticity and consistency of the schemes, the uniform convergence follows as discussed in Section 3.1. The first convergence result is for the elliptic problem, where there is no need for the regularized scheme. Theorem 4.9 Let $$u(x)$$ be the unique viscosity solution of $$F^{1D}[u] = f$$ in $${\it {\Omega}}$$, along with suitable boundary conditions. For each $$h>0$$, let $$u^{1D,e,h}$$ be the uniformly bounded solution of $$F^{1D,e}[u] = f$$. Then $$u^{1D,e,h} \to u$$ locally uniformly, as $$h \to 0$$. Proof. Convergence for the elliptic discretization $$-F^{1D,e}$$ follows from the Barles–Souganidis theorem. □ Unlike in the elliptic problem, in the parabolic problem we need to use the regularized scheme, with the time discretization being given by a forward Euler step. This leads us to $$\label{discretizationtime1D} u(\cdot,t+{\rm d}t) = u(\cdot,t) + {\rm d}t \:F^{1D,e,\delta}[u(\cdot,t)].$$ (4.3) Theorem 4.10 Let $$u(x,t)$$ be the unique viscosity solution of $$u_t = {F}^{1D}[u]$$ in $${\it {\Omega}} \times [0,\infty)$$, along with $$u(x,0) = u_0(x)$$ and suitable boundary conditions. Assume as well that $$K = h^{-1/3}$$ and $$L = h^{-4/3}$$. For each $$h >0$$, let $$u^{1D,e,h}$$ be the uniformly bounded solution of the monotone time discretization (4.3) with $${\rm d}t \leq 1/{C^h}$$ given by (4.2). Then $$u^{1D,e,h} \to u$$ locally uniformly, as $${\rm d}t,h \to 0$$. Proof. The elliptic scheme $$F^{1D,e,\delta}$$ leads to a monotone time discretization, in other words, the solution map satisfies a comparison principle if $${\rm d}t \le 1/{C^h}$$, where $$C^h$$ is the Lipschitz constant of the elliptic scheme by Oberman (2006). Then the Barles–Souganidis theorem (Barles & Souganidis, 1991) applies. □ Remark 4.11 It is also true that, for fixed values of $$K,L$$, there is a unique viscosity solution, $$u^\delta$$, of the regularized PDE. Then, fixing $$K,L$$ and using the discretization above with $${\rm d}t = \mathcal{O}(h^2)$$, the forward Euler method converges uniformly to $$u^\delta$$ as $$h\to 0$$. 5. Nonconvergence of standard finite differences In this section, we show that standard finite differences are unstable for both the one-dimensional and the two-dimensional operators that we study. This instability can be understood from the one-dimensional model, which, if we take $$|u_x| = 1$$, reduces to the operator $$u_{xx}^{1/3}$$. It is certainly plausible that the singularity near $$u_{xx} =0$$ could lead to instabilities. This motivates the regularization introduced in the previous section, which replaces the singularity of the cube root with a linear function with large slope. It also motives the convergent finite difference schemes, which have an explicit time step that ensures the convergence of the time-dependent schemes. We begin with an example where the standard finite scheme does not converge for the one-dimensional model. Next, we consider the time-dependent equation and the associated scheme obtained with the (unregularized) elliptic scheme and a forward Euler step in time. A scaling argument suggest that the choice $${\rm d}t = \mathcal{O}(h^{4/3})$$ should provide a stable scheme. In fact, this scaling argument can be improved to a proof of the maximum principle for the homogeneous equation with the same time step. However, the maximum principle is not enough for convergence (we need the comparison principle) and we demonstrate divergence with that time step. Using a smaller time step $${\rm d}t = \mathcal{O}(h^2)$$ appears to fix the problem. (The standard finite difference scheme diverges for the example we present no matter how small the time step.) The example is then generalized to the two-dimensional operator. 5.1 Nonconvergence of standard finite differences in one dimension Consider the discretization given by inserting the standard centered differences, given by (3.1), $F^{1D,a}[u] = \left(\left(u_x^h\right)^2\left(u_{xx}^h\right)\right)^{1/3}.$ We considered an exact solution $$u_0(x) = \sin(2\pi x)$$ of (AC-1D). Then we solved the time-dependent problem, using the forward Euler time discretization, with initial data given by the solution $$u(x,0) = u_0(x)$$. We found that this solution was unstable for the standard finite difference scheme. The results are displayed in Fig. 1, which illustrates that the numerical solution diverges from the exact solution. This effect persists over different grid sizes and over smaller time steps. Fig. 1. View largeDownload slide Solution of the one-dimensional model equation using standard finite differences at times $$t \in \{0,1,2,5\}$$. Set $${\rm d}t = h^2/2$$ on a $$256$$-point grid. Fig. 1. View largeDownload slide Solution of the one-dimensional model equation using standard finite differences at times $$t \in \{0,1,2,5\}$$. Set $${\rm d}t = h^2/2$$ on a $$256$$-point grid. Next, we consider the discretization of the time-dependent equation $$\label{aff1dt} u_t = F^{1D}[u] - f.$$ (5.1) Using a forward Euler method in time and the elliptic method in space leads to $$\label{aff1dth} u(\cdot,t+dt) = {\it {\Phi}}_t^{1D,e}(u) := u(\cdot,t) + {\rm d}t (F^{1D,e}[u(\cdot,t)]-f).$$ (5.2) A scaling argument suggests $${\rm d}t = \mathcal{O}(h^{\frac{4}{3}})$$ might be stable, because the operator $$F^{1D}$$ is homogeneous to this order in $$h$$. In fact, we are able to prove the following. Lemma 5.1 When $$f=0$$ in (5.1), the solution map $${\it {\Phi}}^{1D,e}$$ satisfies the maximum principle $\min {\it {\Phi}}^{1D,e}_t(u) \leq {\it {\Phi}}^{1D,e}_{t+{\rm d}t}(u) \leq \max {\it {\Phi}}^{1D,e}_t(u),$ provided $${\rm d}t \leq (h^4/2)^{1/3}$$. The proof can be found in Oberman & Salvador (2016). However, as the following example shows, the maximum principle by itself is not enough to guarantee convergence and so the above choice for the time step is not guaranteed to produce a convergent scheme. In practice, the above choice for the time step leads to a divergent, but bounded scheme with nonlinear instabilities. Example 5.2 We solved (5.1) using (5.2). We took $$u(x,0) = u_0$$ where $$u_0(x) = C x^{4/3}$$ and $$f =F^{1D}[u_0]$$. In Fig. 2, we plot the numerical solution obtained at different times with $${\rm d}t = (h^4/4)^{1/3}$$ (The conservative choice of the time step is to account for the fact the equation is not homogeneous.) The exact solution is clearly unstable. For larger times, the solution has the same shape, but with small high-frequency oscillations in time. On the other hand, choosing $${\rm d}t = h^2/2$$ leads to convergence (data are not presented to save space). Fig. 2. View largeDownload slide Plot of the solution obtained described in Example 5.2 at time $$t \in \{0,1,5,20\}$$ on a $$128$$-point grid. Fig. 2. View largeDownload slide Plot of the solution obtained described in Example 5.2 at time $$t \in \{0,1,5,20\}$$ on a $$128$$-point grid. 5.2 Two dimensions In this section, we define and study the standard finite difference scheme for $$F[u]$$. We show that the accuracy degenerates near points where the affine curvature is zero. We give an example of a steady solution where the standard finite difference scheme breaks down. Using (1.1), we can write \begin{align*} F[u] & = \left(u_{xx} u_y^2 - 2u_x u_y u_{xy} + u_{yy} u_x^2\right)^{1/3}. \end{align*} Definition 5.3 Let $$u:\mathbb{R}^2\to\mathbb{R}$$. We define the scheme \begin{align} F^a[u] = \left(u_{xx}^h (u_y^h)^2 - 2u_x^h u_y^h u_{xy}^h + u_{yy}^h (u_x^h)^2\right)^{1/3} \end{align} (AC)a that approximates $$F[u]$$. Remark 5.4 The $$u_{xy}^h$$ term is not elliptic, and consequently $$-F^a$$ is not elliptic. Lemma 5.5 Let $$(x,y) \in {\it {\Omega}}$$ be a reference point on the grid and $$\phi$$ be a smooth function that is defined in a neighborhood of the grid. Then the scheme $$F^a[\phi]$$ defined by (AC)a approximates $$F[\phi]$$ with accuracy $F^a[\phi](x,y) - F[\phi](x,y) = \begin{cases} \mathcal{O}(h^2), & F[\phi](x,y) \neq 0,\\ \mathcal{O}(h^{\frac{2}{3}}), &F[\phi](x,y) = 0. \end{cases}$ Proof. All the standard finite differences used are second-order accurate. Therefore, $\phi_{xx}^h \left(\phi_y^h\right)^2 - 2\phi_x^h \phi_y^h \phi_{xy}^h + \phi_{yy}^h \left(\phi_x^h\right)^2 = \phi_{xx} \phi_y^2 - 2\phi_x \phi_y \phi_{xy} + \phi_{yy} \phi_x^2 + \mathcal{O}(h^2),$ in other words, $$(F^a[\phi])^3 = (F[\phi])^3 + \mathcal{O}(h^2)$$. In addition, the Taylor expansion tells us that $(r+t)^{1/3} = r^{1/3} + \frac{t}{3r^{2/3}} + \mathcal{O}(t^2)$ when $$r \neq 0$$. Hence, when $$F[\phi](x,y) \neq 0$$, it follows that $$F^a[\phi] = F[\phi]+\mathcal{O}(h^2)$$. Otherwise, when $$F[\phi](x,y) = 0$$, we can only conclude that $$F^a[\phi] = F[\phi]+\mathcal{O}(h^{\frac{2}{3}})$$. □ Now we present an example that shows that the scheme $$F^a$$ does not converge. We choose a level set function and a right-hand side so that the equation is a steady state (see Example 7.9(d), for more details). Starting from initial data corresponding to the exact solution, and evolving in time with a forward Euler step, the solution changes to order one. Indeed, the solution does not appear to reach a steady state, even after running for a long time. See Fig. 3 for snapshots in time of the solution. Similar behavior was observed on finer grids (although it took a longer time to reach a comparable change in the solution). Fig. 3. View largeDownload slide Example 7.9 $$(d)$$ with the standard finite difference solver: level sets of the solution at times $$t = 0$$ (upper left), $$t = 15$$ (upper center), $$t = 17$$ (upper right), $$t = 20$$ (lower left), $$t = 40$$ (lower center) and $$t = 50$$ (lower right) with $${\rm d}t = h^2/2$$ on a $$32 \times 32$$ grid. Fig. 3. View largeDownload slide Example 7.9 $$(d)$$ with the standard finite difference solver: level sets of the solution at times $$t = 0$$ (upper left), $$t = 15$$ (upper center), $$t = 17$$ (upper right), $$t = 20$$ (lower left), $$t = 40$$ (lower center) and $$t = 50$$ (lower right) with $${\rm d}t = h^2/2$$ on a $$32 \times 32$$ grid. 6. Convergent finite difference methods Based on the convergence theory discussed above, and on the elliptic discretization of the one-dimensional model equation, we now build a discretization of the affine curvature operator. This discretization is based on the median scheme for the mean curvature operator from Oberman (2004). We could also use the morphological operator (Cattè et al., 1995), which results in a very similar discretization of the operator $$F[u]$$. We establish the accuracy of the discretization and show that it is elliptic. The scheme is augmented to a more accurate, but is still a convergent, filtered scheme, which interpolates between the standard scheme and the elliptic scheme. We regularize the operator, which allows us to build a convergent monotone time discretization. 6.1 Median for $${\it {\Delta}}_1 u$$ In Oberman (2004), an elliptic scheme for $${\it {\Delta}}_1 u(x)$$ is presented based on taking the median of the values of $$u$$ sampled in a small approximately circular neighborhood of $$x$$. The motivation follows from observing, using (1.1), that ${\it {\Delta}}_1 u = u_{\mathbf{t}\mathbf{t}}, \qquad \mathbf{t} = \frac{(-u_y,u_x)}{(u_x^2+u_y^2)^{1/2}},$ where $$\mathbf{t}$$ is the (Euclidean) unit tangent. The median captures an approximation to $$u_{\mathbf{t} \mathbf{t}}$$, the second tangential derivative of $$u$$, because the larger values point in the direction of the gradient and the smaller values point in the opposite directions. We outline the scheme here. We will define it at the reference point $$(x,y)$$. Let $$(x_1,y_1),\ldots,(x_{n_S},y_{n_S})$$ be the $$n_S$$ neighbors and set $${\rm d}\theta = \frac{2\pi}{n_S}$$. We refer to $${\rm d}\theta$$ as the directional resolution. Denote the value of the solution at the point $$(x_i,y_i)$$ by $$u_i$$. Set $$\mathbf{v_i} = (x_i,y_i)-(x,y)$$ and $$(x_{-i},y_{-i}) = (x,y)-\mathbf{v_i}$$. We choose the neighbors such that $(x_{i+1},y_{i+1}) = h \: n_\theta\left(\cos(i{\rm d}\theta),\sin(i {\rm d}\theta)\right) + (e_i,f_i),$ where $$|e_i|,|\,f_i| \leq h$$. Thus, $$h$$ is the spatial resolution and $$n_\theta$$ denotes the width of the stencil. In fact, for $$n_\theta \leq 5$$, the neighbors in the first quadrant are given by $(x_{i+1},y_{i+1}) = h\left(\lfloor n_\theta \cos(i{\rm d}\theta) \rceil, \lfloor n_\theta \sin(i{\rm d}\theta) \rceil\right)\!,$ where $$i = 0,\ldots,\frac{n_S}{4}-1$$, with the points on the remaining three quadrants being obtained by $$\frac{\pi}{2}$$, $$\pi$$ and $$\frac{3\pi}{2}$$ rotations. Here $$\lfloor x \rceil$$ denotes the integer closest to $$x$$. See Fig. 4 for the stencils with $$n_\theta = 3$$ and $$n_\theta = 7$$. Fig. 4. View largeDownload slide Neighbor points of the stencil for $$n_\theta = 3$$ (smaller circle) and $$n_\theta = 7$$ (larger circle). Fig. 4. View largeDownload slide Neighbor points of the stencil for $$n_\theta = 3$$ (smaller circle) and $$n_\theta = 7$$ (larger circle). Definition 6.1 Let $$u:\mathbb{R}^2 \to \mathbb{R}$$. Define the scheme $${\it{\Delta}} _{1}^{e}u(x,y)=2\frac{\underset{i=1,\ldots, {{n}_{S}}}{\mathop{\text{median}}}\,{{u}_{i}}-u(x,y)}{{{(h\,{{n}_{\theta}})}^{2}}}$$ (MC)e that approximates $${\it {\Delta}}_1 u$$. In general, consistency of finite difference schemes follows by Taylor expansions. However, additional care is needed when the PDE is singular, as when $$\nabla u(x) = 0$$ in (MC). We then recall here the definition of consistency for the mean curvature operator. Definition 6.2 The numerical scheme $$F^{h,d\theta}$$ is consistent with $${\it {\Delta}}_1$$ if for every smooth function $$\phi$$ and every $$(x,y) \in \mathbb{R}^2$$ $\lim_{h,{\rm d}\theta \to 0} F^{h,{\rm d}\theta}[\phi] = {\it {\Delta}}_1\phi$ at $$(x,t)$$ if $$\nabla \phi(x,y) \neq 0$$ and $\lambda \leq \liminf_{h,{\rm d}\theta \to 0} F^{h,{\rm d}\theta}[\phi] \leq \limsup_{h,{\rm d}\theta \to 0} F^{h,{\rm d}\theta}[\phi] \leq {\it {\Lambda}}$ at $$(x,t)$$ where $$\lambda,{\it {\Lambda}}$$ are the least and greatest eigenvalues of $$D^2\phi(x,y)$$, otherwise. Lemma 6.3 The finite difference scheme $$-{\it {\Delta}}_1^e u$$, given by (MC)e, is elliptic. Proof. We can rewrite the scheme as $-{\it{\Delta}}_{1}^{e}u=2\frac{\underset{i=1,\ldots, {{n}_{S}}}{\mathop{\text{median}}}\,(u(x,y)-{{u}_{i}})}{{{(h\,{{n}_{\theta}})}^{2}}}.$ Since $$\text{median}$$ is a nondecreasing function, we can immediately conclude that $${\it {\Delta}}_1^e u$$ is elliptic. □ Lemma 6.4 Let $$(x,y) \in {\it {\Omega}}$$ be a reference point on the grid and $$\phi$$ be a smooth function that is defined in a neighborhood of the grid. Then the scheme $${\it {\Delta}}_1^e \phi$$ defined by (MC)e is consistent. Further, it approximates $${\it {\Delta}}_1$$ with accuracy ${\it {\Delta}}_1^e \phi(x,y) = {\it {\Delta}}_1 \phi(x,y) + \mathcal{O}\left((h \: n_\theta)^2 + {\rm d}\theta\right)\!,$ when $$\left\vert{\nabla u(x,y)}\right\vert \neq 0$$. Remark 6.5 Since $${\rm d}\theta = \mathcal{O}\left(\frac{1}{n_\theta}\right)$$, the ‘optimal’ choice is $$h = \mathcal{O}\left(n_\theta^{-3/2}\right)$$, which results in accuracy of $$\mathcal{O}(h^{2/3})$$. However, this also requires a dense stencil. In practice, we use a fairly narrow stencil and combine with a filtered scheme for more accuracy. The proof can be found in Oberman (2004). 6.2 Elliptic scheme for $$-F[u]$$ We now construct an elliptic scheme for $$-F[u]$$ following the same procedure as in Section 4.1. This is accomplished by writing $$F[u]$$ in terms of $$|\nabla u|$$ and $${\it {\Delta}}_1 u$$ as \begin{align*} F[u] & = \left\vert{\nabla u}\right\vert k[u]^{1/3} = (\left\vert{\nabla u}\right\vert^2 {\it {\Delta}}_1 u)^{1/3} = A(\left\vert{\nabla u}\right\vert,{\it {\Delta}}_1 u) \end{align*} and using the elliptic schemes $$\vert{\nabla u^h}\vert^+$$ and $$\vert{\nabla u^h}\vert^-$$ presented in Section 3.3 and $${\it {\Delta}}_1^e u$$ described in Section 6.1. Remark 6.6 We chose to discretize $${\it {\Delta}}_1u$$ with the median scheme (MC)e. However, other schemes could have been, for instance, the morphological scheme (Cattè et al., 1995). Definition 6.7 Let $$u:\mathbb{R}^2\to\mathbb{R}$$. We define the scheme \begin{align} -F^e[u] = A^+\left(\left\vert{\nabla u^h}\right\vert^+, -{\it {\Delta}}_1^e u\right) + A^-\left(-\left\vert{\nabla u^h}\right\vert^-, -{\it {\Delta}}_1^e u \right) \end{align} (AC)e that approximates $$-F[u]$$. Remark 6.8 The above approach can be generalized to obtain an elliptic scheme for the elliptic operator $$-\left\vert{\nabla u}\right\vert G(k[u])$$, where $$G$$ is nondecreasing and homogeneous of order $$\alpha \le 1$$, $$G(t r) = t^\alpha G(r)$$. (This PDE represents curves evolving with normal speed $$G(k)$$.) In such cases, write \begin{align*} \left\vert{\nabla u}\right\vert G(k[u]) = \left\vert{\nabla u}\right\vert^{1-\alpha} G(\left\vert{\nabla u}\right\vert k[u]) = \left\vert{\nabla u}\right\vert^{1-\alpha} G({\it {\Delta}}_1 u). \end{align*} The elliptic scheme is then given by $\left(|\nabla u^h|^+\right)^{1-\alpha} G\left(\left(-{\it {\Delta}}_1^e u\right)^+\right) +\left(|\nabla u^h|^-\right)^{1-\alpha} G\left(\left(-{\it {\Delta}}_1^e u\right)^-\right).$ Now we prove the ellipticity and consistency of the scheme $$-F^e[u]$$. Lemma 6.9 The finite difference scheme $$-F^e[u]$$, given by (AC)e, is elliptic. Proof. The proof is similar to Lemma 4.6. □ Lemma 6.10 Let $$(x,y) \in {\it {\Omega}}$$ be a reference point on the grid and $$\phi$$ be a smooth function that is defined in a neighborhood of the grid. Then the scheme $$F^e[\phi]$$ defined by (AC)e is consistent with $$F[\phi]$$ and has accuracy $F^e[\phi](x,y) - F[\phi](x,y) = \begin{cases} \mathcal{O}\left(h + (h \: n_\theta)^2+ h \: {\rm d}\theta +{\rm d}\theta\right)\!, & \text{if}~F[\phi](x,y) \neq 0,\\ \mathcal{O}\left(((h \: n_\theta)^2 + h \: {\rm d}\theta + {\rm d}\theta)^{1/3}\right)\!, & \text{if}~F[\phi](x,y) = 0 \text{and} \left\vert{\nabla \phi}\right\vert(x,y) \neq 0,\\ \mathcal{O}\left(h^{2/3}\right) & \text{if}~\left\vert{\nabla \phi}\right\vert(x,y) = 0. \end{cases}$ Proof. Suppose first that $$\left\vert{\nabla \phi}\right\vert(x,y) \neq 0$$. We have $\left\vert{\nabla \phi(x,y)}\right\vert^2 = \left(|\nabla \phi^h(x,y)|^\pm\right)^2 + \mathcal{O}(h) \quad \text{and} \quad {\it {\Delta}}_1^e \phi(x,y) = {\it {\Delta}}_1 \phi(x,y) + \mathcal{O}\left({\rm d}\theta + (h \: n_\theta)^2\right)\!,$ (see Lemma 6.4). Hence, when $$F[\phi](x,y) \neq 0$$, $F^e[\phi]^3 = F[\phi]^3 + \mathcal{O}\left(h + (h \: n_\theta)^2 + h \: {\rm d}\theta +{\rm d}\theta\right)$ and, by the Taylor expansion like in Lemma 5.5, we have $$F^e[\phi](x,y) = F[\phi](x,y) + \mathcal{O}(h + (h \: n_\theta)^2 + h \: {\rm d}\theta +{\rm d}\theta)$$. Otherwise, when $$F[\phi](x,y) = 0$$, we necessarily have $${\it {\Delta}}_1 \phi = 0$$ and so $F^e[\phi]^3 = F[\phi]^3 + \mathcal{O}\left((h \: n_\theta)^2 + h \: {\rm d}\theta +{\rm d}\theta\right)\!.$ Therefore, $$F^e[\phi](x,y) = F[\phi](x,y) + \mathcal{O}\left(((h \: n_\theta)^2 + h \: {\rm d}\theta + {\rm d}\theta)^{1/3}\right)$$. When $$\left\vert{\nabla \phi}\right\vert(x,y) = 0$$, $\left\vert{\nabla \phi(x,y)}\right\vert^2 = \left(|\nabla \phi^h(x,y)|^\pm\right)^2 + \mathcal{O}(h^2)$ with $${\it {\Delta}}_1^e \phi$$ being bounded. Hence, $$F^e[\phi] = F[\phi] + \mathcal{O}(h^{2/3})$$. □ Remark 6.11 Since $${\rm d}\theta = \mathcal{O}\left(\frac{1}{n_\theta}\right)$$, the ‘optimal’ choice is $$h = \mathcal{O}\left(n_\theta^{-3/2}\right)$$, which results in accuracy of $$\mathcal{O}(h^{2/9})$$ in the worst case. However, this also requires a dense stencil. In practice, we use a fairly narrow stencil and combine with a filtered scheme for more accuracy. 6.3 Filtered scheme for $$F[u]$$ The filtered scheme, $$F^f[u]$$, first introduced in Froese & Oberman (2013), is defined to be a continuous interpolation between the accurate and the elliptic scheme, which equals the accurate scheme when the two schemes are within $$\epsilon$$ of each other. Here, we define filtered schemes by making use of the auxiliary function $$S^\epsilon:\mathbb{R}^2 \to \mathbb{R}$$, defined to be a continuous function which for $$(a,b) \in R^2$$, is equal to $$a$$ near the diagonal and $$b$$ off the diagonal. Definition 6.12 Define for $$\epsilon \,{>}\, 0$$, $$A_\epsilon \,{=}\, \left\{(a,b) \,{\in}\, \mathbb{R}^2 \mid |a-b| \,{<}\, \epsilon\right\}$$. Set $$\rho \,{=}\, 10\epsilon$$. Define $$d \,{=}\, {\text{dist}}\left((a,b),A_\epsilon\right)$$. $S^\epsilon(a,b) = \begin{cases} a & \text{if}~ (a,b) \in A_\epsilon,\\ \frac{\rho-d}{\rho}a+\frac{d}{\rho}b & \text{if}~ d \leq \rho,\\ b & \text{otherwise}. \end{cases}$ Define the filtered scheme $$\label{filterAC} F^f[u] = S^\epsilon\left(F^a[u],F^e[u]\right),$$ (AC)f where $$\epsilon = \epsilon(h,{\rm d}\theta)$$. Although, theoretically, the only requirement on $$\epsilon$$ to ensure the convergence of the filtered schemes is that $$\epsilon \to 0$$ as $$h,{\rm d}\theta\to0$$, in practice $$\epsilon$$ must be chosen carefully. Intuitively, it should be large enough to permit the accurate scheme to be active where the solution is smooth (as shown in the next lemma) and small enough to force the monotone scheme to be active whenever needed for convergence (for instance, when the solution is singular). Lemma 6.13 Suppose that the formal discretization errors of the schemes $$F^a$$ and $$F^e$$ are $$\mathcal{O}(r^{\beta_a})$$ and $$\mathcal{O}(r^{\beta_e})$$, respectively. Choose $$\alpha$$ such that $$\beta_a > \beta_e > \alpha > 0$$. Then if $$\phi$$ is smooth and $$\epsilon = r^\alpha$$, $$F^f[\phi] = F^a[\phi]$$. Proof. If $$\phi$$ is smooth then ${F^a[\phi]-F^e[\phi]} = {\mathcal{O}(r^{\beta_a})+\mathcal{O}(r^{\beta_e})} = \mathcal{O}(r^{\beta_e}) \leq \mathcal{O}(\epsilon),$ so using the definition of $$F^f$$ it follows that $$F^f[\phi] = F^a[\phi]$$, as desired. □ Remark 6.14 The lemma tells us that heuristically we could choose $$\epsilon = \sqrt{\text{Acc}[F^e]}$$, where $$\text{Acc}[F^e]$$ denotes the accuracy of $$F^e$$ (which corresponds to the choice $$\alpha = \beta_e/2$$ in the above lemma). In the numerical results presented here, we defined $$\epsilon$$ based on the accuracy of the scheme away from the singularities of $$F[u]$$. In practice, the filtered scheme allows us to neglect the error coming from the wide stencil, whereas in theory we still need to send $${\rm d}\theta \to 0$$ for convergence of the filtered scheme. In our numerical examples, we obtain the accuracy of the accurate scheme in most cases. 6.4 Regularization and Forward Euler method As in the for the one-dimensional model equation (see Section 4.2), to build a provably convergent scheme for the time dependent equation (AC) we need to regularize the operator. Write $$F[u] = A(\left\vert{\nabla u}\right\vert,{\it {\Delta}}_1 u)$$, where $$A(p,q) = (p^2q)^{1/3}$$, and regularize the cube root function as before. This leads to $F^\delta[u] = A^\delta(\left\vert{\nabla u}\right\vert,{\it {\Delta}}_1 u),$ where $$A^\delta(p,q) = \text{sgn}(q)\min(\vert{A(p,q)}\vert, K \vert{p}\vert, L \vert{q}\vert)$$. Remark 6.15 The regularized operator, $$F^\delta[u]$$, is still a level set operator. To see this, it is enough to take $$K = L = 1/\delta$$: $F^\delta[u] = \left\vert{\nabla u}\right\vert \text{sgn}(k[u])\min\left(\left\vert{k[u]}\right\vert^{1/3},\frac{\left\vert{k[u]}\right\vert}{\delta},\frac{1}{\delta}\right).$ The operator reduces to either a multiple of the mean curvature operator or the eikonal equation and otherwise we obtain $$F[u]$$. Similar results regarding ellipticity and consistency hold as for the one-dimensional model, which we present without proof. Lemma 6.16 For $$K=K(\delta)$$, $$L=L(\delta)$$ such that $$K\sqrt{L}\geq1$$, define the finite difference scheme \begin{align}\label{elliptic_scheme_regular} -F^{e,\delta}[u] = A^{\delta,+}\left(|\nabla u^h|^+, -{\it {\Delta}}_1^e u \right) + A^{\delta,-}\left(-|\nabla u^h|^-, -{\it {\Delta}}_1^e u \right). \end{align} (AC)e, δ Then, $$-F^{e,\delta}[u]$$ is elliptic and consistent with $$-F^\delta[u]$$. As with the one-dimensional case, $$K$$ and $$L$$ need to be chosen appropriately in order for $$-F^{e,\delta}[u]$$ to be consistent with $$-F[u]$$. Assuming $$K = \mathcal{O}(h^{-\alpha})$$, $$L = \mathcal{O}(h^{-\beta})$$ and $${\rm d}\theta = \mathcal{O}(1/n_\theta) = \mathcal{O}(h^\gamma)$$, full consistency is obtained when $$\alpha \in (0,1)$$, $$\beta \in \left(0,\frac{2}{3}\right)$$ and $$\gamma \in \left(\beta,\frac{2-\beta}{2}\right)$$. The optimal choices are $$\alpha \in [1/9,7/9]$$, $$\beta = 4/9$$, $$\gamma = 2/3$$ and, in such case, the accuracy is $$\mathcal{O}(h^{2/9})$$. As in the one-dimensional model, no accuracy is lost with the regularization (see Remark 6.11). The Lipschitz constant of $$F^{e,\delta}[u]$$ is given by $$\label{Lipconst} C^h = \frac{K}{h} + \frac{2 L}{(h \: n_\theta)^2}.$$ (6.1) In practice, we will choose $$K = c_K h^{-1/9}$$ and $$L = c_L h^{-4/9}$$, which leads to $C^h = c_Kh^{-10/9} +2\frac{c_L}{{n_\theta}^2}h^{-22/9}.$ These results are generalizations of the results proved for the one-dimensional model, whose proofs can be found in Oberman & Salvador (2016). We can finally define and prove the convergence of the discretizations for the time-dependent equation (AC). The time derivative is discretized with an explicit forward Euler step, whereas $$F[u]$$ is discretized using either the elliptic scheme $$F^{e,\delta}$$ or the regularized filtered scheme $$F^{f,\delta}[u]$$ (this is easily defined by replacing the elliptic scheme in $$F^f$$ by its regularized version). This leads to monotone (resp. filtered) time discretization with solution map given by $$\label{discretizationtime} u(\cdot,t+{\rm d}t) = u(\cdot,t) + {\rm d}t F^{e,\delta}[u(\cdot,t)], \quad \left(\text{resp.} = u(\cdot,t) + {\rm d}t F^{f,\delta}[u(\cdot,t)\right)\!.$$ (6.2) 6.5 Convergence theorems Having proved the ellipticity and consistency of the schemes, the uniform convergence follows as discussed in Section 3.1, provided there exists unique viscosity solutions to the PDEs (AC) and $$F[u] = f$$ along with the homogeneous Neumann boundary conditions, which is assured by the theory in Giga (2006), as explained above. We also need existence of solutions to the elliptic and filtered schemes. For the affine curvature evolution PDE (AC), this is immediate, because the schemes are explicit. In the case of the stationary equation, additional work is required to prove existence of solutions. These do not need to be unique to apply the convergence theory, but both existence and uniqueness results for the schemes follow from the discrete comparison principle and from finite dimensional fixed point theorems, as in Oberman (2006) and Froese & Oberman (2013). For readability, the existence of solutions to the schemes will be part of the assumptions in the theorem. In practice, we can get discrete comparison from a general theorem by adding a small multiple of $$u$$ to the scheme, Oberman (2006), or we can do extra work to prove it directly. There are two convergence results. The first is for the elliptic problem, where there is no need for regularization. The second is for the parabolic problem, where we need to regularize and use an explicit Euler time step. Theorem 6.17 Let $${\it {\Omega}}$$ be a bounded convex domain with a $$C^2$$ boundary and $$u$$ denote the unique viscosity solution of $$F[u] = f$$ in $${\it {\Omega}}$$, along with homogeneous Neumann boundary conditions on $$\partial {\it {\Omega}}$$. For each $$\epsilon = \epsilon(h,{\rm d}\theta) >0$$, let $$u^{e,\epsilon}$$, (resp. $$u^{f,\epsilon}$$) be the uniformly bounded solution of $$F^e[u] = f$$, (resp. $$F^f[u] = f$$). Then $u^{e,\epsilon} \to u ~\text{and}~ u^{f,\epsilon} \to u, \quad \text{locally uniformly, as}~ \epsilon \to 0.$ Proof. The assumptions on $${\it {\Omega}}$$ guarantee that the PDE satisfies a comparison principle. Convergence for the elliptic discretization $$-F^e$$ then follows from the Barles–Souganidis theorem. The modification of the proof for filtered schemes can be found in Froese & Oberman (2013). □ Theorem 6.18 Let $${\it {\Omega}}$$ be a bounded convex domain with a $$C^2$$ boundary. Assume that $$u(x,t)$$ is the unique viscosity solution of $$u_t = F[u]$$ in $${\it {\Omega}} \times [0,\infty)$$, along with $$u(x,0) = u_0(x)$$ and homogeneous Neumman boundary conditions. Assume as well that $$K$$ and $$L$$ are picked so that $$F^{e,\delta}$$ is consistent with $$F$$. For each $$\epsilon = \epsilon(h,{\rm d}\theta) >0$$, let $$u^{e,\epsilon}$$ (resp. $$u^{f,\epsilon}$$) be the uniformly bounded solution of the monotone (resp. filtered) time discretization (6.2) with $${\rm d}t \leq 1/{K^h}$$, where $$K^h$$ is the Lipschitz constant of $$F^{e,\delta}[u]$$ (6.1). Then $u^{e,\epsilon} \to u \quad \text{and} \quad u^{f,\epsilon} \to u$ locally uniformly, as $${\rm d}t,\epsilon \to 0$$. Proof. The assumptions on $${\it {\Omega}}$$ guarantee that the PDE satisfies a comparison principle. The elliptic scheme leads to a monotone time discretization, in other words, the solution map satisfies a comparison principle if $${\rm d}t \leq 1/{K^h}$$, where $$K^h$$ is the Lipschitz constant of the elliptic scheme by Oberman (2006). Then the Barles–Souganidis theorem (Barles & Souganidis, 1991) applies. For time discretization of the filtered scheme, the convergence of filtered schemes in Froese & Oberman (2013) applies. □ Remark 6.19 Just like in the one-dimensional model, for fixed values of $$K,L$$, there is a unique viscosity solution, $$u^\delta$$, of the regularized PDE. Then, fixing $$K,L$$ and using the discretization above with $${\rm d}t = \mathcal{O}((h\: n_\theta)^2)$$, the forward Euler method converges uniformly to $$u^\delta$$ as $$h\to 0$$. 7. Numerical results In this section, we start with a simple example to compare the affine invariant curvature motion and the regularized model. We present different examples on the evolution of a single curve in Section 7.1 under the affine curvature motion and compare it with the mean curvature motion. We test the accuracy of the schemes by computing solutions to the time-independent PDE in Section 7.2 and to the time-dependent equation in Section 7.3. We mostly considered stationary problems, because it was easier to generate exact solutions by applying the operator $$F$$ to a given function $$u$$, and including $$f = F[u]$$ as a source function. For the time-dependent problem, we took advantage of the exact solution from Lemma 2.12 to compare the accuracy of the solutions after a short time, $$T= 0.1$$. We test as well in Section 7.4 whether the schemes satisfy numerically the morphology and affine invariance properties that (AC) satisfies (see Theorem 2.10). Definition 7.1 (Parameters for the discretization) We used the regularized schemes with $$K = 20 h^{-1/9}$$ and $$L = 20 h^{-4/9}$$. For the elliptic discretization, we used two different stencils: the narrow elliptic scheme used $$n_\theta = 3$$ and the wider elliptic scheme used $$n_\theta = 7$$. The filtered discretization used the wider elliptic scheme with $\epsilon(h,d\theta) = \sqrt{h} + d\theta/10.$ For the forward Euler method with the elliptic and filtered schemes, we used a constant time step of $${\rm d}t = 1/C_h$$ where $$C_h$$ is given by (6.1). For the forward Euler method with the standard finite difference scheme, we used the same time step as the filtered scheme, except when computing the steady-state solutions in Example 7.9, where we used $${\rm d}t = h^2/2$$ (this choice of time step proved to be enough in practice). The stopping criteria was simply that the $$l^\infty$$ norm of the residual $$\left\vert{F^{f,\delta}[u^n]-f}\right\vert$$ (and similarly for $$F^{e,\delta}, F^a$$) was below $$tol = 10^{-5}$$. Example 7.2 [Comparing the regularized to the unregularized operators] We compare the evolution of an ellipse by the affine invariant curvature motion and by the regularized model. The ellipse should remain an ellipse of fixed eccentricity. The results were obtained by numerically solving (AC) and $$u_t = F^\delta[u]$$ respectively, with the initial condition $u_0(x,y) = \min\left\{\left(\frac{x}{2}\right)^2+y^2-1,1\right\}$ and homogeneous Neumann boundary conditions on $${\it {\Omega}} = [-4,4]^2$$. The narrow elliptic scheme was used for the spatial discretization. The difference between the level sets of the solution of the two equations is visually indistinguishable. Measured in the $$l^\infty$$ norm ranges from $$10^{-6}$$ for early times steps and $$10^{-7}$$ for the later time steps. The level sets of the numerical solution obtained by solving the regularized PDE are depicted in Fig. 5. Fig. 5. View largeDownload slide Evolution of an ellipse by (left) affine curvature, (right) mean curvature (top). Evolution by affine curvature of (left) a diamond and (right) a flatter diamond (bottom). Fig. 5. View largeDownload slide Evolution of an ellipse by (left) affine curvature, (right) mean curvature (top). Evolution by affine curvature of (left) a diamond and (right) a flatter diamond (bottom). Remark 7.3 Although the regularized scheme with the more stringent time step is needed for the convergence proof, in practise we found, as reported in the previous example, using the (unregularized) elliptic discretization with the time step $${\rm d}t = h^2/2$$ gave nearly identical results to within two or more significant digits. For the remaining examples, we present results using the regularized schemes (with $$K = 20h^{-1/9}$$ and $$L = 20h^{-4/9}$$). 7.1 Numerical examples showing curve evolution In this section, we present some numerical examples illustrating the geometric properties of the PDEs. Example 7.4 (Ellipse) We compare the evolution of an ellipse by the affine invariant curvature motion and by mean curvature. For the former, the ellipse should remain an ellipse of fixed eccentricity. For the latter, the ellipse asymptotically approaches a circle instead. The results were obtained by numerically solving (AC) and (MC), respectively, with the initial condition $u_0(x,y) = \min\left\{\left(\frac{x}{2}\right)^2+y^2-1,1\right\}$ and homogeneous Neumann boundary conditions. We took $$[-4,4]^2$$ as the computational domain on a $$128\times 128$$ grid. As for the scheme used, we chose the narrow elliptic schemes for both. In Fig. 5, we plot the zero level sets obtained at $$t \in \{0,0.1,0.3,0.5,0.7,0.9\}$$. Example 7.5 (Diamond) We also compute the solution of (AC) with initial condition $(a) \ u_0(x,y) = \min\left\{|x|+|y|-1,1\right\}, \qquad (b) \ u_0(x,y) = \min\left\{|x|+2|y|-1,1\right\}$ and homogeneous Neumann boundary conditions. The exact solutions are not known. However, we do know that smooth convex curves evolving under affine curvature converge to ellipses until collapsing to a point and that is the behavior we observed here (see Fig. 5). We took $$[-2,2]^2$$ as the computational domain on a $$128 \times 128$$ grid. As for the scheme, we chose again the narrow elliptic scheme. In Fig. 5, we plot the zero level sets of the numerical solution from time $$t = 0$$ to $$t = 0.5$$ in increments of $$0.1$$ for example $$(a)$$ and from $$t=0$$ to $$t=0.3$$ in increments of $$0.05$$ for example $$(b)$$. Example 7.6 (Fan-shape curve) We also compute the solution of (AC) and (MC) with initial condition $u_0(x,y) = \min\left\{c_+^{(1)}(x,y),c_-^{(1)}(x,y),c_+^{(2)}(x,y),c_-^{(2)}(x,y),1\right\}\!,$ where $c_\pm^{(1)}(x,y) = \left(x\pm\frac{1}{2}\right)^2+5\left(y\pm\frac{1}{4}\right)^2-\frac{1}{2} \quad \text{and} \quad c_\pm^{(2)}(x,y) = 5\left(x\pm\frac{1}{4}\right)^2+\left(y\mp\frac{1}{4}\right)^2-\frac{1}{2}$ and homogeneous Neumann boundary conditions. The exact solution is not known. As in the previous example, we took $$[-2,2]^2$$ as the computational domain on a $$128 \times 128$$ grid together with the narrow elliptic scheme. Under the both curvature motions, the curve should initially become convex. At later times, under the affine curvature motion, the curve should evolve to an ellipse as opposed to a circle in the mean curvature case. In Fig. 6, we plot the zero level set of the numerical solutions at time $$t \in \{0,0.05,0.1,0.2\}$$ and observe the exact behavior described. Fig. 6. View largeDownload slide Evolution of a fan-shape like curve under affine curvature motion (top) and mean curvature (bottom) at time $$t~\in~\{0,0.05,0.1,0.2\}$$ (see Example 7.6 for more details). Fig. 6. View largeDownload slide Evolution of a fan-shape like curve under affine curvature motion (top) and mean curvature (bottom) at time $$t~\in~\{0,0.05,0.1,0.2\}$$ (see Example 7.6 for more details). 7.2 Accuracy of stationary solutions We test the accuracy for the following Dirichlet problem $\begin{cases} F[u] = f, &~ \text{in}~ {\it {\Omega}},\\ u = g, & ~\text{on}~ \partial {\it {\Omega}}. \end{cases}$ Solutions are obtained by computing the steady state solution of $$u_t = F[u] - f$$, with $$u(\cdot,t) = g$$ on $$\partial{\it {\Omega}}$$ and $$u(x,0) = u_0(x)$$. We set $${\rm d}t = 1/C_h$$, where $$C_h$$ is given by (6.1) for the elliptic and filtered schemes and $${\rm d}t = h^2/2$$ for the standard finite difference scheme scheme. These examples also demonstrate stability of the time dependent problem for elliptic and filtered scheme, as well as convergence to the steady solution, because we used the time-dependent problem to obtain the solution. Remark 7.7 The unregularized schemes were also used. For these we set $${\rm d}t = h^2/2$$ for all examples, except for the filtered scheme in example (d) where we set $${\rm d}t = h^2/8$$. The results obtained were virtually the same. We set $$u_0(x)$$ to be the exact solution in a layer of seven grid points adjacent to the boundary. As a result, each discretization is initialized at the same set of grid points and, therefore, we can make a fair comparison of their accuracy. On the coarsest grid, we set $$u_0(x) = 0$$ on the interior grid points. To speed up calculations, on finer grids we set $$u_0(x)$$ to be interpolated solutions from the coarser grids at interior grid points. Remark 7.8 When considering Dirichlet boundary conditions, we lose the comparison principle. However, we can always cap off the solution and impose homogeneous boundary conditions. Example 7.9 We consider the following exact solutions \begin{align*} (a) & \quad u(x,y) = x^2+y^2, \quad f(x) = 2 \left(x^2 + y^2\right)^{\frac{1}{3}},\\ (b) & \quad u(x,y) = {\rm e}^{x^2+y^2}, \quad f(x,y) = 2 \left({\rm e}^{3 (x^2 + y^2)}(x^2 + y^2)\right)^{\frac{1}{3}}\!,\\ (c) & \quad u(x,y) = \left(x^2 + y^2\right)^{\frac{1}{3}}, \quad f(x) = \frac{4}{3},\\ (d) & \quad u(x,y) = \frac{\sin(2\pi x)\sin(2\pi y)}{4}, \quad f(x) = \frac{\pi^{\frac{4}{3}}}{2} \left(-\left(2+\cos(4\pi x)+\cos(4\pi y)\right) \sin(2\pi x)\sin(2\pi y)\right)^{\frac{1}{3}}\!, \end{align*} with $${\it {\Omega}} = [-1,1]^2$$. The solutions in (a),(b) and (d) are smooth, but the functions $$f$$ are only Holder continuous, $$C^{0,2/3}$$, with singularities at the origin for (a) and (b), and at several points in example (d). The solution in (c) is only $$C^{0,2/3}$$ with a singularity at the origin, but in this case the function $$f$$ is constant. The results are presented in Table 1. In Example 7.9$$(a)$$ the solution is a quadratic polynomial. The accurate scheme $$F^a$$ gives essentially machine precision, which is not surprising, because this scheme is second-order accurate. The filtered scheme $$F^f$$ gives nearly the same accuracy (with a small discrepancy which could be eliminated by tuning the parameter). On the other hand, both the narrow and the wider elliptic schemes are much less accurate. In contrast, in Example 7.9(b), the solution is smooth, but not quadratic. In this case, we see that the accurate scheme appears to be converging to $$\mathcal{O}(h)$$. The elliptic schemes are less accurate and the filtered scheme is in between. In Example 7.9$$(c)$$, the solution is only Holder continuous. In this case, the elliptic schemes have almost constant error near $$0.01$$ over the range of parameters used. Despite the singular solution, the accurate scheme gives accuracy $$\mathcal{O}(h)$$ and the filtered scheme does just as well. Example 7.9$$(d)$$ shows that the standard finite difference scheme does not converge as discussed in Section 5 (see Fig. 3). The narrow elliptic scheme has almost a constant error $$0.03,$$ and the wider elliptic scheme has an error $$0.1$$ for the smallest grid, decreasing by a factor of two as the grid is refined. The filtered scheme has the best accuracy, achieving an error of $$0.001$$ at the finest grid. In general, we expect better accuracy from the filtered schemes by making a more careful choice of the parameter $$\epsilon$$ (see Remark 6.14). Table 1 Accuracy in the $$l^\infty$$ norm and order of convergence of the schemes for Example 7.9 with regularized schemes Errors and order, Example 7.9 (a) N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.565e$$-$$02 3.978e$$-$$02 2.922e$$-$$07 5.104e$$-$$07 64 2.555e$$-$$02 2.696e$$-$$02 3.019e$$-$$07 2.827e$$-$$07 128 1.623e$$-$$02 1.678e$$-$$02 1.982e$$-$$07 1.984e$$-$$07 256 1.050e$$-$$02 1.055e$$-$$02 9.293e$$-$$08 5.678e$$-$$05 Errors and order, Example 7.9 (b) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 5.462e$$-$$02 7.666e$$-$$02 1.922e$$-$$03 1.985e$$-$$03 64 5.207e$$-$$02 5.872e$$-$$02 9.875e$$-$$04 8.904e$$-$$03 128 4.105e$$-$$02 3.725e$$-$$02 3.385e$$-$$04 7.262e$$-$$03 256 3.173e$$-$$02 2.240e$$-$$02 9.798e$$-$$05 8.065e$$-$$03 Errors and order, Example 7.9 (c) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 1.387e$$-$$02 2.386e$$-$$02 1.129e$$-$$02 1.129e$$-$$02 64 1.310e$$-$$02 7.683e$$-$$03 4.625e$$-$$03 4.625e$$-$$03 128 9.302e$$-$$03 8.202e$$-$$03 1.859e$$-$$03 1.872e$$-$$03 256 6.445e$$-$$03 7.156e$$-$$03 7.426e$$-$$04 7.570e$$-$$04 Errors and order, Example 7.9 (d) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.964e$$-$$02 9.813e$$-$$02 — 1.926e$$-$$02 64 3.788e$$-$$02 4.679e$$-$$02 — 8.309e$$-$$03 128 3.688e$$-$$02 2.367e$$-$$02 — 2.482e$$-$$03 256 3.003e$$-$$02 1.798e$$-$$02 — 9.697e$$-$$04 Errors and order, Example 7.9 (a) N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.565e$$-$$02 3.978e$$-$$02 2.922e$$-$$07 5.104e$$-$$07 64 2.555e$$-$$02 2.696e$$-$$02 3.019e$$-$$07 2.827e$$-$$07 128 1.623e$$-$$02 1.678e$$-$$02 1.982e$$-$$07 1.984e$$-$$07 256 1.050e$$-$$02 1.055e$$-$$02 9.293e$$-$$08 5.678e$$-$$05 Errors and order, Example 7.9 (b) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 5.462e$$-$$02 7.666e$$-$$02 1.922e$$-$$03 1.985e$$-$$03 64 5.207e$$-$$02 5.872e$$-$$02 9.875e$$-$$04 8.904e$$-$$03 128 4.105e$$-$$02 3.725e$$-$$02 3.385e$$-$$04 7.262e$$-$$03 256 3.173e$$-$$02 2.240e$$-$$02 9.798e$$-$$05 8.065e$$-$$03 Errors and order, Example 7.9 (c) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 1.387e$$-$$02 2.386e$$-$$02 1.129e$$-$$02 1.129e$$-$$02 64 1.310e$$-$$02 7.683e$$-$$03 4.625e$$-$$03 4.625e$$-$$03 128 9.302e$$-$$03 8.202e$$-$$03 1.859e$$-$$03 1.872e$$-$$03 256 6.445e$$-$$03 7.156e$$-$$03 7.426e$$-$$04 7.570e$$-$$04 Errors and order, Example 7.9 (d) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.964e$$-$$02 9.813e$$-$$02 — 1.926e$$-$$02 64 3.788e$$-$$02 4.679e$$-$$02 — 8.309e$$-$$03 128 3.688e$$-$$02 2.367e$$-$$02 — 2.482e$$-$$03 256 3.003e$$-$$02 1.798e$$-$$02 — 9.697e$$-$$04 Table 1 Accuracy in the $$l^\infty$$ norm and order of convergence of the schemes for Example 7.9 with regularized schemes Errors and order, Example 7.9 (a) N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.565e$$-$$02 3.978e$$-$$02 2.922e$$-$$07 5.104e$$-$$07 64 2.555e$$-$$02 2.696e$$-$$02 3.019e$$-$$07 2.827e$$-$$07 128 1.623e$$-$$02 1.678e$$-$$02 1.982e$$-$$07 1.984e$$-$$07 256 1.050e$$-$$02 1.055e$$-$$02 9.293e$$-$$08 5.678e$$-$$05 Errors and order, Example 7.9 (b) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 5.462e$$-$$02 7.666e$$-$$02 1.922e$$-$$03 1.985e$$-$$03 64 5.207e$$-$$02 5.872e$$-$$02 9.875e$$-$$04 8.904e$$-$$03 128 4.105e$$-$$02 3.725e$$-$$02 3.385e$$-$$04 7.262e$$-$$03 256 3.173e$$-$$02 2.240e$$-$$02 9.798e$$-$$05 8.065e$$-$$03 Errors and order, Example 7.9 (c) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 1.387e$$-$$02 2.386e$$-$$02 1.129e$$-$$02 1.129e$$-$$02 64 1.310e$$-$$02 7.683e$$-$$03 4.625e$$-$$03 4.625e$$-$$03 128 9.302e$$-$$03 8.202e$$-$$03 1.859e$$-$$03 1.872e$$-$$03 256 6.445e$$-$$03 7.156e$$-$$03 7.426e$$-$$04 7.570e$$-$$04 Errors and order, Example 7.9 (d) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.964e$$-$$02 9.813e$$-$$02 — 1.926e$$-$$02 64 3.788e$$-$$02 4.679e$$-$$02 — 8.309e$$-$$03 128 3.688e$$-$$02 2.367e$$-$$02 — 2.482e$$-$$03 256 3.003e$$-$$02 1.798e$$-$$02 — 9.697e$$-$$04 Errors and order, Example 7.9 (a) N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.565e$$-$$02 3.978e$$-$$02 2.922e$$-$$07 5.104e$$-$$07 64 2.555e$$-$$02 2.696e$$-$$02 3.019e$$-$$07 2.827e$$-$$07 128 1.623e$$-$$02 1.678e$$-$$02 1.982e$$-$$07 1.984e$$-$$07 256 1.050e$$-$$02 1.055e$$-$$02 9.293e$$-$$08 5.678e$$-$$05 Errors and order, Example 7.9 (b) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 5.462e$$-$$02 7.666e$$-$$02 1.922e$$-$$03 1.985e$$-$$03 64 5.207e$$-$$02 5.872e$$-$$02 9.875e$$-$$04 8.904e$$-$$03 128 4.105e$$-$$02 3.725e$$-$$02 3.385e$$-$$04 7.262e$$-$$03 256 3.173e$$-$$02 2.240e$$-$$02 9.798e$$-$$05 8.065e$$-$$03 Errors and order, Example 7.9 (c) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 1.387e$$-$$02 2.386e$$-$$02 1.129e$$-$$02 1.129e$$-$$02 64 1.310e$$-$$02 7.683e$$-$$03 4.625e$$-$$03 4.625e$$-$$03 128 9.302e$$-$$03 8.202e$$-$$03 1.859e$$-$$03 1.872e$$-$$03 256 6.445e$$-$$03 7.156e$$-$$03 7.426e$$-$$04 7.570e$$-$$04 Errors and order, Example 7.9 (d) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.964e$$-$$02 9.813e$$-$$02 — 1.926e$$-$$02 64 3.788e$$-$$02 4.679e$$-$$02 — 8.309e$$-$$03 128 3.688e$$-$$02 2.367e$$-$$02 — 2.482e$$-$$03 256 3.003e$$-$$02 1.798e$$-$$02 — 9.697e$$-$$04 When comparing the different schemes, we have to account for the width of the stencil, because for the elliptic schemes the wider schemes also have a larger spatial discretization error. In general, the accuracy improved with the use of the wider stencil. Moreover, the filtered scheme performed as expected by providing better accuracy than the elliptic scheme and almost as good accuracy as the accurate scheme. The final example shows that the standard finite difference scheme may not converge. This may be due to the fact that, despite the solution being smooth, there were multiple points where $$f$$ was singular. Geometrically, this solution has several points where curves shrunk to zero. We also considered for Example 7.9 the elliptic schemes with $$n_\theta = 1$$ and $$2$$. These provided poor accuracy with the directional resolution error easily dominating the spatial resolution error. The errors do not decrease to zero as we decrease the grid size. This property is consistent with the theoretical results, because convergence is only proved as both $$h,{\rm d}\theta \to 0$$, which is indeed observed in the numerical results. On the other hand, using the narrow and wider schemes, the accuracy of the elliptic scheme was good enough for the filtered scheme to give accuracy comparable to the accurate scheme in many examples. In principle, we still need to send $${\rm d}\theta \to 0$$, but in practice, the rate at which it needs to go to zero is much slower than $$h$$ when the filtered scheme is used. Finally, we point that the choice of $$\epsilon$$ is not an easy one and it is hard to pick an $$\epsilon$$ that yields optimal results in each example presented. By choosing $$\epsilon$$ larger, it is possible to achieve with the filtered schemes the accuracy of the standard schemes in Examples 7.9 (a)–(c). However, such choice is too permissive for Example 7.9 (d). As pointed out in Section 6.3, $$\epsilon$$ needs to be chosen small enough for the monotone scheme to be active to ensure convergence. (This is comparable to the CFL condition in time-dependent equations: methods are convergent as $${\rm d}t,h \to 0$$, with $${\rm d}t$$ satisfying the CFL condition.) 7.3 Accuracy for the time dependent problem Recall here that for general boundary conditions, the time-dependent PDE (AC) requires Neumann boundary conditions for uniqueness of viscosity solutions to hold. We have already established (in the previous section) the stability of the numerical method. In the following examples, we test the accuracy of solutions, comparing two different wide stencil discretizations, along with regularized filtered discretization and the (generally unstable) standard finite difference method. Consequently, we test for the accuracy using Dirichlet boundary conditions coming from the exact solution of Lemma 2.12. Example 7.10 (Neumann BC) We consider the exact solution $u(x,y,t) = \min\left\{t + \frac{3}{4}\left(\frac{b}{a}x^2+\frac{a}{b}y^2\right)^{2/3}-1,0\right\},$ with $$a = 2$$ and $$b = 1$$ (see Lemma 2.12). By taking the minimum with $$0$$, we are imposing homogeneous Neumann boundary conditions, thus avoiding having to deal with boundary of the computation domain. We set $${\it {\Omega}} = [-3,3]^2$$. We display the numerical error in the $$l^\infty$$ norm at time $$T = 0.1$$ in Table 2. Table 2 Error in the $$l^\infty$$ norm of the whole computational domain at time $$t=0.1$$ for the time-dependent Examples 7.10 and 7.11 Errors and order, Example 7.10 N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 4.845e$$-$$02 6.691e$$-$$02 4.894e$$-$$02 4.888e$$-$$02 64 4.432e$$-$$02 4.607e$$-$$02 2.977e$$-$$02 2.975e$$-$$02 128 3.544e$$-$$02 2.823e$$-$$02 2.457e$$-$$02 2.438e$$-$$02 256 2.971e$$-$$02 2.080e$$-$$02 1.747e$$-$$02 1.724e$$-$$02 512 2.764e$$-$$02 1.652e$$-$$02 1.205e$$-$$02 1.182e$$-$$02 Errors and order, Example 7.11 N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 2.182e$$-$$02 1.449e$$-$$02 1.985e$$-$$02 1.985e$$-$$02 64 1.435e$$-$$02 1.160e$$-$$02 1.279e$$-$$02 1.279e$$-$$02 128 9.580e$$-$$03 7.517e$$-$$03 5.566e$$-$$03 5.567e$$-$$03 256 6.404e$$-$$03 4.854e$$-$$03 2.442e$$-$$03 2.409e$$-$$03 512 6.090e$$-$$03 4.288e$$-$$03 1.036e$$-$$03 1.002e$$-$$03 Errors and order, Example 7.10 N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 4.845e$$-$$02 6.691e$$-$$02 4.894e$$-$$02 4.888e$$-$$02 64 4.432e$$-$$02 4.607e$$-$$02 2.977e$$-$$02 2.975e$$-$$02 128 3.544e$$-$$02 2.823e$$-$$02 2.457e$$-$$02 2.438e$$-$$02 256 2.971e$$-$$02 2.080e$$-$$02 1.747e$$-$$02 1.724e$$-$$02 512 2.764e$$-$$02 1.652e$$-$$02 1.205e$$-$$02 1.182e$$-$$02 Errors and order, Example 7.11 N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 2.182e$$-$$02 1.449e$$-$$02 1.985e$$-$$02 1.985e$$-$$02 64 1.435e$$-$$02 1.160e$$-$$02 1.279e$$-$$02 1.279e$$-$$02 128 9.580e$$-$$03 7.517e$$-$$03 5.566e$$-$$03 5.567e$$-$$03 256 6.404e$$-$$03 4.854e$$-$$03 2.442e$$-$$03 2.409e$$-$$03 512 6.090e$$-$$03 4.288e$$-$$03 1.036e$$-$$03 1.002e$$-$$03 Table 2 Error in the $$l^\infty$$ norm of the whole computational domain at time $$t=0.1$$ for the time-dependent Examples 7.10 and 7.11 Errors and order, Example 7.10 N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 4.845e$$-$$02 6.691e$$-$$02 4.894e$$-$$02 4.888e$$-$$02 64 4.432e$$-$$02 4.607e$$-$$02 2.977e$$-$$02 2.975e$$-$$02 128 3.544e$$-$$02 2.823e$$-$$02 2.457e$$-$$02 2.438e$$-$$02 256 2.971e$$-$$02 2.080e$$-$$02 1.747e$$-$$02 1.724e$$-$$02 512 2.764e$$-$$02 1.652e$$-$$02 1.205e$$-$$02 1.182e$$-$$02 Errors and order, Example 7.11 N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 2.182e$$-$$02 1.449e$$-$$02 1.985e$$-$$02 1.985e$$-$$02 64 1.435e$$-$$02 1.160e$$-$$02 1.279e$$-$$02 1.279e$$-$$02 128 9.580e$$-$$03 7.517e$$-$$03 5.566e$$-$$03 5.567e$$-$$03 256 6.404e$$-$$03 4.854e$$-$$03 2.442e$$-$$03 2.409e$$-$$03 512 6.090e$$-$$03 4.288e$$-$$03 1.036e$$-$$03 1.002e$$-$$03 Errors and order, Example 7.10 N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 4.845e$$-$$02 6.691e$$-$$02 4.894e$$-$$02 4.888e$$-$$02 64 4.432e$$-$$02 4.607e$$-$$02 2.977e$$-$$02 2.975e$$-$$02 128 3.544e$$-$$02 2.823e$$-$$02 2.457e$$-$$02 2.438e$$-$$02 256 2.971e$$-$$02 2.080e$$-$$02 1.747e$$-$$02 1.724e$$-$$02 512 2.764e$$-$$02 1.652e$$-$$02 1.205e$$-$$02 1.182e$$-$$02 Errors and order, Example 7.11 N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 2.182e$$-$$02 1.449e$$-$$02 1.985e$$-$$02 1.985e$$-$$02 64 1.435e$$-$$02 1.160e$$-$$02 1.279e$$-$$02 1.279e$$-$$02 128 9.580e$$-$$03 7.517e$$-$$03 5.566e$$-$$03 5.567e$$-$$03 256 6.404e$$-$$03 4.854e$$-$$03 2.442e$$-$$03 2.409e$$-$$03 512 6.090e$$-$$03 4.288e$$-$$03 1.036e$$-$$03 1.002e$$-$$03 Example 7.11 (Dirichlet BC) We consider the exact solution $u(x,y,t) = t + \frac{3}{4}\left(\frac{b}{a}x^2+\frac{a}{b}y^2\right)^{2/3},$ with $$a = 2$$ and $$b = 1$$ (see Lemma 2.12). This is the same example as in Example 7.10, but we consider Dirichlet boundary conditions instead. Therefore, we prescribe the exact solution at a seven-point layer at the boundary for all time $$t$$. This way all schemes are initialized at the same set of grid points, and thus we can compare their accuracy. The error in the $$l^\infty$$ norm at time $$T = 0.1$$ is presented in Table 2. When using Neumman boundary conditions, we observed slow convergence, with errors near $$0.01$$, slowly decreasing as the grid size improved. The accuracy improved as we went from the narrow to the wider elliptic scheme, and further improved as we went to the accurate and then the filtered scheme. In the case of Dirichlet boundary conditions, the accuracy is better overall and the error decrease is slightly faster. The difference is explained by the cap off done in the Neumman boundary conditions that introduces an additional error in the solution. However, this error does not propagate to the whole domain as the level sets of the solution shrink to its interior, and so, when we look at the error away from the cap off, we recover results very similar to the Dirichlet boundary conditions. 7.4 Numerical study of the morphology and affine invariance properties In this section, we test whether our proposed schemes satisfy numerically the morphology and affine invariance properties that (AC) satisfies (see Theorem 2.10). Example 7.12 In this example, we test numerically whether the schemes presented here satisfy the morphology property of the affine curvature evolution ((2) in Theorem 2.10). We consider two examples: $$(a)$$$$g_1(x) = {\rm e}^x$$, $$(b)$$$$g_2(x) = x^3$$. We take $u_0(x,y) = \min\left\{\left(\frac{x}{2}\right)^2+y^2-1,0\right\}$ and compare $${\it {\Phi}}_t(g_v\circ u_0)$$ with $$g_v \circ {\it {\Phi}}_t(u_0)$$ at $$t=1$$ for $$v = 1,2$$. We took $$[-3,3]^2$$ as the computational domain with homogeneous Neumann boundary conditions. The results are displayed in Table 3. The difference in the $$l^\infty$$ norm is one order of magnitude smaller than the observed accuracy for the schemes in Section 7.2. Based on these examples, the morphology property seems to hold numerically. Table 3 Difference in the $$l^\infty$$ norm between $${\it {\Phi}}_t(g_v \circ u_0)$$ and $$g_v\circ {\it {\Phi}}_t(u_0)$$ for $$v=1,2$$ for Example 7.12 Difference in the $$l^\infty$$ norm, Example 7.12 (a) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 8.943e$$-$$03 1.037e$$-$$02 3.419e$$-$$03 3.446e$$-$$03 64 5.709e$$-$$03 6.111e$$-$$03 1.954e$$-$$03 1.970e$$-$$03 128 4.061e$$-$$03 3.257e$$-$$03 1.061e$$-$$03 1.109e$$-$$03 256 3.115e$$-$$03 1.871e$$-$$03 5.907e$$-$$04 6.109e$$-$$04 512 2.604e$$-$$03 1.161e$$-$$03 3.207e$$-$$04 3.397e$$-$$04 Difference in the $$l^\infty$$ norm, Example 7.12 (b) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.450e$$-$$02 7.023e$$-$$02 6.947e$$-$$03 6.947e$$-$$03 64 1.730e$$-$$02 2.842e$$-$$02 1.881e$$-$$03 1.877e$$-$$03 128 1.032e$$-$$02 9.355e$$-$$03 6.254e$$-$$04 6.333e$$-$$04 256 6.894e$$-$$03 4.307e$$-$$03 2.283e$$-$$04 2.307e$$-$$04 512 5.372e$$-$$03 2.316e$$-$$03 8.343e$$-$$05 8.506e$$-$$05 Difference in the $$l^\infty$$ norm, Example 7.12 (a) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 8.943e$$-$$03 1.037e$$-$$02 3.419e$$-$$03 3.446e$$-$$03 64 5.709e$$-$$03 6.111e$$-$$03 1.954e$$-$$03 1.970e$$-$$03 128 4.061e$$-$$03 3.257e$$-$$03 1.061e$$-$$03 1.109e$$-$$03 256 3.115e$$-$$03 1.871e$$-$$03 5.907e$$-$$04 6.109e$$-$$04 512 2.604e$$-$$03 1.161e$$-$$03 3.207e$$-$$04 3.397e$$-$$04 Difference in the $$l^\infty$$ norm, Example 7.12 (b) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.450e$$-$$02 7.023e$$-$$02 6.947e$$-$$03 6.947e$$-$$03 64 1.730e$$-$$02 2.842e$$-$$02 1.881e$$-$$03 1.877e$$-$$03 128 1.032e$$-$$02 9.355e$$-$$03 6.254e$$-$$04 6.333e$$-$$04 256 6.894e$$-$$03 4.307e$$-$$03 2.283e$$-$$04 2.307e$$-$$04 512 5.372e$$-$$03 2.316e$$-$$03 8.343e$$-$$05 8.506e$$-$$05 Table 3 Difference in the $$l^\infty$$ norm between $${\it {\Phi}}_t(g_v \circ u_0)$$ and $$g_v\circ {\it {\Phi}}_t(u_0)$$ for $$v=1,2$$ for Example 7.12 Difference in the $$l^\infty$$ norm, Example 7.12 (a) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 8.943e$$-$$03 1.037e$$-$$02 3.419e$$-$$03 3.446e$$-$$03 64 5.709e$$-$$03 6.111e$$-$$03 1.954e$$-$$03 1.970e$$-$$03 128 4.061e$$-$$03 3.257e$$-$$03 1.061e$$-$$03 1.109e$$-$$03 256 3.115e$$-$$03 1.871e$$-$$03 5.907e$$-$$04 6.109e$$-$$04 512 2.604e$$-$$03 1.161e$$-$$03 3.207e$$-$$04 3.397e$$-$$04 Difference in the $$l^\infty$$ norm, Example 7.12 (b) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.450e$$-$$02 7.023e$$-$$02 6.947e$$-$$03 6.947e$$-$$03 64 1.730e$$-$$02 2.842e$$-$$02 1.881e$$-$$03 1.877e$$-$$03 128 1.032e$$-$$02 9.355e$$-$$03 6.254e$$-$$04 6.333e$$-$$04 256 6.894e$$-$$03 4.307e$$-$$03 2.283e$$-$$04 2.307e$$-$$04 512 5.372e$$-$$03 2.316e$$-$$03 8.343e$$-$$05 8.506e$$-$$05 Difference in the $$l^\infty$$ norm, Example 7.12 (a) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 8.943e$$-$$03 1.037e$$-$$02 3.419e$$-$$03 3.446e$$-$$03 64 5.709e$$-$$03 6.111e$$-$$03 1.954e$$-$$03 1.970e$$-$$03 128 4.061e$$-$$03 3.257e$$-$$03 1.061e$$-$$03 1.109e$$-$$03 256 3.115e$$-$$03 1.871e$$-$$03 5.907e$$-$$04 6.109e$$-$$04 512 2.604e$$-$$03 1.161e$$-$$03 3.207e$$-$$04 3.397e$$-$$04 Difference in the $$l^\infty$$ norm, Example 7.12 (b) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.450e$$-$$02 7.023e$$-$$02 6.947e$$-$$03 6.947e$$-$$03 64 1.730e$$-$$02 2.842e$$-$$02 1.881e$$-$$03 1.877e$$-$$03 128 1.032e$$-$$02 9.355e$$-$$03 6.254e$$-$$04 6.333e$$-$$04 256 6.894e$$-$$03 4.307e$$-$$03 2.283e$$-$$04 2.307e$$-$$04 512 5.372e$$-$$03 2.316e$$-$$03 8.343e$$-$$05 8.506e$$-$$05 Example 7.13 In this example, we do a qualitative test of the affine invariance property, which in practice is what one needs for applications in image analysis. To do so, we plot the level sets of the affine invariant motion by curvature (AC) with $u(x,y) = \left(\frac{x}{2}\right)^2+y^2-1$ and $$u \circ \phi$$ as the initial solutions. For the affine transformations $$\phi(\mathbf{x}) = A\mathbf{x}$$, we consider $(a) \text{(rotation by \pi/4)} A = \begin{bmatrix}\frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2}\\[4pt] -\frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2}\end{bmatrix}, \quad (b) \: A = \begin{bmatrix}\frac{1}{2} & 1\\[4pt] 1 & \frac{1}{2}\end{bmatrix}.$ We take $$[-5,5]^2$$ as the computational domain on a $$256 \times 256$$ grid. In Fig. 7, we plot the zero level set of $${\it {\Phi}}_t(u\circ \phi)$$ and $$\left({\it {\Phi}}_{t(\det \phi)^{2/3}}(u)\right) \circ \phi$$ at $$t=1$$. The filtered scheme provided the best results, being indistinguishable to the naked eye. For a long time, the elliptic scheme did not provide as good results, a consequence of its lower accuracy (see Section 7.2). As for the standard finite difference scheme only in (a) the difference is indistinguishable to the naked eye like the filtered scheme. For (b), where $$A$$ is not a special affine transformation, the difference is significant, but can be removed by taking time steps of half the size (For shorter times all the curves were very close.) We point out that when the affine transformation is a rotation by a multiple of $$\pi/2$$ or a reflection over a line $$L$$ that makes an angle multiple of $$\pi/4$$ with the x-axis, we obtain essentially machine precision. (The difference in the $$l^\infty$$ norm was of the order $$10^{-10}$$.) This is expected, because our stencil is invariant under these transformations. Fig. 7. View largeDownload slide Plot of the zero level sets in Example 7.13 of $${\it {\Phi}}_t(u\circ \phi)$$ and $$\left({\it {\Phi}}_{t(\det \phi)^{2/3}}(u)\right) \circ \phi$$ for regular elliptic scheme (left), standard scheme (center) and regular filtered scheme (right) at time $$t=1$$ with $$\phi$$ given by $$(a)$$ (top) and $$(b)$$ (bottom). Fig. 7. View largeDownload slide Plot of the zero level sets in Example 7.13 of $${\it {\Phi}}_t(u\circ \phi)$$ and $$\left({\it {\Phi}}_{t(\det \phi)^{2/3}}(u)\right) \circ \phi$$ for regular elliptic scheme (left), standard scheme (center) and regular filtered scheme (right) at time $$t=1$$ with $$\phi$$ given by $$(a)$$ (top) and $$(b)$$ (bottom). 8. Conclusions We presented a convergent finite difference discretization of the PDE for motion of level sets by affine curvature in two dimensions. Computational examples demonstrated that the standard finite difference method is unstable, which motivates the need for a convergent method. The foundation of the scheme used an existing wide stencil discretization of the mean curvature operator, combined with an elliptic discretization of the positive and negative eikonal operators, $$\pm \vert{\nabla u}\vert$$. However, explicit time discretizations require Lipschitz continuous operators, which the affine curvature operator fails to be. Thus, we approximated it by a Lipschitz continuous regularization. In theory, the explicit Euler discretization is stable using a time step $$\text{d}\textit{t} \leq \mathcal{O}(h^{22/9})$$, with the constant determined by the width of the stencil. In practice, we achieved numerically equivalent results using $${\rm d}t = h^2/2$$ and without the regularization, although there is no proof of stability with the less restrictive time step. A careful choice of the regularization parameters allowed for the regularized elliptic scheme to maintain the same order of accuracy as the unregularized scheme while being provably convergent. The lower accuracy of both schemes, which results from the singularity of the operator, is overcome by the use of the filtered schemes, which essentially attain the accuracy of the standard finite difference schemes while being provably convergent. Simulations demonstrated the geometric properties of the PDE were nearly preserved by the numerical solutions, including affine invariance, morphological properties and the accurate representation of the shrinking ellipses. Simulations validated the convergence of the elliptic scheme and the improved accuracy of the filtered scheme. Funding FCT doctoral grant (SFRH/BD/84041/2012 to T.S.), in part. References Angenent S. , Sapiro G. & Tannenbaum A. ( 1998 ) On the affine heat equation for non-convex curves. J. Amer. Math. Soc. , 11 , 601 – 634 . Google Scholar CrossRef Search ADS Barles G. & Souganidis P. E. ( 1991 ) Convergence of approximation schemes for fully nonlinear second order equations. Asymptot. Anal. , 4 , 271 – 283 . 92d:35137 Barnett V. ( 1976 ) The ordering of multivariate data. J. R. Stat. Soc. , 139 , 318 – 355 . Calder J. & Smart C. K. ( 2016 ) The limit shape of convex peeling (in preparation). Calder J. , Esedoglu S. & Hero A. O. ( 2014 ) A Hamilton–Jacobi equation for the continuum limit of non-dominated sorting. SIAM J. Math. Anal. , 46 , 603 – 638 . Google Scholar CrossRef Search ADS Cao F. , ( 2003 ) Geometric Curve Evolution and Image Processing . Lecture Notes in Mathematics , vol. 1805 . Berlin : Springer . Google Scholar CrossRef Search ADS Catté F. Dibos F. & Koepfler G. ( 1995 ) A morphological scheme for mean curvature motion and applications to anisotropic diffusion and motion of level sets. SIAM J. Numer. Anal. , 32 , 1895 – 1909 . Google Scholar CrossRef Search ADS Chazelle B. ( 1985 ) On the convex layers of a planar set. IEEE Trans. Inform. Theory , 31 , 509 – 517 . Google Scholar CrossRef Search ADS Ciomaga A. , Monasse P. & Morel J.-M. ( 2011 ) Image visualization and restoration by curvature motions. SIAM Multiscale Model. Simul. , 9 , 834 – 871 . Google Scholar CrossRef Search ADS Crandall M. G. , Ishii H. & Lions P.-L. ( 1992 ) User’s guide to viscosity solutions of second order partial differential equations. Bull. Amer. Math. Soc. (N.S.) , 27 , 1 – 67 . Google Scholar CrossRef Search ADS Deckelnick K. , Dziuk G. & Elliott C. M. ( 2005 ) Computation of geometric partial differential equations and mean curvature flow. Acta Numer. 14 , 139 – 232 . Google Scholar CrossRef Search ADS Esedoglu S. , Ruuth S. & Tsai R. ( 2010 ) Diffusion generated motion using signed distance functions. J. Comput. Phys. , 229 , 1017 – 1042 . Google Scholar CrossRef Search ADS Evans L. C. & Spruck J. ( 1999 ) Motion of level sets by mean curvature, I. Fundamental Contributions to the Continuum Theory of Evolving Phase Interfaces in Solids ( Ball J. M. , Kinderlehrer D. , Podio-Guidugli P. & Slemrod M. ). Berlin : Springer , pp. 328 – 374 . Froese B. D. & Oberman A. M. ( 2013 ) Convergent filtered schemes for the Monge–Ampère partial differential equation. SIAM J. Numer. Anal. , 51 , 423 – 444 . Google Scholar CrossRef Search ADS Gage M. E. ( 1984 ) Curve shortening makes convex curves circular. Invent. Math. , 76 , 357 – 364 . Google Scholar CrossRef Search ADS Giga Y. ( 2006 ) Surface Evolution Equations . Birkhäuser Basel : Springer . Guichard F. ( 1994 ) Axiomatisation des analyses multi-echelles d’images et de films. Ph.D. Thesis , Université Paris Dauphine . Guichard F. & Morel J.-M. ( 1997 ) Partial Differential Equations and Image Iterative Filtering . Institute of Mathematics and its Applications Conference Series , vol. 63 . New York : Oxford University Press , pp. 525 – 562 . Liu R. Y. , Parelius J. M. , Singh K. ( 1999 ) Multivariate analysis by data depth: descriptive statistics, graphics and inference (with discussion and a rejoinder by liu and singh). Ann. Stat. , 27 , 783 – 858 . Merriman B. , Bence J. K. & Osher S. J. ( 1994 ) Motion of multiple junctions: a level set approach. J. Comput. Phys. , 112 , 334 – 363 . Google Scholar CrossRef Search ADS Moisan L. ( 1998 ) Affine plane curve evolution: a fully consistent scheme. IEEE Trans. Image Process. , 7 , 411 – 420 . Google Scholar CrossRef Search ADS PubMed Oberman A. M. ( 2004 ) A convergent monotone difference scheme for motion of level sets by mean curvature. Numer. Math. , 99 , 365 – 379 . Google Scholar CrossRef Search ADS Oberman A. M. ( 2006 ) Convergent difference schemes for degenerate elliptic and parabolic equations: Hamilton–Jacobi equations and free boundary problems. SIAM J. Numer. Anal. , 44 , 879 – 895 (electronic). Google Scholar CrossRef Search ADS Oberman A. M. & Salvador T. ( 2015 ) Filtered schemes for Hamilton–Jacobi equations: a simple construction of convergent accurate difference schemes. J. Comput. Phys. , 284 , 367 – 388 . Google Scholar CrossRef Search ADS Oberman A. M. & Salvador T. ( 2016 ) Numerical method for motion by affine curvature . arXiv preprint arXiv:1610.08831 . Osher S. & Fedkiw R. P. ( 2003 ) Level Set Methods and Dynamic Implicit Surfaces , vol. 153. New York: Springer. Osher S. & Sethian J. A. ( 1988 ) Fronts propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations. J. Comput. Phys. , 79 , 12 – 49 . Google Scholar CrossRef Search ADS Sapiro G. ( 2006 ) Geometric Partial Differential Equations and Image Analysis . Cambridge : Cambridge University Press . Sapiro G. & Tannenbaum A. ( 1994 ) On affine plane curve evolution. J. Funct. Anal. , 119 , 79 – 120 . Google Scholar CrossRef Search ADS Sethian J. A. ( 1999 ) Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science , vol. 3 . Cambridge : Cambridge University Press . Su B. ( 1983 ) Affine Differential Geometry . Beijing : Publisher Science Press ; New York : Gordon & Breach Science Publishers . © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

# Numerical methods for motion of level sets by affine curvature

, Volume Advance Article – Aug 19, 2017
33 pages

/lp/ou_press/numerical-methods-for-motion-of-level-sets-by-affine-curvature-EaHP124kGo
Publisher
Oxford University Press
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drx045
Publisher site
See Article on Publisher Site

### Abstract

Abstract We study numerical methods for the nonlinear partial differential equation that governs the motion of level sets by affine curvature. We show that standard finite difference schemes are nonlinearly unstable. We build convergent finite difference schemes using the theory of viscosity solutions. We demonstrate that our approximate solutions capture the affine invariance and morphological properties of the evolution. Numerical experiments demonstrate the accuracy and stability of the discretization. 1. Introduction 1.1 The affine curvature partial differential equation The affine curvature evolution is one of the most fundamental geometric evolution equations after the mean curvature evolution. It was introduced by Sapiro & Tannenbaum (1994) and Angenent et al. (1998), and has applications in mathematical morphology, edge detection, image smoothing and image enhancement (see Sapiro, 2006). We are motivated by the recent work of Calder & Smart (2016), which provides an application of the affine curvature partial differential equation (PDE) to the statistics of large data sets. The convex hull peeling algorithm (Chazelle, 1985) provides an affine invariant notion of the median and the quantiles of multidimensional probability distributions. Although there is more than one way to measure data depth (Barnett, 1976), affine invariance is an important property for such measures (Liu et al., 1999). The level sets of the solution of the affine curvature PDE, with the right-hand side given by the probability density $$\rho$$, also gives an affine invariant notion of the depth of $$\rho$$. According to Calder & Smart (2016), these two notions of depth are equivalent: the rescaled data depth layers of $$N$$ data points sampled from the density, $$\rho$$, given by convex hull peeling algorithm converge, in the limit $$N \to \infty$$, to the levels given by the solution of the PDE. Compared with convex hull peeling, the PDE characterization is efficient in terms of the number of data points $$N$$: an efficient density estimation method can be used to approximate $$\rho$$, and afterward the PDE solve does not depend on $$N$$. This kind of limiting PDE approach has already been shown to be effective for nondominated sorting (Calder et al., 2014). The planar motion of level sets by affine curvature is governed by the nonlinear PDE $$\label{ACPDE} u_t = F[u] := \left\vert{\nabla u}\right\vert \left (k[u]\right)^{1/3}.$$ (AC) Here $$u = u(x,y):\mathbb{R}^2 \to \mathbb{R}$$, $$\nabla u = (u_x,u_y)$$ denotes the gradient of $$u$$, and $$k[u]$$ denotes the curvature of the level set of $$u$$ $$\label{curvature} k[u] = \text{div}\left(\frac{\nabla u}{\left\vert{\nabla u}\right\vert}\right) = \frac{u_{xx}u_y^2-2u_x u_y u_{xy}+u_{yy}u_x^2}{(u_x^2+u_y^2)^{3/2}}.$$ (1.1) The affine curvature PDE is closely related to the well-known PDE for motion of level sets by mean curvature $$u_t = {\it {\Delta}}_1 u := \left\vert{\nabla u}\right\vert k[u]$$ (MC) studied in the seminal article Osher & Sethian (1988). However, the PDE (AC) exhibits instabilities not found in the mean curvature PDE (MC), as demonstrated below. To resolve these instabilities, we propose a Lipschitz regularization of the operator. The regularized operator is also a geometric PDE, and viscosity solutions converge to solutions of the affine curvature PDE in the limit as the regularization parameter goes to zero [Giga, 2006, Theorem 4.6.1]. The advantage of the regularized operator is that it allows us to build stable, convergent explicit solvers that are not otherwise available. Moreover, the resulting discretization can be combined into a filtered scheme that achieves the higher accuracy of the otherwise unstable standard finite difference scheme. In addition, the numerical solutions exhibit the affine invariance and morphology properties of the evolution. Our approach to the convergent discretization In this paragraph, we present an overview of our approach to building an elliptic discretization for clarity. The details and supporting theory can be found in the sections that follow. Elliptic discretizations are available for $${\it {\Delta}}_1 u$$, rather than for $$k[u]$$, so we rewrite \begin{align}\label{AffPQ} F[u] & = A(\left\vert{\nabla u}\right\vert,{\it {\Delta}}_1 u) = (\left\vert{\nabla u}\right\vert^2 {\it {\Delta}}_1 u)^{1/3}, && \text{where} A(p,q) = \left(p^2 q \right)^{1/3}. \end{align} (1.2) The goal is to make use of available elliptic discretizations of $$\pm \left\vert{\nabla u}\right\vert$$ and $${\it {\Delta}}_1 u$$ to build an elliptic discretization of the full operator $$F[u]$$. However, simply inserting these operators into the function $$A(p,q)$$ is not sufficient: elliptic schemes are built by composing nondecreasing maps with elliptic operators, and $$A(p,q)$$ fails to be nondecreasing. Furthermore, explicit time discretizations require Lipschitz continuous operators, and $$A(p,q)$$ also fails to be Lipschitz continuous. Because the properties of the nonlinear and singular function $$A(p,q)$$ are so important, we first study a model equation in one dimension. Define $$\label{AffPQ1D} F^{1D}[u] := A(u_x, u_{xx}) = \left(\left\vert{u_x}\right\vert^2 \:u_{xx}\right)^{1/3},$$ (1.3) so that the $$F^{1D}[u]$$ operator has the same homogeneity in the first and second derivatives as the $$F[u]$$ operator. Similar to the higher dimensional PDE, the model equation exhibits the instability of standard finite differences. To build an elliptic scheme for the one-dimensional equation, what is needed is a nondecreasing representation of the function $$A(p,q)$$, which is consistent with $$F^{1D}[u]$$. Furthermore, we need a Lipschitz continuous approximation of $$A(p,q)$$, with Lipchitz constant $$K^h$$, to build a monotone discretization of the time-dependent PDE, using a time step $${\rm d}t \le 1/K^h$$. Once this modified function is available, we proceed by inserting the discretization of the two-dimensional operators into the modified function, which results in a convergent scheme. Related numerical work There are by now a large number of numerical methods for the level set mean curvature PDE (MC) (see review papers Deckelnick et al., 2005; Ciomaga et al., 2011). In particular, Catté and Dibos proposed a morphological scheme that satisfies the comparison principle (Cattè et al., 1995) and the first author presented a convergent wide stencil finite difference scheme (Oberman, 2004) using a median formula. For the affine curvature evolution, the recent article Esedoglu et al. (2010) gives a Bence–Merriman–Osher (Merriman et al., 1994) thresholding scheme. Similar to the scheme proposed here, it is monotone, but introduced a different regularization of the cube root, which was needed for theoretical purposes, but not in practice. Thresholding methods are effective for moving a given curve by the evolution and allow for large time steps to be taken. Alvarez and Guichard proposed a local scheme that lacks the affine invariance property (Guichard, 1994). A morphological scheme that generalized Cattè et al. (1995) was proposed by Guichard and Morel for affine curvature (Guichard & Morel, 1997). This inf–sup scheme, although morphologically invariant, has some limitations on the speed at which the level set curves move. In Moisan (1998), a nonlocal geometric morphological scheme is presented for the evolution of a single curve. These are approximated by polygons, and the scheme is shown to be monotone and affine invariant for polygons. Novelty of the work The scheme presented here uses the level set representation (AC). Thus, it moves every level set by the evolution, unlike thresholding methods, but requires a much smaller time step. It is monotone, a desirable property for the convergence, but also for its applications. Moreover, it is shown numerically that it captures the affine invariance and morphological properties. 1.2 Background mathematics Euclidean and Affine Curvature We begin with a parametric description of the affine curvature evolution and make a comparison with the more familiar mean curvature evolution. We refer to Sapiro (2006, Chapter 2) for more details. Consider a curve described parametrically $$\mathcal{C}(s):[a,b] \in \mathbb{R} \to \mathbb{R}^2$$, where $$s$$ parameterizes the curve. If the curve is parameterized by the Euclidean arc length, $$\left\vert{\frac{\text{d} \mathcal{C}} {\text{d} s}}\right\vert = 1,$$ then the curvature, $$k$$, is defined, up to a sign, by the equation $$\mathcal{C}_{ss} = k\mathbf{n}$$, where $$\mathbf{n}$$ denotes the Euclidean normal of the curve. Under mean curvature, the curve evolves with velocity $$k\mathbf{n}$$. The affine curvature arises from a different parameterization of the curve. Define the parameter, $${r}$$ by the condition that the vectors $$C_{r}$$ and $$C_{{r} {r}}$$ form a parallelogram of area $$1$$, $$[\mathcal{C}_{r},\mathcal{C}_{{r} {r}}] = 1$$. (Here the brackets denote the determinant of the matrix whose columns are given by those vectors.) Then the affine curvature of the curve is defined by $$\mu = [\mathcal{C}_{{r} {r}},\mathcal{C}_{{r} {r} {r}}]$$. Affine differential geometry is not defined for nonconvex curves. However, we can still define the affine curvature evolution using the Euclidean curvature by taking the velocity to be $$k^{1/3}\mathbf{n}$$. The affine curvature is the simplest nontrivial affine invariant of the curve (Su, 1983). Under the affine curvature evolution, any convex curve remains convex; any smooth curve remains smooth; any convex smooth curve evolves to an ellipse until it collapses to a single point and any smooth curve becomes convex after a certain time. Moreover, the affine curvature evolution is invariant under the class of special affine transformations, which are defined by matrices with determinant $$1$$ (see Theorem 2.10 below). When compared with mean curvature evolution, this shrinks curves to circles (Gage, 1984) and is invariant under the smaller class of orthogonal transformations. Level Set PDE formulation The level set method (Sethian, 1999; Osher & Fedkiw, 2003) for the affine curvature evolution results in the PDE (AC). The level set method has the following advantages compared with parameterized curve evolution: (i) it provides a natural generalization of the flows when the curve becomes singular and notions such as normals are not well defined; (ii) there is no need to track topological changes, because they are discovered when the corresponding level set is computed; and (iii) it can be discretized on a uniform grid, which is convenient for many applications. In the level set method, a curve is represented implicitly as the level set of the auxiliary function $$u(x,y,t):\mathbb{R}^2\times \mathbb{R} \to\mathbb{R}$$, that is, $\mathcal{C}(t) = \left\{(x,y)\in \mathbb{R}^2 \mid u(x,y,t) = c\right\}$ for some arbitrary constant $$c\in \mathbb{R}$$. If $$u$$ satisfies $$u_t = \left\vert{\nabla u}\right\vert \beta[u(\cdot,t)]$$, for some function $$\beta$$, which depends on the level set of $$u$$, then all its level sets move in the normal direction with speed $$\beta$$. For example, choosing $$\beta = 1$$, we obtain the time-dependent eikonal equation, $$u_t = \left\vert{\nabla u}\right\vert$$. Taking $$\beta = k[u]$$ or $$\left(k[u]\right)^{1/3}$$, we obtain (MC) and (AC), respectively. 2. The affine curvature PDE 2.1 Definition of viscosity solutions Viscosity solutions (Crandall et al., 1992) provide the correct notion of weak solution to a class of degenerate elliptic PDEs, which includes (AC) and (MC). The article Evans & Spruck (1999) focuses on the mean curvature equation. The book Giga (2006) established existence and uniqueness of viscosity solutions in a bounded domain, with Neumann boundary conditions, see Section 3.6 for the Comparison Principle, and Theorem 4.6.1 for the convergence of approximations. It is also shown that the Lipschitz constant of the solution is preserved by the evolution (Section 3.5). See also the book Cao (2003) for viscosity solutions for geometric evolution equations and numerical methods. Definition 2.1 The function $$F: \mathcal{S}^{n} \times \mathbb{R}^n \times \mathbb{R} \times {\it {\Omega}} : \to \mathbb{R}$$ is proper and degenerate elliptic (in the sense of Crandall et al., 1992) if $F(M,p,r,x) \le F(N,p,s,x), \quad \text {for all} M \succeq N, r \le s, \text{and all} x \in {\it {\Omega}}, p \in \mathbb{R}^n,$ where $$Y \preceq X$$ if $$d^\intercal Yd \le d^\intercal X d$$ for all $$d\in \mathbb{R}^n$$. Remark 2.2 By this convention, the Laplacian operator is elliptic when written $$F(M) = -\text{tr}(M)$$. Some authors use the other convention, without the minus sign. Lemma 2.3 The operator $$-F[u]$$ is degenerate elliptic. Proof. We start with the observation that using (1.1), we can write ${\it {\Delta}}_1 u = u_{\mathbf{t}\mathbf{t}}, \qquad \mathbf{t} = \frac{(-u_y,u_x)}{(u_x^2+u_y^2)^{1/2}},$ where $$\mathbf{t}$$ is the (Euclidean) unit tangent. Then, it is clear that $$-{\it {\Delta}}_1$$ is degenerate elliptic. Using the representation (1.2) the degenerate ellipticity of $$-F$$ follows. □ We now give the definition of viscosity solutions on a domain $${\it {\Omega}} \subset \mathbb{R}^2$$. We follow Giga (2006). Let $$T > 0$$. We are interested in the Cauchy problem $$\label{cauchyproblem} \begin{cases} u_t = F[u] & \text{in} (x,t) \in {\it {\Omega}} \times (0,T),\\ B(x,p) := \nu(x)^\intercal p = 0 & \text{on} (x,t) \in \partial {\it {\Omega}} \times (0,T),\\ u(x,0) = u_0(x) & \text{in} x \in \overline{{\it {\Omega}}}, \end{cases}$$ (2.1) where $$\nu$$ is the outward unit normal of $$\partial {\it {\Omega}}$$. For completeness, we recall the definition of upper and lower semicontinuous functions. Definition 2.4 A function $$u:U\to\mathbb{R}$$ is upper (resp. lower) semicontinuous at $$x\in U$$ provided $\limsup_{y \in U\to x} u(y) \leq u(x) \quad \left(\text{resp.} \liminf_{y \in U\to x} u(y) \geq u(x)\right)\!.$ We denote by $$USC(U)$$ (resp. $$LSC(U)$$) the collection of functions that are upper (resp. lower) semicontinuous at all points of $$U$$. Definition 2.5 A function $$u \in USC(\overline{{\it {\Omega}}} \times [0,T])$$ is called a viscosity subsolution of (2.1) if $$u(x,0) \leq u_0(x)$$ for $$x \in \overline{{\it {\Omega}}}$$ and for any $$\varphi \in C^2(\overline{{\it {\Omega}}}\times[0,T])$$ such that $$u-\varphi$$ has a local maximum at $$(x,t) \in \overline{{\it {\Omega}}} \times (0,T)$$ then $\begin{cases} \varphi_t(x,t) - F[\varphi](x,t) \leq 0 & \text{if}~ x \in {\it {\Omega}},\\ \min\{\varphi_t(x,t) - F[\varphi](x,t), B(x,\nabla \phi(x,t))\} \leq 0 & \text{if}~ \in \partial {\it {\Omega}}. \end{cases}$ Similarly, $$u \in LSC(\overline{{\it {\Omega}}} \times [0,T])$$ is called a viscosity supersolution of (2.1) if $$u(x,0) \geq u_0(x)$$ for $$x \in \overline{{\it {\Omega}}}$$ and for any $$\varphi \in C^2(\overline{{\it {\Omega}}}\times[0,T])$$ such that $$u-\varphi$$ has a local minimum at $$(x,t) \in \overline{{\it {\Omega}}} \times (0,T)$$ then $\begin{cases} \varphi_t(x,t) - F[\varphi](x,t) \geq 0 & \text{if}~ x \in {\it {\Omega}},\\ \max\{\varphi_t(x,t) - F[\varphi](x,t), B(x,\nabla \phi(x,t))\} \geq 0 & \text{if}~ \in \partial {\it {\Omega}}. \end{cases}$ Finally, we call $$u$$ a viscosity solution of (2.1) if $$u$$ is both a viscosity subsolution and a viscosity supersolution. We have the following uniqueness result, from Giga (2006). Theorem 2.6 (Comparison Principle) Suppose that $${\it {\Omega}}$$ is convex with $$C^2$$ boundary $$\partial {\it {\Omega}}$$. Let $$u$$ and $$v$$ be a viscosity subsolution and supersolution of (2.1). Then $$u \leq v$$ in $${\it {\Omega}} \times (0,T)$$. Remark 2.7 We will also be interested in the static equation $$F[u] = f$$. The definition of viscosity solution and comparison principle presented here for the parabolic equation generalize to the static equation. 2.2 Invariance properties of the PDE We present some properties of (AC), starting with some properties of the affine curvature operator $$F[u]$$. The proofs can be found in Oberman & Salvador (2016). Lemma 2.8 Let $$u \in C^2(\mathbb{R}^2)$$. Then the operator $$F$$ has the following properties: (1) Rescaling: for $$h > 0$$, define $$v(x,y) := u(x/h,y/h)$$ $F[v] = h^{-4/3}F[u].$ (2) Morphology: for $$g \in C^1(\mathbb{R})$$, $F[g \circ u] = g^\prime(u) F[u].$ (3) Affine invariance: for any affine map $$\phi(\mathbf{x}) = A\mathbf{x}+\textbf{b}$$, $F[u \circ \phi] = (\det A)^{2/3} F[u] \circ \phi.$ Remark 2.9 For $${\it {\Delta}}_1 u$$, rescaling gives $${\it {\Delta}}_1 v = h^{-2} {\it {\Delta}}_1 u$$ and invariance only holds for orthogonal transformations. Using the Lemma, the properties can be generalized for the affine curvature evolution (AC). Theorem 2.10 Consider $${\it {\Phi}}_t:u_0\mapsto u(\cdot,t)$$ the solution map of (AC). Then $${\it {\Phi}}_t$$ satisfies the following properties: (1) Monotonicity$$:$$$$u \leq v \Rightarrow {\it {\Phi}}_t(u) \leq {\it {\Phi}}_t(v)$$. (2) Morphology/Relabelling: for any monotone scalar function $$g$$, ${\it {\Phi}}_t(g \circ u) = g \circ {\it {\Phi}}_t(u).$ (3) Affine invariance: for any affine map $$\phi(\mathbf{x}) = A\mathbf{x}+\textbf{b}$$, ${\it {\Phi}}_t(u\circ \phi) = \left({\it {\Phi}}_{t(\det A)^{2/3}}(u)\right) \circ \phi.$ Remark 2.11 Note that our rescaling factor property $$(iii)$$ differs from the one in Moisan (1998); however, both formulas agree (and give unity) for special affine transformations that have determinant $$1$$. 2.3 An exact solution We now give an example of an exact solution for (AC). A proof can be found in Oberman & Salvador (2016). See also Giga (2006, Section 1.7.4). Lemma 2.12 Given $$a,b > 0$$, let $$u:\mathbb{R}^2\times [0,\infty)$$ be given by the shifted ellipse $u(x,y,t) = t + \frac{3}{4}(ab)^{2/3}\left(\left(\frac{x}{a}\right)^2+\left(\frac{y}{b}\right)^2\right)^{2/3} = t + \frac{3}{4}\left(\frac{b}{a}x^2+\frac{a}{b}y^2\right)^{2/3}.$ Then $$u$$ is a classical solution of (AC). Moreover, we conclude that ellipses remain ellipses of fixed eccentricity under the motion by affine curvature. 3. Finite difference equations and convergence 3.1 Elliptic finite difference schemes The framework developed in Barles & Souganidis (1991) provides conditions for when approximation schemes converge to the unique viscosity solution of a PDE. The article Oberman (2006) defines elliptic finite difference schemes, which are monotone and stable. Theorem 3.1 Consider an elliptic PDE that satisfies a comparison principle for viscosity solutions. A consistent, stable and monotone approximation scheme converges locally uniformly to the (unique) viscosity solution. Consider the domain $${\it {\Omega}} \subset D$$ where $$D$$ is an $$n$$-cube. We use a uniform grid of spacing, $$h$$, in $$D$$, and define the finite difference grid, $$G^h = \{x \in h \mathbb{Z}^n ~\mid ~ x\in D \}$$. Definition 3.2 (Elliptic finite difference equation) Let $$C(G^h)$$ denote the set of grid functions, $$u: G^h \to \mathbb{R}$$. A finite difference operator is a map $$F^h: C(G^h) \to C(G^h)$$, which has the following form, $F^h[u](x) = F^h(x, u(x), u(x) - u(\cdot)),$ where $$u(\cdot)$$ indicates the values of the grid function $$u$$. It has stencil width$$W$$ if $$F^h(x, u(x), u(x) - u(\cdot))$$ depends only on values $$u(y)$$ for $$\Vert{y-x}\Vert_\infty/h \le W$$. A solution of the finite difference scheme is a grid function that satisfies the equation $$F^h[u](x) = 0$$ for all $$x \in G^h$$. The finite difference operator is elliptic if $r \le s, ~ v(\cdot) \le w(\cdot) \implies {F^h}\left(x,r,v(\cdot)\right) \leq {F^h}(x,s,w(\cdot)).$ Definition 3.3 (Consistent) The scheme $${F^h}$$ is consistent with the equation $$F$$ if for any smooth function $$\phi$$ and $$x\in{{\it {\Omega}}}$$, $\lim_{h \to 0,y\to x} {F^h}[\phi](y) = F(D^2 \phi(x), \nabla \phi(x),\phi(x),x).$ The scheme is accurate to order $$k$$ if for any smooth function $$\phi$$ in a neighborhood of $$G^h$$ and $$x \in G^h \cap {\it {\Omega}}$$ ${F^h}[\phi](x) = F(D^2 \phi(x), \nabla \phi(x),\phi(x),x) + \mathcal{O}(h^k).$ Remark 3.4 Consistency and accuracy are usually verified by a Taylor series argument. Definition 3.5 (Discrete Comparison Principle) Given the finite difference operator $${F^h}: C(G^h) \to C(G^h)$$, the comparison principle holds for $${F^h}$$ if $$\label{comparison} {F^h}[u](x) \le {F^h}[v](x)~ \text{for all} ~x \implies u(x) \le v(x) ~\text{for all}~ x.$$ (Comp) Remark 3.6 In the discrete comparison principle, the boundary conditions are encoded in $${F^h}$$, the assumption $${F^h}[u] \leq {F^h}[v]$$ means $$u \le v$$ at Dirichlet boundary points. Uniqueness of solutions follows from the discrete comparison principle, because if $$u,v$$ are solutions, then $${F^h}(u) = {F^h}(v) = 0$$, so $$u \le v$$ and $$u\ge v$$, and thus $$u=v$$. Definition 3.7 (Lipschitz constant of the scheme) The finite difference operator $${F^h}: C(G^h) \to C(G^h)$$ is Lipschitz continuous with constant $$K^h$$ if $$K^h$$ is the smallest constant such that $\left\vert{{F^h}\left(x,r,v(\cdot)\right) - {F^h}\left(x,s,w(\cdot)\right)}\right\vert \le K^h \max(\vert{r-s}\vert, \Vert{v - w}\Vert_\infty), ~\text {for all}~ x \in G^h.$ Remark 3.8 In the definition above, the maximum on the right-hand side can be replaced with a maximum over the neighboring grid values without changing the Lipschitz constant. In Oberman (2006), it is shown that the Euler map $$u \to u - \rho F[u]$$ with constant $$\rho \le 1/K^h$$ is monotone and a (nonstrict) contraction. By adding an arbitrarily small multiple of $$u$$ to the scheme, the comparison principle holds, and the Euler map is a strict contraction. The resulting scheme provides a convergent discretization of the solution of the parabolic PDE $$u_t + F[u]=0$$, because the discretization is monotone and stable. Filtered schemes were first proposed in Froese & Oberman (2013) in the context of the Monge–Ampère equation to overcome the limited accuracy of the wide-stencil elliptic schemes. The idea is to blend a stable elliptic convergent scheme with an accurate scheme and retain the advantages of both: stability and convergence of the former and higher accuracy of the latter. The only requirement on the accurate scheme is consistency. Filtered schemes have also been used for Hamilton–Jacobi equations in Oberman & Salvador (2015). We require that the difference between the filtered scheme and the elliptic scheme is uniformly bounded. Under this assumption, the proof of the Barles–Souganidis theorem can be modified to obtain convergence results for filtered schemes (Froese & Oberman, 2013). 3.2 Example finite difference schemes Define the standard finite difference operators for $$u_x, u_y$$ and $$u_{xy}$$. Definition 3.9 (Standard finite differences for $$u_x$$, $$u_y$$ and $$u_{xy}$$) \label{standard_fd} \begin{aligned} & u_x^h := \frac{u(x+h,y)-u(x-h,y)}{2h} = u_x(x,y) + \mathcal{O}(h^2),\\ & u_{xx}^h := \frac{u(x+h,y)-2u(x,y)+u(x-h,y)}{h^2} = u_{xx}(x,y) + \mathcal{O}(h^2),\\ & u_{xy}^h := \frac{u(x+h,y+h)+u(x-h,y-h)-u(x+h,y-h)-u(x-h,y+h)}{4h^2} = u_{xy}(x,y) + \mathcal{O}(h^2), \end{aligned} (3.1) and, similarly for $$u_y^h$$, and $$u_{yy}^h$$ in the $$y$$ coordinate. The operator $$-u_{xx}^h$$ is elliptic, the others are not. Definition 3.10 Define the backward and forward first-order derivatives \begin{align*} D^-_x[u](x,y) &:= \frac{u(x,y)-u(x-h,y)}{h} = u_x(x,y) + \mathcal{O}(h), \\ -D^+_x[u](x,y) &:= \frac{u(x,y)-u(x+h,y)}{h} = -u_x(x,y) + \mathcal{O}(h), \end{align*} and, similarly for $$D^-_y[u]$$ and $$D^-_y[u]$$. The operators $$D^-_x$$ and $$-D^+_x$$ are elliptic. 3.3 Nondecreasing functions and elliptic discretizations To build elliptic difference schemes, we study the properties of nondecreasing maps. This is required, because, in general, elliptic schemes are built composing nondecreasing maps with elliptic terms. Moreover, for some nonlinear elliptic PDEs, the domain of ellipticity is restricted, and thus we need to build nondecreasing extensions of the underlying functions. In this section, we define elliptic schemes for $$\left\vert{\nabla u}\right\vert$$ and $$-\left\vert{\nabla u}\right\vert$$. This is done both as an illustration of the general principle and because these schemes will be needed for our discretization of the two-dimensional affine curvature operator. Definition 3.11 (Nondecreasing functions) For $$x,y \in \mathbb{R}^N$$ we say $$x \leq y$$ if $$x_i \leq y_i$$ for all $$i =1,\dots, N$$. The function $$F:\mathbb{R}^N \to \mathbb{R}$$ is nondecreasing, $$F \in ND(\mathbb{R}^N)$$, if $x \leq y \implies F(x) \leq F(y).$ Write $$\mathbb{R}^+ = \{x \in \mathbb{R} \mid x \geq 0 \}$$ and $$\mathbb{R}^- = \{x \in \mathbb{R} \mid x \leq 0 \}$$. Furthermore, if $$F \in ND(\mathbb{R}^N)$$ and $$F: \mathbb{R}^N \to \mathbb{R}^+$$ (resp. $$F: \mathbb{R}^N \to \mathbb{R}^-)$$ we write $$F \in ND^+(\mathbb{R}^N)$$ (resp. $$F \in ND^-(\mathbb{R}^N)$$). Remark 3.12 When $$f \in C^1(\mathbb{R}^N)$$, if $$f$$ is nondecreasing in each variable, i.e., $$f_{x_i} \geq 0$$ for all $$i =1,\dots, N$$, then $$f \in ND(\mathbb{R}^N)$$. Example 3.13 The function $$x^+ = \max(x,0) \in ND^+(\mathbb{R})$$ and $$x^- = \min(x,0)$$ is in $$ND^-(\mathbb{R})$$. Write $$N(x,y) = \sqrt{x^2+y^2}$$, and define $$N^+(x,y) := N(x^+,y^+)$$ and $$N^-(x,y) := -N(x^-,y^-)$$. Then $$N^+ \in ND^+(\mathbb{R}^2)$$, and $$N^- \in ND^-(\mathbb{R}^2)$$. Furthermore, $$N^+ = N$$ on $$\{x, y \geq 0 \}$$, $$N^- = -N$$ on $$\{x, y \leq 0 \}$$. Example 3.14 (Upwinding) More generally, if we consider the operator $$a(x)u_x$$, then the upwind discretization $a^+ D^-_x[u] + a^-D^+_x[u]$ is first-order accurate and elliptic. The first example of using elliptic discretizations as building blocks is given by the following. Example 3.15 Define $\left\vert{u_x^h}\right\vert^+ = \max\left\{-D^+_xu,D^-_xu,0\right\}, \quad -\left\vert{u_x^h}\right\vert^- = \min\left\{-D^+_xu,D^-_xu,0\right\}$ that approximate $$\vert{u_x}\vert$$ and $$-\vert{u_x}\vert$$to first order. The operators $$\left\vert{u_x^h}\right\vert^+$$ and $$-\left\vert{u_x^h}\right\vert^-$$ are elliptic, the former is non-negative and the latter is nonpositive. Example 3.16 Define $\vert{\nabla u^h}\vert^+ = N \left(|u_x^h|^+,|u_y^h|^+\right), \qquad -\vert{\nabla u^h}\vert^- = -N\left(-|u_x^h|^-,-|u_y^h|^-\right)$ which are elliptic, consistent with $$\left\vert{\nabla u}\right\vert$$, and $$-\left\vert{\nabla u}\right\vert$$, respectively, and first-order accurate. 4. Numerical methods for the model equation In this section, we build convergent discretizations to the one-dimensional model of our equation defined by $$\label{1Dmodel} F^{1D}[u] := A(u_x,u_{xx}) = \left(u_x^2 \:u_{xx}\right)^{1/3} = f, \quad x \in (-1,1),$$ (AC-1D) as well as the parabolic PDE, $u_t = F^{1D}[u]-f, \quad \text{for} (x,t) \in (-1,1) \times (0,\infty),$ along with initial and boundary conditions. To do so, we need to build elliptic discretizations for $$F^{1D}$$, which is achieved by studying the nonlinear function $$A(p,q)$$. A naive approach would suggest to simply substitute the elliptic discretizations for $$\left\vert{u_x}\right\vert$$ and $$u_{xx}$$ into $$A(p,q)$$. However, this is not possible since $$A(p,q) \not\in ND(\mathbb{R}^2)$$: $$dA/dp = 2/3(q/p)^{1/3}$$ and so $$A$$ fails to be nondecreasing when $$pq < 0$$. Our approach is the following. First, in Section 4.1, we write $$A(p,q)$$ as a sum of two nondecreasing functions in terms of $$\left\vert{p}\right\vert$$, $$-\left\vert{p}\right\vert$$, $$q$$. These terms are replaced by $$\vert{u^h_x}\vert$$, $$-\vert{u^h_x}\vert$$, $$u_{xx}^h$$, which are elliptic and consistent, and thus a consistent and elliptic discretization for $$F^{1D}$$ is build. Then, in Section 4.2, we present a Lipschitz regularization of the function $$A(p,q)$$ and proceed similarly. The Lipschitz regularization is needed to build a provably convergent explicit scheme for the parabolic PDE as discussed in Section 3.1. Finally, in Section 4.3, we present the convergence results. To save space, we omit some of the proofs. They can be found in Oberman & Salvador (2016). 4.1 An elliptic discretization of the one-dimensional operator In the next lemma, we decompose $$A(p,q)$$ into the sum of two nondecreasing functions. Lemma 4.1 Define $$A^+(p,q) \,{=}\, A(p^+,q^+)$$ and $$A^-(p,q) \,{=}\, A(p^-,q^-)$$. Then $$A^+ \,{\in}\, ND^+(\mathbb{R}^2)$$, $$A^- \,{\in}\, ND^-(\mathbb{R}^2)$$ and $-A(p,q) = A(p,-q) = A^+(\vert{p}\vert,-q) + A^-(-\vert{p}\vert,-q), \quad \text{for all}~ p,q.$ We can now build an elliptic scheme for $$-F^{1D}[u]$$. Lemma 4.2 Define the finite difference operator \begin{align} -F^{1D,e}[u] = A^+\left(\left\vert{u_x^h}\right\vert^+, -u_{xx}^h \right) + A^-\left(-\left\vert{u_x^h}\right\vert^-,-u_{xx}^h \right), \end{align} (AC)1D,e where the finite difference operators involved are defined above. Then $$-F^{1D,e}$$ is elliptic and consistent with the $$-F^{1D}$$. Proof. The operators $$|u_x^h|^+$$, $$-|u_x^h|^-$$, $$-u_{xx}^h$$ are elliptic. Then, due to Lemma 4.1, the operator $$-F^{1D,e}[u]$$ consists of the composition of the nondecreasing functions $$A^+$$ and $$A^-$$ with the three preceding operators, with the first two taking the place of $$\vert{p}\vert$$ and $$-\vert{p}\vert$$ and the last one taking the place of $$-q$$. Because the operators are elliptic and the functions are nondecreasing, $$-F^{1D,e}[u]$$ is elliptic. Consistency follows from the consistency of each of the schemes used. □ Remark 4.3 (Accuracy) One challenge in forming a discretization of (AC-1D) is that the accuracy breaks down near $$u_{xx} = 0$$. For example, if we use second-order accurate approximations of $$u_{xx}$$, the accuracy of finite differences for the operator is only $$\mathcal{O}(h^{2/3})$$. Consider the expression $A(p+h^k,q+h^2) = ((p + h^k)^2(q+ h^2))^{1/3},$ where $$k = 1$$ or $$2$$. If $$p\not=0, q=0$$ we get $A(p+h^k,q+h^2) = h^{2/3} (p + h^k)^{2/3} = \mathcal{O}(h^{2/3}).$ So regardless of the accuracy of the discretization of the $$u_x$$ term, the overall accuracy of the scheme cannot be better than $$\mathcal{O}(h^{2/3})$$. In fact, a similar argument shows that the order of accuracy is $$2k/3$$ near $$p=0, q\not=0$$. Our elliptic discretization, which uses $$k=1$$, will have an order of accuracy $$2/3$$. 4.2 Lipschitz regularization of the one-dimensional operator The function $$A(p,q)$$ fails to be Lipschitz continuous near the axis $$p=0$$, $$q=0$$. Hence, the elliptic scheme presented in the previous section is not Lipschitz, a property that is required to build a monotone convergent scheme for the time-dependent problem. Thus, using the notation for the sign function $$\text{sgn}(q) = q/\vert{q}\vert$$ for $$q\not=0$$ and $$\text{sgn}(0) = 0$$ otherwise, we regularize $$A(p,q)$$ as follows. Definition 4.4 Define for $$K= K(\delta), L = L(\delta) > 0$$, the regularized function $$\label{Adelta} A^\delta(p,q) = \text{sgn}(q)\min(\vert{A(p,q)}\vert,K\vert{p}\vert,L\vert{q}\vert).$$ (4.1) Naturally, we can then define the regularized PDE operator as $F^{1D,\delta}[u] = A^\delta(u_x,u_{xx}).$ Defining an elliptic and consistent scheme for $$F^{1D,\delta}[u]$$ is accomplished by noticing that the discretization (AC)1D,e generalizes when we replace $$A$$ with the regularized version $$A^\delta$$ in each term: just like for $$A(p,q)$$ in Lemma 4.1, we can decompose $$A^\delta(p,q)$$ into the sum of two nondecreasing functions. This is achieved in the following lemma. Lemma 4.5 Define $$A^{\delta,+}(p,q) = A^\delta(p^+,q^+)$$ and $$A^{\delta,-}(p,q) = A^\delta(p^-,q^-)$$. Then $-A^\delta(p,q) = A^{\delta,+}(\vert{p}\vert,-q)+A^{\delta,-}(-\vert{p}\vert,-q),$ where $$A^{\delta,+} \in ND^+(\mathbb{R}^2)$$ and $$A^{\delta,-} \in ND^-(\mathbb{R}^2)$$. Ellipticity and consistency of the regularized scheme with respect to $$-F^{1D,e,\delta}[u]$$ follows now easily, just like in Lemma 4.2. For this reason, we omit the proof. Lemma 4.6 For $$K = K(\delta)$$, $$L=L(\delta) > 0$$, define the finite difference scheme $$\label{elliptic_scheme_regular1D} -F^{1D,e,\delta}[u] = A^{\delta,+}\left(\left\vert{u_x^h}\right\vert^+, -u_{xx}^h \right) + A^{\delta,-}\left(-\left\vert{u_x^h}\right\vert^-, -u_{xx}^h \right).$$ (AC)1,D,e,δ Then $$-F^{1D,e,\delta}[u]$$ is elliptic and consistent with $$-F^{1D,\delta}[u]$$. Proving consistency of the regularized scheme with respect to $$-F^{1D}[u]$$ requires extra work as the parameters $$K$$, $$L$$ need to be chosen carefully. The next theorem accomplishes precisely that. Theorem 4.7 Asssume $$K = h^{-1/3}$$ and $$L = h^{-4/3}$$. Let $$x \in {\it {\Omega}}$$ be a reference point on the grid and $$\phi$$ be a smooth function that is defined in a neighborhood of the grid. Then the scheme $$F^{1D,e,\delta}$$ defined by (AC)1D,e,δ is consistent with $$F^{1D}$$ and has accuracy $F^{1D,e,\delta}[\phi](x) - F^{1D}[\phi](x) = \mathcal{O}(h^{2/3}).$ Moreover, $$F^{1D,e,\delta}$$ is Lipschitz continuous with constant $$C^h$$ given by $$\label{Lipconst1D} C^h = \frac{K}{h}+\frac{2L}{h^2} = h^{-4/3} + 2h^{-10/3}.$$ (4.2) Remark 4.8 As a result of the proof, we show as well that $$F^{1D,e,\delta}$$ is consistent with $$F^{1D,\delta}$$ with accuracy $F^{1D,e,\delta}[\phi](x) - F^{1D,\delta}[\phi](x) = \mathcal{O}\left(h^{1-\alpha} + h^{2-\beta}\right),$ if $$K = \mathcal{O}(h^{-\alpha})$$ and $$L = \mathcal{O}(h^{-\beta})$$. The optimal choice is $$\alpha = 1/3$$ and $$\beta = 4/3$$, in which case the overall accuracy of $$F^{1D,e,\delta}$$ and $$F^{1D,e}$$ with respect to $$F^{1D}$$ is the same (see Remark 4.3), meaning that no accuracy is lost due to the regularization. 4.3 Convergence theorems for the one-dimensional model Having proved the ellipticity and consistency of the schemes, the uniform convergence follows as discussed in Section 3.1. The first convergence result is for the elliptic problem, where there is no need for the regularized scheme. Theorem 4.9 Let $$u(x)$$ be the unique viscosity solution of $$F^{1D}[u] = f$$ in $${\it {\Omega}}$$, along with suitable boundary conditions. For each $$h>0$$, let $$u^{1D,e,h}$$ be the uniformly bounded solution of $$F^{1D,e}[u] = f$$. Then $$u^{1D,e,h} \to u$$ locally uniformly, as $$h \to 0$$. Proof. Convergence for the elliptic discretization $$-F^{1D,e}$$ follows from the Barles–Souganidis theorem. □ Unlike in the elliptic problem, in the parabolic problem we need to use the regularized scheme, with the time discretization being given by a forward Euler step. This leads us to $$\label{discretizationtime1D} u(\cdot,t+{\rm d}t) = u(\cdot,t) + {\rm d}t \:F^{1D,e,\delta}[u(\cdot,t)].$$ (4.3) Theorem 4.10 Let $$u(x,t)$$ be the unique viscosity solution of $$u_t = {F}^{1D}[u]$$ in $${\it {\Omega}} \times [0,\infty)$$, along with $$u(x,0) = u_0(x)$$ and suitable boundary conditions. Assume as well that $$K = h^{-1/3}$$ and $$L = h^{-4/3}$$. For each $$h >0$$, let $$u^{1D,e,h}$$ be the uniformly bounded solution of the monotone time discretization (4.3) with $${\rm d}t \leq 1/{C^h}$$ given by (4.2). Then $$u^{1D,e,h} \to u$$ locally uniformly, as $${\rm d}t,h \to 0$$. Proof. The elliptic scheme $$F^{1D,e,\delta}$$ leads to a monotone time discretization, in other words, the solution map satisfies a comparison principle if $${\rm d}t \le 1/{C^h}$$, where $$C^h$$ is the Lipschitz constant of the elliptic scheme by Oberman (2006). Then the Barles–Souganidis theorem (Barles & Souganidis, 1991) applies. □ Remark 4.11 It is also true that, for fixed values of $$K,L$$, there is a unique viscosity solution, $$u^\delta$$, of the regularized PDE. Then, fixing $$K,L$$ and using the discretization above with $${\rm d}t = \mathcal{O}(h^2)$$, the forward Euler method converges uniformly to $$u^\delta$$ as $$h\to 0$$. 5. Nonconvergence of standard finite differences In this section, we show that standard finite differences are unstable for both the one-dimensional and the two-dimensional operators that we study. This instability can be understood from the one-dimensional model, which, if we take $$|u_x| = 1$$, reduces to the operator $$u_{xx}^{1/3}$$. It is certainly plausible that the singularity near $$u_{xx} =0$$ could lead to instabilities. This motivates the regularization introduced in the previous section, which replaces the singularity of the cube root with a linear function with large slope. It also motives the convergent finite difference schemes, which have an explicit time step that ensures the convergence of the time-dependent schemes. We begin with an example where the standard finite scheme does not converge for the one-dimensional model. Next, we consider the time-dependent equation and the associated scheme obtained with the (unregularized) elliptic scheme and a forward Euler step in time. A scaling argument suggest that the choice $${\rm d}t = \mathcal{O}(h^{4/3})$$ should provide a stable scheme. In fact, this scaling argument can be improved to a proof of the maximum principle for the homogeneous equation with the same time step. However, the maximum principle is not enough for convergence (we need the comparison principle) and we demonstrate divergence with that time step. Using a smaller time step $${\rm d}t = \mathcal{O}(h^2)$$ appears to fix the problem. (The standard finite difference scheme diverges for the example we present no matter how small the time step.) The example is then generalized to the two-dimensional operator. 5.1 Nonconvergence of standard finite differences in one dimension Consider the discretization given by inserting the standard centered differences, given by (3.1), $F^{1D,a}[u] = \left(\left(u_x^h\right)^2\left(u_{xx}^h\right)\right)^{1/3}.$ We considered an exact solution $$u_0(x) = \sin(2\pi x)$$ of (AC-1D). Then we solved the time-dependent problem, using the forward Euler time discretization, with initial data given by the solution $$u(x,0) = u_0(x)$$. We found that this solution was unstable for the standard finite difference scheme. The results are displayed in Fig. 1, which illustrates that the numerical solution diverges from the exact solution. This effect persists over different grid sizes and over smaller time steps. Fig. 1. View largeDownload slide Solution of the one-dimensional model equation using standard finite differences at times $$t \in \{0,1,2,5\}$$. Set $${\rm d}t = h^2/2$$ on a $$256$$-point grid. Fig. 1. View largeDownload slide Solution of the one-dimensional model equation using standard finite differences at times $$t \in \{0,1,2,5\}$$. Set $${\rm d}t = h^2/2$$ on a $$256$$-point grid. Next, we consider the discretization of the time-dependent equation $$\label{aff1dt} u_t = F^{1D}[u] - f.$$ (5.1) Using a forward Euler method in time and the elliptic method in space leads to $$\label{aff1dth} u(\cdot,t+dt) = {\it {\Phi}}_t^{1D,e}(u) := u(\cdot,t) + {\rm d}t (F^{1D,e}[u(\cdot,t)]-f).$$ (5.2) A scaling argument suggests $${\rm d}t = \mathcal{O}(h^{\frac{4}{3}})$$ might be stable, because the operator $$F^{1D}$$ is homogeneous to this order in $$h$$. In fact, we are able to prove the following. Lemma 5.1 When $$f=0$$ in (5.1), the solution map $${\it {\Phi}}^{1D,e}$$ satisfies the maximum principle $\min {\it {\Phi}}^{1D,e}_t(u) \leq {\it {\Phi}}^{1D,e}_{t+{\rm d}t}(u) \leq \max {\it {\Phi}}^{1D,e}_t(u),$ provided $${\rm d}t \leq (h^4/2)^{1/3}$$. The proof can be found in Oberman & Salvador (2016). However, as the following example shows, the maximum principle by itself is not enough to guarantee convergence and so the above choice for the time step is not guaranteed to produce a convergent scheme. In practice, the above choice for the time step leads to a divergent, but bounded scheme with nonlinear instabilities. Example 5.2 We solved (5.1) using (5.2). We took $$u(x,0) = u_0$$ where $$u_0(x) = C x^{4/3}$$ and $$f =F^{1D}[u_0]$$. In Fig. 2, we plot the numerical solution obtained at different times with $${\rm d}t = (h^4/4)^{1/3}$$ (The conservative choice of the time step is to account for the fact the equation is not homogeneous.) The exact solution is clearly unstable. For larger times, the solution has the same shape, but with small high-frequency oscillations in time. On the other hand, choosing $${\rm d}t = h^2/2$$ leads to convergence (data are not presented to save space). Fig. 2. View largeDownload slide Plot of the solution obtained described in Example 5.2 at time $$t \in \{0,1,5,20\}$$ on a $$128$$-point grid. Fig. 2. View largeDownload slide Plot of the solution obtained described in Example 5.2 at time $$t \in \{0,1,5,20\}$$ on a $$128$$-point grid. 5.2 Two dimensions In this section, we define and study the standard finite difference scheme for $$F[u]$$. We show that the accuracy degenerates near points where the affine curvature is zero. We give an example of a steady solution where the standard finite difference scheme breaks down. Using (1.1), we can write \begin{align*} F[u] & = \left(u_{xx} u_y^2 - 2u_x u_y u_{xy} + u_{yy} u_x^2\right)^{1/3}. \end{align*} Definition 5.3 Let $$u:\mathbb{R}^2\to\mathbb{R}$$. We define the scheme \begin{align} F^a[u] = \left(u_{xx}^h (u_y^h)^2 - 2u_x^h u_y^h u_{xy}^h + u_{yy}^h (u_x^h)^2\right)^{1/3} \end{align} (AC)a that approximates $$F[u]$$. Remark 5.4 The $$u_{xy}^h$$ term is not elliptic, and consequently $$-F^a$$ is not elliptic. Lemma 5.5 Let $$(x,y) \in {\it {\Omega}}$$ be a reference point on the grid and $$\phi$$ be a smooth function that is defined in a neighborhood of the grid. Then the scheme $$F^a[\phi]$$ defined by (AC)a approximates $$F[\phi]$$ with accuracy $F^a[\phi](x,y) - F[\phi](x,y) = \begin{cases} \mathcal{O}(h^2), & F[\phi](x,y) \neq 0,\\ \mathcal{O}(h^{\frac{2}{3}}), &F[\phi](x,y) = 0. \end{cases}$ Proof. All the standard finite differences used are second-order accurate. Therefore, $\phi_{xx}^h \left(\phi_y^h\right)^2 - 2\phi_x^h \phi_y^h \phi_{xy}^h + \phi_{yy}^h \left(\phi_x^h\right)^2 = \phi_{xx} \phi_y^2 - 2\phi_x \phi_y \phi_{xy} + \phi_{yy} \phi_x^2 + \mathcal{O}(h^2),$ in other words, $$(F^a[\phi])^3 = (F[\phi])^3 + \mathcal{O}(h^2)$$. In addition, the Taylor expansion tells us that $(r+t)^{1/3} = r^{1/3} + \frac{t}{3r^{2/3}} + \mathcal{O}(t^2)$ when $$r \neq 0$$. Hence, when $$F[\phi](x,y) \neq 0$$, it follows that $$F^a[\phi] = F[\phi]+\mathcal{O}(h^2)$$. Otherwise, when $$F[\phi](x,y) = 0$$, we can only conclude that $$F^a[\phi] = F[\phi]+\mathcal{O}(h^{\frac{2}{3}})$$. □ Now we present an example that shows that the scheme $$F^a$$ does not converge. We choose a level set function and a right-hand side so that the equation is a steady state (see Example 7.9(d), for more details). Starting from initial data corresponding to the exact solution, and evolving in time with a forward Euler step, the solution changes to order one. Indeed, the solution does not appear to reach a steady state, even after running for a long time. See Fig. 3 for snapshots in time of the solution. Similar behavior was observed on finer grids (although it took a longer time to reach a comparable change in the solution). Fig. 3. View largeDownload slide Example 7.9 $$(d)$$ with the standard finite difference solver: level sets of the solution at times $$t = 0$$ (upper left), $$t = 15$$ (upper center), $$t = 17$$ (upper right), $$t = 20$$ (lower left), $$t = 40$$ (lower center) and $$t = 50$$ (lower right) with $${\rm d}t = h^2/2$$ on a $$32 \times 32$$ grid. Fig. 3. View largeDownload slide Example 7.9 $$(d)$$ with the standard finite difference solver: level sets of the solution at times $$t = 0$$ (upper left), $$t = 15$$ (upper center), $$t = 17$$ (upper right), $$t = 20$$ (lower left), $$t = 40$$ (lower center) and $$t = 50$$ (lower right) with $${\rm d}t = h^2/2$$ on a $$32 \times 32$$ grid. 6. Convergent finite difference methods Based on the convergence theory discussed above, and on the elliptic discretization of the one-dimensional model equation, we now build a discretization of the affine curvature operator. This discretization is based on the median scheme for the mean curvature operator from Oberman (2004). We could also use the morphological operator (Cattè et al., 1995), which results in a very similar discretization of the operator $$F[u]$$. We establish the accuracy of the discretization and show that it is elliptic. The scheme is augmented to a more accurate, but is still a convergent, filtered scheme, which interpolates between the standard scheme and the elliptic scheme. We regularize the operator, which allows us to build a convergent monotone time discretization. 6.1 Median for $${\it {\Delta}}_1 u$$ In Oberman (2004), an elliptic scheme for $${\it {\Delta}}_1 u(x)$$ is presented based on taking the median of the values of $$u$$ sampled in a small approximately circular neighborhood of $$x$$. The motivation follows from observing, using (1.1), that ${\it {\Delta}}_1 u = u_{\mathbf{t}\mathbf{t}}, \qquad \mathbf{t} = \frac{(-u_y,u_x)}{(u_x^2+u_y^2)^{1/2}},$ where $$\mathbf{t}$$ is the (Euclidean) unit tangent. The median captures an approximation to $$u_{\mathbf{t} \mathbf{t}}$$, the second tangential derivative of $$u$$, because the larger values point in the direction of the gradient and the smaller values point in the opposite directions. We outline the scheme here. We will define it at the reference point $$(x,y)$$. Let $$(x_1,y_1),\ldots,(x_{n_S},y_{n_S})$$ be the $$n_S$$ neighbors and set $${\rm d}\theta = \frac{2\pi}{n_S}$$. We refer to $${\rm d}\theta$$ as the directional resolution. Denote the value of the solution at the point $$(x_i,y_i)$$ by $$u_i$$. Set $$\mathbf{v_i} = (x_i,y_i)-(x,y)$$ and $$(x_{-i},y_{-i}) = (x,y)-\mathbf{v_i}$$. We choose the neighbors such that $(x_{i+1},y_{i+1}) = h \: n_\theta\left(\cos(i{\rm d}\theta),\sin(i {\rm d}\theta)\right) + (e_i,f_i),$ where $$|e_i|,|\,f_i| \leq h$$. Thus, $$h$$ is the spatial resolution and $$n_\theta$$ denotes the width of the stencil. In fact, for $$n_\theta \leq 5$$, the neighbors in the first quadrant are given by $(x_{i+1},y_{i+1}) = h\left(\lfloor n_\theta \cos(i{\rm d}\theta) \rceil, \lfloor n_\theta \sin(i{\rm d}\theta) \rceil\right)\!,$ where $$i = 0,\ldots,\frac{n_S}{4}-1$$, with the points on the remaining three quadrants being obtained by $$\frac{\pi}{2}$$, $$\pi$$ and $$\frac{3\pi}{2}$$ rotations. Here $$\lfloor x \rceil$$ denotes the integer closest to $$x$$. See Fig. 4 for the stencils with $$n_\theta = 3$$ and $$n_\theta = 7$$. Fig. 4. View largeDownload slide Neighbor points of the stencil for $$n_\theta = 3$$ (smaller circle) and $$n_\theta = 7$$ (larger circle). Fig. 4. View largeDownload slide Neighbor points of the stencil for $$n_\theta = 3$$ (smaller circle) and $$n_\theta = 7$$ (larger circle). Definition 6.1 Let $$u:\mathbb{R}^2 \to \mathbb{R}$$. Define the scheme $${\it{\Delta}} _{1}^{e}u(x,y)=2\frac{\underset{i=1,\ldots, {{n}_{S}}}{\mathop{\text{median}}}\,{{u}_{i}}-u(x,y)}{{{(h\,{{n}_{\theta}})}^{2}}}$$ (MC)e that approximates $${\it {\Delta}}_1 u$$. In general, consistency of finite difference schemes follows by Taylor expansions. However, additional care is needed when the PDE is singular, as when $$\nabla u(x) = 0$$ in (MC). We then recall here the definition of consistency for the mean curvature operator. Definition 6.2 The numerical scheme $$F^{h,d\theta}$$ is consistent with $${\it {\Delta}}_1$$ if for every smooth function $$\phi$$ and every $$(x,y) \in \mathbb{R}^2$$ $\lim_{h,{\rm d}\theta \to 0} F^{h,{\rm d}\theta}[\phi] = {\it {\Delta}}_1\phi$ at $$(x,t)$$ if $$\nabla \phi(x,y) \neq 0$$ and $\lambda \leq \liminf_{h,{\rm d}\theta \to 0} F^{h,{\rm d}\theta}[\phi] \leq \limsup_{h,{\rm d}\theta \to 0} F^{h,{\rm d}\theta}[\phi] \leq {\it {\Lambda}}$ at $$(x,t)$$ where $$\lambda,{\it {\Lambda}}$$ are the least and greatest eigenvalues of $$D^2\phi(x,y)$$, otherwise. Lemma 6.3 The finite difference scheme $$-{\it {\Delta}}_1^e u$$, given by (MC)e, is elliptic. Proof. We can rewrite the scheme as $-{\it{\Delta}}_{1}^{e}u=2\frac{\underset{i=1,\ldots, {{n}_{S}}}{\mathop{\text{median}}}\,(u(x,y)-{{u}_{i}})}{{{(h\,{{n}_{\theta}})}^{2}}}.$ Since $$\text{median}$$ is a nondecreasing function, we can immediately conclude that $${\it {\Delta}}_1^e u$$ is elliptic. □ Lemma 6.4 Let $$(x,y) \in {\it {\Omega}}$$ be a reference point on the grid and $$\phi$$ be a smooth function that is defined in a neighborhood of the grid. Then the scheme $${\it {\Delta}}_1^e \phi$$ defined by (MC)e is consistent. Further, it approximates $${\it {\Delta}}_1$$ with accuracy ${\it {\Delta}}_1^e \phi(x,y) = {\it {\Delta}}_1 \phi(x,y) + \mathcal{O}\left((h \: n_\theta)^2 + {\rm d}\theta\right)\!,$ when $$\left\vert{\nabla u(x,y)}\right\vert \neq 0$$. Remark 6.5 Since $${\rm d}\theta = \mathcal{O}\left(\frac{1}{n_\theta}\right)$$, the ‘optimal’ choice is $$h = \mathcal{O}\left(n_\theta^{-3/2}\right)$$, which results in accuracy of $$\mathcal{O}(h^{2/3})$$. However, this also requires a dense stencil. In practice, we use a fairly narrow stencil and combine with a filtered scheme for more accuracy. The proof can be found in Oberman (2004). 6.2 Elliptic scheme for $$-F[u]$$ We now construct an elliptic scheme for $$-F[u]$$ following the same procedure as in Section 4.1. This is accomplished by writing $$F[u]$$ in terms of $$|\nabla u|$$ and $${\it {\Delta}}_1 u$$ as \begin{align*} F[u] & = \left\vert{\nabla u}\right\vert k[u]^{1/3} = (\left\vert{\nabla u}\right\vert^2 {\it {\Delta}}_1 u)^{1/3} = A(\left\vert{\nabla u}\right\vert,{\it {\Delta}}_1 u) \end{align*} and using the elliptic schemes $$\vert{\nabla u^h}\vert^+$$ and $$\vert{\nabla u^h}\vert^-$$ presented in Section 3.3 and $${\it {\Delta}}_1^e u$$ described in Section 6.1. Remark 6.6 We chose to discretize $${\it {\Delta}}_1u$$ with the median scheme (MC)e. However, other schemes could have been, for instance, the morphological scheme (Cattè et al., 1995). Definition 6.7 Let $$u:\mathbb{R}^2\to\mathbb{R}$$. We define the scheme \begin{align} -F^e[u] = A^+\left(\left\vert{\nabla u^h}\right\vert^+, -{\it {\Delta}}_1^e u\right) + A^-\left(-\left\vert{\nabla u^h}\right\vert^-, -{\it {\Delta}}_1^e u \right) \end{align} (AC)e that approximates $$-F[u]$$. Remark 6.8 The above approach can be generalized to obtain an elliptic scheme for the elliptic operator $$-\left\vert{\nabla u}\right\vert G(k[u])$$, where $$G$$ is nondecreasing and homogeneous of order $$\alpha \le 1$$, $$G(t r) = t^\alpha G(r)$$. (This PDE represents curves evolving with normal speed $$G(k)$$.) In such cases, write \begin{align*} \left\vert{\nabla u}\right\vert G(k[u]) = \left\vert{\nabla u}\right\vert^{1-\alpha} G(\left\vert{\nabla u}\right\vert k[u]) = \left\vert{\nabla u}\right\vert^{1-\alpha} G({\it {\Delta}}_1 u). \end{align*} The elliptic scheme is then given by $\left(|\nabla u^h|^+\right)^{1-\alpha} G\left(\left(-{\it {\Delta}}_1^e u\right)^+\right) +\left(|\nabla u^h|^-\right)^{1-\alpha} G\left(\left(-{\it {\Delta}}_1^e u\right)^-\right).$ Now we prove the ellipticity and consistency of the scheme $$-F^e[u]$$. Lemma 6.9 The finite difference scheme $$-F^e[u]$$, given by (AC)e, is elliptic. Proof. The proof is similar to Lemma 4.6. □ Lemma 6.10 Let $$(x,y) \in {\it {\Omega}}$$ be a reference point on the grid and $$\phi$$ be a smooth function that is defined in a neighborhood of the grid. Then the scheme $$F^e[\phi]$$ defined by (AC)e is consistent with $$F[\phi]$$ and has accuracy $F^e[\phi](x,y) - F[\phi](x,y) = \begin{cases} \mathcal{O}\left(h + (h \: n_\theta)^2+ h \: {\rm d}\theta +{\rm d}\theta\right)\!, & \text{if}~F[\phi](x,y) \neq 0,\\ \mathcal{O}\left(((h \: n_\theta)^2 + h \: {\rm d}\theta + {\rm d}\theta)^{1/3}\right)\!, & \text{if}~F[\phi](x,y) = 0 \text{and} \left\vert{\nabla \phi}\right\vert(x,y) \neq 0,\\ \mathcal{O}\left(h^{2/3}\right) & \text{if}~\left\vert{\nabla \phi}\right\vert(x,y) = 0. \end{cases}$ Proof. Suppose first that $$\left\vert{\nabla \phi}\right\vert(x,y) \neq 0$$. We have $\left\vert{\nabla \phi(x,y)}\right\vert^2 = \left(|\nabla \phi^h(x,y)|^\pm\right)^2 + \mathcal{O}(h) \quad \text{and} \quad {\it {\Delta}}_1^e \phi(x,y) = {\it {\Delta}}_1 \phi(x,y) + \mathcal{O}\left({\rm d}\theta + (h \: n_\theta)^2\right)\!,$ (see Lemma 6.4). Hence, when $$F[\phi](x,y) \neq 0$$, $F^e[\phi]^3 = F[\phi]^3 + \mathcal{O}\left(h + (h \: n_\theta)^2 + h \: {\rm d}\theta +{\rm d}\theta\right)$ and, by the Taylor expansion like in Lemma 5.5, we have $$F^e[\phi](x,y) = F[\phi](x,y) + \mathcal{O}(h + (h \: n_\theta)^2 + h \: {\rm d}\theta +{\rm d}\theta)$$. Otherwise, when $$F[\phi](x,y) = 0$$, we necessarily have $${\it {\Delta}}_1 \phi = 0$$ and so $F^e[\phi]^3 = F[\phi]^3 + \mathcal{O}\left((h \: n_\theta)^2 + h \: {\rm d}\theta +{\rm d}\theta\right)\!.$ Therefore, $$F^e[\phi](x,y) = F[\phi](x,y) + \mathcal{O}\left(((h \: n_\theta)^2 + h \: {\rm d}\theta + {\rm d}\theta)^{1/3}\right)$$. When $$\left\vert{\nabla \phi}\right\vert(x,y) = 0$$, $\left\vert{\nabla \phi(x,y)}\right\vert^2 = \left(|\nabla \phi^h(x,y)|^\pm\right)^2 + \mathcal{O}(h^2)$ with $${\it {\Delta}}_1^e \phi$$ being bounded. Hence, $$F^e[\phi] = F[\phi] + \mathcal{O}(h^{2/3})$$. □ Remark 6.11 Since $${\rm d}\theta = \mathcal{O}\left(\frac{1}{n_\theta}\right)$$, the ‘optimal’ choice is $$h = \mathcal{O}\left(n_\theta^{-3/2}\right)$$, which results in accuracy of $$\mathcal{O}(h^{2/9})$$ in the worst case. However, this also requires a dense stencil. In practice, we use a fairly narrow stencil and combine with a filtered scheme for more accuracy. 6.3 Filtered scheme for $$F[u]$$ The filtered scheme, $$F^f[u]$$, first introduced in Froese & Oberman (2013), is defined to be a continuous interpolation between the accurate and the elliptic scheme, which equals the accurate scheme when the two schemes are within $$\epsilon$$ of each other. Here, we define filtered schemes by making use of the auxiliary function $$S^\epsilon:\mathbb{R}^2 \to \mathbb{R}$$, defined to be a continuous function which for $$(a,b) \in R^2$$, is equal to $$a$$ near the diagonal and $$b$$ off the diagonal. Definition 6.12 Define for $$\epsilon \,{>}\, 0$$, $$A_\epsilon \,{=}\, \left\{(a,b) \,{\in}\, \mathbb{R}^2 \mid |a-b| \,{<}\, \epsilon\right\}$$. Set $$\rho \,{=}\, 10\epsilon$$. Define $$d \,{=}\, {\text{dist}}\left((a,b),A_\epsilon\right)$$. $S^\epsilon(a,b) = \begin{cases} a & \text{if}~ (a,b) \in A_\epsilon,\\ \frac{\rho-d}{\rho}a+\frac{d}{\rho}b & \text{if}~ d \leq \rho,\\ b & \text{otherwise}. \end{cases}$ Define the filtered scheme $$\label{filterAC} F^f[u] = S^\epsilon\left(F^a[u],F^e[u]\right),$$ (AC)f where $$\epsilon = \epsilon(h,{\rm d}\theta)$$. Although, theoretically, the only requirement on $$\epsilon$$ to ensure the convergence of the filtered schemes is that $$\epsilon \to 0$$ as $$h,{\rm d}\theta\to0$$, in practice $$\epsilon$$ must be chosen carefully. Intuitively, it should be large enough to permit the accurate scheme to be active where the solution is smooth (as shown in the next lemma) and small enough to force the monotone scheme to be active whenever needed for convergence (for instance, when the solution is singular). Lemma 6.13 Suppose that the formal discretization errors of the schemes $$F^a$$ and $$F^e$$ are $$\mathcal{O}(r^{\beta_a})$$ and $$\mathcal{O}(r^{\beta_e})$$, respectively. Choose $$\alpha$$ such that $$\beta_a > \beta_e > \alpha > 0$$. Then if $$\phi$$ is smooth and $$\epsilon = r^\alpha$$, $$F^f[\phi] = F^a[\phi]$$. Proof. If $$\phi$$ is smooth then ${F^a[\phi]-F^e[\phi]} = {\mathcal{O}(r^{\beta_a})+\mathcal{O}(r^{\beta_e})} = \mathcal{O}(r^{\beta_e}) \leq \mathcal{O}(\epsilon),$ so using the definition of $$F^f$$ it follows that $$F^f[\phi] = F^a[\phi]$$, as desired. □ Remark 6.14 The lemma tells us that heuristically we could choose $$\epsilon = \sqrt{\text{Acc}[F^e]}$$, where $$\text{Acc}[F^e]$$ denotes the accuracy of $$F^e$$ (which corresponds to the choice $$\alpha = \beta_e/2$$ in the above lemma). In the numerical results presented here, we defined $$\epsilon$$ based on the accuracy of the scheme away from the singularities of $$F[u]$$. In practice, the filtered scheme allows us to neglect the error coming from the wide stencil, whereas in theory we still need to send $${\rm d}\theta \to 0$$ for convergence of the filtered scheme. In our numerical examples, we obtain the accuracy of the accurate scheme in most cases. 6.4 Regularization and Forward Euler method As in the for the one-dimensional model equation (see Section 4.2), to build a provably convergent scheme for the time dependent equation (AC) we need to regularize the operator. Write $$F[u] = A(\left\vert{\nabla u}\right\vert,{\it {\Delta}}_1 u)$$, where $$A(p,q) = (p^2q)^{1/3}$$, and regularize the cube root function as before. This leads to $F^\delta[u] = A^\delta(\left\vert{\nabla u}\right\vert,{\it {\Delta}}_1 u),$ where $$A^\delta(p,q) = \text{sgn}(q)\min(\vert{A(p,q)}\vert, K \vert{p}\vert, L \vert{q}\vert)$$. Remark 6.15 The regularized operator, $$F^\delta[u]$$, is still a level set operator. To see this, it is enough to take $$K = L = 1/\delta$$: $F^\delta[u] = \left\vert{\nabla u}\right\vert \text{sgn}(k[u])\min\left(\left\vert{k[u]}\right\vert^{1/3},\frac{\left\vert{k[u]}\right\vert}{\delta},\frac{1}{\delta}\right).$ The operator reduces to either a multiple of the mean curvature operator or the eikonal equation and otherwise we obtain $$F[u]$$. Similar results regarding ellipticity and consistency hold as for the one-dimensional model, which we present without proof. Lemma 6.16 For $$K=K(\delta)$$, $$L=L(\delta)$$ such that $$K\sqrt{L}\geq1$$, define the finite difference scheme \begin{align}\label{elliptic_scheme_regular} -F^{e,\delta}[u] = A^{\delta,+}\left(|\nabla u^h|^+, -{\it {\Delta}}_1^e u \right) + A^{\delta,-}\left(-|\nabla u^h|^-, -{\it {\Delta}}_1^e u \right). \end{align} (AC)e, δ Then, $$-F^{e,\delta}[u]$$ is elliptic and consistent with $$-F^\delta[u]$$. As with the one-dimensional case, $$K$$ and $$L$$ need to be chosen appropriately in order for $$-F^{e,\delta}[u]$$ to be consistent with $$-F[u]$$. Assuming $$K = \mathcal{O}(h^{-\alpha})$$, $$L = \mathcal{O}(h^{-\beta})$$ and $${\rm d}\theta = \mathcal{O}(1/n_\theta) = \mathcal{O}(h^\gamma)$$, full consistency is obtained when $$\alpha \in (0,1)$$, $$\beta \in \left(0,\frac{2}{3}\right)$$ and $$\gamma \in \left(\beta,\frac{2-\beta}{2}\right)$$. The optimal choices are $$\alpha \in [1/9,7/9]$$, $$\beta = 4/9$$, $$\gamma = 2/3$$ and, in such case, the accuracy is $$\mathcal{O}(h^{2/9})$$. As in the one-dimensional model, no accuracy is lost with the regularization (see Remark 6.11). The Lipschitz constant of $$F^{e,\delta}[u]$$ is given by $$\label{Lipconst} C^h = \frac{K}{h} + \frac{2 L}{(h \: n_\theta)^2}.$$ (6.1) In practice, we will choose $$K = c_K h^{-1/9}$$ and $$L = c_L h^{-4/9}$$, which leads to $C^h = c_Kh^{-10/9} +2\frac{c_L}{{n_\theta}^2}h^{-22/9}.$ These results are generalizations of the results proved for the one-dimensional model, whose proofs can be found in Oberman & Salvador (2016). We can finally define and prove the convergence of the discretizations for the time-dependent equation (AC). The time derivative is discretized with an explicit forward Euler step, whereas $$F[u]$$ is discretized using either the elliptic scheme $$F^{e,\delta}$$ or the regularized filtered scheme $$F^{f,\delta}[u]$$ (this is easily defined by replacing the elliptic scheme in $$F^f$$ by its regularized version). This leads to monotone (resp. filtered) time discretization with solution map given by $$\label{discretizationtime} u(\cdot,t+{\rm d}t) = u(\cdot,t) + {\rm d}t F^{e,\delta}[u(\cdot,t)], \quad \left(\text{resp.} = u(\cdot,t) + {\rm d}t F^{f,\delta}[u(\cdot,t)\right)\!.$$ (6.2) 6.5 Convergence theorems Having proved the ellipticity and consistency of the schemes, the uniform convergence follows as discussed in Section 3.1, provided there exists unique viscosity solutions to the PDEs (AC) and $$F[u] = f$$ along with the homogeneous Neumann boundary conditions, which is assured by the theory in Giga (2006), as explained above. We also need existence of solutions to the elliptic and filtered schemes. For the affine curvature evolution PDE (AC), this is immediate, because the schemes are explicit. In the case of the stationary equation, additional work is required to prove existence of solutions. These do not need to be unique to apply the convergence theory, but both existence and uniqueness results for the schemes follow from the discrete comparison principle and from finite dimensional fixed point theorems, as in Oberman (2006) and Froese & Oberman (2013). For readability, the existence of solutions to the schemes will be part of the assumptions in the theorem. In practice, we can get discrete comparison from a general theorem by adding a small multiple of $$u$$ to the scheme, Oberman (2006), or we can do extra work to prove it directly. There are two convergence results. The first is for the elliptic problem, where there is no need for regularization. The second is for the parabolic problem, where we need to regularize and use an explicit Euler time step. Theorem 6.17 Let $${\it {\Omega}}$$ be a bounded convex domain with a $$C^2$$ boundary and $$u$$ denote the unique viscosity solution of $$F[u] = f$$ in $${\it {\Omega}}$$, along with homogeneous Neumann boundary conditions on $$\partial {\it {\Omega}}$$. For each $$\epsilon = \epsilon(h,{\rm d}\theta) >0$$, let $$u^{e,\epsilon}$$, (resp. $$u^{f,\epsilon}$$) be the uniformly bounded solution of $$F^e[u] = f$$, (resp. $$F^f[u] = f$$). Then $u^{e,\epsilon} \to u ~\text{and}~ u^{f,\epsilon} \to u, \quad \text{locally uniformly, as}~ \epsilon \to 0.$ Proof. The assumptions on $${\it {\Omega}}$$ guarantee that the PDE satisfies a comparison principle. Convergence for the elliptic discretization $$-F^e$$ then follows from the Barles–Souganidis theorem. The modification of the proof for filtered schemes can be found in Froese & Oberman (2013). □ Theorem 6.18 Let $${\it {\Omega}}$$ be a bounded convex domain with a $$C^2$$ boundary. Assume that $$u(x,t)$$ is the unique viscosity solution of $$u_t = F[u]$$ in $${\it {\Omega}} \times [0,\infty)$$, along with $$u(x,0) = u_0(x)$$ and homogeneous Neumman boundary conditions. Assume as well that $$K$$ and $$L$$ are picked so that $$F^{e,\delta}$$ is consistent with $$F$$. For each $$\epsilon = \epsilon(h,{\rm d}\theta) >0$$, let $$u^{e,\epsilon}$$ (resp. $$u^{f,\epsilon}$$) be the uniformly bounded solution of the monotone (resp. filtered) time discretization (6.2) with $${\rm d}t \leq 1/{K^h}$$, where $$K^h$$ is the Lipschitz constant of $$F^{e,\delta}[u]$$ (6.1). Then $u^{e,\epsilon} \to u \quad \text{and} \quad u^{f,\epsilon} \to u$ locally uniformly, as $${\rm d}t,\epsilon \to 0$$. Proof. The assumptions on $${\it {\Omega}}$$ guarantee that the PDE satisfies a comparison principle. The elliptic scheme leads to a monotone time discretization, in other words, the solution map satisfies a comparison principle if $${\rm d}t \leq 1/{K^h}$$, where $$K^h$$ is the Lipschitz constant of the elliptic scheme by Oberman (2006). Then the Barles–Souganidis theorem (Barles & Souganidis, 1991) applies. For time discretization of the filtered scheme, the convergence of filtered schemes in Froese & Oberman (2013) applies. □ Remark 6.19 Just like in the one-dimensional model, for fixed values of $$K,L$$, there is a unique viscosity solution, $$u^\delta$$, of the regularized PDE. Then, fixing $$K,L$$ and using the discretization above with $${\rm d}t = \mathcal{O}((h\: n_\theta)^2)$$, the forward Euler method converges uniformly to $$u^\delta$$ as $$h\to 0$$. 7. Numerical results In this section, we start with a simple example to compare the affine invariant curvature motion and the regularized model. We present different examples on the evolution of a single curve in Section 7.1 under the affine curvature motion and compare it with the mean curvature motion. We test the accuracy of the schemes by computing solutions to the time-independent PDE in Section 7.2 and to the time-dependent equation in Section 7.3. We mostly considered stationary problems, because it was easier to generate exact solutions by applying the operator $$F$$ to a given function $$u$$, and including $$f = F[u]$$ as a source function. For the time-dependent problem, we took advantage of the exact solution from Lemma 2.12 to compare the accuracy of the solutions after a short time, $$T= 0.1$$. We test as well in Section 7.4 whether the schemes satisfy numerically the morphology and affine invariance properties that (AC) satisfies (see Theorem 2.10). Definition 7.1 (Parameters for the discretization) We used the regularized schemes with $$K = 20 h^{-1/9}$$ and $$L = 20 h^{-4/9}$$. For the elliptic discretization, we used two different stencils: the narrow elliptic scheme used $$n_\theta = 3$$ and the wider elliptic scheme used $$n_\theta = 7$$. The filtered discretization used the wider elliptic scheme with $\epsilon(h,d\theta) = \sqrt{h} + d\theta/10.$ For the forward Euler method with the elliptic and filtered schemes, we used a constant time step of $${\rm d}t = 1/C_h$$ where $$C_h$$ is given by (6.1). For the forward Euler method with the standard finite difference scheme, we used the same time step as the filtered scheme, except when computing the steady-state solutions in Example 7.9, where we used $${\rm d}t = h^2/2$$ (this choice of time step proved to be enough in practice). The stopping criteria was simply that the $$l^\infty$$ norm of the residual $$\left\vert{F^{f,\delta}[u^n]-f}\right\vert$$ (and similarly for $$F^{e,\delta}, F^a$$) was below $$tol = 10^{-5}$$. Example 7.2 [Comparing the regularized to the unregularized operators] We compare the evolution of an ellipse by the affine invariant curvature motion and by the regularized model. The ellipse should remain an ellipse of fixed eccentricity. The results were obtained by numerically solving (AC) and $$u_t = F^\delta[u]$$ respectively, with the initial condition $u_0(x,y) = \min\left\{\left(\frac{x}{2}\right)^2+y^2-1,1\right\}$ and homogeneous Neumann boundary conditions on $${\it {\Omega}} = [-4,4]^2$$. The narrow elliptic scheme was used for the spatial discretization. The difference between the level sets of the solution of the two equations is visually indistinguishable. Measured in the $$l^\infty$$ norm ranges from $$10^{-6}$$ for early times steps and $$10^{-7}$$ for the later time steps. The level sets of the numerical solution obtained by solving the regularized PDE are depicted in Fig. 5. Fig. 5. View largeDownload slide Evolution of an ellipse by (left) affine curvature, (right) mean curvature (top). Evolution by affine curvature of (left) a diamond and (right) a flatter diamond (bottom). Fig. 5. View largeDownload slide Evolution of an ellipse by (left) affine curvature, (right) mean curvature (top). Evolution by affine curvature of (left) a diamond and (right) a flatter diamond (bottom). Remark 7.3 Although the regularized scheme with the more stringent time step is needed for the convergence proof, in practise we found, as reported in the previous example, using the (unregularized) elliptic discretization with the time step $${\rm d}t = h^2/2$$ gave nearly identical results to within two or more significant digits. For the remaining examples, we present results using the regularized schemes (with $$K = 20h^{-1/9}$$ and $$L = 20h^{-4/9}$$). 7.1 Numerical examples showing curve evolution In this section, we present some numerical examples illustrating the geometric properties of the PDEs. Example 7.4 (Ellipse) We compare the evolution of an ellipse by the affine invariant curvature motion and by mean curvature. For the former, the ellipse should remain an ellipse of fixed eccentricity. For the latter, the ellipse asymptotically approaches a circle instead. The results were obtained by numerically solving (AC) and (MC), respectively, with the initial condition $u_0(x,y) = \min\left\{\left(\frac{x}{2}\right)^2+y^2-1,1\right\}$ and homogeneous Neumann boundary conditions. We took $$[-4,4]^2$$ as the computational domain on a $$128\times 128$$ grid. As for the scheme used, we chose the narrow elliptic schemes for both. In Fig. 5, we plot the zero level sets obtained at $$t \in \{0,0.1,0.3,0.5,0.7,0.9\}$$. Example 7.5 (Diamond) We also compute the solution of (AC) with initial condition $(a) \ u_0(x,y) = \min\left\{|x|+|y|-1,1\right\}, \qquad (b) \ u_0(x,y) = \min\left\{|x|+2|y|-1,1\right\}$ and homogeneous Neumann boundary conditions. The exact solutions are not known. However, we do know that smooth convex curves evolving under affine curvature converge to ellipses until collapsing to a point and that is the behavior we observed here (see Fig. 5). We took $$[-2,2]^2$$ as the computational domain on a $$128 \times 128$$ grid. As for the scheme, we chose again the narrow elliptic scheme. In Fig. 5, we plot the zero level sets of the numerical solution from time $$t = 0$$ to $$t = 0.5$$ in increments of $$0.1$$ for example $$(a)$$ and from $$t=0$$ to $$t=0.3$$ in increments of $$0.05$$ for example $$(b)$$. Example 7.6 (Fan-shape curve) We also compute the solution of (AC) and (MC) with initial condition $u_0(x,y) = \min\left\{c_+^{(1)}(x,y),c_-^{(1)}(x,y),c_+^{(2)}(x,y),c_-^{(2)}(x,y),1\right\}\!,$ where $c_\pm^{(1)}(x,y) = \left(x\pm\frac{1}{2}\right)^2+5\left(y\pm\frac{1}{4}\right)^2-\frac{1}{2} \quad \text{and} \quad c_\pm^{(2)}(x,y) = 5\left(x\pm\frac{1}{4}\right)^2+\left(y\mp\frac{1}{4}\right)^2-\frac{1}{2}$ and homogeneous Neumann boundary conditions. The exact solution is not known. As in the previous example, we took $$[-2,2]^2$$ as the computational domain on a $$128 \times 128$$ grid together with the narrow elliptic scheme. Under the both curvature motions, the curve should initially become convex. At later times, under the affine curvature motion, the curve should evolve to an ellipse as opposed to a circle in the mean curvature case. In Fig. 6, we plot the zero level set of the numerical solutions at time $$t \in \{0,0.05,0.1,0.2\}$$ and observe the exact behavior described. Fig. 6. View largeDownload slide Evolution of a fan-shape like curve under affine curvature motion (top) and mean curvature (bottom) at time $$t~\in~\{0,0.05,0.1,0.2\}$$ (see Example 7.6 for more details). Fig. 6. View largeDownload slide Evolution of a fan-shape like curve under affine curvature motion (top) and mean curvature (bottom) at time $$t~\in~\{0,0.05,0.1,0.2\}$$ (see Example 7.6 for more details). 7.2 Accuracy of stationary solutions We test the accuracy for the following Dirichlet problem $\begin{cases} F[u] = f, &~ \text{in}~ {\it {\Omega}},\\ u = g, & ~\text{on}~ \partial {\it {\Omega}}. \end{cases}$ Solutions are obtained by computing the steady state solution of $$u_t = F[u] - f$$, with $$u(\cdot,t) = g$$ on $$\partial{\it {\Omega}}$$ and $$u(x,0) = u_0(x)$$. We set $${\rm d}t = 1/C_h$$, where $$C_h$$ is given by (6.1) for the elliptic and filtered schemes and $${\rm d}t = h^2/2$$ for the standard finite difference scheme scheme. These examples also demonstrate stability of the time dependent problem for elliptic and filtered scheme, as well as convergence to the steady solution, because we used the time-dependent problem to obtain the solution. Remark 7.7 The unregularized schemes were also used. For these we set $${\rm d}t = h^2/2$$ for all examples, except for the filtered scheme in example (d) where we set $${\rm d}t = h^2/8$$. The results obtained were virtually the same. We set $$u_0(x)$$ to be the exact solution in a layer of seven grid points adjacent to the boundary. As a result, each discretization is initialized at the same set of grid points and, therefore, we can make a fair comparison of their accuracy. On the coarsest grid, we set $$u_0(x) = 0$$ on the interior grid points. To speed up calculations, on finer grids we set $$u_0(x)$$ to be interpolated solutions from the coarser grids at interior grid points. Remark 7.8 When considering Dirichlet boundary conditions, we lose the comparison principle. However, we can always cap off the solution and impose homogeneous boundary conditions. Example 7.9 We consider the following exact solutions \begin{align*} (a) & \quad u(x,y) = x^2+y^2, \quad f(x) = 2 \left(x^2 + y^2\right)^{\frac{1}{3}},\\ (b) & \quad u(x,y) = {\rm e}^{x^2+y^2}, \quad f(x,y) = 2 \left({\rm e}^{3 (x^2 + y^2)}(x^2 + y^2)\right)^{\frac{1}{3}}\!,\\ (c) & \quad u(x,y) = \left(x^2 + y^2\right)^{\frac{1}{3}}, \quad f(x) = \frac{4}{3},\\ (d) & \quad u(x,y) = \frac{\sin(2\pi x)\sin(2\pi y)}{4}, \quad f(x) = \frac{\pi^{\frac{4}{3}}}{2} \left(-\left(2+\cos(4\pi x)+\cos(4\pi y)\right) \sin(2\pi x)\sin(2\pi y)\right)^{\frac{1}{3}}\!, \end{align*} with $${\it {\Omega}} = [-1,1]^2$$. The solutions in (a),(b) and (d) are smooth, but the functions $$f$$ are only Holder continuous, $$C^{0,2/3}$$, with singularities at the origin for (a) and (b), and at several points in example (d). The solution in (c) is only $$C^{0,2/3}$$ with a singularity at the origin, but in this case the function $$f$$ is constant. The results are presented in Table 1. In Example 7.9$$(a)$$ the solution is a quadratic polynomial. The accurate scheme $$F^a$$ gives essentially machine precision, which is not surprising, because this scheme is second-order accurate. The filtered scheme $$F^f$$ gives nearly the same accuracy (with a small discrepancy which could be eliminated by tuning the parameter). On the other hand, both the narrow and the wider elliptic schemes are much less accurate. In contrast, in Example 7.9(b), the solution is smooth, but not quadratic. In this case, we see that the accurate scheme appears to be converging to $$\mathcal{O}(h)$$. The elliptic schemes are less accurate and the filtered scheme is in between. In Example 7.9$$(c)$$, the solution is only Holder continuous. In this case, the elliptic schemes have almost constant error near $$0.01$$ over the range of parameters used. Despite the singular solution, the accurate scheme gives accuracy $$\mathcal{O}(h)$$ and the filtered scheme does just as well. Example 7.9$$(d)$$ shows that the standard finite difference scheme does not converge as discussed in Section 5 (see Fig. 3). The narrow elliptic scheme has almost a constant error $$0.03,$$ and the wider elliptic scheme has an error $$0.1$$ for the smallest grid, decreasing by a factor of two as the grid is refined. The filtered scheme has the best accuracy, achieving an error of $$0.001$$ at the finest grid. In general, we expect better accuracy from the filtered schemes by making a more careful choice of the parameter $$\epsilon$$ (see Remark 6.14). Table 1 Accuracy in the $$l^\infty$$ norm and order of convergence of the schemes for Example 7.9 with regularized schemes Errors and order, Example 7.9 (a) N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.565e$$-$$02 3.978e$$-$$02 2.922e$$-$$07 5.104e$$-$$07 64 2.555e$$-$$02 2.696e$$-$$02 3.019e$$-$$07 2.827e$$-$$07 128 1.623e$$-$$02 1.678e$$-$$02 1.982e$$-$$07 1.984e$$-$$07 256 1.050e$$-$$02 1.055e$$-$$02 9.293e$$-$$08 5.678e$$-$$05 Errors and order, Example 7.9 (b) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 5.462e$$-$$02 7.666e$$-$$02 1.922e$$-$$03 1.985e$$-$$03 64 5.207e$$-$$02 5.872e$$-$$02 9.875e$$-$$04 8.904e$$-$$03 128 4.105e$$-$$02 3.725e$$-$$02 3.385e$$-$$04 7.262e$$-$$03 256 3.173e$$-$$02 2.240e$$-$$02 9.798e$$-$$05 8.065e$$-$$03 Errors and order, Example 7.9 (c) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 1.387e$$-$$02 2.386e$$-$$02 1.129e$$-$$02 1.129e$$-$$02 64 1.310e$$-$$02 7.683e$$-$$03 4.625e$$-$$03 4.625e$$-$$03 128 9.302e$$-$$03 8.202e$$-$$03 1.859e$$-$$03 1.872e$$-$$03 256 6.445e$$-$$03 7.156e$$-$$03 7.426e$$-$$04 7.570e$$-$$04 Errors and order, Example 7.9 (d) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.964e$$-$$02 9.813e$$-$$02 — 1.926e$$-$$02 64 3.788e$$-$$02 4.679e$$-$$02 — 8.309e$$-$$03 128 3.688e$$-$$02 2.367e$$-$$02 — 2.482e$$-$$03 256 3.003e$$-$$02 1.798e$$-$$02 — 9.697e$$-$$04 Errors and order, Example 7.9 (a) N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.565e$$-$$02 3.978e$$-$$02 2.922e$$-$$07 5.104e$$-$$07 64 2.555e$$-$$02 2.696e$$-$$02 3.019e$$-$$07 2.827e$$-$$07 128 1.623e$$-$$02 1.678e$$-$$02 1.982e$$-$$07 1.984e$$-$$07 256 1.050e$$-$$02 1.055e$$-$$02 9.293e$$-$$08 5.678e$$-$$05 Errors and order, Example 7.9 (b) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 5.462e$$-$$02 7.666e$$-$$02 1.922e$$-$$03 1.985e$$-$$03 64 5.207e$$-$$02 5.872e$$-$$02 9.875e$$-$$04 8.904e$$-$$03 128 4.105e$$-$$02 3.725e$$-$$02 3.385e$$-$$04 7.262e$$-$$03 256 3.173e$$-$$02 2.240e$$-$$02 9.798e$$-$$05 8.065e$$-$$03 Errors and order, Example 7.9 (c) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 1.387e$$-$$02 2.386e$$-$$02 1.129e$$-$$02 1.129e$$-$$02 64 1.310e$$-$$02 7.683e$$-$$03 4.625e$$-$$03 4.625e$$-$$03 128 9.302e$$-$$03 8.202e$$-$$03 1.859e$$-$$03 1.872e$$-$$03 256 6.445e$$-$$03 7.156e$$-$$03 7.426e$$-$$04 7.570e$$-$$04 Errors and order, Example 7.9 (d) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.964e$$-$$02 9.813e$$-$$02 — 1.926e$$-$$02 64 3.788e$$-$$02 4.679e$$-$$02 — 8.309e$$-$$03 128 3.688e$$-$$02 2.367e$$-$$02 — 2.482e$$-$$03 256 3.003e$$-$$02 1.798e$$-$$02 — 9.697e$$-$$04 Table 1 Accuracy in the $$l^\infty$$ norm and order of convergence of the schemes for Example 7.9 with regularized schemes Errors and order, Example 7.9 (a) N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.565e$$-$$02 3.978e$$-$$02 2.922e$$-$$07 5.104e$$-$$07 64 2.555e$$-$$02 2.696e$$-$$02 3.019e$$-$$07 2.827e$$-$$07 128 1.623e$$-$$02 1.678e$$-$$02 1.982e$$-$$07 1.984e$$-$$07 256 1.050e$$-$$02 1.055e$$-$$02 9.293e$$-$$08 5.678e$$-$$05 Errors and order, Example 7.9 (b) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 5.462e$$-$$02 7.666e$$-$$02 1.922e$$-$$03 1.985e$$-$$03 64 5.207e$$-$$02 5.872e$$-$$02 9.875e$$-$$04 8.904e$$-$$03 128 4.105e$$-$$02 3.725e$$-$$02 3.385e$$-$$04 7.262e$$-$$03 256 3.173e$$-$$02 2.240e$$-$$02 9.798e$$-$$05 8.065e$$-$$03 Errors and order, Example 7.9 (c) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 1.387e$$-$$02 2.386e$$-$$02 1.129e$$-$$02 1.129e$$-$$02 64 1.310e$$-$$02 7.683e$$-$$03 4.625e$$-$$03 4.625e$$-$$03 128 9.302e$$-$$03 8.202e$$-$$03 1.859e$$-$$03 1.872e$$-$$03 256 6.445e$$-$$03 7.156e$$-$$03 7.426e$$-$$04 7.570e$$-$$04 Errors and order, Example 7.9 (d) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.964e$$-$$02 9.813e$$-$$02 — 1.926e$$-$$02 64 3.788e$$-$$02 4.679e$$-$$02 — 8.309e$$-$$03 128 3.688e$$-$$02 2.367e$$-$$02 — 2.482e$$-$$03 256 3.003e$$-$$02 1.798e$$-$$02 — 9.697e$$-$$04 Errors and order, Example 7.9 (a) N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.565e$$-$$02 3.978e$$-$$02 2.922e$$-$$07 5.104e$$-$$07 64 2.555e$$-$$02 2.696e$$-$$02 3.019e$$-$$07 2.827e$$-$$07 128 1.623e$$-$$02 1.678e$$-$$02 1.982e$$-$$07 1.984e$$-$$07 256 1.050e$$-$$02 1.055e$$-$$02 9.293e$$-$$08 5.678e$$-$$05 Errors and order, Example 7.9 (b) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 5.462e$$-$$02 7.666e$$-$$02 1.922e$$-$$03 1.985e$$-$$03 64 5.207e$$-$$02 5.872e$$-$$02 9.875e$$-$$04 8.904e$$-$$03 128 4.105e$$-$$02 3.725e$$-$$02 3.385e$$-$$04 7.262e$$-$$03 256 3.173e$$-$$02 2.240e$$-$$02 9.798e$$-$$05 8.065e$$-$$03 Errors and order, Example 7.9 (c) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 1.387e$$-$$02 2.386e$$-$$02 1.129e$$-$$02 1.129e$$-$$02 64 1.310e$$-$$02 7.683e$$-$$03 4.625e$$-$$03 4.625e$$-$$03 128 9.302e$$-$$03 8.202e$$-$$03 1.859e$$-$$03 1.872e$$-$$03 256 6.445e$$-$$03 7.156e$$-$$03 7.426e$$-$$04 7.570e$$-$$04 Errors and order, Example 7.9 (d) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.964e$$-$$02 9.813e$$-$$02 — 1.926e$$-$$02 64 3.788e$$-$$02 4.679e$$-$$02 — 8.309e$$-$$03 128 3.688e$$-$$02 2.367e$$-$$02 — 2.482e$$-$$03 256 3.003e$$-$$02 1.798e$$-$$02 — 9.697e$$-$$04 When comparing the different schemes, we have to account for the width of the stencil, because for the elliptic schemes the wider schemes also have a larger spatial discretization error. In general, the accuracy improved with the use of the wider stencil. Moreover, the filtered scheme performed as expected by providing better accuracy than the elliptic scheme and almost as good accuracy as the accurate scheme. The final example shows that the standard finite difference scheme may not converge. This may be due to the fact that, despite the solution being smooth, there were multiple points where $$f$$ was singular. Geometrically, this solution has several points where curves shrunk to zero. We also considered for Example 7.9 the elliptic schemes with $$n_\theta = 1$$ and $$2$$. These provided poor accuracy with the directional resolution error easily dominating the spatial resolution error. The errors do not decrease to zero as we decrease the grid size. This property is consistent with the theoretical results, because convergence is only proved as both $$h,{\rm d}\theta \to 0$$, which is indeed observed in the numerical results. On the other hand, using the narrow and wider schemes, the accuracy of the elliptic scheme was good enough for the filtered scheme to give accuracy comparable to the accurate scheme in many examples. In principle, we still need to send $${\rm d}\theta \to 0$$, but in practice, the rate at which it needs to go to zero is much slower than $$h$$ when the filtered scheme is used. Finally, we point that the choice of $$\epsilon$$ is not an easy one and it is hard to pick an $$\epsilon$$ that yields optimal results in each example presented. By choosing $$\epsilon$$ larger, it is possible to achieve with the filtered schemes the accuracy of the standard schemes in Examples 7.9 (a)–(c). However, such choice is too permissive for Example 7.9 (d). As pointed out in Section 6.3, $$\epsilon$$ needs to be chosen small enough for the monotone scheme to be active to ensure convergence. (This is comparable to the CFL condition in time-dependent equations: methods are convergent as $${\rm d}t,h \to 0$$, with $${\rm d}t$$ satisfying the CFL condition.) 7.3 Accuracy for the time dependent problem Recall here that for general boundary conditions, the time-dependent PDE (AC) requires Neumann boundary conditions for uniqueness of viscosity solutions to hold. We have already established (in the previous section) the stability of the numerical method. In the following examples, we test the accuracy of solutions, comparing two different wide stencil discretizations, along with regularized filtered discretization and the (generally unstable) standard finite difference method. Consequently, we test for the accuracy using Dirichlet boundary conditions coming from the exact solution of Lemma 2.12. Example 7.10 (Neumann BC) We consider the exact solution $u(x,y,t) = \min\left\{t + \frac{3}{4}\left(\frac{b}{a}x^2+\frac{a}{b}y^2\right)^{2/3}-1,0\right\},$ with $$a = 2$$ and $$b = 1$$ (see Lemma 2.12). By taking the minimum with $$0$$, we are imposing homogeneous Neumann boundary conditions, thus avoiding having to deal with boundary of the computation domain. We set $${\it {\Omega}} = [-3,3]^2$$. We display the numerical error in the $$l^\infty$$ norm at time $$T = 0.1$$ in Table 2. Table 2 Error in the $$l^\infty$$ norm of the whole computational domain at time $$t=0.1$$ for the time-dependent Examples 7.10 and 7.11 Errors and order, Example 7.10 N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 4.845e$$-$$02 6.691e$$-$$02 4.894e$$-$$02 4.888e$$-$$02 64 4.432e$$-$$02 4.607e$$-$$02 2.977e$$-$$02 2.975e$$-$$02 128 3.544e$$-$$02 2.823e$$-$$02 2.457e$$-$$02 2.438e$$-$$02 256 2.971e$$-$$02 2.080e$$-$$02 1.747e$$-$$02 1.724e$$-$$02 512 2.764e$$-$$02 1.652e$$-$$02 1.205e$$-$$02 1.182e$$-$$02 Errors and order, Example 7.11 N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 2.182e$$-$$02 1.449e$$-$$02 1.985e$$-$$02 1.985e$$-$$02 64 1.435e$$-$$02 1.160e$$-$$02 1.279e$$-$$02 1.279e$$-$$02 128 9.580e$$-$$03 7.517e$$-$$03 5.566e$$-$$03 5.567e$$-$$03 256 6.404e$$-$$03 4.854e$$-$$03 2.442e$$-$$03 2.409e$$-$$03 512 6.090e$$-$$03 4.288e$$-$$03 1.036e$$-$$03 1.002e$$-$$03 Errors and order, Example 7.10 N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 4.845e$$-$$02 6.691e$$-$$02 4.894e$$-$$02 4.888e$$-$$02 64 4.432e$$-$$02 4.607e$$-$$02 2.977e$$-$$02 2.975e$$-$$02 128 3.544e$$-$$02 2.823e$$-$$02 2.457e$$-$$02 2.438e$$-$$02 256 2.971e$$-$$02 2.080e$$-$$02 1.747e$$-$$02 1.724e$$-$$02 512 2.764e$$-$$02 1.652e$$-$$02 1.205e$$-$$02 1.182e$$-$$02 Errors and order, Example 7.11 N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 2.182e$$-$$02 1.449e$$-$$02 1.985e$$-$$02 1.985e$$-$$02 64 1.435e$$-$$02 1.160e$$-$$02 1.279e$$-$$02 1.279e$$-$$02 128 9.580e$$-$$03 7.517e$$-$$03 5.566e$$-$$03 5.567e$$-$$03 256 6.404e$$-$$03 4.854e$$-$$03 2.442e$$-$$03 2.409e$$-$$03 512 6.090e$$-$$03 4.288e$$-$$03 1.036e$$-$$03 1.002e$$-$$03 Table 2 Error in the $$l^\infty$$ norm of the whole computational domain at time $$t=0.1$$ for the time-dependent Examples 7.10 and 7.11 Errors and order, Example 7.10 N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 4.845e$$-$$02 6.691e$$-$$02 4.894e$$-$$02 4.888e$$-$$02 64 4.432e$$-$$02 4.607e$$-$$02 2.977e$$-$$02 2.975e$$-$$02 128 3.544e$$-$$02 2.823e$$-$$02 2.457e$$-$$02 2.438e$$-$$02 256 2.971e$$-$$02 2.080e$$-$$02 1.747e$$-$$02 1.724e$$-$$02 512 2.764e$$-$$02 1.652e$$-$$02 1.205e$$-$$02 1.182e$$-$$02 Errors and order, Example 7.11 N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 2.182e$$-$$02 1.449e$$-$$02 1.985e$$-$$02 1.985e$$-$$02 64 1.435e$$-$$02 1.160e$$-$$02 1.279e$$-$$02 1.279e$$-$$02 128 9.580e$$-$$03 7.517e$$-$$03 5.566e$$-$$03 5.567e$$-$$03 256 6.404e$$-$$03 4.854e$$-$$03 2.442e$$-$$03 2.409e$$-$$03 512 6.090e$$-$$03 4.288e$$-$$03 1.036e$$-$$03 1.002e$$-$$03 Errors and order, Example 7.10 N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 4.845e$$-$$02 6.691e$$-$$02 4.894e$$-$$02 4.888e$$-$$02 64 4.432e$$-$$02 4.607e$$-$$02 2.977e$$-$$02 2.975e$$-$$02 128 3.544e$$-$$02 2.823e$$-$$02 2.457e$$-$$02 2.438e$$-$$02 256 2.971e$$-$$02 2.080e$$-$$02 1.747e$$-$$02 1.724e$$-$$02 512 2.764e$$-$$02 1.652e$$-$$02 1.205e$$-$$02 1.182e$$-$$02 Errors and order, Example 7.11 N Narrow elliptic ($$n_\theta = 3$$) Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 2.182e$$-$$02 1.449e$$-$$02 1.985e$$-$$02 1.985e$$-$$02 64 1.435e$$-$$02 1.160e$$-$$02 1.279e$$-$$02 1.279e$$-$$02 128 9.580e$$-$$03 7.517e$$-$$03 5.566e$$-$$03 5.567e$$-$$03 256 6.404e$$-$$03 4.854e$$-$$03 2.442e$$-$$03 2.409e$$-$$03 512 6.090e$$-$$03 4.288e$$-$$03 1.036e$$-$$03 1.002e$$-$$03 Example 7.11 (Dirichlet BC) We consider the exact solution $u(x,y,t) = t + \frac{3}{4}\left(\frac{b}{a}x^2+\frac{a}{b}y^2\right)^{2/3},$ with $$a = 2$$ and $$b = 1$$ (see Lemma 2.12). This is the same example as in Example 7.10, but we consider Dirichlet boundary conditions instead. Therefore, we prescribe the exact solution at a seven-point layer at the boundary for all time $$t$$. This way all schemes are initialized at the same set of grid points, and thus we can compare their accuracy. The error in the $$l^\infty$$ norm at time $$T = 0.1$$ is presented in Table 2. When using Neumman boundary conditions, we observed slow convergence, with errors near $$0.01$$, slowly decreasing as the grid size improved. The accuracy improved as we went from the narrow to the wider elliptic scheme, and further improved as we went to the accurate and then the filtered scheme. In the case of Dirichlet boundary conditions, the accuracy is better overall and the error decrease is slightly faster. The difference is explained by the cap off done in the Neumman boundary conditions that introduces an additional error in the solution. However, this error does not propagate to the whole domain as the level sets of the solution shrink to its interior, and so, when we look at the error away from the cap off, we recover results very similar to the Dirichlet boundary conditions. 7.4 Numerical study of the morphology and affine invariance properties In this section, we test whether our proposed schemes satisfy numerically the morphology and affine invariance properties that (AC) satisfies (see Theorem 2.10). Example 7.12 In this example, we test numerically whether the schemes presented here satisfy the morphology property of the affine curvature evolution ((2) in Theorem 2.10). We consider two examples: $$(a)$$$$g_1(x) = {\rm e}^x$$, $$(b)$$$$g_2(x) = x^3$$. We take $u_0(x,y) = \min\left\{\left(\frac{x}{2}\right)^2+y^2-1,0\right\}$ and compare $${\it {\Phi}}_t(g_v\circ u_0)$$ with $$g_v \circ {\it {\Phi}}_t(u_0)$$ at $$t=1$$ for $$v = 1,2$$. We took $$[-3,3]^2$$ as the computational domain with homogeneous Neumann boundary conditions. The results are displayed in Table 3. The difference in the $$l^\infty$$ norm is one order of magnitude smaller than the observed accuracy for the schemes in Section 7.2. Based on these examples, the morphology property seems to hold numerically. Table 3 Difference in the $$l^\infty$$ norm between $${\it {\Phi}}_t(g_v \circ u_0)$$ and $$g_v\circ {\it {\Phi}}_t(u_0)$$ for $$v=1,2$$ for Example 7.12 Difference in the $$l^\infty$$ norm, Example 7.12 (a) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 8.943e$$-$$03 1.037e$$-$$02 3.419e$$-$$03 3.446e$$-$$03 64 5.709e$$-$$03 6.111e$$-$$03 1.954e$$-$$03 1.970e$$-$$03 128 4.061e$$-$$03 3.257e$$-$$03 1.061e$$-$$03 1.109e$$-$$03 256 3.115e$$-$$03 1.871e$$-$$03 5.907e$$-$$04 6.109e$$-$$04 512 2.604e$$-$$03 1.161e$$-$$03 3.207e$$-$$04 3.397e$$-$$04 Difference in the $$l^\infty$$ norm, Example 7.12 (b) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.450e$$-$$02 7.023e$$-$$02 6.947e$$-$$03 6.947e$$-$$03 64 1.730e$$-$$02 2.842e$$-$$02 1.881e$$-$$03 1.877e$$-$$03 128 1.032e$$-$$02 9.355e$$-$$03 6.254e$$-$$04 6.333e$$-$$04 256 6.894e$$-$$03 4.307e$$-$$03 2.283e$$-$$04 2.307e$$-$$04 512 5.372e$$-$$03 2.316e$$-$$03 8.343e$$-$$05 8.506e$$-$$05 Difference in the $$l^\infty$$ norm, Example 7.12 (a) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 8.943e$$-$$03 1.037e$$-$$02 3.419e$$-$$03 3.446e$$-$$03 64 5.709e$$-$$03 6.111e$$-$$03 1.954e$$-$$03 1.970e$$-$$03 128 4.061e$$-$$03 3.257e$$-$$03 1.061e$$-$$03 1.109e$$-$$03 256 3.115e$$-$$03 1.871e$$-$$03 5.907e$$-$$04 6.109e$$-$$04 512 2.604e$$-$$03 1.161e$$-$$03 3.207e$$-$$04 3.397e$$-$$04 Difference in the $$l^\infty$$ norm, Example 7.12 (b) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.450e$$-$$02 7.023e$$-$$02 6.947e$$-$$03 6.947e$$-$$03 64 1.730e$$-$$02 2.842e$$-$$02 1.881e$$-$$03 1.877e$$-$$03 128 1.032e$$-$$02 9.355e$$-$$03 6.254e$$-$$04 6.333e$$-$$04 256 6.894e$$-$$03 4.307e$$-$$03 2.283e$$-$$04 2.307e$$-$$04 512 5.372e$$-$$03 2.316e$$-$$03 8.343e$$-$$05 8.506e$$-$$05 Table 3 Difference in the $$l^\infty$$ norm between $${\it {\Phi}}_t(g_v \circ u_0)$$ and $$g_v\circ {\it {\Phi}}_t(u_0)$$ for $$v=1,2$$ for Example 7.12 Difference in the $$l^\infty$$ norm, Example 7.12 (a) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 8.943e$$-$$03 1.037e$$-$$02 3.419e$$-$$03 3.446e$$-$$03 64 5.709e$$-$$03 6.111e$$-$$03 1.954e$$-$$03 1.970e$$-$$03 128 4.061e$$-$$03 3.257e$$-$$03 1.061e$$-$$03 1.109e$$-$$03 256 3.115e$$-$$03 1.871e$$-$$03 5.907e$$-$$04 6.109e$$-$$04 512 2.604e$$-$$03 1.161e$$-$$03 3.207e$$-$$04 3.397e$$-$$04 Difference in the $$l^\infty$$ norm, Example 7.12 (b) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.450e$$-$$02 7.023e$$-$$02 6.947e$$-$$03 6.947e$$-$$03 64 1.730e$$-$$02 2.842e$$-$$02 1.881e$$-$$03 1.877e$$-$$03 128 1.032e$$-$$02 9.355e$$-$$03 6.254e$$-$$04 6.333e$$-$$04 256 6.894e$$-$$03 4.307e$$-$$03 2.283e$$-$$04 2.307e$$-$$04 512 5.372e$$-$$03 2.316e$$-$$03 8.343e$$-$$05 8.506e$$-$$05 Difference in the $$l^\infty$$ norm, Example 7.12 (a) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 8.943e$$-$$03 1.037e$$-$$02 3.419e$$-$$03 3.446e$$-$$03 64 5.709e$$-$$03 6.111e$$-$$03 1.954e$$-$$03 1.970e$$-$$03 128 4.061e$$-$$03 3.257e$$-$$03 1.061e$$-$$03 1.109e$$-$$03 256 3.115e$$-$$03 1.871e$$-$$03 5.907e$$-$$04 6.109e$$-$$04 512 2.604e$$-$$03 1.161e$$-$$03 3.207e$$-$$04 3.397e$$-$$04 Difference in the $$l^\infty$$ norm, Example 7.12 (b) N Narrow elliptic $$(n_\theta = 3)$$ Wide elliptic ($$n_\theta = 7$$) Standard Filter 32 3.450e$$-$$02 7.023e$$-$$02 6.947e$$-$$03 6.947e$$-$$03 64 1.730e$$-$$02 2.842e$$-$$02 1.881e$$-$$03 1.877e$$-$$03 128 1.032e$$-$$02 9.355e$$-$$03 6.254e$$-$$04 6.333e$$-$$04 256 6.894e$$-$$03 4.307e$$-$$03 2.283e$$-$$04 2.307e$$-$$04 512 5.372e$$-$$03 2.316e$$-$$03 8.343e$$-$$05 8.506e$$-$$05 Example 7.13 In this example, we do a qualitative test of the affine invariance property, which in practice is what one needs for applications in image analysis. To do so, we plot the level sets of the affine invariant motion by curvature (AC) with $u(x,y) = \left(\frac{x}{2}\right)^2+y^2-1$ and $$u \circ \phi$$ as the initial solutions. For the affine transformations $$\phi(\mathbf{x}) = A\mathbf{x}$$, we consider $(a) \text{(rotation by \pi/4)} A = \begin{bmatrix}\frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2}\\[4pt] -\frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2}\end{bmatrix}, \quad (b) \: A = \begin{bmatrix}\frac{1}{2} & 1\\[4pt] 1 & \frac{1}{2}\end{bmatrix}.$ We take $$[-5,5]^2$$ as the computational domain on a $$256 \times 256$$ grid. In Fig. 7, we plot the zero level set of $${\it {\Phi}}_t(u\circ \phi)$$ and $$\left({\it {\Phi}}_{t(\det \phi)^{2/3}}(u)\right) \circ \phi$$ at $$t=1$$. The filtered scheme provided the best results, being indistinguishable to the naked eye. For a long time, the elliptic scheme did not provide as good results, a consequence of its lower accuracy (see Section 7.2). As for the standard finite difference scheme only in (a) the difference is indistinguishable to the naked eye like the filtered scheme. For (b), where $$A$$ is not a special affine transformation, the difference is significant, but can be removed by taking time steps of half the size (For shorter times all the curves were very close.) We point out that when the affine transformation is a rotation by a multiple of $$\pi/2$$ or a reflection over a line $$L$$ that makes an angle multiple of $$\pi/4$$ with the x-axis, we obtain essentially machine precision. (The difference in the $$l^\infty$$ norm was of the order $$10^{-10}$$.) This is expected, because our stencil is invariant under these transformations. Fig. 7. View largeDownload slide Plot of the zero level sets in Example 7.13 of $${\it {\Phi}}_t(u\circ \phi)$$ and $$\left({\it {\Phi}}_{t(\det \phi)^{2/3}}(u)\right) \circ \phi$$ for regular elliptic scheme (left), standard scheme (center) and regular filtered scheme (right) at time $$t=1$$ with $$\phi$$ given by $$(a)$$ (top) and $$(b)$$ (bottom). Fig. 7. View largeDownload slide Plot of the zero level sets in Example 7.13 of $${\it {\Phi}}_t(u\circ \phi)$$ and $$\left({\it {\Phi}}_{t(\det \phi)^{2/3}}(u)\right) \circ \phi$$ for regular elliptic scheme (left), standard scheme (center) and regular filtered scheme (right) at time $$t=1$$ with $$\phi$$ given by $$(a)$$ (top) and $$(b)$$ (bottom). 8. Conclusions We presented a convergent finite difference discretization of the PDE for motion of level sets by affine curvature in two dimensions. Computational examples demonstrated that the standard finite difference method is unstable, which motivates the need for a convergent method. The foundation of the scheme used an existing wide stencil discretization of the mean curvature operator, combined with an elliptic discretization of the positive and negative eikonal operators, $$\pm \vert{\nabla u}\vert$$. However, explicit time discretizations require Lipschitz continuous operators, which the affine curvature operator fails to be. Thus, we approximated it by a Lipschitz continuous regularization. In theory, the explicit Euler discretization is stable using a time step $$\text{d}\textit{t} \leq \mathcal{O}(h^{22/9})$$, with the constant determined by the width of the stencil. In practice, we achieved numerically equivalent results using $${\rm d}t = h^2/2$$ and without the regularization, although there is no proof of stability with the less restrictive time step. A careful choice of the regularization parameters allowed for the regularized elliptic scheme to maintain the same order of accuracy as the unregularized scheme while being provably convergent. The lower accuracy of both schemes, which results from the singularity of the operator, is overcome by the use of the filtered schemes, which essentially attain the accuracy of the standard finite difference schemes while being provably convergent. Simulations demonstrated the geometric properties of the PDE were nearly preserved by the numerical solutions, including affine invariance, morphological properties and the accurate representation of the shrinking ellipses. Simulations validated the convergence of the elliptic scheme and the improved accuracy of the filtered scheme. Funding FCT doctoral grant (SFRH/BD/84041/2012 to T.S.), in part. References Angenent S. , Sapiro G. & Tannenbaum A. ( 1998 ) On the affine heat equation for non-convex curves. J. Amer. Math. Soc. , 11 , 601 – 634 . Google Scholar CrossRef Search ADS Barles G. & Souganidis P. E. ( 1991 ) Convergence of approximation schemes for fully nonlinear second order equations. Asymptot. Anal. , 4 , 271 – 283 . 92d:35137 Barnett V. ( 1976 ) The ordering of multivariate data. J. R. Stat. Soc. , 139 , 318 – 355 . Calder J. & Smart C. K. ( 2016 ) The limit shape of convex peeling (in preparation). Calder J. , Esedoglu S. & Hero A. O. ( 2014 ) A Hamilton–Jacobi equation for the continuum limit of non-dominated sorting. SIAM J. Math. Anal. , 46 , 603 – 638 . Google Scholar CrossRef Search ADS Cao F. , ( 2003 ) Geometric Curve Evolution and Image Processing . Lecture Notes in Mathematics , vol. 1805 . Berlin : Springer . Google Scholar CrossRef Search ADS Catté F. Dibos F. & Koepfler G. ( 1995 ) A morphological scheme for mean curvature motion and applications to anisotropic diffusion and motion of level sets. SIAM J. Numer. Anal. , 32 , 1895 – 1909 . Google Scholar CrossRef Search ADS Chazelle B. ( 1985 ) On the convex layers of a planar set. IEEE Trans. Inform. Theory , 31 , 509 – 517 . Google Scholar CrossRef Search ADS Ciomaga A. , Monasse P. & Morel J.-M. ( 2011 ) Image visualization and restoration by curvature motions. SIAM Multiscale Model. Simul. , 9 , 834 – 871 . Google Scholar CrossRef Search ADS Crandall M. G. , Ishii H. & Lions P.-L. ( 1992 ) User’s guide to viscosity solutions of second order partial differential equations. Bull. Amer. Math. Soc. (N.S.) , 27 , 1 – 67 . Google Scholar CrossRef Search ADS Deckelnick K. , Dziuk G. & Elliott C. M. ( 2005 ) Computation of geometric partial differential equations and mean curvature flow. Acta Numer. 14 , 139 – 232 . Google Scholar CrossRef Search ADS Esedoglu S. , Ruuth S. & Tsai R. ( 2010 ) Diffusion generated motion using signed distance functions. J. Comput. Phys. , 229 , 1017 – 1042 . Google Scholar CrossRef Search ADS Evans L. C. & Spruck J. ( 1999 ) Motion of level sets by mean curvature, I. Fundamental Contributions to the Continuum Theory of Evolving Phase Interfaces in Solids ( Ball J. M. , Kinderlehrer D. , Podio-Guidugli P. & Slemrod M. ). Berlin : Springer , pp. 328 – 374 . Froese B. D. & Oberman A. M. ( 2013 ) Convergent filtered schemes for the Monge–Ampère partial differential equation. SIAM J. Numer. Anal. , 51 , 423 – 444 . Google Scholar CrossRef Search ADS Gage M. E. ( 1984 ) Curve shortening makes convex curves circular. Invent. Math. , 76 , 357 – 364 . Google Scholar CrossRef Search ADS Giga Y. ( 2006 ) Surface Evolution Equations . Birkhäuser Basel : Springer . Guichard F. ( 1994 ) Axiomatisation des analyses multi-echelles d’images et de films. Ph.D. Thesis , Université Paris Dauphine . Guichard F. & Morel J.-M. ( 1997 ) Partial Differential Equations and Image Iterative Filtering . Institute of Mathematics and its Applications Conference Series , vol. 63 . New York : Oxford University Press , pp. 525 – 562 . Liu R. Y. , Parelius J. M. , Singh K. ( 1999 ) Multivariate analysis by data depth: descriptive statistics, graphics and inference (with discussion and a rejoinder by liu and singh). Ann. Stat. , 27 , 783 – 858 . Merriman B. , Bence J. K. & Osher S. J. ( 1994 ) Motion of multiple junctions: a level set approach. J. Comput. Phys. , 112 , 334 – 363 . Google Scholar CrossRef Search ADS Moisan L. ( 1998 ) Affine plane curve evolution: a fully consistent scheme. IEEE Trans. Image Process. , 7 , 411 – 420 . Google Scholar CrossRef Search ADS PubMed Oberman A. M. ( 2004 ) A convergent monotone difference scheme for motion of level sets by mean curvature. Numer. Math. , 99 , 365 – 379 . Google Scholar CrossRef Search ADS Oberman A. M. ( 2006 ) Convergent difference schemes for degenerate elliptic and parabolic equations: Hamilton–Jacobi equations and free boundary problems. SIAM J. Numer. Anal. , 44 , 879 – 895 (electronic). Google Scholar CrossRef Search ADS Oberman A. M. & Salvador T. ( 2015 ) Filtered schemes for Hamilton–Jacobi equations: a simple construction of convergent accurate difference schemes. J. Comput. Phys. , 284 , 367 – 388 . Google Scholar CrossRef Search ADS Oberman A. M. & Salvador T. ( 2016 ) Numerical method for motion by affine curvature . arXiv preprint arXiv:1610.08831 . Osher S. & Fedkiw R. P. ( 2003 ) Level Set Methods and Dynamic Implicit Surfaces , vol. 153. New York: Springer. Osher S. & Sethian J. A. ( 1988 ) Fronts propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations. J. Comput. Phys. , 79 , 12 – 49 . Google Scholar CrossRef Search ADS Sapiro G. ( 2006 ) Geometric Partial Differential Equations and Image Analysis . Cambridge : Cambridge University Press . Sapiro G. & Tannenbaum A. ( 1994 ) On affine plane curve evolution. J. Funct. Anal. , 119 , 79 – 120 . Google Scholar CrossRef Search ADS Sethian J. A. ( 1999 ) Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science , vol. 3 . Cambridge : Cambridge University Press . Su B. ( 1983 ) Affine Differential Geometry . Beijing : Publisher Science Press ; New York : Gordon & Breach Science Publishers . © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

### Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Aug 19, 2017

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations