Preconditioned Pseudo-Spectral Gradient Flow for Computing the Steady-State of Space Fractional Cahn-Allen Equations With Variable Coefficients

The aim of this paper is to propose some efficient and accurate numerical methods to compute the steady-state of variable coefficients space fractional Cahn-Allen equations. The approach combines an adaptive time stepping semi-implicit gradient flow method to minimize the fractional energy functional and pseudo-spectral approximation schemes. Based on the use of a preconditioned GMRES, the space fractional Cahn-Allen equation is then solved efficiently. The full methodology is supported by the numerical solution of a one-dimensional problem.

1 Introduction

During the last 2 decades, the mathematical analysis and numerical simulation of fractional partial differential equations became one of the most important hot topics in applied and computational mathematics [1–4]. One of the reasons of this growing interest is related to the fact that these new models provide a better description of some physical phenomena. Indeed, for example, for the heat equation, considering space fractional effects allows to handle anomalous diffusion [5] when the stochastic process driving the medium is a Lévy α-stable flight. Fractional PDEs are also of interest when modeling e.g., the behavior of turbulent flows [6, 7], the chaotic effects in the dynamics of conservative systems [8] but also for the transport of contaminant in groundwater flow [9, 10] or finally in mathematical finance [11]. In [12], a fast algorithm based on time two-mesh finite element scheme, which aims at solving nonlinear problems quickly, is considered to numerically solve the nonlinear space fractional Cahn-Allen equations with smooth and non-smooth solutions. In [13], different schemes has been effectively employed for finding the exact solutions and dynamics of the Cahn-Allen model and the dispersive predator-prey model. The goal of this paper is to contribute to the numerical solution to such systems, and most particularly to a space fractional Cahn-Allen equation [14] with variable coefficients. The integer order equation is usually used to describe the complex phase separation and coarsening phenomena in a solid. Here, the fractional order appears as an exponent α/2 (1 ≤ α ≤ 2) in the laplacian operator, α = 2 corresponding to the standard situation. The fractional Cahn-Allen equation can then be seen as a more general situation than the standard one. It has been introduced (see e.g., [15–17]) for studying the competing stable phases having an identical Lyapunov functional density for Lévy processes of order α/2. This allows to include the situation where the random walk has correlations, non-Gaussian or non-Markovian memory effects which cannot be described by the standard Laplacian.

A recent interesting application [18, 19] considers image segmentation techniques based on the fractional Cahn-Allen equation, with examples in medical imaging. Here, the objective of the paper is to show how the widely used pseudo-spectral approximation schemes which are well-adapted to solve constant coefficients PDEs (see e.g., [20–23]) can be adapted to the variable coefficients case for computing the steady-state of fractional Cahn-Allen equations based on the gradient flow formulation and its suitable discretization.

To this end, we introduce in Section 2 the problem setting for an example of space fractional Cahn-Allen equation (the approach could be probably extended to other kinds of equations), and we provide a few useful definitions for the pseudo-spectral approximation, in particular concerning the spectral definition of the fractional laplacian. The gradient flow formulation is also given. Section 3 is devoted to the development of numerical schemes for the gradient flow as well as an adaptive time stepping method. We also explain how to get an efficient implementation of the pseudo-spectral approximation schemes, most particularly regarding the use of a Krylov subspace iterative solver with a well-designed preconditioner. In Section 4, we analyze the full methodology in the case of a one-dimensional example. This shows that the newly designed method performs very well for various fractional orders. Finally, we conclude in Section 5.

2 Problem Setting and Definition

The problem that we propose to solve is the following: find $u:{\mathbb{R}}^{*}\times \mathrm{\Omega }\to \mathbb{R}$ solution to the variable coefficients space fractional Cahn-Allen (SFCA) equation

$\left\{\begin{array}{c}{\partial }_{t}u=-\kappa \left(\mathbf{x}\right)\left(\epsilon {\left(-\mathrm{\Delta }\right)}^{\alpha /2}u+f\left(u\right)\right),\left(t,\mathbf{x}\right)\in \mathbb{R}\times \mathrm{\Omega },\hfill \\ u\left(t=0,\mathbf{x}\right)=u_0\left(\mathbf{x}\right),\mathbf{x}\in \mathrm{\Omega },\hfill \\ {\partial }_{\mathbf{n}}u\left(\mathbf{x}\right)=0,t>0,\mathbf{x}\in \mathrm{\Gamma }≔\partial \mathrm{\Omega }.\hfill \end{array}\right.(1)$

In the above system, u≔u(t, x) represents the concentration of one of the mixtures at time t and point x appearing in the modeling of the phase of a binary mixture. The domain of computation is a d-dimensional rectangular box: $\mathrm{\Omega }≔{\prod }_{k=1}^d\left.\right]a_k;b_k\left[\right.$, with − ∞ < a_k < b_k < ∞, for k = 1, … , d, with d = 1, 2, 3. Vector n denotes the outwardly directed unit normal vector to Ω. The strictly positive function κ is the variable mobility. For simplifications, it is standard to set this function as equal to 1, but we consider also here the nontrivial situation for pseudo-spectral methods where it is x-dependent. In the standard constant coefficients CA equation with integer order α = 2, the (usually small) constant ɛ > 0 stands for the interfacial width, which captures the dominating effect of reaction kinetics and stays for effective diffusivity.

In this paper, we consider the SFCA equation for a fractional laplacian (−Δ)^α/2, 1 < α ≤ 2, with homogeneous Neumann boundary conditions [24]. Nevertheless, homogeneous Dirichlet or periodic boundary conditions can also be included. In the one-dimensional case, we introduce λ_j as the eigenvalues and φ_j, $j\in \mathbb{N}$, as the associated eigenfunctions of the operator $-\mathrm{\Delta }≔-{\partial }_x^2$ with homogeneous Neumann boundary condition:

$\left\{\begin{array}{c}-\mathrm{\Delta }{\phi }_j={\lambda }_j{\phi }_j,\mathrm{i}\mathrm{n}\mathrm{\Omega }≔\left.\right]\left(a,b\right)\left[\right.,\hfill \\ {\partial }_{\mathbf{n}}{\phi }_j=0,\mathrm{a}\mathrm{t}x=a\mathrm{a}\mathrm{n}\mathrm{d}b.\hfill \end{array}\right.(2)$

These functions, which are explicitly defined by

${\lambda }_j={\left(\frac{j\pi }{b-a}\right)}^2,{\phi }_j\left(x\right)≔\sqrt{\frac{2}{b-a}}\cos \left(\frac{j\pi \left(x-a\right)}{b-a}\right),(3)$

form a complete set of orthogonal eigenfunctions. In the case of a homogeneous Dirichlet boundary condition, then one would get

${\lambda }_j={\left(\frac{\left(j+1\right)\pi }{b-a}\right)}^2,{\phi }_j\left(x\right)≔\sqrt{\frac{2}{b-a}}\sin \left(\frac{\left(j+1\right)\pi \left(x-a\right)}{b-a}\right).(4)$

Let us define $u\in {\mathcal{U}}_{\alpha /2}$ as

${\mathcal{U}}_{\alpha /2}≔\left\{u=\sum _{j=0}^{\infty }{\widehat{u}}_j{\phi }_j,{\widehat{u}}_j=<u,{\phi }_j>,\sum _{j=0}^{\infty }{\mu }_j^{\alpha }|{\widehat{u}}_j{|}^2<\infty \right\}.$

Then, if we assume that $u\in {\mathcal{U}}_{\alpha /2}$, the spectral definition of the fractional laplacian is

${\left(-{\partial }_x^2\right)}^{\alpha /2}u=\sum _{j=0}^{\infty }{\mu }_j^{\alpha }{\widehat{u}}_j{\phi }_j,$

with ${\mu }_j≔\sqrt{{\lambda }_j}$, taking the principal determination of the square-root.

The nonlinear term f is given by f′(u) = F(u). The choice of the potential function F, also called Helmholtz free-energy density, is related to the physical situation, e.g., the case of a convex quartic double-well potential or a non-convex logarithmic potential [25, 26]. Here, to fix the ideas, we consider the first choice with a variable space function, for x ∈ Ω, 0 < β(x) ≤ 1, generalizing the standard constant situation β = 1, and setting

$F\left(u\right)=\frac{1}{4\kappa \left(\mathbf{x}\right)\beta \left(\mathbf{x}\right)}{\left(\beta \left(\mathbf{x}\right)u^2-1\right)}^2,(5)$

leading to

$f\left(u\right)=F{\left(u\right)}^{\prime }=\frac{1}{\kappa \left(\mathbf{x}\right)}\left(\beta \left(\mathbf{x}\right)u^2-1\right)u.(6)$

Let us introduce the generalized energy functional

$E_{\alpha /2}\left(u\right)≔\frac{1}{2}{\left(\epsilon {\left(-\mathrm{\Delta }\right)}^{\alpha /2}u,u\right)}_{0,\mathrm{\Omega }}+{\int }_{\mathrm{\Omega }}F\left(u\right)d\mathbf{x},(7)$

with (⋅,⋅)_0,Ω as the L²(Ω)-norm. As in [21, 24], we define $\bar{u}={(-\mathrm{\Delta })}^{\alpha /2-1}u$, with 1 < α ≤ 2, which leads to the fractional gradient as ∇^α/2≔∇(−Δ)^α/2–1, and then ${\nabla }^{\alpha }u=\nabla \bar{u}$. Therefore, one gets

$E_{\alpha /2}\left(u\right)\left(t\right)≔\frac{\epsilon }{2}\Vert {\nabla }^{\alpha /2}u{\Vert }_{0,\mathrm{\Omega }}^2+{\int }_{\mathrm{\Omega }}F\left(u\right)d\mathbf{x},(8)$

setting $\Vert u{\Vert }_{0,\mathrm{\Omega }}^2≔{(u,u)}_{0,\mathrm{\Omega }}$. Based on our notations, then solving Equation 1 is equivalent to the gradient flow formulation

$\left\{\begin{array}{c}{\partial }_{t}u=-\kappa \left(\mathbf{x}\right){\nabla }_u{E}_{\alpha /2}\left(u\right)=0,\left(t,\mathbf{x}\right)\in \mathbb{R}\times \mathrm{\Omega },\hfill \\ u\left(t=0,\mathbf{x}\right)=u_0\left(\mathbf{x}\right),\mathbf{x}\in \mathrm{\Omega },\hfill \\ {\partial }_{\mathbf{n}}u\left(\mathbf{x}\right)=0,t>0,\mathbf{x}\in \mathrm{\Gamma }≔\partial \mathrm{\Omega }.\hfill \end{array}\right.(9)$

One also can prove that the energy is nonincreasing, after integration over the domain Ω, taking a test-function v = (ɛ(−Δ)^α/2u + f(u)), and since κ is a positive function

$E_{\alpha /2}\left(t\right)\le E_{\alpha /2}\left(s\right),\mathrm{f}\mathrm{o}\mathrm{r}s\le t,$

$$

3 Numerical Schemes

3.1 Minimization Algorithms

When one wants to reach the long time solution to Eq. 1, i.e., to compute the steady-state, we can directly discretize the initial boundary-value system Eq. 9 as a time-dependent problem which is a standard approach. An alternative solution consists in looking at the long-time solution to the gradient flow system Eq. 9 which cancels the energy gradient (Euler equation), based on an optimization algorithm related to the minimization of the energy functional. A standard way is to use a gradient-type method by building a sequence of minimizers u⁽ⁿ⁾ which is supposed to converge towards the steady-state, hence corresponding to the long-time dynamics solution of the gradient flow, i.e., for large values of n.

Based on this point of view, we can for example build the following simple scheme

$\left\{\begin{array}{c}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{o}\mathrm{n}n\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{l}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\hfill \\ u^{\left(n+1\right)}=u^{\left(n\right)}-{\delta }^{\left(n\right)}\kappa \nabla E_{\alpha /2}\left(u^{\left(n\right)}\right),\mathbf{x}\in \mathrm{\Omega },\hfill \\ u^{\left(0\right)}=u_0,\mathbf{x}\in \mathrm{\Omega },\hfill \\ {\partial }_{\mathbf{n}}{u}^{\left(n\right)}\left(\mathbf{x}\right)=0,\mathbf{x}\in \mathrm{\Gamma }≔\partial \mathrm{\Omega },\hfill \\ \mathrm{e}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}n\hfill \end{array}\right.(10)$

with

$\nabla E_{\alpha /2}\left(u^{\left(n\right)}\right)≔\epsilon {\left(-\mathrm{\Delta }\right)}^{\alpha /2}{u}^{\left(n\right)}+f\left(u^{\left(n\right)}\right),$

and where ρ⁽ⁿ⁾ > 0 is the local descent step of the gradient method which is computed in such a way that the energy is indeed decaying between two successive steps n and n + 1.

Let us now make a few remarks. If one considers that the step δ⁽ⁿ⁾ is a constant step, which means that we use a constant-step minimization gradient algorithm, then scheme Eq. 10 would correspond to an explicit Euler discretization in time of Eq. 1 with constant time step Δt = ρ. In the case of the computation of an optimal step, or at least an adapted gradient step ρ⁽ⁿ⁾, then one gets an explicit Euler scheme with time stepping. It is well-known that the steepest descent method with constant step is not optimal and may suffer from convergence problems, therefore requiring the computation of an optimal step. To this end, we use an optimal linear search method following (the starting value ρ = 1 can be changed)

$\left\{\begin{array}{c}\rho =1\hfill \\ \mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{u}\mathrm{t}\mathrm{e}{\varphi }^1=u^{\left(n\right)}-\rho \kappa \nabla E_{\alpha /2}\left(u^{\left(n\right)}\right)\hfill \\ \mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{u}\mathrm{t}\mathrm{e}{\varphi }^{1/2}=u^{\left(n\right)}-\frac{\rho }{2}\kappa \nabla E_{\alpha /2}\left(u^{\left(n\right)}\right)\hfill \\ \mathrm{w}\mathrm{h}\mathrm{i}\mathrm{l}\mathrm{e}{E}_{\alpha /2}\left({\varphi }^1\right)>E_{\alpha /2}\left({\varphi }^{1/2}\right)\mathrm{d}\mathrm{o}\hfill \\ \rho =\rho /2\hfill \\ {\varphi }^1={\varphi }^{1/2}\hfill \\ \mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{u}\mathrm{t}\mathrm{e}{\varphi }^{1/2}=u^{\left(n\right)}-\frac{\rho }{2}\kappa \nabla E_{\alpha /2}\left(u^{\left(n\right)}\right)\hfill \\ \mathrm{e}\mathrm{n}\mathrm{d}\mathrm{w}\mathrm{h}\mathrm{i}\mathrm{l}\mathrm{e}\hfill \\ {\rho }^{\left(k\right)}=\rho \hfill \\ u^{\left(n+1\right)}={\varphi }^{1/2}\hfill \end{array}\right.(11)$

which will also be used for all the other alternative schemes below.

Concerning the stopping criterion related to Eq. 10, one possibility consists in forcing a maximal value of n. Based on the interpretation of the steepest descent step ρ⁽ⁿ⁾ as a local time step Δt_n, this would correspond to fixing a maximal time of computation t_N for a certain a priori known value of N. Nevertheless, N is never a priori known and should be determined dynamically into the loop on n thanks to a stopping criterion. Here, we propose the following strong convergence test

$|E_{\alpha /2}^{\left(n+1\right)}-E_{\alpha /2}^{\left(n\right)}|\le {ϵ}_1,(12)$

for a very small value of ϵ₁, where $E_{\alpha /2}^{(n)}≔E_{\alpha /2}(u^{(n)})$. Another possibility is to fix that the uniform norm on Ω between two successive solutions u⁽ⁿ⁾ and u⁽ⁿ⁺¹⁾ is very small (i.e. ‖u⁽ⁿ⁺¹⁾ − u⁽ⁿ⁾‖_∞,Ω ≤ ϵ₁) but this may lead to a criterion which does not provide the minimum.

The stability and convergence properties of the steepest descent algorithm depend on the fact that it is also fully explicit. Here we consider the following alternative scheme where we implicit the energy gradient term

$\left\{\begin{array}{c}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{o}\mathrm{n}n\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{l}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\hfill \\ u^{\left(n+1\right)}=u^{\left(n\right)}-{\delta }^{\left(n\right)}\kappa \nabla E_{\alpha /2}\left(u^{\left(n+1\right)}\right),\mathbf{x}\in \mathrm{\Omega },\hfill \\ u^{\left(0\right)}=u_0,\mathbf{x}\in \mathrm{\Omega },\hfill \\ {\partial }_{\mathbf{n}}{u}^{\left(n+1\right)}\left(\mathbf{x}\right)=0,t>0,\mathbf{x}\in \mathrm{\Gamma }≔\partial \mathrm{\Omega },\hfill \\ \mathrm{e}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}n\hfill \end{array}\right.(13)$

which, in terms of PDEs can be rewritten as

$\left\{\begin{array}{c}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{o}\mathrm{n}n\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{l}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\hfill \\ u^{\left(n+1\right)}=u^{\left(n\right)}-{\delta }^{\left(n\right)}\kappa \left(\mathbf{x}\right)\left(\epsilon {\left(-\mathrm{\Delta }\right)}^{\alpha /2}{u}^{\left(n+1\right)}+f\left(u^{\left(n+1\right)}\right)\right),\mathbf{x}\in \mathrm{\Omega },\hfill \\ u^{\left(0\right)}=u_0,\mathbf{x}\in \mathrm{\Omega },\hfill \\ {\partial }_{\mathbf{n}}{u}^{\left(n+1\right)}\left(\mathbf{x}\right)=0,t>0,\mathbf{x}\in \mathrm{\Gamma }≔\partial \mathrm{\Omega },\hfill \\ \mathrm{e}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}n\hfill \end{array}\right.(14)$

The problem related to the above fully implicit (Euler) scheme is that we have to solve a nonlinear equation at each gradient step, leading to a higher computational cost if for example a fixed-point or Newton-type algorithm is used, most particularly in the framework of pseudo-spectral approximation schemes. Here, we consider the semi-implicit (Euler) (E, for short) scheme

$\left\{\begin{array}{c}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{o}\mathrm{n}n\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{l}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\hfill \\ u^{\left(n+1\right)}=u^{\left(n\right)}-{\delta }^{\left(n\right)}\kappa \left(\mathbf{x}\right)\left(\epsilon {\left(-\mathrm{\Delta }\right)}^{\alpha /2}{u}^{\left(n+1\right)}+f\left(u^{\left(n\right)}\right)\right),\mathbf{x}\in \mathrm{\Omega },\hfill \\ u^{\left(0\right)}=u_0,\mathbf{x}\in \mathrm{\Omega },\hfill \\ {\partial }_{\mathbf{n}}{u}^{\left(n+1\right)}\left(\mathbf{x}\right)=0,t>0,\mathbf{x}\in \mathrm{\Gamma }≔\partial \mathrm{\Omega },\hfill \\ \mathrm{e}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}n\hfill \end{array}\right.(15)$

which avoids any nonlinear solution and leads to solving the linear PDE formulation

$\left\{\begin{array}{c}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{o}\mathrm{n}n\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{l}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\hfill \\ \left(\frac{1}{{\delta }^{\left(n\right)}}I+\kappa \left(\mathbf{x}\right)\epsilon {\left(-\mathrm{\Delta }\right)}^{\alpha /2}\right)u^{\left(n+1\right)}=\frac{1}{{\delta }^{\left(n\right)}}{u}^{\left(n\right)}-\kappa \left(\mathbf{x}\right)f\left(u^{\left(n\right)}\right),\mathbf{x}\in \mathrm{\Omega },\hfill \\ u^{\left(0\right)}=u_0,\mathbf{x}\in \mathrm{\Omega },\hfill \\ {\partial }_{\mathbf{n}}{u}^{\left(n+1\right)}\left(\mathbf{x}\right)=0,t>0,\mathbf{x}\in \mathrm{\Gamma }≔\partial \mathrm{\Omega },\hfill \\ \mathrm{e}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}n\hfill \end{array}\right.(16)$

where I is the identity operator. In fact, it can be seen as the formulation of the fully-implicit (Euler) scheme Eq. 15 with one iteration of the fixed-point algorithm. If one uses the decomposition of the Ginzburg–Landau free energy into a linear kinetic part and a nonlinear part related to the potential function F

$E_{\alpha /2}\left(u\right)≔E_{\alpha /2}^{\text{kin}}\left(u\right)+E_{\alpha /2}^F\left(u\right),(17)$

setting

$\left\{\begin{array}{c}{E}_{\alpha /2}^{\text{kin}}\left(u\right)≔\frac{\epsilon }{2}\Vert {\nabla }^{\alpha /2}u{\Vert }_{0,\mathrm{\Omega }}^2,\hfill \\ E_{\alpha /2}^F\left(u\right)≔{\int }_{\mathrm{\Omega }}F\left(u\right)d\mathbf{x},\hfill \end{array}\right.(18)$

then Eq. 15 can also be seen as an algorithm where the kinetic energy is implicitly approximated while the potentiel energy is explicit resolved

$\left\{\begin{array}{c}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{o}\mathrm{n}n\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{l}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\hfill \\ u^{\left(n+1\right)}=u^{\left(n\right)}-{\delta }^{\left(n\right)}\kappa \left(\mathbf{x}\right)\nabla E_{\alpha /2}^{\text{kin}}\left(u^{\left(n+1\right)}\right)-{\delta }^{\left(n\right)}\kappa \left(\mathbf{x}\right)\nabla E_{\alpha /2}^F\left(u^{\left(n\right)}\right),\mathbf{x}\in \mathrm{\Omega },\hfill \\ u^{\left(0\right)}=u_0,\mathbf{x}\in \mathrm{\Omega },\hfill \\ {\partial }_{\mathbf{n}}{u}^{\left(n+1\right)}\left(\mathbf{x}\right)=0,t>0,\mathbf{x}\in \mathrm{\Gamma }≔\partial \mathrm{\Omega },\hfill \\ \mathrm{e}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}n\hfill \end{array}\right.(19)$

An alternative to the semi-implicit Euler scheme is to use a semi-implicit Crank-Nicolson (CN, for brevity) scheme which writes

$\left\{\begin{array}{c}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{o}\mathrm{n}n\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{l}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\hfill \\ \left(\frac{2}{{\delta }^{\left(n\right)}}I+\kappa \left(\mathbf{x}\right)\epsilon {\left(-\mathrm{\Delta }\right)}^{\alpha /2}\right)u^{\left(n+1/2\right)}=\frac{2}{{\delta }^{\left(n\right)}}{u}^{\left(n\right)}-\frac{\kappa }{2}f\left(u^{\left(n\right)}\right),\mathbf{x}\in \mathrm{\Omega },\hfill \\ u^{\left(0\right)}=u_0,\mathbf{x}\in \mathrm{\Omega },\hfill \\ {\partial }_{\mathbf{n}}{u}^{\left(n+1/2\right)}\left(\mathbf{x}\right)=0,t>0,\mathbf{x}\in \mathrm{\Gamma }≔\partial \mathrm{\Omega },\hfill \\ \mathrm{e}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}n\hfill \end{array}\right.(20)$

and is similar to Eq. 19 but evaluating the kinetic energy at the mid-point u^n+1/2≔(u⁽ⁿ⁺¹⁾ + u⁽ⁿ⁾)/2 as well as the Neumann boundary condition. For both the E and CN schemes, the line search algorithm Eq. 11 is adapted thanks to the formulation.

3.2 Pseudo-Spectral Discretization in Space and Linear System Solution

In a practical computation, we approximate the solution u by a J-terms finite-dimensional approximation $u_J≔{\sum }_{j=0}^J{\widehat{u}}_j{\phi }_j$, $J\in \mathbb{N}$. To this end, the fast cosine/sine transform can be used for computing the fractional laplacian when Neumann/Dirichlet boundary conditions are used. For the case of periodic boundary conditions, FFT-based methods are applied.

When one uses a fast transform (cosine, sine, FFT), we can proceed to the efficient computation of the application of an operator to a given function. For example, let us consider that one wants to numerically apply the fractional laplacian of order α to a given function v, i.e., we wish to compute w≔(−Δ)^α/2v, where v is a given field known by its values v_j at the nodes x_j≔ − L/2 + (2j + 1)h/2, for 0 ≤ j ≤ J − 1 and where h = L/J denotes the uniform spatial meshsize (L = b − a). Then, we get the practical efficient calculation of w by the cosine (dct) and inverse cosine (idct) transforms in one-dimension through the Matlab code

$\begin{array}{c}\mathtt{L}=\mathtt{b}-\mathtt{a};\hfill \\ \mathtt{J}=\mathtt{2}\mathtt{^}\mathtt{8};\hfill \\ \mathtt{h}=\mathtt{L}/\mathtt{J};\hfill \\ \mathtt{l}\mathtt{a}\mathtt{m}\mathtt{b}\mathtt{d}\mathtt{a}\mathtt{j}=\left(\left(\mathtt{0}:\mathtt{J}-\mathtt{1}\right)\ast \mathtt{p}\mathtt{i}/\mathtt{L}\right)\mathtt{^}\mathtt{2}{.}^{\prime };\hfill \\ \mathtt{w}=\mathtt{i}\mathtt{d}\mathtt{c}\mathtt{t}\left(\mathtt{l}\mathtt{a}\mathtt{m}\mathtt{b}\mathtt{d}\mathtt{a}\mathtt{j}.\mathtt{^}\left(\mathtt{a}\mathtt{l}\mathtt{p}\mathtt{h}\mathtt{a}/\mathtt{2}\right).\ast \mathtt{d}\mathtt{c}\mathtt{t}\left(\mathtt{v}\right)\right);\hfill \end{array}(21)$

For all the previous schemes presented in section 3.1, and for a constant function κ, the finite-dimensional solution in space can be obtained directly thanks to the above pseudo-spectral discretization. Indeed, even for the semi-implicit schemes, the operators to invert, i.e. $({({\delta }^{(n)})}^{-1}I+\kappa \epsilon {(-\mathrm{\Delta })}^{\alpha /2})$ for the E schemes, can be computed directly through their symbols. For the Neumann boundary conditions, this means that we can apply the discrete symbols of ${({({\delta }^{(n)})}^{-1}I+\kappa \epsilon {(-\mathrm{\Delta })}^{\alpha /2})}^{-1}$ in the dual space thanks to its spectral decomposition involving the eigenvalues ${({({\delta }^{(n)})}^{-1}I+\kappa \epsilon {\mu }_j^{\alpha })}^{-1}$, 0 ≤ j ≤ J − 1, and in a similar way to Eq. 21.

However, for the more complex situation that we consider in this paper, this is no longer possible. Indeed, a variable coefficient operator cannot be inverted explicitly in the Fourier space. Let us for example consider the situation of a full space representation (no boundary condition). Then, we can use the inverse Fourier representation of the operator $({({\delta }^{(n)})}^{-1}I+\kappa (\mathbf{x})\epsilon {(-\mathrm{\Delta })}^{\alpha /2})$ following

$\begin{array}{c}w=A\left(\mathbf{x},D_{\mathbf{x}}\right)v≔\left({\left({\delta }^{\left(n\right)}\right)}^{-1}I+\kappa \left(\mathbf{x}\right)\epsilon {\left(-\mathrm{\Delta }\right)}^{\alpha /2}\right)v\hfill \\ ={\mathcal{F}}^{-1}\left(a\left(\mathbf{x},\mathit{\xi }\right)\mathcal{F}\left(u\right)\left(\mathit{\xi }\right)\right),\hfill \end{array}(22)$

where $\mathit{\xi }\in {\mathbb{R}}^d$ is the dual Fourier variable, the direct and inverse Fourier transform being denoted by $\mathcal{F}$ and ${\mathcal{F}}^{-1}$, respectively, and defined by

$\mathcal{F}\left(f\right)\left(\mathit{\xi }\right)=\widehat{f}\left(\mathit{\xi }\right)={\int }_{{\mathbb{R}}^d}{e}^{-i\mathit{x}\cdot \mathit{\xi }}f\left(\mathit{x}\right)d\mathit{x},$

and

${\mathcal{F}}^{-1}\left(f\right)\left(\mathit{x}\right)=\frac{1}{{\left(2\pi \right)}^d}{\int }_{{\mathbb{R}}^d}{e}^{i\mathit{x}\cdot \mathit{\xi }}\widehat{f}\left(\mathit{\xi }\right)d\mathit{\xi }.$

The symbol of the pseudodifferential operator [27] is given here by the expression

$a\left(\mathbf{x},\mathit{\xi }\right)≔\left({\left({\delta }^{\left(n\right)}\right)}^{-1}+\kappa \left(\mathbf{x}\right)\epsilon \Vert \mathit{\xi }{\Vert }^{\alpha }\right).(23)$

The inner product of two vectors p and q in ${\mathbb{R}}^d$ is such that: $\mathbf{p}\cdot \mathbf{q}≔{\sum }_{k=1}^d{p}_k{q}_k$, and the associated norm is $\Vert \mathbf{p}\Vert ≔\sqrt{\mathbf{p}\cdot \mathbf{p}}$. The symbol a of the operator A, which is of order α, then generates an operator of order α acting between a Sobolev space of order m onto m − α. Let us define B as the (left) pseudo-inverse operator of A, i.e., BA = I modulo a smoothing operator. The operator B can be computed recursively through an asymptotic expansion of its total symbol b≔b(x, ξ) by the formula [27].

${\sigma }_{BA}=b♯a\left(\mathbf{x},\mathit{\xi }\right)\sim \sum _{\mathit{\eta }\in {\mathbb{N}}^d}\frac{1}{\mathit{\eta }!}{\partial }_{\mathit{\xi }}^{\mathit{\eta }}b\left(\mathbf{x},\mathit{\xi }\right)D_{\mathbf{x}}^{\mathit{\eta }}a\left(\mathbf{x},\mathit{\xi }\right)=1,(24)$

where σ_BA is the symbol of the operator BA. The operator D^η is: $D^{\mathit{\eta }}=D_1^{{\eta }_1}\dots ,D_d^{{\eta }_d}$, D = (D₁, … , D_d), setting $D_k≔-i{\partial }_k=-i{\partial }_{x_k}$, k = 1, … , d, and $i≔\sqrt{-1}$. The multi-index η is defined as a vector η≔(η₁, … , η_d) with d positive integer components, i.e. $\mathit{\eta }\in {\mathbb{N}}^d$. The factorial is: $\mathit{\eta }!={\mathrm{\Pi }}_{k=1}^d{\eta }_k!$. Thanks to these notations, one gets: ${\sigma }_{{(-\mathrm{\Delta })}^{\alpha /2}}=\Vert \mathit{\xi }{\Vert }^{\alpha }$. If κ is a positive constant, then the composition Eq. 24 leads to

$b\left(\mathbf{x},\mathit{\xi }\right)=b\left(\mathit{\xi }\right)=a{\left(\mathbf{x},\mathit{\xi }\right)}^{-1}=a{\left(\mathit{\xi }\right)}^{-1},$

since a does not depend on x and then the infinite sum Eq. 24 reduces to only one term for η = 0. When a is x-dependent, then the asymptotic expansion is infinite and b includes an infinite number of terms. The operator B_−α/2 with symbol b_−α/2≔a(x,ξ)⁻¹ is a regularizing pseudodifferential operator of order − α which is not the inverse of A. Indeed, we have: B_−α/2A = I + R, where R is a pseudodifferential operator of negative order which is however usually not an infinitely smoothing operator.

Now, if we come back to Eq. 22, we see that we can rewrite it as

$\begin{array}{c}w=A\left(\mathbf{x},D_{\mathbf{x}}\right)v≔\left({\left({\delta }^{\left(n\right)}\right)}^{-1}I+\kappa \left(\mathbf{x}\right)\epsilon {\left(-\mathrm{\Delta }\right)}^{\alpha /2}\right)v\hfill \\ ={\left({\delta }^{\left(n\right)}\right)}^{-1}v+\kappa \left(\mathbf{x}\right)\epsilon {\mathcal{F}}^{-1}\left(\Vert \mathit{\xi }{\Vert }^{\alpha }\widehat{v}\left(\mathit{\xi }\right)\right),\hfill \end{array}(25)$

which can be discretized directly thanks to a FFT, fast cosine or sine transform according to the boundary conditions. Let us write the corresponding discrete operation as

$\mathbf{w}≔\left[\left[A_n\right]\right]\mathbf{v}≔\left({\left({\delta }^{\left(n\right)}\right)}^{-1}I+\epsilon \mathit{\kappa }\left[\left[{\left(-\mathrm{\Delta }\right)}^{\alpha /2}\right]\right]\right)\mathbf{v},(26)$

where v and w denotes the given and unknown vectors, respectively, at the discretization points, I is the identity matrix of size J, κ is the diagonal matrix with the entries κ(x_j), for 0 ≤ j ≤ J − 1, and [[(−Δ)^α/2]]v denotes the application of the fast algorithm for evaluating Av at the discrete level. Therefore, [[A]] is a matrix-free operator given through a fast evaluation procedure (with e.g., a function call $\mathtt{@}(\mathtt{v})\mathtt{A}(\mathtt{v})$ in Matlab) requiring $\mathcal{O}(J\log J)$ operations and a memory cost of $\mathcal{O}(J)$ entries.

As previously seen, using the E or CN approximations needs the resolution of a linear system. When κ is a constant function, this can be trivially done thanks to the explicit inversion of the operator through the inverse of its symbol. For a function κ depending on x, then the process for solving one of the linear systems is to use an iterative solver instead of a direct method. Among them, fixed-point methods [28] (Jacobi, Gauss-Seidel, successive approximation) can be used. Nevertheless, it is well-known that they are not robust since their convergence rate strongly depends on the spectral radius of the iteration matrix. More efficient methods are based on Krylov subspace iterative solvers [28, 29]. We choose here the Generalized Minimal Residual (GMRES) method [29] which is extremely robust (other methods, like e.g., BiCGStab, have been tested but are proved to be less efficient within the framework of this paper; we do not report the results for the sake of conciseness). This method can solve matrix-free linear systems by calling an argument which is the function $\mathtt{@}(\mathtt{v})\mathtt{A}(\mathtt{v})$ related to the matrix-vector product w≔[[A_n]]v. The convergence rate of the GMRES is related to the spectral distribution of the matrix of the linear system and the cost of applying it to our problem is therefore $\mathcal{O}(n_n^{iter}J\log J)$, where $n_n^{iter}$ is the number of iterations required for solving the linear system with a tolerance ϵ₂. To accelerate the computations, one usually includes a preconditioner [[B_n]] which is an approximation of the inverse of [[A_n]]. From the above discussion, a natural choice is

$\mathbf{v}≔\left[\left[B_n\right]\right]\mathbf{w},(27)$

where [[B_n]] is the discretization of the operator with symbol $a_n^{-1}$. This operator is the inverse of [[A_n]] if κ is a constant but is only an approximation for a variable function κ. Its application has a cost of $\mathcal{O}(J\log J)$ operations and the memory cost if $\mathcal{O}(J)$. As it will be seen later during the numerical examples, the use of [[B_n]] may strongly accelerate the convergence rate. Concretely, for the one-dimensional case and the cosine transform (Neumann boundary condition), the Matlab procedure is simply

$\mathtt{v}=\mathtt{i}\mathtt{d}\mathtt{c}\mathtt{t}\left(\mathtt{1}./\left(\mathtt{1}/\mathtt{d}\mathtt{e}\mathtt{l}\mathtt{t}\mathtt{a}\mathtt{n}+\mathtt{k}\mathtt{a}\mathtt{p}\mathtt{p}\mathtt{a}.\ast \mathtt{e}\mathtt{p}\mathtt{s}\mathtt{i}\mathtt{l}\mathtt{o}\mathtt{n}.\ast \mathtt{l}\mathtt{a}\mathtt{m}\mathtt{b}\mathtt{d}\mathtt{a}\mathtt{j}.\mathtt{^}\left(\mathtt{a}\mathtt{l}\mathtt{p}\mathtt{h}\mathtt{a}/\mathtt{2}\right)\right).\ast \mathtt{d}\mathtt{c}\mathtt{t}\left(\mathtt{w}\right)\right);(28)$

where $\mathtt{k}\mathtt{a}\mathtt{p}\mathtt{p}\mathtt{a}$ is encoded as a vector representing the diagonal matrix. In the case of a constant parameter κ, the GMRES equipped with this preconditioner solves the linear system in one iteration (as a direct solver would do). A useful alternative preconditioner is $[[B_n^m]]$, which is based on the mean ${\kappa }^m=J^{-1}{\sum }_{j=0}^{J-1}{\kappa }_j$. In practice, $[[B_n^m]]$ is much more robust than [[B_n]].

4 Computational Aspects and Numerical Example

4.1 Computational Aspects

In the following, we call E and CN the semi-implicit Euler and Crank-Nicolson schemes Eqs 16, 20, respectively, without any preconditioner for the GMRES. When the preconditioner $[[B_n^m]]$ is used in the GMRES method, the algorithms are called PE and PCN. The GMRES is considered without restarting and with a tolerance ɛ₂ = 10^–12. For each iteration n of the minimization method, the number of iterations of the GMRES is then called $N_{\alpha /2,J}^{\text{m},n}$ (m = E, CN, PE, or PCN) accordingly to the method, fractional order α/2 and number of spatial discretization points J. The total computational cost at iteration n is essentially related to the number of Matrix-Vector (MV) operations based on [[A_n]], and is twice when the preconditioner $[[B_n^m]]$ (or [[B_n]]) is applied since it requires an additional function call. Concerning the minimization algorithms, the stopping criterion Eq. 12 is fixed through ϵ₁ = 10^–12 and used with the linear line search algorithm. If we call ${\bar{N}}_{\alpha /2,J}^{m,n}$ the total number of MV products for the whole procedure and for a method m at step n, then we have

${\bar{N}}_{\alpha /2,J}^{\text{m},n}≔\sum _{\ell =1}^n{N}_{\alpha /2,J}^{\text{m},\ell },{\bar{N}}_{\alpha /2,J}^{\text{m}}≔{\bar{N}}_{\alpha /2,J}^{\text{m},n_{\alpha /2,J}^{\text{m}}}=\sum _{n=1}^{n_{\alpha /2,J}^{\text{m}}}{N}_{\alpha /2,J}^{\text{m},n},$

where $n_{\alpha /2,J}^{\text{m}}$ is the iteration number n to reach the tolerance in the minimization algorithm for the method m. As a consequence, the algorithm has a computational cost of $\mathcal{O}({\bar{N}}_{\alpha /2,J}^{\text{m}}J\log J)$ operations and needs a storage of $\mathcal{O}(J)$ coefficients. In the case of a higher-dimensional problem, J denotes the total number of unknowns on the d-grid.

Let us introduce the discrete L²-norm as: $\forall \mathbf{v}\in {\mathbb{R}}^J,\Vert \mathbf{v}{\Vert }_0≔h^{1/2}\Vert \mathbf{v}\Vert$, with the associated discrete inner product: (u,v)₀≔hu ⋅v. Then, the converged discrete energy at level J is given by

$E^{\left(n_{\alpha /2,J}^{\text{m}}\right)}≔E_{\alpha /2,J}\left({\mathbf{u}}^{\left(n_{\alpha /2,J}^{\text{m}}\right)}\right)≔\frac{\epsilon }{2}\Vert \left[\left[{\nabla }^{\alpha /2}\right]\right]{\mathbf{u}}^{\left(n_{\alpha /2,J}^{\text{m}}\right)}{\Vert }_0^2+{\left(F\left({\mathbf{u}}^{\left(n_{\alpha /2,J}^{\text{m}}\right)}\right),\mathbf{1}\right)}_0.$

The converged solution for a method m, the fractional PDE of order α/2, a number of spatial grid points J and a number of iterations $n_{\alpha /2,J}^{\text{m}}$ is called ${\mathbf{u}}^{(n_{\alpha /2,J}^{\text{m}})}$. Finally, the residual energy of a method m at iteration n for a fractional order α/2 and spatial discretization J is given by

$\mathrm{\Delta }{E}_{\alpha /2,J}^{\text{m},n}≔|E_{\alpha /2,J}\left({\mathbf{u}}^{m,\left(n\right)}\right)-E_{\alpha /2,J}\left({\mathbf{u}}^{m,\left(n-1\right)}\right)|.(29)$

Therefore, representing $\mathrm{\Delta }{E}_{\alpha /2,J}^{\text{m},n}$ vs. ${\bar{N}}_{\alpha /2,J}^{\text{m},n}$ provides the correct way of measuring the convergence speed.

4.2 A One-Dimensional Example

To illustrate the methods, we consider an example in dimension one. The computational domain is $\mathrm{\Omega }=\left.\right]-\mathrm{1,1}\left[\right.$. The initial data is given by

$u_0\left(x\right)=0.5\sin \left(3\pi x/2\right)\left(\cos \left(\pi x\right)-1\right).$

The parameter ɛ is fixed to 10^–2 and β is the function

$\beta \left(x\right)≔0.6+0.4\cos \left(20\pi x\right).$

Concerning the variable function κ, we consider

$\kappa \left(x\right)≔1+0.5\cos \left(10\pi x\right).(30)$

Therefore, κ^m = 1 while κ oscillates. The number of grid points is set to J = 2^p, for $p\in \mathbb{N}$, leading to h = 2^−p+1.

Let us first report on Figure 1 (left) the computed solutions ${\mathbf{u}}^{(n_{1,J}^{\text{E}})}$ when considering κ and $\bar{\kappa }$ based on the Euler formulation Eq. 16, and the solution obtained by the direct inversion of the operator with the inverse symbol of a, i.e. a⁻¹, and the variable function κ. Here, the situation considers the case α = 2 which is the integer order CA equation. The number of grid points is J = 2¹⁰. The solution based on the inverted symbol does not provide a spatially converging solution, which is expected from the previous explanations. As seen on the figure, considering κ^m is not enough to get the true solution, but it provides an idea of the solution. This means that a preconditioner based on κ^m is expected to be robust when included in the GMRES. We now report the converged solution for various fractional orders α = 2, 1.75, 1.5, 1.25, 1.01 on Figure 1 (right). The effect of the oscillating function κ can be clearly observed, and the solution tends to be more singular as α goes to 1. The respective converged energies are: 76.56, 64.90, 53.01, 38.64 and 25.18.

FIGURE 1

FIGURE 1. 1d case: (A) (α = 2): converged solutions for 1) the variable coefficients case with the unpreconditioned Euler scheme and variable function κ, 2) with the direct inversion of the operator by using the variable function κ and 3) with the Euler scheme by using the average constant κ^m. (B): converged solutions for various values of the fractional order α.

We analyze now the behavior of the global iterative method on Figure 2. The spatial discretization parameter is J = 2¹⁰ for the same function cases as before, but with α = 1.5 and α = 1.9. On Figure 2, we report the residual history for the E, PE, CN and PCN methods, the preconditioner being based on the inverse of the operator with average symbol κ^m. As it can be seen, the semi-implicit Euler iterative scheme converges faster than the CN method, while preconditioning helps a lot for both methods by dividing the number of iterations by a factor 5–6 (α = 1.5) to 20 (α = 1.9). When α is close to 1, the preconditioner effect is less visible while it is more important for α close to 2. This can be understood by the fact that the condition number is then larger as α increases, leading to a need for preconditioning (−Δ)^α/2. In addition, the preconditioner helps a lot for the convergence when refining the spatial discretization thanks to J, i.e. to get a higher accuracy. As a conclusion, and from other simulations in different configurations (not reported here for the sake of conciseness), it appears that the PE method always provides the best convergence rate.

FIGURE 2

FIGURE 2. 1d case: (A) (α = 1.5) and (B) (α = 1.9): energy residual history of the E, PE, CN and PCN methods.

5 Conclusion

In this paper, we introduced some full discretizations of the gradient flow formulation of space fractional Cahn-Allen equations with variable coefficients in space. To this end, we use an Euler or Crank-Nicolson time discretization of the gradient flow formulation, with adaptive time stepping. In addition, a Fourier-based pseudo-spectral discretization scheme in space allows for an efficient and robust computation of each time step, in conjunction with a GMRES solver accelerated by a spectral preconditioner. Some simulations in dimension one show that the method based on the Euler scheme with the spectral preconditioner performs the best. This conclusion should extend to 2D/3D problems since all the ingredients here can be adapted and higher-dimensional grids are cartesian extensions of the one-dimensional case. In some future works, we will address more numerical examples in 2D and an extension of the method to the fractional Cahn-Hillard equation [30, 31] that is used for example in [32] as fractional inpainting model.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding authors.

Author Contributions

SA and CC developed the theoretical formalism, performed the analytic calculations and performed the numerical simulations. Authors contributed equally to the final version of the manuscript.

Funding

The authors would like to thank the Deanship of Scientific Research at Umm Al-Qura University for supporting this work by Grand code 18-SCI-1-01-0017.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Karniadakis GE, Hesthaven JS, Podlubny I. Special Issue on "Fractional PDEs: Theory, Numerics, and Applications". J Comput Phys (2015) 293:1–3. doi:10.1016/j.jcp.2015.04.007

CrossRef Full Text | Google Scholar

2. Podlubny I. Fractional Differential Equations. San Diego: Academic Press (1999).

Google Scholar

3. Baleanu D, Diethelm K, Scalas E, Trujillo JJ. Fractional Calculus, Models and Numerical Methods. Series on Complexity, Nonlinearity and Chaos, 3. Boston, Mass: World Scientific (2012).

Google Scholar

4. Herrmann R. Fractional Calculus, An Introduction for Physicists. GigaHedron, Germany: World Scientific (2011).

Google Scholar

5. Metzler R, Klafter J. The Random Walk's Guide to Anomalous Diffusion: a Fractional Dynamics Approach. Phys Rep (2000) 339(1):1–77. doi:10.1016/s0370-1573(00)00070-3

CrossRef Full Text | Google Scholar

6. Carreras BA, Lynch VE, Zaslavsky GM. Anomalous Diffusion and Exit Time Distribution of Particle Tracers in Plasma Turbulence Model. Phys Plasmas (2001) 8:5096–103. doi:10.1063/1.1416180

CrossRef Full Text | Google Scholar

7. Shlesinger MF, West BJ, Klafter J. Lévy Dynamics of Enhanced Diffusion: Application to Turbulence. Phys Rev Lett (1987) 58:1100–3. doi:10.1103/physrevlett.58.1100

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Zaslavsky GM, Stevens D, Weitzner H. Self-similar Transport in Incomplete Chaos. Phys Rev E (1993) 48:1683–94. doi:10.1103/physreve.48.1683

CrossRef Full Text | Google Scholar

9. Meerschaert MM, Benson DA, Wheatcraft SW. Subordinated Advection-Dispersion Equation for Contaminant Transport. Water Resour Res (2001) 37:1543–50. doi:10.1029/2000WR900409

CrossRef Full Text | Google Scholar

10. Benson DA, Wheatcraft SW, Meerschaert MM. The Fractional-Order Governing Equation of Lévy Motion. Water Resour Res (2000) 36:1413–23. doi:10.1029/2000wr900032

CrossRef Full Text | Google Scholar

11. Scalas E, Gorenflo R, Mainardi F. Fractional Calculus and Continuous-Time Finance. Phys A: Stat Mech Appl (2000) 284:376–84. doi:10.1016/s0378-4371(00)00255-7

CrossRef Full Text | Google Scholar

12. Yin B, Liu Y, Li H, He S. Fast Algorithm Based on TT-M FE System for Space Fractional Allen-Cahn Equations with Smooth and Non-smooth Solutions. J Comput Phys (2019) 379:351–72. doi:10.1016/j.jcp.2018.12.004

CrossRef Full Text | Google Scholar

13. Ullah MS, Harun-Or-Roshid , Ali MZ, Noor NFM. Novel Dynamics of Wave Solutions for Cahn-Allen and Diffusive Predator-Prey Models Using MSE Scheme. Partial Differ Equations Appl Maths (2021) 3:100017. doi:10.1016/j.padiff.2020.100017

CrossRef Full Text | Google Scholar

14. Allen SM, Cahn JW. A Microscopic Theory for Antiphase Boundary Motion and its Application to Antiphase Domain Coarsening. Acta Metal (1979) 27:1085–95. doi:10.1016/0001-6160(79)90196-2

CrossRef Full Text | Google Scholar

15. Algahtani OJJ. Comparing the Atangana-Baleanu and Caputo-Fabrizio Derivative with Fractional Order: Allen Cahn Model. Chaos Solitons Fractals (2016) 89:552–9. doi:10.1016/j.chaos.2016.03.026

CrossRef Full Text | Google Scholar

16. Nec Y, Nepomnyashchy AA, Golovin AA. Front-type Solutions of Fractional Allen-Cahn Equation. Phys D: Nonlinear Phenomena (2008) 237(24):3237–51. doi:10.1016/j.physd.2008.08.002

CrossRef Full Text | Google Scholar

17. Lee S, Lee D. The Fractional Allen-Cahn Equation with the Sextic Potential. Appl Maths Comput (2019) 351:176–92. doi:10.1016/j.amc.2019.01.037

CrossRef Full Text | Google Scholar

18. Benes M, Chalupecky V, Mikula K. Geometrical Image Segmentation by the Allen-Cahn Equation. Appl Numer Maths (2004) 51(2-3):187–205.

Google Scholar

19. Lee D, Lee S. Image Segmentation Based on Modified Fractional Allen-Cahn Equation. Math Probl Eng (2019) 2019:3980181. doi:10.1155/2019/3980181

CrossRef Full Text | Google Scholar

20. Zhai S, Gui D, Zhao J, Feng X. High Accuracy Spectral Method for the Space-Fractional Diffusion Equations. JMS (2014) 47(3):274–86. doi:10.4208/jms.v47n3.14.03

CrossRef Full Text | Google Scholar

21. Bueno-Orovio A, Kay D, Burrage K. Fourier Spectral Methods for Fractional-In-Space Reaction-Diffusion Equations. Bit Numer Math (2014) 54(4):937–54. doi:10.1007/s10543-014-0484-2

CrossRef Full Text | Google Scholar

22. Zhai S, Weng Z, Feng X. Fast Explicit Operator Splitting Method and Time-step Adaptivity for Fractional Non-local Allen-Cahn Model. Appl Math Model (2016) 40(2):1315–24. doi:10.1016/j.apm.2015.07.021

CrossRef Full Text | Google Scholar

23. Alzahrani SS, Khaliq AQM. Fourier Spectral Exponential Time Differencing Methods for Multi-Dimensional Space-Fractional Reaction-Diffusion Equations. J Comput Appl Maths (2019) 361:157–75. doi:10.1016/j.cam.2019.04.001

CrossRef Full Text | Google Scholar

24. Lischke A, Pang G, Gulian M, Song F, Glusa C, Zheng X, et al. What Is the Fractional Laplacian? A Comparative Review with New Results. J Comput Phys (2020) 404:109009. doi:10.1016/j.jcp.2019.109009

CrossRef Full Text | Google Scholar

25. Choi J-W, Lee HG, Jeong D, Kim J. An Unconditionally Gradient Stable Numerical Method for Solving the Allen-Cahn Equation. Physica A: Stat Mech its Appl (2009) 388:1791–803. doi:10.1016/j.physa.2009.01.026

CrossRef Full Text | Google Scholar

26. Dai S, Du Q. Motion of Interfaces Governed by the Cahn--Hilliard Equation with Highly Disparate Diffusion Mobility. SIAM J Appl Math (2012) 72(6):1818–41. doi:10.1137/120862582

CrossRef Full Text | Google Scholar

27. Taylor ME. Pseudodifferential Operators. Princeton, NJ: Princeton University Press (1981).

Google Scholar

28. Saad Y. Iterative Methods for Sparse Linear Systems. Boston: PWS Pub. Co. (1996).

Google Scholar

29. Saad Y, Schultz MH. GMRES: a Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM J Sci Stat Comput (1986) 7(3):856–69. doi:10.1137/0907058

CrossRef Full Text | Google Scholar

30. Ainsworth M, Mao Z. Analysis and Approximation of a Fractional Cahn--Hilliard Equation. SIAM J Numer Anal (2017) 55(4):1689–718. doi:10.1137/16m1075302

CrossRef Full Text | Google Scholar

31. Weng Z, Zhai S, Feng X. A Fourier Spectral Method for Fractional-In-Space Cahn-Hilliard Equation. Appl Math Model (2017) 42:462–77. doi:10.1016/j.apm.2016.10.035

CrossRef Full Text | Google Scholar

32. Bosch J, Stoll M. A Fractional Inpainting Model Based on the Vector-Valued Cahn--Hilliard Equation. SIAM J Imaging Sci (2015) 8:2352–82. doi:10.1137/15m101405x

CrossRef Full Text | Google Scholar