Spectral approximation of fractional PDEs in image processing and phase field modeling

TL;DR

This paper develops a spectral approximation method for fractional PDEs, demonstrating its effectiveness in image processing and phase field modeling through theoretical analysis and numerical experiments.

Contribution

It introduces a spectral numerical scheme for fractional PDEs and analyzes its efficiency, specifically applied to image processing and phase field models.

Findings

Spectral method effectively approximates fractional PDEs.

Numerical experiments confirm the method's efficiency.

Applications include image processing and phase field modeling.

Abstract

Fractional differential operators provide an attractive mathematical tool to model effects with limited regularity properties. Particular examples are image processing and phase field models in which jumps across lower dimensional subsets and sharp transitions across interfaces are of interest. The numerical solution of corresponding model problems via a spectral method is analyzed. Its efficiency and features of the model problems are illustrated by numerical experiments.

Equations

$$(-\Delta {)}^s\equiv (-{\Delta }_{{\mathbb{T}}^d}{)}^s$$

$$E(u)=|Du|({\mathbb{T}}^d)+\frac{\alpha }{2}\Vert u-g{\Vert }^2.$$

$$E(u)=\frac{1}{2}\Vert (-\Delta {)}^{s/2}u{\Vert }^2+\frac{\alpha }{2}\Vert (-\Delta {)}^{-\beta /2}(u-g){\Vert }^2,$$

$$(-\Delta {)}^{s}u+\alpha (-\Delta {)}^{-\beta }(u-g)=0\text{mbox}in{\mathbb{T}}^d,$$

$$BV({\mathbb{T}}^d)\cap L^{\infty }({\mathbb{T}}^d)\subset B_{2,\infty }^{1/2}({\mathbb{T}}^d)$$

$$B_{2,\infty }^{1/2}({\mathbb{T}}^d)\subset H^s({\mathbb{T}}^d)$$

$$BV({\mathbb{T}}^d)\cap L^{\infty }({\mathbb{T}}^d)\subset H^s({\mathbb{T}}^d),$$

$$E_{\epsilon }(u)=\frac{1}{2}\Vert (-\Delta {)}^{s/2}u{\Vert }^2+{\epsilon }^{-2}{\int }_{{\mathbb{T}}^d}F(u)dx$$

$$\begin{array}{rl}{\partial }_{t}u+(-\Delta {)}^{s}u+{\epsilon }^{-2}f(u)& =0\text{mbox}in(0,T)\times {\mathbb{T}}^d,\\ u(0,\cdot )& =u_0\text{mbox}in{\mathbb{T}}^d,\end{array}$$

$$\varepsilon^{-2}\mathrel{\mathop{\mathchar 58\relax}}=\left\{\begin{array}[]{ll}\widetilde{\varepsilon}^{-2s}&\mbox{if }s\in(0,1/2),\\ |\log\widetilde{\varepsilon}|&\mbox{if }s=1/2,\\ \widetilde{\varepsilon}^{1-2s}&\mbox{if }s\in(1/2,1)\end{array}\right.$$

$$\begin{array}{rl}(-\Delta {)}^{-\alpha }{\partial }_{t}u+(-\Delta {)}^{s}u+{\epsilon }^{-2}f(u)& =0\text{mbox}in(0,T)\times {\mathbb{T}}^d,\\ u(0,\cdot )& =u_0\text{mbox}in{\mathbb{T}}^d,\end{array}$$

$$\Phi^{k}=\big{(}e^{{\mathrm{i}k\cdot x}_{j}}\mathrel{\mathop{\mathchar 58\relax}}j\in\mathbb{N}_{n}^{d}\big{)},$$

$$(V,W{)}_n=\frac{(2\pi {)}^d}{n^d}\sum _{j\in {\mathbb{N}}_n^d}{v}_j{\overline{w}}_j.$$

$${\tilde{v}}_k=(V,{\Phi }^k{)}_n.$$

$$\mathcal{S}_{n}=\big{\{}v_{n}\in C(\mathbb{T}^{d};\mathbb{C})\mathrel{\mathop{\mathchar 58\relax}}v_{n}(x)=\sum_{k\in\mathbb{Z}_{n}^{d}}c_{k}\varphi^{k}(x),\,c_{k}\in\mathbb{C}\big{\}},$$

$$(v_n,w_n)=(V,W{)}_n.$$

$$I_{n}v=\frac{1}{(2\pi {)}^d}\sum _{k\in {\mathbb{Z}}_n^d}{\tilde{v}}_k{\phi }^k.$$

$${\hat{v}}_k=(v,{\phi }^k).$$

$$(P_{n}v,w_n)=(v,w_n).$$

$$|v{|}_{\mu }^2=\sum _{k\in {\mathbb{Z}}^d}|k{|}^{2\mu }|{\hat{v}}_k{|}^2<\infty .$$

$$|\psi {|}_{-\mu }^2=\sum _{k\in {\mathbb{Z}}^d\setminus \{0\}}|k{|}^{-2\mu }|{\hat{\psi }}_k{|}^2<\infty ,$$

$$|v-P_{n}v|_{\lambda}\leq\Big{(}\frac{n}{2}\Big{)}^{-(\mu-\lambda)}|v|_{\mu}.$$

$$|v-I_{n}v|_{\lambda}\leq c_{d,\lambda,\mu}\Big{(}\frac{n}{2}\Big{)}^{-(\mu-\lambda)}|v|_{\mu}$$

$$|v_{n}|_{r}\leq\max_{{k\in\mathbb{Z}_{n}^{d}}}|k|^{r-s}|v_{n}|_{s}\leq\Big{(}\frac{n}{2}\Big{)}^{r-s}|v_{n}|_{s}.$$

$$\accentset{\circ}{H}^{r}(\mathbb{T}^{d};\mathbb{C})=\big{\{}v\in H^{r}(\mathbb{T}^{d};\mathbb{C})\mathrel{\mathop{\mathchar 58\relax}}\widehat{v}_{0}=0\big{\}}.$$

$$(-\Delta {)}^{s}v=\frac{1}{(2\pi {)}^d}\sum _{k\in {\mathbb{Z}}^d\setminus \{0\}}|k{|}^{2s}{\hat{v}}_k{\phi }^k.$$

$$(-\Delta {)}^{s}u=f.$$

$$u=\frac{1}{(2\pi {)}^d}\sum _{k\in {\mathbb{Z}}^d\setminus \{0\}}|k{|}^{-2s}{\hat{f}}_k{\phi }^k,$$

$$|u{|}_{\mu +2s}=|f{|}_{\mu },$$

$$(-\Delta {)}^s\text{mathchar}58\text{accentset}\circ H^r({\mathbb{T}}^d;\mathbb{C})\to \text{accentset}\circ H^{r-2s}({\mathbb{T}}^d;\mathbb{C}).$$

Full text

**AMS subject classification **35S15, 26A33, 65R20, 65N12, 65N35, 49M05, 49K20

1 Introduction

Let $\mathbb{T}^{d}$, $d\geq 1$, be the $d$-dimensional torus. The purpose of this paper is to study the approximation of problems involving the fractional Laplace operator of order $2s$

[TABLE]

using the Fourier spectral method and to illustrate the importance of fractional differential operators. Such operators appear in various models with periodic boundary conditions, see [18, 28, 32]. The approach discussed here extends to problems with other boundary conditions such as Dirichlet or Neumann boundary conditions. For a discussion on nonhomogeneous boundary conditions we refer to [2].

Motivated by applications including fracture mechanics and turbulence, see [24, 16, 21], problems with fractional derivatives have recently gained a lot of interest. Several experiments suggest the presence of fractional derivatives, for instance, the electrical signal propagation in a cardiac tissue [10]. The appearance of fractional derivatives there is attributed to the heterogeneity of the underlying medium. A question arises: can one, for example, tailor the diffusion coefficient in [10] to get the same effects as the fractional model? Indeed to arrive at a direct evidence justifying the presence of fractional derivatives is a difficult question to address. This paper is an attempt to partially address this question by considering two specific applications where the presence of spectral fractional operators makes a significant difference. In particular, we illustrate the effect and advantages of fractional derivatives on two, by now, classical problems: image denoising and phase field modeling.

A well-known total variation based image denoising model is the so-called Rudin–Osher–Fatemi (ROF) model [33] which seeks a minimizer $u\in BV(\mathbb{T}^{d})\cap L^{2}(\mathbb{T}^{d})$ for

[TABLE]

Here $\mathbb{T}^{d}$ denotes the image domain, $\|\cdot\|$ is the norm in $L^{2}(\mathbb{T}^{d};\mathbb{C})$ with corresponding inner product $(\cdot,\cdot)$, and $\alpha>0$ is a regularization parameter. The function $g\mathrel{\mathop{\mathchar 58\relax}}\mathbb{T}^{d}\rightarrow\mathbb{C}$ represents the given observed possibly noisy image. The first term in $E$ is the total variation $|Du|(\mathbb{T}^{d})$ which has a regularizing effect but at the same time allows for discontinuities which may represent edges in the image. The second term is the fidelity term which measures the distance to the given image. Often, weaker norms such as the $H^{-1}$ norm are considered to define the latter term. While the existence and uniqueness of a minimizer can be established via the direct method of calculus of variations, the non-differentiability of the total variation is challenging from a computational point of view. In fact, a non-exhaustive list of papers that have attempted to resolve this are [5, 7, 15, 27, 19, 25, 30, 6]. Another question to ask is, whether natural images belong to $BV(\mathbb{T}^{d})\cap L^{2}(\mathbb{T}^{d})$. The paper [26] shows that natural images are incompletely represented by $BV(\mathbb{T}^{d})$ functions. We will handle both these shortcomings by replacing the total variation term in (1.1) by a squared fractional Sobolev norm. In other words, we propose to minimize

[TABLE]

with $0<s<1$ and $\beta\in[0,1]$. The first order necessary and sufficient optimality condition determines the unique minimizer $u$ via

[TABLE]

which is a linear elliptic partial differential equation (PDE) that can be efficiently solved using, for instance, the Fourier spectral method (which is the focus of this paper) or the so-called Caffarelli-Silvestre extension (in $\mathbb{R}^{n}$) [11] and the Stinga-Torrea extension (in bounded domains) [32, 37], see also [31]. Our experiments reveal that the fractional model (1.3) leads to results which are comparable to those provided by the ROF model but at a significantly reduced computational effort (see Section 4). We remark that fractional derivatives have been used in image processing before; see [13, 14], where the authors use spectral methods, and [23], where the authors use the finite element method. However, in both these cases the authors solve a fractional dynamical system with initial condition given by $g$. For completeness, we also refer to [4] where the authors consider a fractional norm equivalent regularization in the context of optimal control problems and parameter identification problems – this equivalent norm was realized using a mutilevel approach.

A mathematical justification of our choice of (1.2) as a substitute to (1.1) is given next. We seek $u$ solving (1.3) in a fractional Sobolev space $H^{s}(\mathbb{T}^{d})$. Moreover, we notice that if $g\in L^{\infty}(\mathbb{T}^{d})$ then following Theorem 3.5 part (1)(b) of [38] it is possible to show that $u\in L^{\infty}(\mathbb{T}^{d})$, see also [3]. We will next see that $BV(\mathbb{T}^{d})\cap L^{\infty}(\mathbb{T}^{d})$ is contained in $H^{s}(\mathbb{T}^{d})$ for $s<1/2$. Indeed by Lemma 38.1 of [39] we have the following continuous embedding

[TABLE]

where $B^{1/2}_{2,\infty}(\mathbb{T}^{d})$ is a Besov space. In addition, using Proposition 1.2 of [20], see also [29, pg. 1222] and [41, Section 3], we have the following continuous embedding

[TABLE]

provided that $s<1/2$. Finally, combining the inclusions we arrive at

[TABLE]

which justifies our energy functional (1.2). We remark that the regularizing quadratic term in (1.2) does not have the gradient sparsity property of the total variation norm. This effect however cannot be proven for the ROF model in general due the presence of the quadratic fidelity term but is certainly visible in experiments.

As a second example we consider gradient flows of the energy functional

[TABLE]

with initial condition $u(0)=u_{0}$. The $L^{2}$-gradient flow of (1.4) leads to the fractional Allen–Cahn equation

[TABLE]

where $0<s<1$, $T>0$, and $f=F^{\prime}$ is typically nonlinear in $u$. Moreover,

[TABLE]

where $0<\widetilde{\varepsilon}<1$. When $s=1$, we set $\varepsilon^{-2}\mathrel{\mathop{\mathchar 58\relax}}=\widetilde{\varepsilon}^{-2}$. The $H^{-\alpha}$ gradient flow, with $\alpha\in(0,1]$ leads to the fractional Cahn–Hilliard equation

[TABLE]

By testing (1.6) with a constant function it is easy to check that (1.6) is mass conserving. We note that throughout this article we consider $\widetilde{\varepsilon}$ as a fixed small number. The aforementioned scaling of $\widetilde{\varepsilon}$ is the right scaling to obtain a sharp interface limit as $\widetilde{\varepsilon}\to 0$ we refer the reader to [35]. We remark that even though the optimality system in case of image denoising is linear (1.3), the system for the fractional phase field model is nonlinear (1.6) and controlling these nonlinearities in the presence of fractional derivatives turned out to be a nonobvious task.

When $s=1$, a standard numerical method requires a fine mesh resolution around interfaces to capture sharp transitions [8]. The fully discrete scheme proposed in this paper is unconditionally stable and supported by a rigorous error analysis. In our experiments we observe that using the spectral method and choosing small values for $s$, it is possible to obtain sharp interfaces on relatively coarse meshes and moderate values for $\widetilde{\varepsilon}$.

We remark that the use of spectral methods in the context of phase field models (when $s=1$, $\alpha=0$ or $\alpha=1$) has been considered before, see [9, 17]. We further remark that the recent paper [36] also investigates a fractional Allen–Cahn equation and uses the fractional Riemann-Liouville derivative which is different from our definition. In addition, few analytical details are provided. We also refer to [1] which discusses analytical properties of a fractional Cahn–Hilliard equation with fixed $s=1$. The authors report that the dynamics in case $\alpha>0$ and $s=1$ are closer to the classical Cahn-Hilliard equation than to the Allen-Cahn equation. The error analysis provided there is restricted to spatial discretizations while the used fully discrete scheme treats the nonlinearity explicitly and is hence only conditionally stable.

We remark that the goal of this article is to show possible applications of fractional PDEs. The simple image denoising problem serves as a model problem in which the effect of fractional derivatives of different order becomes directly apparent. The nonlinear evolution model defined by the fractional phase field equation combines different effects so that an interpretation of the effect of fractional derivatives of different order requires a more careful interpretation. Our experiments are meant to illustrate these effects for different parameters $s$.

This paper is organized as follows: In Section 2 we recall facts about spectral interpolation estimates in fractional Sobolev spaces. The fractional Laplace operator and its discretization by the spectral method are addressed in Section 3. We present the details on the fractional image denoising problem in Section 4. Section 5 is devoted to a general error analysis for a numerical scheme covering both the fractional Allen–Cahn and Cahn–Hilliard equations. We conclude with several illustrative numerical examples in Section 6.

2 Spectral approximation

In this section we specify notation needed to define the discrete Fourier transformation and recall elementary approximation results in fractional Sobolev spaces.

2.1 Discrete Fourier transformation

We consider the $2\pi$-periodic torus $\mathbb{T}^{d}$ and the set of grid points $(x_{j}\mathrel{\mathop{\mathchar 58\relax}}j\in\mathbb{N}_{n}^{d})$ on $\mathbb{T}^{d}$ defined by $x_{j}=(j_{1},\dots,j_{d})\frac{2\pi}{n},$ where \mathbb{N}_{n}^{d}=\big{\{}j=(j_{1},\dots,j_{d})\in\mathbb{Z}^{d}\mathrel{\mathop{\mathchar 58\relax}}0\leq j_{i}\leq n-1\big{\}}. A family of grid functions $(\Phi^{k}\mathrel{\mathop{\mathchar 58\relax}}k\in\mathbb{Z}_{n}^{d})$ is defined by

[TABLE]

where \mathbb{Z}_{n}^{d}=\big{\{}k=(k_{1},\dots,k_{d})\in\mathbb{Z}^{d}\mathrel{\mathop{\mathchar 58\relax}}-n/2\leq k_{i}\leq n/2-1\big{\}} and $\mathrm{i}^{2}=-1$. For grid functions $V=(v_{j}\mathrel{\mathop{\mathchar 58\relax}}j\in\mathbb{N}_{n}^{d})$ and $W=(w_{j}\mathrel{\mathop{\mathchar 58\relax}}j\in\mathbb{N}_{n}^{d})$ we define the discrete scalar product

[TABLE]

The associated norm is denoted $\|\cdot\|_{n}$. Notice that the family $(\Phi^{k}\mathrel{\mathop{\mathchar 58\relax}}k\in\mathbb{Z}_{n}^{d})$ defines an orthogonal basis for the space of grid functions with $\|\Phi^{k}\|_{n}=(2\pi)^{d/2}$. The discrete Fourier transform of a grid function $V=(v_{j}\mathrel{\mathop{\mathchar 58\relax}}j\in\mathbb{N}_{n}^{d})$ is the coefficient vector $\widetilde{V}=(\widetilde{v}_{k}\mathrel{\mathop{\mathchar 58\relax}}{k\in\mathbb{Z}_{n}^{d}})$ with

[TABLE]

With these coefficients we have $V=(2\pi)^{-d}\sum_{k\in\mathbb{Z}_{n}^{d}}\widetilde{v}_{k}\Phi^{k}.$

2.2 Trigonometric interpolation

We consider the space of trigonometric polynomials defined via

[TABLE]

with the functions $\varphi^{k}(x)=e^{\mathrm{i}k\cdot x}$ which define an orthogonal basis for $\mathcal{S}_{n}$ with respect to the inner product in $L^{2}(\mathbb{T}^{d};\mathbb{C})$. With $v\in C(\mathbb{T}^{d};\mathbb{C})$ we associate a grid function $V=(v_{j}\mathrel{\mathop{\mathchar 58\relax}}{j\in\mathbb{N}_{n}^{d}})$ via $v_{j}=v(x_{j})$, $j\in\mathbb{N}_{n}^{d}$. Notice that for $v_{n},w_{n}\in\mathcal{S}_{n}$ with associated grid functions $V,W$ we have

[TABLE]

The discrete Fourier transformation gives rise to a nodal interpolation operator.

Definition 1**.**

Given $v\in C(\mathbb{T}^{d};\mathbb{C})$ with nodal values $V=(v_{j}\mathrel{\mathop{\mathchar 58\relax}}{j\in\mathbb{N}_{n}^{d}})$ and discrete Fourier coefficients $\widetilde{V}=(\widetilde{v}_{k}\mathrel{\mathop{\mathchar 58\relax}}{k\in\mathbb{Z}_{n}^{d}})$, the trigonometric interpolant $I_{n}v\in\mathcal{S}_{n}$ of $v$ is defined via

[TABLE]

Remark 2**.**

*(i) Note that $I_{n}v(x_{j})=v(x_{j})$ for all ${j\in\mathbb{N}_{n}^{d}}$.

(ii) We have $\widetilde{v}_{k}=(v,\varphi^{k})_{n}$ for all ${k\in\mathbb{Z}_{n}^{d}}$.

(iii) We have $(I_{n}v,w_{n})_{n}=(v,w_{n})_{n}$ for all $w_{n}\in\mathcal{S}_{n}$ and $v\in C(\mathbb{T}^{d};\mathbb{C})$.*

The (continuous) Fourier transform of a function $v\in L^{2}(\mathbb{T}^{d};\mathbb{C})$ is the coefficient vector $\widehat{V}=(\widehat{v}_{k}\mathrel{\mathop{\mathchar 58\relax}}k\in\mathbb{Z}^{d})$ defined by

[TABLE]

Note that here the $L^{2}$ inner product is used instead of its discrete approximation. With respect to convergence in $L^{2}(\mathbb{T}^{d};\mathbb{C})$ we have that $v=(2\pi)^{-d}\sum_{k\in\mathbb{Z}^{d}}\widehat{v}_{k}\varphi^{k},$ and, in particular, Plancherel’s formula $(v,w)=(2\pi)^{-d}(\widehat{v},\widehat{w})_{\ell^{2}(\mathbb{Z}^{d})}$.

2.3 Approximation in Sobolev spaces

We analyze the approximation properties of the interpolation operator $I_{n}$ in terms of Sobolev norms and with the help of the $L^{2}$ projection onto $\mathcal{S}_{n}$ which is obtained by truncation of the Fourier series of a function.

Definition 3**.**

The $L^{2}$ projection $P_{n}\mathrel{\mathop{\mathchar 58\relax}}L^{2}(\mathbb{T}^{d};\mathbb{C})\to\mathcal{S}_{n}$ is for $v\in L^{2}(\mathbb{T}^{d};\mathbb{C})$ defined by the condition that for all $w_{n}\in\mathcal{S}_{n}$ we have

[TABLE]

Note that for every $v\in L^{2}(\mathbb{T}^{d};\mathbb{C})$ we have $P_{n}v=\sum_{k\in\mathbb{Z}_{n}^{d}}\widehat{v}_{k}\varphi^{k}$. The following definition is motivated by the fact that $(\widehat{\partial^{\alpha}v})_{k}=\mathrm{i}^{|\alpha|}k^{\alpha}\widehat{v}_{k}$ for every $\alpha\in\mathbb{N}_{0}^{d}$.

Definition 4**.**

Given $\mu\geq 0$ the Sobolev space $H^{\mu}(\mathbb{T}^{d};\mathbb{C})$ consists of all functions $v\in L^{2}(\mathbb{T}^{d};\mathbb{C})$ with

[TABLE]

Its dual $H^{-\mu}(\mathbb{T}^{d};\mathbb{C})$ consists of all linear functionals $\psi\mathrel{\mathop{\mathchar 58\relax}}H^{\mu}(\mathbb{T}^{d};\mathbb{C})\to\mathbb{C}$ with

[TABLE]

where $\widehat{\psi}_{k}=\psi(\phi_{k})$.

The Sobolev spaces allow us to quantify approximation properties of the operators $P_{n}$ and $I_{n}$. We refer the reader to Chapter 8 in [34] for details.

Lemma 5** (Projection error).**

For $\lambda,\mu\in\mathbb{R}$ with $\lambda\leq\mu$ and $v\in H^{\mu}(\mathbb{T}^{d};\mathbb{C})$ we have

[TABLE]

By comparing $P_{n}$ and $I_{n}$ we obtain a trigonometric interpolation estimate. It is shown in Remark 8.3.1 of [34] that the conditions of the following result cannot be improved in general.

Lemma 6** (Interpolation error).**

If $\mu>d/2$, $0\leq\lambda\leq\mu$, and $v\in H^{\mu}(\mathbb{T}^{d};\mathbb{C})$ we have

[TABLE]

with a constant $c_{d,\lambda,\mu}>0$ that is independent of $v$ and $n$.

We conclude the section with an inverse estimate. Particularly, for every function $v_{n}\in\mathcal{S}_{n}$ and $r\geq s$ we have that

[TABLE]

3 Fractional Laplace operator

We define subspaces of Sobolev spaces via

[TABLE]

For $r\geq 0$ the subspaces consist of Sobolev functions with vanishing mean. On the subspaces $\accentset{\circ}{H}^{r}(\mathbb{T}^{d};\mathbb{C})$ the corresponding seminorms $|\cdot|_{r}$ are norms.

Definition 7**.**

For $s,\mu\geq 0$ and $v\in H^{\mu}(\mathbb{T}^{d};\mathbb{C})$ the fractional Laplacian of $v$ is the (generalized) function $(-\Delta)^{s}v\in\accentset{\circ}{H}^{\mu-2s}(\mathbb{T}^{d};\mathbb{C})$ defined by

[TABLE]

Given $f\in L^{2}(\mathbb{T}^{d};\mathbb{C})$ with vanishing mean the fractional Poisson problem seeks $u\in\accentset{\circ}{H}^{s}(\mathbb{T}^{d};\mathbb{C})$ with

[TABLE]

The unique solution to (3.1) is given by

[TABLE]

and in fact satisfies $u\in\accentset{\circ}{H}^{2s}(\mathbb{T}^{d};\mathbb{C})$. More generally, for $f\in\accentset{\circ}{H}^{\mu}(\mathbb{T}^{d};\mathbb{C})$ we have

[TABLE]

i.e., the fractional Laplacian defines an isometric isomorphism

[TABLE]

We define the fractional Laplace operator of negative order as the inverse of $(-\Delta)^{s}$, i.e.,

[TABLE]

Note that for $r,s\in\mathbb{R}$ and $v\in\accentset{\circ}{H}^{r}(\mathbb{T}^{d};\mathbb{C})$ with $(-\Delta)^{s}v\in L^{2}(\mathbb{T}^{d};\mathbb{C})$ we have

[TABLE]

If $s\leq r$ we have the continuous embedding property

[TABLE]

The discretized fractional Poisson problem seeks for a given $f_{n}\in\mathcal{S}_{n}$ with vanishing mean a function $u_{n}\in\mathcal{S}_{n}$ with

[TABLE]

The uniquely defined solution is given by

[TABLE]

There are two noticeable differences between the continuous (3.2) and the discrete solutions (3.5). Besides the finite and the infinite sums, $u_{n}$ contains the discrete Fourier coefficients $\widetilde{f}_{k}$ and $u$ contains the continuous Fourier coefficients $\widehat{f}_{k}$. The following a priori error estimates hold.

Theorem 8**.**

Let $u$ and $u_{n}$ solve the continuous (3.1) and the discrete (3.4) problems, respectively. We have

[TABLE]

In particular, if $f\in\accentset{\circ}{H}^{\mu}(\mathbb{T}^{d};\mathbb{C})$ and $f_{n}=P_{n}f$ we have

[TABLE]

while if $f\in\accentset{\circ}{H}^{\nu}(\mathbb{T}^{d};\mathbb{C})$ with $\nu>d/2$ and $f_{n}=I_{n}f$ we have

[TABLE]

Proof.

In view of (3.1) and (3.4) we have

[TABLE]

This implies the general estimate and in combination with Lemma 5 the estimate in case $f_{n}=P_{n}f$. With (3.3) and Lemma 6 we deduce that

[TABLE]

which implies the estimate. ∎

4 Fractional image denoising

Our second problem is a replacement of the ROF image denoising model (1.1). Given an image $g\in L^{2}(\mathbb{T}^{d};\mathbb{C})$ we propose to minimize

[TABLE]

with $0<s<1$ and $\beta\in[0,1]$. The minimization is carried out over $H^{s}(\mathbb{T}^{d};\mathbb{C})$ when $\beta=0$ and over $\accentset{\circ}{H}^{s}(\mathbb{T}^{d};\mathbb{C})$ when $\beta$ is positive. In the latter case we assume that $g$ has a vanishing mean. The existence and uniqueness of a minimizer follows by using the direct method in the calculus of variations. The first order necessary and sufficient optimality condition determines the unique minimizer $u$ via

[TABLE]

We note that since $g\in L^{2}(\mathbb{T}^{d};\mathbb{C})$, we have $u\in H^{2(s+\beta)}(\mathbb{T}^{d};\mathbb{C})$. In particular, the solution to (4.2) is

[TABLE]

The discretized problem seeks for a given $g_{n}\in\mathcal{S}_{n}$ a function $u_{n}\in\mathcal{S}_{n}$ with

[TABLE]

The uniquely defined solution is given by

[TABLE]

Theorem 9**.**

Let $u$ and $u_{n}$ solve the continuous and the discrete problems (4.3) and (4.2), respectively. We have that

[TABLE]

In particular, if $g\in H^{\mu}(\mathbb{T}^{d};\mathbb{C})$ and $g_{n}=P_{n}g$ we have

[TABLE]

while if $g\in H^{\nu}(\mathbb{T}^{d};\mathbb{C})$ with $\nu>d/2$ and $g_{n}=I_{n}g$ we have

[TABLE]

Proof.

Testing the difference of (4.2) and (4.3) by $u-u_{n}$ implies that

[TABLE]

The estimates follow from using $(a+b)^{2}\leq 2(a^{2}+b^{2})$ and arguing as in the proof of Theorem 8. ∎

5 Fractional phase field equations

Given parameters $\alpha,s\geq 0$ we recall the fractional Cahn-Hilliard equation (1.6)

[TABLE]

on a $d$-dimensional torus $\mathbb{T}^{d}$ and with initial condition $u(0)=u_{0}$. We recall that $\alpha=0$ gives rise to the fractional Allen–Cahn equation (1.5). Below we will impose the restrictions $s>0$ and $s\geq\alpha$.

We assume a splitting of the nonnegative potential $F$ into convex and concave parts $F^{cx}$ and $F^{cv}$ which induces a decomposition of $f=F^{\prime}$ into a monotone and an antimonotone part

[TABLE]

We assume for simplicity that $f^{cx}$ and $f^{cv}$ are smooth and Lipschitz continuous. The latter condition is justified by a maximum principle in the case for the Allen–Cahn equation and $L^{\infty}$ bounds for solutions of the Cahn–Hilliard equation [12] corresponding to $(\alpha,s)=(0,1)$ and $(\alpha,s)=(1,1)$, respectively.

5.1 Numerical scheme and error analysis

The numerical scheme computes iterates $(u_{n}^{k})_{k=0,\dots,K}\subset\mathcal{S}_{n}$ via

[TABLE]

where $d_{t}w^{k}=(w^{k}-w^{k-1})/\tau$ with $\tau>0$ being the time step-size and $u_{n}^{0}$ is a suitable approximation of $u^{0}$. By applying the operator $(-\Delta)^{\alpha}$ and testing the resulting identity with constant functions we observe the mass conservation property $(d_{t}u_{n}^{k},1)=0$ if $\alpha>0$. Existence of the iterates is established via convex minimization problems; if the convex part of $F$ is quadratic then $f^{cx}$ is linear and the scheme (5.2) defines a linear system of equations. The scheme is unconditionally energy stable in the sense that we have

[TABLE]

with the discrete energy functional

[TABLE]

This follows from testing (5.2) with $d_{t}u_{n}^{k}$, using

[TABLE]

and noting that as a consequence of convexity and concavity we have

[TABLE]

Assuming initial data with $E_{\varepsilon}^{n}(u_{n}^{0})\leq c$ as $n\to\infty$, the energy estimate provides a priori bounds on interpolants of the approximations which allows us to select an accumulation point

[TABLE]

as $\tau\to 0$ and $n\to\infty$. Its identification as a solution for the fractional Cahn–Hilliard equation follows from the Aubin–Lions lemma provided that $s>0$. Uniqueness of solutions is a consequence of the assumed Lipschitz continuity of $f$. For an error analysis we note that $u\in C([0,T];L^{2}(\mathbb{T}^{d}))$, let $u^{k}=u(t_{k})$ with $t_{k}=k\tau$, and define

[TABLE]

where $P_{n}$ is the orthogonal projection given in Definition 3. Note that we have $(P_{n}v,w_{n})_{n}=(P_{n}v,w_{n})$ but $(v,w_{n})_{n}\neq(v,w_{n})$ unless $v$ belongs to $\mathcal{S}_{n}$. We omit the subscript $n$ whenever the scalar product is applied to two functions belonging to $\mathcal{S}_{n}$. For ease of readability we abbreviate

[TABLE]

The sequences $(u_{n}^{k})$ and $(u^{k})$ satisfy the discrete equations

[TABLE]

where we used that $P_{n}$ commutes with $(-\Delta)^{r}$ for every $r\in\mathbb{R}$. Subtracting the identities leads to the error equation

[TABLE]

with the discretization errors

[TABLE]

Testing the error equation with $e_{n}^{k}$ shows that we have

[TABLE]

To bound the first term on the right-hand side we insert $f_{\varepsilon}(u_{n}^{k})$, use Lemma 6, the inverse estimate (2.1), and insert $P_{n}u^{k}$ to deduce with the Lipschitz continuity of $f_{\varepsilon}$ that

[TABLE]

For the second term we have that

[TABLE]

where we used the inverse estimate $\|w_{n}\|\leq n^{\alpha}|w_{n}|_{-\alpha}$ if $\alpha>0$. For the third term we assume that $u\in C^{2}([0,T];H^{-\alpha}(\mathbb{T}^{d}))$ and estimate

[TABLE]

A combination of the estimates, multiplication by $\tau$, and summation over $k=1,2,\dots,K$, yield that

[TABLE]

If $\alpha=0$ we may directly apply the discrete Gronwall lemma to obtain an error estimate. If $\alpha>0$ we assume $\alpha\leq s$, require that $(u_{n}^{0},1)=(u^{0},1)$ so that $(e_{n}^{k},1)=0$, and use the bound

[TABLE]

to deduce with Young’s inequality the estimate

[TABLE]

Here, $D_{0}$ is the sum of the first four terms on the right-hand side of the above estimate and $D_{1}=2(3|f_{\varepsilon}|_{\text{Lip}})^{2}$. The discrete Gronwall lemma leads to the estimate

[TABLE]

for all $K$ with $K\tau\leq T$ provided that $\tau D_{1}\leq 1/2$. With the triangle inequality and approximation estimates for $P_{n}$ we obtain the following error estimate.

Theorem 10**.**

Let $u\in C([0,T];H^{s}(\mathbb{T}^{d}))\cap C^{2}([0,T];H^{-\alpha}(\mathbb{T}^{d}))$ solve (1.6) and let the sequence $(u_{n}^{k})_{k=0,\dots,K}\subset\mathcal{S}_{n}$ be defined via (5.2). There exists a constant $c_{\varepsilon}>0$ such that we have

[TABLE]

for all $\tau>0$ and $n\in\mathbb{N}$.

The constant $c_{\varepsilon}$, in general, depends exponentially on $\varepsilon^{-1}$. For the derivation of this estimate we used the indicated regularity assumption. By standard arguments, see [40] the regularity assumption on the exact solution can be weakened to the conditions

[TABLE]

The suboptimal term $\tau n^{\alpha}$ corresponds to the semi-implicit treatment of the nonlinearity which makes the scheme (5.2) fully practical.

5.2 Improved error estimate via spectral bounds

A significantly improved error estimate can be obtained if additional analytical knowledge about the evolution is available, e.g., in the form of lower bounds for the principal eigenvalue

[TABLE]

For ease of presentation we only consider the fractional Allen–Cahn equation with $\alpha=0$ and outline the main arguments following [22, 8]. We focus on the contribution to the error equation resulting from the nonlinearity and write it with abstract consistency functionals $\mathcal{C}_{n}^{k}$ as

[TABLE]

A precise formula for $\mathcal{C}_{n}^{k}$ is obtained from subtracting the projected partial differential equation evaluated at $t_{k}$ onto $\mathcal{S}_{n}$ from the numerical scheme as above, e.g., in case of a fully implicit numerical scheme with a nodal interpolation of the nonlinear term we have

[TABLE]

To relate the error equation to the principal eigenvalue we require a controlled failure of monotonicity for $f$ in the sense that there exists a constant $c_{f}>0$ such that for all $a,b\in\mathbb{R}$ we have

[TABLE]

With this estimate we deduce with $c_{f}^{\prime}=\|f^{\prime}\|_{L^{\infty}(\mathbb{R})}$ that

[TABLE]

We incorporate the eigenvalue $\lambda^{k}=\lambda(t_{k})$ via

[TABLE]

Choosing $\theta=\varepsilon^{2}$ and letting $\mu^{k}=\max\{-\lambda^{k},0\}$ thus leads to the estimate

[TABLE]

Rearranging terms gives

[TABLE]

In an inductive argument we may assume that $\varepsilon^{-2}\|e_{n}^{k}\|_{L^{3}(\mathbb{T}^{d})}^{3}\leq c\|e_{n}^{k}\|^{2}$ and use the discrete Gronwall lemma to obtain an error estimate that depends exponentially on the principal eigenvalue $\lambda^{k}$. Hence, if $\lambda^{k}$ is uniformly bounded from below the resulting error estimate depends only algebraically on $\varepsilon^{-1}$. More generally, it suffices to assume that a discrete time integral of $\lambda^{k}$ is uniformly bounded from below. This allows us to cover large classes of evolutions including topological changes.

6 Numerical Examples

In this Section, we present several numerical examples. In Section 6.1 we discuss the approximation of the fractional Poisson problem. Section 6.2 is devoted to image denoising problem. In Section 6.3 we study features of the fractional Allen–Cahn equation. We conclude with experiments for the fractional Cahn–Hilliard equation in Section 6.4.

6.1 Approximation of the Poisson problem

To construct a nonsmooth solution for the fractional Poisson problem we first let $w\in C(\mathbb{T})$ be defined via

[TABLE]

Since $w(0)=w(2\pi)$ we find for $k\neq 0$ with an integration-by-parts that

[TABLE]

i.e., $\widehat{w}_{k}=-4/k^{2}$ if $k$ is odd and $\widehat{w}_{k}=0$ if $k$ is even. We have $\widehat{w}_{0}=\pi^{2}$. We then let $u\in C(\mathbb{T}^{d})$ be for $x\in\mathbb{T}^{d}$ defined via

[TABLE]

We have $\widehat{u}_{k}=\widehat{w}_{k_{1}}\cdots\widehat{w}_{k_{d}}$ if $k\neq 0$ and $\widehat{u}_{0}=0$. We set $f=(-\Delta)^{s}u$, i.e., for ${k\in\mathbb{Z}_{n}^{d}}$ let $\widehat{f}_{k}=|k|^{2s}\widehat{u}_{k}$. Note that $f\in L^{2}(\mathbb{T}^{d})$ if and only if $s<1/2$. We choose the approximation $f_{n}=P_{n}f$ which is explicitly available here. The output for $s=1/2$ and $n=16$ is displayed in Figure 1. In contrast to solutions of the classical Poisson problem we observe here the occurrence of kinks in the solution.

6.2 Fractional image denoising

The fractional Laplacian with $s<1/2$ is closely related to the total variation norm (see Section 1) which motivates its application in image processing. Given a noisy image $g\in L^{2}(\mathbb{T}^{d})$ we define a regularized image $u\in H^{2s}(\mathbb{T}^{d})$ via

[TABLE]

The fidelity parameter $\alpha$ penalizes the deviation of $u$ from $g$ in the $L^{2}$ metric. Alternatively, this distance can be taken in the weaker metric of $H^{-1}(\mathbb{T}^{d})$ which leads to the equation

[TABLE]

where we assume that $g$ has vanishing mean and look for $u$ with vanishing mean. The results of two experiments are displayed in the rows of Figure 2. In the first experiment we set $s=0.42$, $\alpha=10$, and $n=1566$. In the second experiment we used $s=0.35$, $\alpha=5\times 10^{3}$ and $n=256$. The first and second columns display the original and noisy images, respectively. The third and fourth columns show the results of $L^{2}$ and $H^{-1}$ fidelity. We note that in the first example, where the additive noise is normally distributed with mean zero and standard deviation $0.15$, the $L^{2}$-fidelity almost perfectly recovers the original image reflecting the fact that for Gaussian noise this is statistically the optimal choice. In the second example where the noise is given by the nodal interpolant of the sinusoidal function

[TABLE]

we obtain better recovery using the $H^{-1}$-fidelity. We remark that it took 0.2 sec and 0.006 sec to solve the first and the second problem in Matlab on a MacBook Pro with an 2.8 GHz Intel Core i7 processor (16 GB 1600 MHz DDR3 RAM).

Comparison between fractional and total variation

We next compare the fractional and the total variation (TV) based models (1.1) with the help of three examples. This comparison was carried out in Python which was specifically chosen due to the availability of SciKit-Image toolbox [42]. In all the examples we assume that original image denoted by $o$ (without noise) is known. For the fractional case we compute the optimal parameters $(s,\alpha)$ by solving a minimization problem: $\min_{u}$ $\|u-o\|^{2}$, subject to $u$ solving $(-\Delta)^{s}u+\alpha(u-g)=0.$ We further assume that $(s,\alpha)$ lies in a closed convex set, i.e., $0.05\leq s\leq 0.5$ and $1\leq\alpha\leq 50$. We solve this optimization problem using an in-built optimization algorithm in Python. The corresponding optimization problem for TV is solved for $\alpha$ using a genetic algorithm.

Our first example uses a picture of Gauss (cf. Figure 3, top row). The left image is the original image with $n=n_{x}=n_{y}=1566$. Our second example (cf. Figure 3, middle row) uses a synthetic image with $n=n_{x}=n_{y}=1500$. Our final example (cf. Figure 3, bottom row) is based on an in-built image from SciKit with different number of points in the $x$ and $y$ directions, i.e., $n_{x}=400$, $n_{y}=600$. We note that even though the approach discussed in Section 4 assumes $n=n_{x}=n_{y}$ it is directly extended to this case where $n_{x}\neq n_{y}$. In the second column (from the left), in all the examples, we have added a normally distributed noise with mean zero and standard deviation $0.15$, we denote the resulting noisy image by $g$.

For the fractional case the optimal parameters are $(s,\alpha)=(0.49,50)$ for the first example and $(s,\alpha)=(0.5,50)$ for the second and third examples. On the other hand, for TV case the optimal parameters are $\alpha=3.58$, $\alpha=1.0$, and $\alpha=8.74$, respectively. Using these parameters we solve the corresponding image denoising problems using the fractional approach (third column from the left) and compare the results with the in-built TV algorithm from [42]. We observe that the two approaches give comparable results in all three examples. The fractional approach is significantly cheaper as only two discrete Fourier transformations have to be computed. The CPU times for implementations of the fractional and TV approach in Matlab and Python, respectively, are provided in the caption of Figure 3. They show a reduction of the computing times by factors 10-100. We remark that certain aspects in the Chambolle–Pock algorithm implemented in Scikit such as the specification of a suitable stopping criterion and choice of step-sizes may lead to different results.

6.3 Fractional Allen–Cahn equation

We consider the fractional Allen-Cahn equation (1.5) with a given initial function $u_{0}\in L^{2}(\mathbb{T}^{d})$. The function $f$ is the derivative of a double well potential $F$ with quadratic growth, i.e.,

[TABLE]

which leads to linear systems of equations in our semi-implicit time discretization. Snapshots of the evolutions with $\widetilde{\varepsilon}=1/8$, $n=512$, and $\Delta t=1/100$ at $t=1$, $t=4$, $t=12$, and $t=20$ (rowwise) are shown in Figure 4. The first column corresponds to $s=1$, second to $s=0.45$, third to $s=0.30$, and finally fourth to $s=0.15$. Clearly the interface in case of a fractional model is sharper, however the dynamics are slower.

6.4 Fractional Cahn-Hilliard equation

We next study the fractional Cahn–Hilliard equation specified in (1.6) defining $F$ and $f$ as in Subsection 6.3. In a first experiment we define the initial condition $u_{0}$ as

[TABLE]

Snapshots of the evolutions with $\widetilde{\varepsilon}=1/8$, $n=512$, and $\Delta t=1/100$ at $t=0.25$, $t=0.50$, $t=2$, and $t=3$ (rowwise) are shown in Figure 5 for $s=1$, $s=0.45$, $s=0.30$, and $s=0.15$ (columnwise) with $\alpha=1$ in all cases.

In a second experiment we focus on the coarsening dynamics of the fractional phase field equation. As in [1] the initial condition is given by $u_{0}=2\phi-1+\delta$ where $\delta$ is a random perturbation uniformly distributed in [-0.2,0.2]. Snapshots of the evolutions with $\phi=0.5$, $\widetilde{\varepsilon}=1/8$, $n=512$, and $\Delta t=1/100$ at $t=0.25$, $t=0.5$, $t=1$, and $t=1.5$ (rowwise) are shown in Figure 6. The first two columns correspond to $\alpha=1$ with $s=1$ and $s=0.20$, respectively. The last two columns are obtained with $\alpha=1/2$ with $s=1$ and $s=0.20$, respectively.

Acknowledgement

We thank Pablo Stinga for stimulating discussions. We also thank Diego Torrejon for help with Python.

Funding

The work of the first author is partially supported by NSF grant DMS-1521590.