Numerical computations of split Bregman method for fourth order total variation flow

TL;DR

This paper develops numerical methods based on the split Bregman framework for fourth order total variation flow, including discretization techniques and a new shrinkage operator, with experiments demonstrating their effectiveness.

Contribution

It introduces a novel numerical approach for fourth order total variation flow using split Bregman, including new discretization and shrinkage operators.

Findings

Effective discretization of the problem using piecewise constant functions

Development of a new shrinkage operator for Spohn's model

Numerical experiments validate the proposed methods

Abstract

The split Bregman framework for Osher-Sol\'e-Vese (OSV) model and fourth order total variation flow are studied. We discretize the problem by piecewise constant function and compute $\nabla(-\Delta_{\mathrm{av}})^{-1}$ approximately and exactly. Furthermore, we provide a new shrinkage operator for Spohn's fourth order model. Numerical experiments are demonstrated for fourth order problems under periodic boundary condition.

Tables

Table 1. Table 1: Time step k 𝑘 k which satisfies ‖ u k ‖ ∞ < 10 − 4 subscript norm superscript u 𝑘 superscript 10 4 \|\textbf{u}^{k}\|_{\infty}<10^{-4} , 10 − 6 superscript 10 6 10^{-6} and 10 − 8 superscript 10 8 10^{-8} .

parameters	$\tau$	$T^{\ast }\left(u^0\right)/\tau ⪅$
$N=100$, $\lambda =h^{-3}$	${10}^{-6}$	17911	4032	41769	135755
$N=100$, $\lambda =10h^{-3}$	${10}^{-7}$	179112	40311	60579	333015
$N=200$, $\lambda =10h^{-3}$	$1.25\times {10}^{-8}$	1432898	322491	592634	1267927

Equations

$$\frac{\partial u}{\partial t}\in -{\partial }_{H}E(u)\text{mbox}fort>0,$$

$${\partial }_{H}E(u)=\left\{p\in H:E(v)-E(u)\ge (p,v-u{)}_H\text{mbox}forallv\in H\right\}.$$

$$u_t=-\Delta \left(\mathrm{div}\left(\beta \frac{\nabla u}{|\nabla u|}+|\nabla u{|}^{p-2}\nabla u\right)\right),$$

$$E(u)=\beta {\int }_{\Omega }|Du|+\frac{1}{p}{\int }_{\Omega }|Du{|}^p.$$

$$u={\mathrm{argmin}}_{u\in L^2(\Omega )}\left\{{\int }_{\Omega }|Du|+\frac{\lambda }{2}\Vert u-f{\Vert }_{L^2(\Omega )}^2\right\},$$

$$u_t=\mathrm{div}\left(\frac{\nabla u}{|\nabla u|}\right)+\lambda (u-f).$$

$$u={\mathrm{argmin}}_{u\in H^{-1}(\Omega )}\left\{{\int }_{\Omega }|Du|+\frac{\lambda }{2}\Vert u-f{\Vert }_{H^{-1}(\Omega )}^2\right\},$$

$$u_t=-\Delta \left(\mathrm{div}\left(\frac{\nabla u}{|\nabla u|}\right)\right)+\lambda (u-f).$$

$$u_t=\mathrm{div}\left(\frac{\nabla u}{|\nabla u|}\right)$$

$$u^{k+1}={\mathrm{argmin}}_{u\in L^2(\Omega )}\left\{{\int }_{\Omega }|Du|+\frac{1}{2\tau }\Vert u-u^k{\Vert }_{L^2(\Omega )}^2\right\},$$

$$u_t=-\Delta \left(\mathrm{div}\left(\frac{\nabla u}{|\nabla u|}\right)\right)$$

$$\frac{\partial u}{\partial t}-|\nabla u|\mathrm{div}\left(\frac{\nabla u}{|\nabla u|}\right)=0.$$

$${\hat{f}}_T(\xi )=⟨f,e^{-2\pi i\xi x}{⟩}_{{\mathcal{D}}^{\prime }(\mathbb{T}),\mathcal{D}(\mathbb{T})}.$$

$${\widehat{df/dx}}_T(\xi )={⟨f,\frac{d}{dx}{e}^{-2\pi i\xi x}⟩}_{{\mathcal{D}}^{\prime }(\mathbb{T}),\mathcal{D}(\mathbb{T})}=-2\pi i\xi {\hat{f}}_T(\xi )$$

$$f(x)=\sum _{\xi \in \mathbb{Z}}{\hat{f}}_T(\xi )e^{2\pi i\xi x}$$

$$L^2(\mathbb{T})=\left\{f\in {\mathcal{D}}^{\prime }(\mathbb{T}):\sum _{\xi \in \mathbb{Z}}|{\hat{f}}_T(\xi ){|}^2<\infty \right\},$$

$$H^1(\mathbb{T})=\left\{f\in L^2(\mathbb{T}):\sum _{\xi \in \mathbb{Z}}{\xi }^2|{\hat{f}}_T(\xi ){|}^2<\infty \right\}=\left\{f\in {\mathcal{D}}^{\prime }(\mathbb{T}):\sum _{\xi \in \mathbb{Z}}(1+{\xi }^2)|{\hat{f}}_T(\xi ){|}^2<\infty \right\},$$

$$H^{-1}(\mathbb{T})=\left\{f\in {\mathcal{D}}^{\prime }(\mathbb{T}):\sum _{\xi \in \mathbb{Z}}(1+{\xi }^2{)}^{-1}|{\hat{f}}_T(\xi ){|}^2<\infty \right\}.$$

$$⟨f,g{⟩}_{H^{-1}(\mathbb{T}),H^1(\mathbb{T})}=\sum _{\xi \in \mathbb{Z}}{\hat{f}}_T(\xi )\overline{{\hat{g}}_T(\xi )}={\int }_{\mathbb{T}}f(x)\overline{g(x)}dx$$

$$L_{\mathrm{av}}^2(\mathbb{T})$$

$$H_{\mathrm{av}}^1(\mathbb{T})$$

$$H_{\mathrm{av}}^{-1}(\mathbb{T})$$

$$L_{\mathrm{av}}^2(\mathbb{T})=\left\{f\in {\mathcal{D}}^{\prime }(\mathbb{T}):\sum _{\xi \ne 0}|{\hat{f}}_T(\xi ){|}^2<\infty \text{mbox}and{\hat{f}}_T(0)=0\right\},$$

$$H_{\mathrm{av}}^1(\mathbb{T})=\left\{f\in {\mathcal{D}}^{\prime }(\mathbb{T}):\sum _{\xi \ne 0}{\xi }^2|{\hat{f}}_T(\xi ){|}^2<\infty \text{mbox}and{\hat{f}}_T(0)=0\right\},$$

$$H_{\mathrm{av}}^{-1}(\mathbb{T})=\left\{f\in {\mathcal{D}}^{\prime }(\mathbb{T}):\sum _{\xi \ne 0}{\xi }^{-2}|{\hat{f}}_T(\xi ){|}^2<\infty \text{mbox}and{\hat{f}}_T(0)=0\right\}.$$

$$(f,g{)}_{L_{\mathrm{av}}^2(\mathbb{T})}$$

$$(f,g{)}_{H_{\mathrm{av}}^1(\mathbb{T})}$$

$$(f,g{)}_{H_{\mathrm{av}}^{-1}(\mathbb{T})}$$

$$\Vert f{\Vert }_{L_{\mathrm{av}}^2(\mathbb{T})}$$

$$\Vert f{\Vert }_{H_{\mathrm{av}}^1(\mathbb{T})}$$

$$-\Delta u(x)=\sum _{\xi \in \mathbb{Z}}{\widehat{(-\Delta u)}}_T(\xi )e^{2\pi i\xi x}=\sum _{\xi \in \mathbb{Z}}4{\pi }^2{\xi }^2{\hat{u}}_T(\xi )e^{2\pi i\xi x}=f(x),$$

$$\Vert f{\Vert }_{H_{\mathrm{av}}^{-1}(\mathbb{T})}^2$$

$$=\sum _{\xi \ne 0}\frac{1}{4{\pi }^2}{\xi }^{-2}{\widehat{(-\Delta u)}}_T(\xi )\overline{{\widehat{(-\Delta u)}}_T(\xi )}$$

Full text

Key words: Total variation flow, Singular diffusion, Osher-Solé-Vese model, surface relaxation, fourth order model, split Bregman,

2010 Mathematics Subject Classification: 35K25, 35K67, 65M06

1 Introduction

A gradient flow has been of great interest in mathematics and mathematical physics because several evolution equations can be regarded as the gradient flows. For example, mathematical models in materials sciences including the Allen-Cahn equation and mean curvature flow can be regarded as second order $L^{2}$-gradient flows. The Cahn-Hilliard equation is interpreted as a fourth-order $H^{-1}$-gradient flow.

We are interested in several important examples of gradient flows which are of the form

[TABLE]

where $H$ is a Hilbert space, $E:H\to\mathbb{R}\cup\{\infty\}$ is a convex, lower semi-continuous functional and the subdifferential $\partial_{H}$ is defined as

[TABLE]

In this paper, we consider gradient flows (1) with convex energy $E$ but may be very singular. We give a few examples. Spohn [36] has proposed a mathematical model for the relaxation of a crystalline surface below the roughening temperature;

[TABLE]

where $\beta>0$ and $p>1$. Kashima [24] has presented this model as a fourth order $H^{-1}$-gradient flow for energy functional

[TABLE]

Furthermore, the total variation flow, which is the gradient flow for total variation energy, has been studied well in image processing. In 1992, Rudin, Osher and Fatemi [35] have introduced the total variation to image denoising and reconstruction. Their model, which is known as the Rudin-Osher-Fatemi (ROF) model, is described as

[TABLE]

where $\Omega\subset\mathbb{R}^{2}$ is bounded domain and $f:\Omega\to\mathbb{R}$ is a given noisy image. This introduces the second order total variation flow

[TABLE]

On the other hand, Osher, Solé and Vese [33] have introduced the $H^{-1}$ fidelity and provided Osher-Solé-Vese (OSV) model

[TABLE]

where $H^{-1}(\Omega)=(H^{1}_{0}(\Omega))^{*}$. Their model performs better on textured or oscillatory images. Equation (7) gives the fourth order total variation flow

[TABLE]

Performing numerical computations for the ROF model, the OSV model and total variation flow have difficulties because of its singularity. Several studies have suggested numerical schemes for the ROF model and second order total variation flow. Especially, the split Bregman framework is well-known as an efficient solver for the ROF model. The aim of this paper is to provide a new numerical scheme, which is based on the backward Euler method and split Bregman framework, for fourth order total variation flow and Spohn’s fourth order model. Numerical experiments are demonstrated for fourth order problems under periodic boundary condition.

The split Bregman method, which is based on the Bregman iterative scheme [2], has been studied and performed in image processing (for example, see [31]). Goldstein and Osher [20] have proposed the alternating split Bregman method. Their method separates the “$L^{1}$” minimization part and “$L^{2}$” part. The alternating split Bregman method has several advantages. They have mentioned that the “$L^{2}$” part is differentiable, and the shrinking method can be applied to the “$L^{1}$” part for the ROF model. Therefore it is extremely efficient solver and easy to code. The split Bregman framework can be performed for second order total variation flow easily. For example,

[TABLE]

introduces the subdifferential formulation $u_{t}\in-\partial F(u)$, where $F(u)=\int_{\Omega}|Du|$. We let $u_{t}\approx(u^{k+1}-u^{k})/\tau$ and apply the backward Euler method to the subdifferential formulation, then we obtain

[TABLE]

where $\tau$ is time step size. This is essentially the same problem as the ROF model, therefore the split Bregman framework can be applied to second order total variation flow.

In this paper, we propose the split Bregman framework for the OSV model (7), fourth order total variation flow

[TABLE]

and Spohn’s fourth order model (3). For simplicity, we consider one-dimensional torus $\mathbb{T}$. We introduce spatial discretization by piecewise constant functions, then we compute $\nabla(-\Delta_{\mathrm{av}})^{-1}v_{h}$ for piecewise constant function $v_{h}$ approximately or exactly. We apply the discrete gradient and discrete inverse Laplacian in our first scheme. In our second scheme, we calculate the inverse Laplacian for piecewise constant functions directly by using the second degree B-spline, which is continuously differentiable piecewise polynomials. The problem can be reduced into a minimization problem on the Euclidean space, which is included in earlier studies for the ROF model. Therefore we can apply the split Bregman framework to fourth order problems. Several theoretical results such as the convergence [5] can be applied to our scheme directly. Both of our two schemes are demonstrated for fourth order problems, and we can check that they perform quite well. Furthermore, we introduce a new shrinkage operator for Spohn’s fourth order model. This enables to perform the numerical experiment for the relaxation of a crystalline surface below the roughening temperature quite effectively. Our scheme can be extended to fourth order problems on the two-dimensional torus. We also suggest a shrinkage operator for two-dimensional Spohn’s model.

Let us quickly overview some earlier results. There are many mathematical studies for the second and fourth order total variation flow. The well-posedness for fourth order total variation flow can be proved by considering the right hand side in (11) as a subdifferential of a total variation in $H^{-1}(\Omega)$ (see [24]). This enables us to use the theory of maximal monotone operators [27, 3]. On the other hand, Elliott and Smitheman [9] have proved the well-posedness for fourth order total variation flow by using the Galerkin method. Adjusting the methods in [15], Giga, Kuroda and Matsuoka [16] have established the extinction time estimate under Dirichlet boundary condition. Numerical computations which include anisotropic diffusion are performed in [29] for second order models. Note that even for the second order total variation flow (9), because of singularity at $\nabla u=0$, the speed of evolution is determined by nonlocal quantities. Therefore the definition of the solution itself is nontrivial. For the second order model, the comparison principle holds, and the theory of viscosity solutions is applicable to show well-posedness for a wide class of total variation type equations [12, 13]. However, for the fourth order problem, the comparison principle does not hold in general (see [11, Theorem 3.7]), and the theory of viscosity solutions is not applicable. For more details of mathematical analysis, we refer the reader to [11] and references therein.

Several studies have considered the fourth order problem under periodic boundary condition. Kashima [25] has studied the characterization of subdifferential in $H_{\mathrm{av}}^{-1}(\mathbb{T}^{d})$. The exact profile of the fourth order total variation flow has been studied in [11]. The extinction time estimate under periodic boundary condition has been established in [15]. A duality based numerical scheme which applies the forward-backward splitting has been proposed in [17]. Kohn and Versieux [26] have performed the numerical computation for Spohn’s model. Their numerical scheme is based on the backward Euler method, mixed finite element method and regularization for singularity. They have proved the convergence by combining the regularization error estimate with temporal and spatial semidiscretization error estimates.

The application of the split Bregman framework to crystalline flow has also been studied through what is called a level-set method. A crystalline mean curvature flow has been proposed independently in [1] and [37]. Oberman, Osher, Takei and Tsai [30] have proposed applying the split Bregman framework to the level-set equation of mean curvature flow. Požár [34] has studied self-similar solutions of three dimensional crystalline mean curvature flow and presented a numerical scheme which is based on the finite element method and split Bregman framework. However, all calculations given there are for the second order model.

A level-set method for mean curvature flow interprets the motion of curvature flow by a level-set of a solution of

[TABLE]

It is a powerful tool to calculate evolution which experiences topological change. It was introduced by Osher and Sathian [32] as a numerical way to calculate the mean curvature flow. Note that the level-set mean curvature equation (12) looks similar to (9). However, the singularity of (12) at $\nabla u=0$ is weaker than one of (9) because of the multiplier $|\nabla u|$. Therefore it is not necessary to study nonlocal quantities for the level-set mean curvature equation. Its analytic foundation such as well-posedness and comparison principle has been established in [8, 10]. For more details, we refer the readers to [14]. Very recently, the analytic foundation of the level-set method is extended to crystalline flow by Požár and the first author [18, 19] and Chambolle, Morini and Ponsiglione [7] and with Novaga [6] independently. Their methods are quite different.

This paper is organized as follows. Section 2 states the definition of $H^{-1}_{\mathrm{av}}(\mathbb{T})$ and the total variation. We introduce the discretization by piecewise constant functions in Section 3. Furthermore, we propose two schemes for reducing $\|\cdot\|_{H^{-1}_{\mathrm{av}}(\mathbb{T})}$ into Euclidean norm. Section 4 presents the split Bregman framework for the OSV model and fourth order total variation flow problem. In Section 5, we describe the shrinking method for Spohn’s model. This report presents numerical examples on one-dimensional torus in Section 6. Finally, we extend our scheme to two-dimensional fourth order problems under periodic boundary condition in Section 7.

2 Preliminary

2.1 Fourier analysis on the torus $\mathbb{T}$

First, we review some of the standard facts on the Fourier analysis for one-dimensional torus $\mathbb{T}=\mathbb{R}/\mathbb{Z}$. In this paper, we regard $\mathbb{T}$ as an interval $[0,1]$ with periodic boundary condition. The Fourier transform for $f\in L^{1}(\mathbb{T})$ and definition of Sobolev space on $\mathbb{R}$ are explained in [21, Chapter 3] and [22, Chapter 1.3], respectively. Let $\mathcal{D}(\mathbb{T})$ be the complex-valued function space $C^{\infty}(\mathbb{T})$ endowed with the usual test function topology and $\mathcal{D}^{\prime}(\mathbb{T})$ be its dual. The Fourier coefficient $\widehat{f}_{T}(\xi)\in\mathbb{C}$ for $f\in\mathcal{D}^{\prime}(\mathbb{T})$ is defined by the generalized Fourier transform (for example, see [23, Chapter 5]);

[TABLE]

The generalized Fourier transform satisfies similar properties to Fourier transform for $f\in L^{1}(\mathbb{T})$, for example,

[TABLE]

for all $f\in\mathcal{D}^{\prime}(\mathbb{T})$. Furthermore, the Fourier coefficients $\widehat{f}_{T}(\xi)\in\mathbb{C}$ satisfies

[TABLE]

for all $f\in\mathcal{D}^{\prime}(\mathbb{T})$ (see [23, Chapter 8.2]). In this Fourier series, the convergence should be understood in the natural topology of $\mathcal{D}^{\prime}(\mathbb{T})$. It is well-known that $\mathcal{D}(\mathbb{T})$ is dense in $L^{2}(\mathbb{T})$, therefore we have $L^{2}(\mathbb{T})\simeq(L^{2}(\mathbb{T}))^{*}\subset\mathcal{D}^{\prime}(\mathbb{T})$, where $(L^{2}(\mathbb{T}))^{*}$ is the dual space of $L^{2}(\mathbb{T})$ (for example, see [4, Chapter 5.2]). Using the generalized Fourier transform, the Lebesgue space and the Sobolev space on $\mathbb{T}$ are defined as follows:

[TABLE]

[TABLE]

[TABLE]

Note that the duality pairing can be described formally as

[TABLE]

for all $f\in H^{-1}(\mathbb{T})$ and $g\in H^{1}(\mathbb{T})$.

2.2 The inverse Laplacian $(-\Delta_{\mathrm{av}})^{-1}$

We consider the functions on $\mathbb{T}$ whose average are equal to zero. Let

[TABLE]

These definitions agree with the following ones:

[TABLE]

[TABLE]

[TABLE]

It is easy to check that each of these spaces are Hilbert space endowed with the inner products

[TABLE]

respectively. These inner products introduce the norms $\|\cdot\|_{L^{2}_{\mathrm{av}}(\mathbb{T})}$, $\|\cdot\|_{H^{1}_{\mathrm{av}}(\mathbb{T})}$ and $\|\cdot\|_{H^{-1}_{\mathrm{av}}(\mathbb{T})}$. It is easy to check that

[TABLE]

Fix $u\in H^{1}_{\mathrm{av}}(\mathbb{T})$ arbitrarily. Let $c(\xi)=4\pi^{2}\xi^{2}\widehat{u}_{T}(\xi)\in\mathbb{C}$ and $f(x)=\sum_{\xi\in\mathbb{Z}}c(\xi)e^{2\pi i\xi x}$, then we have $c(0)=0$ and $\sum_{\xi\neq 0}\xi^{-2}|c(\xi)|^{2}=16\pi^{4}\sum_{\xi\neq 0}\xi^{2}|\widehat{u}_{T}(\xi)|^{2}<\infty$. This implies $f\in H^{-1}_{\mathrm{av}}(\mathbb{T})$. Moreover,

[TABLE]

where $-\Delta u=-d^{2}u/dx^{2}$. Consequently, $u=(-\Delta_{\mathrm{av}})^{-1}f$ defines $(-\Delta_{\mathrm{av}})^{-1}:H^{1}_{\mathrm{av}}(\mathbb{T})\to H^{-1}_{\mathrm{av}}(\mathbb{T})$. We call this operator the inverse Laplacian. Let $f\in H^{-1}_{\mathrm{av}}(\mathbb{T})$ and $u=(-\Delta)_{\mathrm{av}}^{-1}f\in H^{1}_{\mathrm{av}}(\mathbb{T})$, then we have

[TABLE]

This implies

[TABLE]

for all $f\in H^{-1}_{\mathrm{av}}(\mathbb{T})$, where $\nabla=d/dx$.

2.3 Bounded variation and $H^{-1}$ fidelity for the torus $\mathbb{T}$

We recall the spaces of functions of bounded variation in one-dimensional torus.

Definition 1** (Definition 3.3.13 of [21]).**

For a measurable function $f$ on $\mathbb{T}$ which is defined everywhere, we define the total variation as

[TABLE]

where the supremum is taken over all partition of the interval $[0,1]$. we say $f$ is bounded variation if the total variation of $f$ is bounded. Moreover, we define

[TABLE]

Remark**.**

In the definition, $D$ can be regarded as the distributional derivative, and $Df$ can be identified with a signed Borel measure.

Remark**.**

The total variation on a general open set $\Omega\subset\mathbb{R}^{d}$ is defined as

[TABLE]

and the space of bounded variation is defined as

[TABLE]

It is well-known that if $v\in W^{1,1}(\Omega)$, then

[TABLE]

and therefore $W^{1,1}(\Omega)\subset BV(\Omega)\subset L^{1}(\Omega)$.

We define the functional $\Phi:H^{-1}_{\mathrm{av}}(\mathbb{T})\to\mathbb{R}\cup\{\infty\}$ as

[TABLE]

Note that $\Phi:H^{-1}_{\mathrm{av}}(\mathbb{T})\to\mathbb{R}\cup\{\infty\}$ is nonnegative, proper, lower semi-continuous and convex. In this paper, we consider the gradient flow equation of the form

[TABLE]

where the subdifferential $\partial_{H^{-1}_{\mathrm{av}}(\mathbb{T})}$ is defined as

[TABLE]

for any convex functional $F:H^{-1}_{\mathrm{av}}(\mathbb{T})\to\mathbb{R}\cup\{\infty\}$ and $u\in H^{-1}_{\mathrm{av}}(\mathbb{T})$. It is well-known that the theory of maximal monotone operators shows the existence and uniqueness of solution $u\in C([0,\infty),H^{-1}_{\mathrm{av}}(\mathbb{T}))$ to equation (35) (for example, see [27]).

Let $\tau>0$ be the temporal step size. We consider the backward Euler method for gradient flow equation (35); for given $u^{k}\in H^{-1}_{\mathrm{av}}(\mathbb{T})$, find $u^{k+1}\in H^{-1}_{\mathrm{av}}(\mathbb{T})$ such that

[TABLE]

This can be reduced to solving the following minimization problem:

[TABLE]

Since $\Phi$ is convex, such $u^{k+1}$ is uniquely determined. The convergence of backward Euler method has been proved in [27]. Note that equation (38) is similar to the OSV model [33] which can be described as

[TABLE]

where $f\in H^{-1}_{\mathrm{av}}(\mathbb{T})$ is given data and $\lambda>0$ is an artificial parameter. The existence result in convex analysis (for example, see [4, Cor3.23]) gives that the minimizer $u\in BV(\mathbb{T})\cap H^{-1}_{\mathrm{av}}(\mathbb{T})$ exists.

Hereafter, we consider the following minimization problem: find $u\in H^{-1}_{\mathrm{av}}(\mathbb{T})$ such that

[TABLE]

where $f\in H^{-1}_{\mathrm{av}}(\mathbb{T})$ is a given data or $f=u^{k}$, and $\lambda$ is a given parameter or $\lambda=1/\tau$. This involves both of (OSV) and the backward Euler method for (gradient flow). Furthermore, $(\mbox{P}0)$ introduces the following constrained problem:

[TABLE]

Remark**.**

When we consider the Spohn’s model

[TABLE]

the subdifferential formulation is given as $u_{t}\in-\partial_{H^{-1}_{\mathrm{av}}(\mathbb{T})}\widetilde{\Phi}(u)$, where

[TABLE]

Therefore the backward Euler method yields

[TABLE]

Then we consider the constraint problem

[TABLE]

3 Discretization for total variation flow and OSV model

3.1 Discretization for minimization problem

We introduce the (spatial) discretization for the problem $(\mbox{P}1)$. Let $N\in\mathbb{N}$ be the partition number, $h=1/N$ and $x_{n}=nh$. We regard $x_{0}=x_{N}$, then $\{x_{n}\}_{n=0}^{N}$ gives an uniform partition for $\mathbb{T}$. Furthermore, we let $x_{n+1/2}=(x_{n}+x_{n+1})/2=(n+1/2)h$ for $n=-1,0,\dots,N$, where $x_{-1/2}$ and $x_{N+1/2}$ are identified with $x_{N-1/2}$ and $x_{1/2}$, respectively. Then we define the following spaces of piecewise constant functions:

[TABLE]

where $I=[0,1]$, $I_{n}=[x_{n-1/2},x_{n+1/2})$, $\mathbb{P}_{0}(I_{n})$ is a space of constant functions on interval $I_{n}$ and $\boldsymbol{1}_{I_{n}}$ is its characteristic function. Note that $V_{h0}$ is a finite dimensional subspace of $H^{-1}_{\mathrm{av}}(\mathbb{T})$.

Furthermore, we define $D_{h}:V_{h0}\to\widehat{V}_{h}\cap L^{2}_{\mathrm{av}}(I)$ as

[TABLE]

where $v_{0}$ is identified with $v_{N}$. Let $d_{h}=D_{h}v_{h}\in\widehat{V}_{h}$, $\textbf{d}=(d_{1},\dots,d_{N})^{\mathrm{T}}\in\mathbb{R}^{N}$ for $d_{h}=\sum_{n=1}^{N}d_{n}\boldsymbol{1}_{[x_{n-1},x_{n})}$, $\widetilde{\textbf{v}}=(v_{1},v_{2},\dots,v_{N})^{\mathrm{T}}\in\mathbb{R}^{N}$ for $v_{h}=\sum_{n=1}^{N}v_{n}\boldsymbol{1}_{I_{n}}\in V_{h0}$, then we have

[TABLE]

where $\nabla_{h}:\mathbb{R}^{N}\to\mathbb{R}^{N}$ is the discrete gradient

[TABLE]

Note that $D_{h}v_{h}\in\widehat{V}_{h}\cap L^{2}_{\mathrm{av}}(I)\subset L^{2}_{\mathrm{av}}(I)$ for all $v_{h}\in V_{h0}$; however, $D_{h}v_{h}\not\in H^{-1}_{\mathrm{av}}(\mathbb{T})$ because it does not satisfy the periodic boundary condition (see Figure 1).

Here we introduce the discretized problem for (P$1$);

[TABLE]

where $f_{h}\in V_{h0}$ is given data or $f_{h}=u_{h}^{k}$, and $\textbf{d}=(d_{1},\dots,d_{N})^{\mathrm{T}}$ for $d_{h}=\sum_{n=1}^{N}d_{n}\boldsymbol{1}_{[x_{n-1},x_{n})}$. Furthermore, we introduce the unconstrained problem

[TABLE]

Remark**.**

In this paper, we use $\|d_{h}-D_{h}u_{h}\|_{L^{2}(I)}^{2}$. This enables to apply the shrinking method to minimization problem in the split Bregman framework.

3.2 Corresponding matrix form

We reduce $(\mathrm{P}2)_{h}$ to the matrix formulation. Let $\textbf{d}=(d_{1},\dots,d_{N})^{\mathrm{T}}\in\mathbb{R}^{N}$ for $d_{h}=\sum_{n=1}^{N}d_{n}\boldsymbol{1}_{[x_{n-1},x_{n})}$, $\widetilde{\textbf{u}}=(u_{1},\dots,u_{N})^{\mathrm{T}}\in\mathbb{R}^{N}$ and $\textbf{u}=(u_{1},\dots,u_{N-1})^{\mathrm{T}}\in\mathbb{R}^{N-1}$ for $u_{h}=\sum_{n=1}^{N}u_{n}\boldsymbol{1}_{I_{n}}\in V_{h0}$, then we have

[TABLE]

where $S_{N}=h\nabla_{h}\in\mathbb{R}^{N\times N}$. Furthermore, $u_{h}\in V_{h0}$ implies $u_{N}=-\sum_{n=1}^{N-1}u_{n}$, that is, $\widetilde{\textbf{u}}=R_{N}\textbf{u}$, where

[TABLE]

Therefore

[TABLE]

Next, we consider two expressions of $\|v_{h}\|_{H^{-1}_{\mathrm{av}}(\mathbb{T})}^{2}$ for $v_{h}\in V_{h0}$. Recall that equation (28) implies

[TABLE]

We propose two schemes for considering $\nabla(-\Delta_{\mathrm{av}})^{-1}$. The first scheme is to approximate $\nabla(-\Delta_{\mathrm{av}})^{-1}$ by using the discrete gradient $\nabla_{h}\in\mathbb{R}^{N\times N}$ and the discrete Laplacian

[TABLE]

Let $\widetilde{\textbf{v}}=(v_{1},v_{2},\dots,v_{N})^{\mathrm{T}}\in\mathbb{R}^{N}$ and $\textbf{v}=(v_{1},\dots,v_{N-1})^{\mathrm{T}}\in\mathbb{R}^{N-1}$ for $v_{h}\in V_{h0}$, then $\widetilde{\textbf{v}}=R_{N}\textbf{v}$. We define $\textbf{w}\in\mathbb{R}^{N-1}$ and $\widetilde{\textbf{w}}\in\mathbb{R}^{N}$ for $w_{h}\in V_{h0}$ in the same way. Letting $\widetilde{\textbf{v}}=-\Delta_{h}\widetilde{\textbf{w}}$ implies

[TABLE]

Multiplying the (unique) pseudo-inverse matrix

[TABLE]

yields $\textbf{v}=L_{N}(-\Delta_{h})R_{N}\textbf{w}=h^{-2}L_{N}S_{N}^{\mathrm{T}}S_{N}R_{N}\textbf{w}$. For simplicity of notation, we let

[TABLE]

It is easy to check that

[TABLE]

satisfies $\det A_{N}=N^{2}\neq 0$, therefore we have $\det(-\Delta_{\mathrm{av}})_{h}^{-1}\neq 0$. This implies

[TABLE]

Our first scheme is to approximate $(-\Delta_{\mathrm{av}})^{-1}$ by $R_{N}(-\Delta_{\mathrm{av}})_{h}^{-1}$, instead of $(-\Delta_{h})^{-1}$ which does not exist. This yields

[TABLE]

that is,

[TABLE]

For simplicity, we let $J=SR_{N}A_{N}^{-1}\in\mathbb{R}^{N\times(N-1)}$, then our first scheme can be described as

[TABLE]

for all $v_{h}\in V_{h0}$.

Remark**.**

When we apply $H^{-s}_{\mathrm{av}}(\mathbb{T})$ norm $(0<s<1)$ to $(\mbox{P}1)_{h}$, the discrete inverse Laplacian $(-\Delta_{\mathrm{av}})_{h}^{-s}$ can be introduced by the discrete Fourier transform (for example, see [17]).

Our second scheme is to compute $\nabla(-\Delta_{\mathrm{av}})^{-1}$ directly. It requires the second degree piecewise polynomial which has continuous derivative. We define the second degree periodic B-spline basis functions (see Figure 2)

[TABLE]

We identify $B_{-1}\equiv B_{N-1}$, $B_{0}\equiv B_{N}$ and $B_{1}\equiv B_{N+1}$. The B-spline basis functions have continuous derivative (see Figure 3)

[TABLE]

Therefore we have

[TABLE]

Fix $v_{h}\in V_{h0}$ arbitrarily, then there exits $w_{h}\in\operatorname{span}\{B_{1},\dots,B_{N}\}$ such that $w_{h}=(-\Delta_{\mathrm{av}}^{-1})v_{h}\in H^{1}_{\mathrm{av}}(\mathbb{T})$. It is easy to check that

[TABLE]

Let $\sum_{n=1}^{N}w_{n}=0$, then equation (67) implies $\sum_{n=1}^{N}w_{n}B_{n}\in H^{1}_{\mathrm{av}}(\mathbb{T})$. Furthermore, we let

[TABLE]

Then we have

[TABLE]

This implies

[TABLE]

Multiplying the pseudo-inverse matrix $L_{N}$ yields

[TABLE]

Therefore we have

[TABLE]

The definition, combining with equation (28) gives $\|v_{h}\|_{H^{-1}_{\mathrm{av}}(\mathbb{T})}^{2}=\|\nabla w_{h}\|_{L^{2}(\mathbb{T})}^{2}$, where

[TABLE]

is a piecewise linear function which satisfies

[TABLE]

This implies

[TABLE]

where

[TABLE]

We identify $\phi_{-1/2}=\phi_{N-1/2}$, $\phi_{1/2}=\phi_{N+1/2}$ (see Figure 4). It is easy to check that

[TABLE]

for all $n=1,\dots,N$. Therefore we have

[TABLE]

where

[TABLE]

Let

[TABLE]

where $a=\dfrac{\sqrt{3}+1}{2\sqrt{3}}$ and $b=\dfrac{\sqrt{3}-1}{2\sqrt{3}}$, then $T^{\mathrm{T}}T=M_{N}$. Summarizing the above argument, our second scheme can be described as

[TABLE]

Let $H=TSR_{N}A_{N}^{-1}=TJ\in\mathbb{R}^{N\times(N-1)}$ for simplicity of notation, then we have

[TABLE]

Applying equation (54), (63) and (80) to $(\mbox{P}2)_{h}$ implies the following two discretized problems;

[TABLE]

where $\textbf{f}\in\mathbb{R}^{N-1}$ is given as $f\in V_{h0}$ or $\textbf{f}=\textbf{u}^{k}$. Recall that the matrix $J$ is introduced by the approximation $\nabla(-\Delta_{\mathrm{av}})^{-1}\approx\nabla_{h}R_{N}(-\Delta_{\mathrm{av}})_{h}^{-1}$. On the other hand, we obtain $H$ by using $\nabla(-\Delta_{\mathrm{av}})^{-1}$ exactly. Therefore (81a) can be regarded as an approximation of (81b), which is equivalent to $(\mbox{P}2)_{h}$.

4 Split Bregman framework

In this section, we review the alternating split Bregman framework in [20] for the problem

[TABLE]

where $K\in\mathbb{R}^{N\times(N-1)}$ is equal to $J$ or $H$. Recall that $(\mbox{P}3K)_{h}$ is an approximation of the discrete problem for $(\mbox{P}0)$;

[TABLE]

Let

[TABLE]

The Bregman method replaces $\Psi(\textbf{u},\textbf{d})$ into the Bregman distance and iteratively solves

[TABLE]

where the Bregman distance $D_{\Psi}^{\textbf{p}^{k}}$ is defined as

[TABLE]

and $\textbf{p}^{k}=(\textbf{p}_{u}^{k},\textbf{p}_{d}^{k})\in\mathbb{R}^{N-1}\times\mathbb{R}^{N}$ is defined as

[TABLE]

Thanks to $\Psi:\mathbb{R}^{N-1}\times\mathbb{R}^{N}\to\mathbb{R}$ is convex and lower semi-continuous, the Bregman distance $D_{\Psi}^{\textbf{p}^{k}}(\cdot,(\textbf{u}^{k},\textbf{d}^{k}))$ is also convex and lower semi-continuous. Applying the usual existence result of convex analysis (see [4, Cor3.23]) gives that there exists a minimizer $(\textbf{u}^{k+1},\textbf{d}^{k+1})$. Furthermore, by using induction we can show that

[TABLE]

is independent of $(\textbf{u},\textbf{d})$, where $\boldsymbol{\alpha}^{k+1}\in\mathbb{R}^{N}$ is defined as

[TABLE]

This implies the minimizer $(\textbf{u}^{k+1},\textbf{d}^{k+1})$ of problem (85) satisfies

[TABLE]

This is the split Bregman iteration for the problem $(\mbox{P}3K)_{h}$. Finally, we apply the alternating split Bregman algorithm and obtain

[TABLE]

where $\textbf{f}\in\mathbb{R}^{N-1}$ is given data or $\textbf{f}=\textbf{u}^{k}$, $\boldsymbol{\alpha}^{0}=\textbf{0}$, $\textbf{u}^{0}$ is given as 0 or initial condition, and $\textbf{d}^{0}=S_{N}R_{N}\textbf{u}^{0}$. This satisfies the following convergence result.

Lemma 1** (Theorem 3.2 of [5]).**

Suppose that $(\mbox{P}0K)_{h}$ has a minimizer $\textbf{u}^{*}\in\mathbb{R}^{N-1}$, then $\textbf{u}^{k}$ which determined by $(\mbox{P}4K)_{h}$ satisfies

[TABLE]

Furthermore, if a minimizer $\textbf{u}^{*}$ of $(\mbox{P}0K)_{h}$ is unique, then $\lim_{k\to\infty}\|\textbf{u}^{k+1}-\textbf{u}^{*}\|_{2}=0$.

The functional in (91a) is differentiable with respect to u, and the minimization (91b) can be reduced to the shrinking method

[TABLE]

where $(\textbf{v})_{n}$ is the $n$-entry of vector v and

[TABLE]

Therefore, the problem $(\mbox{P}4K)_{h}$ introduces

[TABLE]

where $\textbf{f}\in\mathbb{R}^{N-1}$ is given data or $\textbf{f}=\textbf{u}^{k}$, $\boldsymbol{\alpha}^{0}=\textbf{0}$, $\textbf{u}^{0}$ is given as 0 or initial condition, and $\textbf{d}^{0}=S_{N}R_{N}\textbf{u}^{0}$.

5 Shrinking method for Spohn’s model

In this section, we consider the split Bregman framework for Spohn’s model

[TABLE]

which can be regarded as the gradient flow problem for energy functional

[TABLE]

where $\beta>0$ and $p>1$. This energy is considered in model for the relaxation of a crystalline surface below the roughening temperature (for example, see [26]). If we replace $w=(Du)^{p}$, the alternating split Bregman method introduces nonlinear problem. In this paper we always assume $p=3$, and we apply the constraint $d=Du$ to $|Du|^{3}$. The alternating split Bregman method implies

[TABLE]

We consider the Euler-Lagrange equation for equation (97b);

[TABLE]

For simplicity of notation, we let $x=(\textbf{d}^{k+1})_{n}$, $a=1/(\mu h)>0$ and $\rho=(SR_{N}\textbf{u}^{k+1}+\boldsymbol{\alpha}^{k})_{n}$. This, combining with $p=3$ gives

[TABLE]

Suppose that $x>0$, then we have $a\beta<\rho$ and

[TABLE]

By the similar way, supposing $x<0$ yields $\rho<-a\beta$ and

[TABLE]

If $-a\beta<\rho<a\beta$, we let $x=0$. These observations provide the shrinking operator of the form

[TABLE]

Applying this to equation (97b) gives

[TABLE]

where $\rho_{n}^{k+1}=(S_{N}R_{N}\textbf{u}^{k+1}+\boldsymbol{\alpha}^{k})_{n}$.

6 Numerical example

6.1 Example 1: Comparison of two schemes

Here we show numerical examples of $(\mbox{P}5K)_{h}$. Note that equation (91b) implies that $\mu$ should satisfy $\mu=O(h^{-1})$. Moreover, this and equation (91a) yield $\lambda=O(h^{-3})$ is necessary for reasonable computation.

In this paper, we always regard $\mathbb{T}$ as an interval $[0,1]$ with periodic boundary condition. Our first numerical example is the gradient flow (35) with the initial condition

[TABLE]

Note that the similar example is computed in [17]. They essentially apply the matrix $J$ and compute the gradient flow problem without split Bregman method. Their scheme requires $\tau=\lambda^{-1}=O(h^{5})$ for $H^{-1}_{\mathrm{av}}$ fidelity.

We check the difference between $K=J$ and $K=H$. Figure 5 shows two numerical results with the same parameters $N=40$, $\lambda=h^{-3}$ and $\mu=5h^{-1}$. Numerical results $\textbf{u}^{k}\in\mathbb{R}^{N-1}$ are represented as piecewise constant functions $u_{h}^{k}\in V_{h0}$. They are different because the matrix $J$ is introduced by discrete gradient and discrete inverse Laplacian. This difference is expected to be small if we consider sufficiently small $h$. Figure 6 shows evolution of numerical solutions for $N=200$, $\lambda=h^{-3}$ and $\mu=5h^{-1}$. We infer from them that (81a) can provide sufficiently accurate result.

6.2 Example 2: Discontinuity and symmetry

Our second numerical example for (35) is

[TABLE]

where $a=450$ and $r=1/15$. In [11], a class of initial data including (104) as an example has been studied analytically. They rigorously proved that the solution becomes discontinuous instantaneously. Their analysis gives an exact profile of the fourth order gradient flow. Note that because of uniqueness of a solution, the symmetry of initial profile is preserved during evolution. We can check that our numerical result shows the discontinuity and symmetry approximately (see Figure 7a, 7b). We use $K=J$, $N=200$, $\lambda=25h^{-3}$ and $\mu=15h^{-1}$. Furthermore, we note that we can compute until $\textbf{u}^{k}\approx\textbf{0}$ easily, because our scheme can be stable for $\tau=\lambda^{-1}=O(h^{3})$.

6.3 Example 3: Extinction time

Our third example for (35) is

[TABLE]

which gives

[TABLE]

Figure 8a shows evolution of numerical solution for third example. We use $N-200$, $\lambda=20h^{-3}$ and $\mu=30h^{-1}$ for .Figure 8a. Recall that our numerical scheme can compute the evolution until $\textbf{u}^{k}\approx\textbf{0}$. easily. Furthermore, applying the extinction time estimate [15, Theorem 3.11] to one-dimensional torus implies

[TABLE]

where $T^{*}(u^{0})$ is the extinction time for the initial condition $u^{0}\in H^{-1}_{\mathrm{av}}(\mathbb{T})$ and the constant $C^{*}$ satisfies $\|f\|_{H^{-1}_{\mathrm{av}}(\mathbb{T})}\leq C^{*}\int_{\mathbb{T}}|Df|$ for all $f\in H^{-1}_{\mathrm{av}}(\mathbb{T})$. It is easy to check that

[TABLE]

for all $f\in H^{-1}_{\mathrm{av}}(\mathbb{T})$. Therefore, the extinction time for $u^{0}(x)=-\cos(2\pi x)$ can be estimated as

[TABLE]

The numerical solution is expected to be “extinct” in

[TABLE]

Table 1 shows the time step number $k$ such that $\|\textbf{u}^{k}\|_{\infty}<10^{-4}$, $10^{-6}$ and $10^{-8}$ for each parameters. This result shows that we can get $\|\textbf{u}^{k}\|_{\infty}\lessapprox\tau$ in reasonable iteration number which is expected in (110), however, it requires more iteration to obtain smaller $\|\textbf{u}^{k}\|_{\infty}$.

6.4 Example 4: Spohn’s model

Our fourth example is split Bregman framework for Spohn’s fourth order model (4), which is described in Section 5. Recall that we suppose that $p=3$ in this paper. Therefore we can apply the shrinkage operator (102) to split Bregman framework for Spohn’s model. Figure 8b shows the numerical example for $u^{0}(x)=-\cos(2\pi x)$, $\beta=0.5$, $N=200$, $\lambda=50h^{-3}$ and $\mu=30h^{-1}$.

7 Two dimensional case

The fourth order total variation flow and Spohn’s model on two dimensional torus $\mathbb{T}^{2}$ can be computed by the similar way to one dimensional case. We can define $L^{2}_{\mathrm{av}}(\mathbb{T}^{2})$, $H^{1}_{\mathrm{av}}(\mathbb{T}^{2})$, $H^{-1}_{\mathrm{av}}(\mathbb{T}^{2})$ and $(-\Delta_{\mathrm{av}})^{-1}:H^{-1}_{\mathrm{av}}(\mathbb{T}^{2})\to H^{1}_{\mathrm{av}}(\mathbb{T}^{2})$ by the generalized Fourier transform.

First, the fourth order isotropic total variation flow introduces the constraint problem

[TABLE]

where $D_{x}$, $D_{y}$ is distributional derivative for each variable. Note that

[TABLE]

where $\nabla_{x}=\partial/\partial x$ and $\nabla_{y}=\partial/\partial y$. Let $N_{x}$, $N_{y}$ be the partition number, $h_{x}=1/N_{x}$, $h_{y}=1/N_{y}$, $x_{n}=nh_{x}$ and $y_{n}=nh_{y}$. Furthermore, we let $Q_{n_{x},n_{y}}=[x_{n_{x}-1/2},x_{n_{x}+1/2})\times[y_{n_{y}-1/2},y_{n_{y}+1/2})$ and $\widehat{Q}_{n_{x},n_{y}}=[x_{n_{x}-1},x_{n_{x}})\times[y_{n_{y}-1},y_{n_{y}})$. Then we consider the space of piecewise constant functions

[TABLE]

where $\Omega=[0,1)^{2}$. Any element $d_{h}\in\widehat{V}_{h}$ is described as $d_{h}=\sum_{n_{x},n_{y}}d_{n_{x},n_{y}}\boldsymbol{1}_{\widehat{Q}_{n_{x},n_{y}}}$. Let

[TABLE]

for $v_{h}\in V_{h0}$ and $d_{h}\in\widehat{V}_{h}$. We define $D_{xh},D_{yh}:V_{h0}\to\widehat{V}_{h}\cap L^{2}_{\mathrm{av}}(\Omega)$ as

[TABLE]

This gives

[TABLE]

where $R_{N_{x}N_{y}}\in\mathbb{R}^{(N_{x}N_{y})\times(N_{x}N_{y}-1)}$ is defined as equation (53) and $\nabla_{xh}$, $\nabla_{yh}$ are the discrete gradient

[TABLE]

where $I_{N}\in\mathbb{R}^{N\times N}$ is the identity matrix and $\otimes$ is the Kronecker product. Then our discretized problem is described as

[TABLE]

where $\textbf{d}_{xy}\in\mathbb{R}^{N_{x}\times N_{y}}$ is defined as

[TABLE]

and $K_{x},K_{y}\in\mathbb{R}^{(N_{x}N_{y})\times(N_{x}N_{y}-1)}$ are deduced from $\nabla_{x}(-\Delta_{\mathrm{av}})^{-1}$ and $\nabla_{y}(-\Delta_{\mathrm{av}})^{-1}$, respectively. For example, we can approximate the inverse Laplacian by using

[TABLE]

This yields that our first scheme for two dimensional case is described as $K_{x}=J_{x}$ and $K_{y}=J_{y}$, where

[TABLE]

If we let $h_{x}=h_{y}=h$, then it is required that $\lambda=O(h^{-4})$ and $\tau=O(h^{-2})$. The split Bregman framework gives

[TABLE]

where $\textbf{f}\in\mathbb{R}^{N_{x}N_{y}-1}$ is given data or $\textbf{f}=\textbf{u}^{k}$, $\boldsymbol{\alpha}_{x}^{0}=\boldsymbol{\alpha}_{x}^{0}=\textbf{0}$, $\textbf{u}^{0}$ is given as 0 or initial condition, and $\textbf{d}_{x}^{0}=h_{x}\nabla_{xh}R_{N_{x}N_{y}}\textbf{u}^{0}$, $\textbf{d}_{y}^{0}=h_{y}\nabla_{yh}R_{N_{x}N_{y}}\textbf{u}^{0}$.

Note that the equation (121b) is essentially the same formulation as the one of split Bregman framework for second order isotropic problem, which is mentioned in [20]. The Euler-Lagrange equation for equation (121b) yields

[TABLE]

for all $n=1,\dots,N_{x}N_{y}$. We consider the approximation

[TABLE]

where

[TABLE]

Applying them into equations (122a) and (122b) give the following shrinkage formula, which are equivalent to ones of [20, Section 4.1]:

[TABLE]

Figure 9a shows the numerical result of fourth order isotropic total variation flow (11) in $\mathbb{T}^{2}$ with initial data $u^{0}(x,y)=x(x-1)y(y-1)-1/36$. We use $N_{x}=N_{y}=40$, $\lambda=5h^{-4}$ and $\mu=20h^{-2}$.

Next, the fourth order anisotropic total variation flow

[TABLE]

Letting $F(u)=\int_{\mathbb{T}^{2}}\left(|D_{x}u|+|D_{y}u|\right)$ implies that formally we have

[TABLE]

therefore $u_{t}\in-\partial_{H^{-1}_{\mathrm{av}}(\mathbb{T}^{2})}F$. We apply the backward Euler method and obtain

[TABLE]

which introduces the constraint problem

[TABLE]

This, combining with the split Bregman framework gives

[TABLE]

We can apply the shrinking method (94) to equations (127b) and (127c). Figure 9b presents the evolution of fourth order anisotropic total variation flow for $u^{0}(x,y)=x(x-1)y(y-1)-1/36$, $N_{x}=N_{y}=40$, $\lambda=5h^{-4}$ and $\mu=20h^{-2}$. For second order anisotropic total variation flow, Łasica, Moll and Mucha [28] have considered rectangular domain $\Omega\subset\mathbb{R}^{2}$ or $\Omega=\mathbb{R}^{2}$ and rigorously proved that if the initial profile is piecewise constant, then the exact solution is piecewise constant. We can infer from our numerical experiment 9b that their theoretical result is true also for fourth order anisotropic total variation flow.

Finally, we consider two dimensional Spohn’s fourth order model. The split Bregman framework provides

[TABLE]

The Euler-Lagrange equation for (128b) can be approximated by equation (123). In this paper, we always suppose that $p=3$. Note that the approximation (123) implies

[TABLE]

We obtain approximated Euler-Lagrange equations

[TABLE]

By the similar way to one dimensional case, we provide the shrinkage operators of the form

[TABLE]

Figure 9c shows the numerical result of split Bregman framework for Spohn’s forth order model. We use $p=3$, $\beta=0.25$, $N_{x}=N_{y}=40$, $\lambda=1.25h^{-4}$ and $\mu=5h^{-2}$. Moreover, we use the initial value $u^{0}(x,y)=x(x-1)y(y-1)-1/36$, which is considered in [26]. We can obtain the similar numerical result quite effectively by split Bregman framework.

8 Conclusion

In this study, we propose a new numerical scheme for the OSV model, fourth order total variation flow and Spohn’s fourth order model. Our scheme is based on the split Bregman framework for the ROF model and second order total variation flow. We demonstrate several numerical examples for one dimensional and two dimensional problems under periodic boundary condition. We use the parameters $\lambda=O(h^{-3})$, $\mu=O(h^{-1})$ for one dimensional case, and $\lambda=O(h^{-4})$, $\mu=O(h^{-2})$ for two dimensional case. For fourth order total variation flow, our numerical results approximately represent the flat facet and discontinuity, which is expected by the theoretical result for the exact profile. Furthermore, we propose new shrinkage operators for Spohn’s model. Numerical results for Spohn’s model show facet and relaxation.

Acknowledgement

A part of the work of the second author was done when he was a postdoc fellow at the University of Tokyo. Its hospitality is gratefully acknowledged. The work of the first author was partly supported by the Japan Society for the Promotion of Science through the grant No. 26220702 (Kiban S), No. 19H00639 (Kiban A), No. 18H05323 (Kaitaku), No. 17H01091 (Kiban A) and No. 16H03948 (Kiban B).