Journal ArticleDOI

# A finite difference method for an initial–boundary value problem with a Riemann–Liouville–Caputo spatial fractional derivative

01 Jan 2021-Journal of Computational and Applied Mathematics (North-Holland)-Vol. 381, pp 113020

TL;DR: A fractional Friedrichs’ inequality is derived and is used to prove that the problem approaches a steady-state solution when the source term is zero, and it is proved that the scheme converges with first order in the maximum norm.

AbstractAn initial–boundary value problem with a Riemann–Liouville-Caputo space fractional derivative of order α ∈ ( 1 , 2 ) is considered, where the boundary conditions are reflecting. A fractional Friedrichs’ inequality is derived and is used to prove that the problem approaches a steady-state solution when the source term is zero. The solution of the general problem is approximated using a finite difference scheme defined on a uniform mesh and the error analysis is given in detail for typical solutions which have a weak singularity near the spatial boundary  x = 0 . It is proved that the scheme converges with first order in the maximum norm. Numerical results are given that corroborate our theoretical results for the order of convergence of the difference scheme, the approach of the solution to steady state, and mass conservation.

Topics: Fractional calculus (60%), Boundary value problem (58%), Finite difference method (54%), Rate of convergence (54%), Singularity (53%)

### 1. Introduction

• The problem considered in this paper is inspired by [1], which gives a lengthy discussion of various types of fractional initial-boundary value problem and the boundary conditions that are appropriate for each type.
• Here the authors assume that the function v is such that the definitions make sense.
• In [1] the quantity Dα−1C,x u is called the Caputo fractional flux.
• Three numerical examples are given in Section 5, to illustrate their theoretical results.
• Note that C can take different values in different places.

### 2. Some properties of the solution

• In this section the authors first derive a Friedrichs’ inequality for Caputo derivatives (Lemma 1).
• Lemma 1 (Friedrichs’ inequality for Caputo derivatives).
• A related result was obtained in [1, but using a very technical argument.
• Corollary 1 (Stability and uniqueness of solution).

### 3. Finite difference scheme

• ,N. To discretise the spatial derivative in (1a) the authors follow [8], where the two-point boundary value problem corresponding to (1a) was considered.
• The time derivative in (1a) is discretised by the backward Euler method.
• Thus, the diagonal entries of A are positive and its off-diagonal entries are nonpositive.
• This completes the inductive step and the proof.

### 4. Error estimate for scheme

• In this section an error bound for their difference scheme is derived.
• To convert these bounds to an error estimate for the computed solution, the authors shall employ a barrier function whose construction is discussed in Section 4.3.
• For the time-dependent problem (1), one expects that the classical time derivative will not affect the behaviour of the spatial derivatives.
• The mesh function {Ψ̃nm}M,Nm=0,n=0 is called a discrete barrier function.

### 5. Numerical experiments

• The authors present numerical results for three examples.
• In the first example, the exact solution is known and satisfies the bounds (17); its numerical results illustrate the error estimates of Theorem 1.
• In the third example, convergence to steady state and mass conservation are discussed.
• The computed orders of converge again agree with Theorem 1.

Did you find this useful? Give us your feedback

Content maybe subject to copyright    Report

A ﬁnite dierence method for an initial-boundary value
problem with a Riemann-Liouville-Caputo spatial fractional
derivative
Jos
´
e Luis Gracia
1
, Martin Stynes
b,
a
IUMA and Department of Applied Mathematics, University of Zaragoza, Spain.
b
Applied and Computational Mathematics Division, Beijing Computational Science Research
Center, Haidian District, Beijing 100193, P.R.China.
Abstract
An initial-boundary value problem with a Riemann-Liouville-Caputo space frac-
tional derivative of order α (1, 2) is considered, where the boundary conditions
are reﬂecting. A fractional Friedrichs’ inequality is derived and is used to prove that
the problem approaches a steady-state solution when the source term is zero. The
solution of the general problem is approximated using a ﬁnite dierence scheme
deﬁned on a uniform mesh and the error analysis is given in detail for typical solu-
tions which have a weak singularity near the spatial boundary x = 0. It is proved
that the scheme converges with ﬁrst order in the maximum norm. Numerical re-
sults are given that corroborate our theoretical results for the order of convergence
of the dierence scheme, the approach of the solution to steady state, and mass
conservation.
Keywords: Fractional dierential equation, time-dependent problem,
Riemann-Liouville-Caputo fractional derivative, weak singularity, discrete
comparison principle, ﬁnite dierence scheme, steady-state problem
1. Introduction
The problem considered in this paper is inspired by [1], which gives a lengthy
discussion of various types of fractional initial-boundary value problem and the
boundary conditions that are appropriate for each type. In particular, we shall focus
on the “Caputo fractional ﬂux” and reﬂecting boundary conditions of [1, Section
6].
Corresponding author: m.stynes@csrc.ac.cn
Preprint submitted to Elsevier October 21, 2020

Set := (0, L) and Q := × (0, T ]. For (x, t) and constant r > 0, deﬁne
the Riemann-Liouville integral operator I
r
x
of order r by
I
r
x
v(x, t) :=
1
Γ(r)
Z
x
s=0
(x s)
r1
v(s, t) ds.
Then for any positive constant β with n 1 < β < n where n is a positive integer,
the Caputo fractional derivative of order β is deﬁned by
D
β
C,x
v(x, t) := I
nβ
x
n
v
x
n
(x, t).
Here we assume that the function v is such that the deﬁnitions make sense.
Let α be constant with 1 < α < 2. In this paper, we examine the initial-
boundary value problem
u
t
D
α
RLC,x
u = f for (x, t) Q, (1a)
u(x, 0) = φ(x) for x , (1b)
D
α1
C,x
u(0, t) = D
α1
C,x
u(L, t) = 0 for t (0, T], (1c)
where D
α
RLC,x
is the Riemann-Liouville-Caputo fractional derivative of order α,
which is deﬁned by
D
α
RLC,x
u(x, t) :=
x
D
α1
C,x
u(x, t) for x > 0 and 0 < t T.
This hybrid fractional derivative D
α
RLC,x
has been suggested by several researchers,
from both modelling and mathematical viewpoints; see [8] for references. In [1]
the quantity D
α1
C,x
u is called the Caputo fractional ﬂux.
The left boundary condition in (1c) is deﬁned by
0 = D
α1
C,x
u(0, t) := lim
x0
+
D
α1
C,x
u(x, t).
It is suitable for certain physical models [3] and removes a troublesome singularity
from the solution u(x, t) at x = 0; see [1, 5, 8] and Remark 1 below. Furthermore,
both boundary conditions in (1c) are reﬂecting and ensure that mass is conserved;
see the discussion in [1].
Remark 1. In [8] it is shown that D
α1
C,x
u(0, t) = 0 in (1c) is equivalent to the
classical Neumann boundary condition u
x
(0, t) = 0 if D
α
RLC,x
u(·, t) C[0, L]. Our
analysis in Sections 3 and 4 assumes that the solution u of problem (1) has this
regularity, so for convenience our numerical method will discretise u
x
(0, t) = 0
2

α1
C,x
u(0, t) = 0; but the other boundary condition D
α1
C,x
u(L, t) = 0
of (1c) cannot be simpliﬁed in the same way and must be handled directly.
Remark 2. For the function x
β
with α 1 < β < 1, one has [4, p. 193]
D
α1
C,x
x
β
=
Γ(β + 1)
Γ(β α + 2)
x
βα+1
and
d
dx
x
β
= βx
β1
.
Here β α + 1 > 0 while β 1 < 0, i.e., for these functions x
β
, the boundary
condition D
α1
C,x
u(0, t) = 0 of (1c) does not imply the Neumann boundary condition
u
x
= 0, unlike the situation described in Remark 1. Thus replacing D
α1
C,x
u(0, t) = 0
by u
x
(x, 0) = 0 means we exclude functions u(x, t) that behave like a multiple
of x
β
as x 0 for some ﬁxed value of t. But for such functions, in (1a) one has
D
α
RLC,x
u(x, t) O(x
βα
) near x = 0 with β α < 0, so we are merely excluding
certain singularities in the terms of the dierential equation.
In [2, Proposition 19] it is proved that problem (1) with f 0 is well-posed
in the Banach space L
1
(0, L) for each value of t. We shall assume a reasonable
amount of smoothness of the solution see equation (17) below. Our aim is to
approximate the solution u of problem (1) by a ﬁnite dierence method whose
analysis requires these bounds on derivatives.
The structure of the paper is as follows. In Section 2 we discuss some proper-
ties of the solution u of problem (1). A fractional Friedrichs’ inequality for Caputo
fractional derivatives is established and is used to prove that u converges to the
steady state solution when f 0. A ﬁnite dierence scheme for solving (1) on
a uniform mesh is deﬁned in Section 3 and it is shown that it satisﬁes a discrete
comparison principle. In Section 4 this principle and an appropriate barrier func-
tion are used to prove that the solution of the ﬁnite dierence scheme converges to
u with ﬁrst order in the discrete maximum norm. Three numerical examples are
given in Section 5, to illustrate our theoretical results.
Notation: Denote by AC[0, L] the set of absolutely continuous functions on [0, L]
and by L
p
(0, L) the usual Lebesgue space with norm k · k
L
P
(0,L)
. Throughout the
paper, C denotes a generic constant that can depend on the data of the problem (1)
but it is independent of the mesh used for its numerical solution. Note that C can
take dierent values in dierent places.
2. Some properties of the solution
In this section we ﬁrst derive a Friedrichs’ inequality for Caputo derivatives
(Lemma 1). Hence, in Lemma 2, we prove convergence of the solution u to the
3

R
L
x=0
φ(x) dx
/L in the special case when f 0.
This lemma implies uniqueness of the solution to (1) and stability of this solution
in terms of perturbations of the initial condition; see Corollary 1.
Lemma 1 (Friedrichs’ inequality for Caputo derivatives). Let β (0, 1). Suppose
that v AC[0, L] with
R
v dx = 0 and kD
β
C,x
vk
L
2
(0,L)
< . Then
kvk
L
2
(0,L)
L
β
Γ(β + 1)
kD
β
C,x
vk
L
2
(0,L)
. (2)
Proof. For all x (0, L], by [4, Theorem 3.8] one has
z(x) := v(x) v(0) = (I
β
x
D
β
C,x
v)(x) = (ω
β
D
β
C,x
v)(x),
where denotes convolution and ω
β
(x) := x
β1
/Γ(β). Hence, similarly to [6,
Lemma 2.6], we get
kzk
L
2
(0,L)
= kω
β
D
β
C,x
vk
L
2
(0,L)
kω
β
k
L
1
(0,L)
kD
β
C,x
vk
L
2
(0,L)
=
L
β
Γ(β + 1)
kD
β
C,x
vk
L
2
(0,L)
, (3)
using Young’s inequality for convolutions. But
kzk
2
L
2
(0,L)
=
Z
L
x=0
v
2
(x) 2v(0)v(x) + v
2
(0)
dx = kvk
2
L
2
(0,L)
+ Lv
2
(0),
since
R
vdx = 0. Thus, (3) implies (2).
Remark 3. If v AC[0, L] and 0 < β < 1, then kD
β
C,x
vk
L
p
(0,L)
< for 1 p < 1
by [4, Lemma 2.12]. Thus in Lemma 1 the hypothesis v AC[0, L] implies the
hypothesis kD
β
C,x
vk
L
2
(0,L)
< when β < 1/2.
The next lemma shows that when f 0, the solution of problem (1) converges
(in the L
2
(0, L) sense) to the steady-state solution as t . A related result
was obtained in [1, Appendix] (see Remark 5 below), but using a very technical
argument. Our simpler proof is based on the Friedrichs’ inequality of Lemma 1.
Lemma 2. Assume that f 0 in (1). Assume also that the solution of problem (1)
satisﬁes u(x, ·) AC[0, T ] for each x and u
x
(·, t) AC[0, L] for each t. Then
u(x, t)
1
L
Z
L
x=0
φ(x) dx
L
2
(0,L)
0 as t .
4

Proof. Set
v(x, t) = u(x, t)
1
L
Z
L
x=0
φ(x) dx .
Then
v
t
D
α
RLC,x
v = 0 for (x, t) Q, (4a)
v(x, 0) = φ(x)
1
L
Z
L
x=0
φ(x) dx for x (0, L), (4b)
D
α1
C,x
v(0, t) = D
α1
C,x
v(L, t) = 0 for t (0, T]. (4c)
Observe that v has the additional property that
Z
L
x=0
v(x, 0) dx = 0.
Note ﬁrst that mass is conserved, i.e., for each t we have
Z
L
x=0
v(x, t) =
Z
L
x=0
v(x, 0) = 0; (5)
to see this, integrate v
t
D
α
RLC,x
v = 0 over [0, L] × [0, t] and use (4c) to eliminate
the x-derivative terms.
Let t [0, T ] be arbitrary. Let k be a constant that is chosen later. Multi-
ply (4a) by e
kt
v(x, t) then integrate over [0, L] × [0, t]. Using e
kt
vv
t
= e
kt
(v
2
)
t
/2 and
integration by parts in time and in space, we get
0 =
1
2
Z
L
x=0
e
kt
v
2
(x, t) v
2
(x, 0)
dx
1
2
Z
t
s=0
ke
ks
Z
L
x=0
v
2
(x, s) ds dx
+
Z
t
s=0
e
ks
Z
L
x=0
v
x
(x, s)D
α1
C,x
v(x, s) dx ds. (6)
These integrations by parts are justiﬁed since by hypothesis v(x, ·) AC[0, T] for
each x and v
x
(·, t) AC[0, L] for each t, so D
α1
C,x
v(·, t) = I
2α
x
v
x
(·, t) AC[0, L]
by [9, Lemma 2.3]. But by [10, Lemma 3.1] one has
Z
L
x=0
v
x
(x, s)D
α1
C,x
v(x, s) dx =
Z
L
x=0
v
x
(x, s)(I
2α
x
v
x
)(x, s) dx
cos
(2 α)π
2
Z
L
x=0
I
1α/2
x
v
x
2
(x, s) dx
5

##### Citations
More filters

Posted Content
Abstract: This paper derives physically meaningful boundary conditions for fractional diffusion equations, using a mass balance approach. Numerical solutions are presented, and theoretical properties are reviewed, including well-posedness and steady state solutions. Absorbing and reflecting boundary conditions are considered, and illustrated through several examples. Reflecting boundary conditions involve fractional derivatives. The Caputo fractional derivative is shown to be unsuitable for modeling fractional diffusion, since the resulting boundary value problem is not positivity preserving.

39 citations

##### References
More filters

Journal ArticleDOI
Abstract: We discuss existence, uniqueness, and structural stability of solutions of nonlinear differential equations of fractional order. The differential operators are taken in the Riemann–Liouville sense and the initial conditions are specified according to Caputo's suggestion, thus allowing for interpretation in a physically meaningful way. We investigate in particular the dependence of the solution on the order of the differential equation and on the initial condition, and we relate our results to the selection of appropriate numerical schemes for the solution of fractional differential equations.

2,775 citations

Journal ArticleDOI
Abstract: In this article a theoretical framework for the Galerkin finite element approximation to the steady state fractional advection dispersion equation is presented. Appropriate fractional derivative spaces are defined and shown to be equivalent to the usual fractional dimension Sobolev spaces Hs. Existence and uniqueness results are proven, and error estimates for the Galerkin approximation derived. Numerical results are included that confirm the theoretical estimates. © 2005 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2005

639 citations

Journal ArticleDOI
TL;DR: The final convergence result shows clearly how the regularity of the solution and the grading of the mesh affect the order of convergence of the difference scheme, so one can choose an optimal mesh grading.
Abstract: A reaction-diffusion problem with a Caputo time derivative of order $\alpha\in (0,1)$ is considered. The solution of such a problem is shown in general to have a weak singularity near the initial time $t=0$, and sharp pointwise bounds on certain derivatives of this solution are derived. A new analysis of a standard finite difference method for the problem is given, taking into account this initial singularity. This analysis encompasses both uniform meshes and meshes that are graded in time, and includes new stability and consistency bounds. The final convergence result shows clearly how the regularity of the solution and the grading of the mesh affect the order of convergence of the difference scheme, so one can choose an optimal mesh grading. Numerical results are presented that confirm the sharpness of the error analysis.

343 citations

Journal ArticleDOI
Abstract: A class of nonlocal models based on the use of fractional derivatives (FDs) is proposed to describe nondiffusive transport in magnetically confined plasmas. FDs are integro-differential operators that incorporate in a unified framework asymmetric non-Fickian transport, non-Markovian (“memory”) effects, and nondiffusive scaling. To overcome the limitations of fractional models in unbounded domains, we use regularized FDs that allow the incorporation of finite-size domain effects, boundary conditions, and variable diffusivities. We present an α-weighted explicit/implicit numerical integration scheme based on the Grunwald-Letnikov representation of the regularized fractional diffusion operator in flux conserving form. In sharp contrast with the standard diffusive model, the strong nonlocality of fractional diffusion leads to a linear in time response for a decaying pulse at short times. In addition, an anomalous fractional pinch is observed, accompanied by the development of an uphill transport region where the “effective” diffusivity becomes negative. The fractional flux is in general asymmetric and, for steady states, it has a negative (toward the core) component that enhances confinement and a positive component that increases toward the edge and leads to poor confinement. The model exhibits the characteristic anomalous scaling of the confinement time, τ, with the system’s size, L, τ∼Lα, of low-confinement mode plasma where 1<α<2 is the order of the FD operator. Numerical solutions of the model with an off-axis source show that the fractional inward transport gives rise to profile peaking reminiscent of what is observed in tokamak discharges with auxiliary off-axis heating. Also, cold-pulse perturbations to steady sates in the model exhibit fast, nondiffusive propagation phenomena that resemble perturbative experiments.

100 citations

Posted Content
Abstract: In this article we investigate the solution of the steady-state fractional diffusion equation on a bounded domain in $\real^{1}$. From an analysis of the underlying model problem, we postulate that the fractional diffusion operator in the modeling equations is neither the Riemann-Liouville nor the Caputo fractional differential operators. We then find a closed form expression for the kernel of the fractional diffusion operator which, in most cases, determines the regularity of the solution. Next we establish that the Jacobi polynomials are pseudo eigenfunctions for the fractional diffusion operator. A spectral type approximation method for the solution of the steady-state fractional diffusion equation is then proposed and studied.

88 citations