Some comments on using fractional derivative operators in modeling non-local diffusion processes

Abstract

We start with a general governing equation for diffusion transport, written in a conserved form, in which the phenomenological flux laws can be constructed in a number of alternative ways. We pay particular attention to flux laws that can account for non-locality through space fractional derivative operators. The available results on the well posedness of the governing equations using such flux laws are discussed. A discrete control volume numerical solution of the general conserved governing equation is developed and a general discrete treatment of boundary conditions, independent of the particular choice of flux law, is presented. The numerical properties of the scheme resulting from the flux laws are analyzed. We use numerical solutions of various test problems to compare the operation and predictive ability of two discrete fractional diffusion flux laws based on the Caputo (C) and Riemann–Liouville (RL) derivatives respectively. When compared with the C flux-law we note that the RL flux law includes an additional term, that, in a phenomenological sense, acts as an apparent advection transport. Through our test solutions we show that, when compared to the performance of the C flux-law, this extra term can lead to RL-flux law predictions that may be physically and mathematically unsound. We conclude, by proposing a parsimonious definition for a fractional derivative based flux law that removes the ambiguities associated with the selection between non-local flux laws based on the RL and C fractional derivatives.

Introduction

A non-local transport phenomenon can be defined as one where the transport flux at given point depends not only on values at that point but also on values at remote locations. For example, a non-local diffusion flux can be defined as a weighted sum of potential gradients evaluated at points throughout the domain. Due to their non-local nature, many works (e.g., see [1], [2], [3] have been directed at using fractional derivative operators to construct non-local transport models in physical settings. Recent work [4] studied the numerical solution of the diffusion equation

$${\partial }_{t}u(x,t)=\kappa L^{1+\alpha }u(x,t)+f,0<\alpha \le 1,x\in \Omega ,t>0,$$

where $\Omega$ is the open interval $(0,1)$, the operator $L^{1+\alpha }$ is a fractional derivative of order $1<1+\alpha \le 2$, $\kappa$ is an appropriately scaled diffusivity, and $f$ is a source/sink term. By numerically predicting the fate of an initial non-negative pulse $u(x,0)\ge 0$ on $\Omega$, the authors of [4] investigate three alternatives for the fractional operator $L^{1+\alpha }$, a Riemann–Liouville (RL${}^{1+\alpha }$) derivative, a Caputo (C${}^{1+\alpha }$) derivative, and a so called Patie–Simon (PS${}^{1+\alpha }$) derivative. A major and important contribution from the work of [4], is the development of a comprehensive and consistent treatment of boundary conditions for numerical solution of the fractional heat conduction equation. The work shows, under a variety of boundary and initial conditions, that numerical solutions using the RL${}^{1+\alpha }$ and PS${}^{1+\alpha }$ derivatives, while not-identical, are always conservative. This, however, is not the case for the C${}^{1+\alpha }$ derivative, under some boundary and initial conditions this operator may not be conservative. The authors of [4] demonstrate this by considering a problem with fixed boundary values, $u(0,t)=u(1,t)=0$, presenting a numerical solution, based on the C${}^{1+\alpha }$ derivative, where an initial non-negative pulse $u(x,0)$ produces later time solutions $u(x,t)$ with negative values.

An alternative form for the fractional diffusion equation is

$${\partial }_{t}u(x,t)={\partial }_x\left(\kappa L^{\alpha }u(x,t)\right)+f(x,t),0<\alpha \le 1,x\in \Omega ,t>0,$$

where we interpret the right-hand side as the conserved statement the divergence of a flux $q=L^{\alpha }u(x,t)$. We note that in using the equation form suggested by [4], i.e., Eq. (1), conservation will only by ensured if ${\partial }_x\left(L^{\alpha }u(x,t)\right)\equiv L^{\alpha +1}u(x,t)$. This is true in the case the operator $L^{\alpha }$ is the RL${}^{\alpha }$ derivative of order $0<\alpha \le 1$. For Caputo derivatives, however, the condition that ${\partial }_x\left(C^{\alpha }u(x,t)\right)\equiv C^{\alpha +1}u(x,t)$, only holds if ${\partial }_{x}u(0,t)=0$ (see Remark 6.5 in [4] and Subsection 3.1. below). This explains, why, in some cases use of the order ($1<1+\alpha \le 2)$ C${}^{\alpha +1}$ fractional derivative in Eq. (1) is not conservative. Indeed, we note that a core advantage of the PS${}^{\alpha +1}$ derivative, adopted by [4], is that, by its construction, ${\partial }_x\left({\text{C}}^{\alpha }u(x,t)\right)\equiv {\text{PS}}^{\alpha +1}u(x,t)$ and thus, when used in Eq. (1) conservation is enforced.

Our task here is to study the space fractional diffusion equation in $\Omega$, under a variety of flux definitions. In a departure from [4], we will base our analysis and solutions on the conserved form of the diffusion transport equation, Eq. (2). Our investigation generates a number of insights and results related to the solution of fractional diffusion equations in bounded domains and also highlights possible advantages in using the Caputo derivative in seeking discrete (numerical) solutions.

We start with a general conserved statement for a transient diffusion transport written in terms of a general flux $q$. On specific definition of the flux we arrive at three alternative diffusion transport models, the classic Fourier second law ($q\sim {\partial }_{x}u$), a Caputo model (${\partial }_{x_C}^{\alpha }u$), and a Riemann–Liouville model (${\partial }_{x_{RL}}^{\alpha }u$).

To lay a foundation for development of numerical solutions for the general conserved model we present some analytical results related to the well posedness of fractional flux models.

We gather here in Section 3 the currently available results, see [5], [6], [7], [8], on the well-posedness of (1) or (2) augmented with an array of boundary conditions and initial data. It turns out that we can consider the data from the following list: Dirichlet conditions, fixed flux, fixed slope conditions and mixed ones.

Another important issue is the validity of the maximum principle. It holds for some problems written in a conservative form. We present not only a rigorous argument, but also discuss the consequence of this principle on the invariance of Eq. (2) under a re-scaling; i.e., regardless of the measurement scale used the actual value of the potential $u$ (e.g., temperature) should always be the same. We show, by example, that Eq. (2) with the flux given in terms of the Riemann–Liouville derivative fails to satisfy this invariance. On the other hand this scale invariance and the maximum principle hold if we use the Caputo derivative in Eq. (2). We offer the explanation in Section 7 and further comments in Section 7.5.

A thorough discussion of the boundary value problem (2) is not complete without addressing the issue of the boundary conditions. This is done in Section 3, where we ask the question whether or not we can evaluate $u$, the flux or the derivative of $u$ at the boundary of $\Omega$ for $u$ a solution to (2). In particular, we note that the condition ${\partial }_{x_C}^{\alpha }u\left(0\right)=0$ is automatically satisfied for sufficiently regular $u$, see Proposition 3.1. We explain it in detail in Section 3.1. This section is a good place to discuss the correctness of the definition of ${\partial }_{x_C}^{\alpha }$ and ${\partial }_{x_{RL}}^{\alpha }$.

Following this theory, we develop an explicit control volume discretization for the general conserved model, identifying a range of possible boundary condition treatments, the operation of which are essentially invariant to the precise definition of the flux used. A numerical analysis investigating the consistency and stability of these schemes is conducted. Then, taking a similar approach to the one adopted in [4] we study the fate of an initial non-negative pulse in $\Omega$ with both the Caputo and Riemann–Liouville models. As expected these tests give identical results to those provided in [4] and illustrate the predictive differences between using Caputo and Riemann–Liouville flux models. However, due to the fact that our analysis is based on a conserved treatment we are able to write down explicit discrete expressions for the Riemann–Liouville and Caputo fluxes. These expressions show that difference between the Caputo and Riemann–Liouville flux models reduce to the addition of a correction term to the Riemann–Liouville flux with the form of an apparent advection transport. This observation provides us with a direct and clear explanation of the predictive differences between Caputo and Riemann–Liouville fluxes. In addition it also allows us to identify situations where the Riemann–Liouville flux definition may lead to predictions which are physically and mathematically unsound. To conclude we suggest an alternative fractional flux law that essentially negates the differences between the Caputo and Riemann–Liouville flux models.