Next Article in Journal
Existence Results for Nonlinear Fractional Differential Inclusions via q-ROF Fixed Point
Next Article in Special Issue
Fourth-Order Numerical Solutions for a Fuzzy Time-Fractional Convection–Diffusion Equation under Caputo Generalized Hukuhara Derivative
Previous Article in Journal
New Results for Homoclinic Fractional Hamiltonian Systems of Order α∈(1/2,1]
Previous Article in Special Issue
Sequential Caputo–Hadamard Fractional Differential Equations with Boundary Conditions in Banach Spaces
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Second-Order Crank-Nicolson-Type Scheme for Nonlinear Space–Time Reaction–Diffusion Equations on Time-Graded Meshes

1
Department of Mathematics, University of Louisiana at Lafayette, Lafayette, LA 70504, USA
2
Department of Mathematical Sciences, Middle Tennessee State University, Murfreesboro, TN 37132, USA
3
Department of Mathematics, Clarkson University, Potsdam, NY 13699, USA
*
Author to whom correspondence should be addressed.
Fractal Fract. 2023, 7(1), 40; https://doi.org/10.3390/fractalfract7010040
Submission received: 11 November 2022 / Revised: 23 December 2022 / Accepted: 26 December 2022 / Published: 30 December 2022

Abstract

:
A weak singularity in the solution of time-fractional differential equations can degrade the accuracy of numerical methods when employing a uniform mesh, especially with schemes involving the Caputo derivative (order α , ), where time accuracy is of the order ( 2 α ) or ( 1 + α ) . To deal with this problem, we present a second-order numerical scheme for nonlinear time–space fractional reaction–diffusion equations. For spatial resolution, we employ a matrix transfer technique. Using graded meshes in time, we improve the convergence rate of the algorithm. Furthermore, some sharp error estimates that give an optimal second-order rate of convergence are presented and proven. We discuss the stability properties of the numerical scheme and elaborate on several empirical examples that corroborate our theoretical observations.

1. Introduction

The last decade has witnessed tremendous developments in practical methods to solve fractional differential equations. These problems are of particular importance because they can provide a better model for understanding complex phenomena such as memory-dependent processes [1,2,3], material properties [4], diffusion in media with memory [5,6], groundwater modeling [7,8], and control theory [9]. Recently, many researchers have adopted fractional-order models to predict and gain insight into the evolution of the COVID-19 pandemic. This is possible due to the memory/hereditary properties inherent in the fractional-order derivatives, cf. [10,11,12,13,14].
We study a nonlinear time–space fractional reaction–diffusion problem in the form
c D 0 , t α u = κ Δ β 2 u ( x , t ) + g ( u ) , in Ω × ( 0 , T ) , u ( x , 0 ) = φ ( x ) , x Ω R , u ( x , t ) | Ω = 0
where Ω is bounded in R , Ω denotes the boundary of Ω , κ is the diffusion coefficient, ( Δ ) β 2 denotes the Laplacian of a fractional order β , 1 < β 2 , and g ( u ) is a sufficiently smooth function. The α -order Caputo derivative, 0 < α 1 , c D 0 , t α u , in variable t, is adopted here and defined as
c D 0 , t α u ( x , t ) = 1 Γ ( 1 α ) 0 t ( t s ) α u ( x , s ) s d s .
The operator ( Δ ) β 2 is taken here as the spectral Laplacian of a fractional order β , 1 < β 2 given by
( Δ ) β 2 u ( x ) = k = 1 c k λ k β 2 ϕ k ( x ) .
where λ k and ϕ k ( x ) are the eigenvalues and eigenfunctions of the Laplace operator Δ on Ω , respectively, and c k are Fourier coefficients of u (in { ϕ k ( x ) } ) (see [15]).
Several numerical methods for solving Problem (1) or some of its counterparts have been developed and investigated by authors based on the uniform discretization of the Caputo derivative (see, for example, [16,17,18,19,20,21,22,23]). In most cases, the derived schemes were either ( 2 α ) or ( 1 + α ) accurate in time. This somewhat reduced order of convergence is to be expected due to a singular kernel ( t s ) α embedded in the time derivative. This, in turn, has motivated the research question of how to improve the order beyond ( 2 α ) and ( 1 + α ) . The natural choice is to use nonuniform time meshes.
Brunner [24] made use of meshes that are graded in order to improve the accuracy of the approximation to a Volterra integral equation of the second kind with a weakly singular kernel employing collocation methods. Zhang et al. [25] developed a numerical method for a linear counterpart of (1) based on the nonuniform discretization of the Caputo derivative and the compact difference method for spatial discretization. Their theoretical analysis and numerical examples showed the efficiency of their methods. Lyu and Vong [26] proposed a high-order method to resolve a time-fractional Benjamin–Bona–Mahony equation over a nonuniform temporal mesh. Stynes et al. [27] investigated the stability and error analysis of a finite difference scheme using a uniform mesh and meshes graded in time. Liao et al. [28] investigated the convergence and stability of an L 1 technique to solve linear reaction–subdiffusion equations with the Caputo derivative. Kopteva [29] discussed the error analysis of the L 1 method for a fractional-order parabolic problem in two and three dimensions using both uniform and graded meshes. Wang and Zhou [30] proved the convergence of the corrected k-step backward difference formula without imposing further regularity assumptions on the solution of the semilinear subdiffusion equation. Mustapha [31] developed an L 1 scheme for subdiffusion equations with Riemann–Liouville time-fractional derivatives on nonuniform time intervals. He used the regularity of the solution and the properties of the nonuniform mesh to obtain a second-order accurate scheme.
In this present study, we propose a predictor–corrector numerical scheme for solving (1) based on time-graded meshes. The scheme is similar to the one given in [18]. We then explore the regularity properties of the solution and some of the properties of the meshes to derive a scheme that is second-order accurate in time.
The remaining sections are organized as follows. Section 2 briefly discusses the spatial discretization method, and thereafter we derive the time-stepping scheme for the solution of (3). In Section 3, we discuss the error analysis and stability of the scheme. In Section 4, we give some numerical examples to illustrate the convergence of the scheme. Finally, in Section 5 there are concluding remarks.

2. Numerical Scheme

2.1. Matrix Transfer Technique for Spatial Discretizations

Let M be given and we denote by x j , for 0 j M , a 1D uniform grid point of size h. It was shown in [32] that
Δ β 2 u ( x ) A β 2 u ( x )
where
A = 1 h 2 2 1 0 0 0 1 2 1 0 0 0 0 1 2 1 0 0 0 1 2 , u ( x ) = u ( x 1 ) u ( x 2 ) u ( x M 2 ) u ( x M 1 ) ,
then, the error in (2) is of order 2. This approximation can be thought of as transferring the matrix approximation of Δ to approximate Δ β 2 . Applying this technique to (1), we obtain a system of nonlinear time-fractional differential equations in the form
c D 0 , t α u + A β 2 u = g ( u ) , u ( 0 ) = u 0 ,
where u and g ( u ) are the vectors that denote the nodal values of u and g , respectively, and we have chosen κ = 1 , without any loss of generality. In addition, we write u and g ( u ) instead of u and g ( u ) to denote the vectors of the node values. Consequently, the major concern of this paper is the accurate time discretization of Problem (3).

2.2. Time Discretizations

In this section, the second-order time-stepping scheme over time-graded meshes is considered for solving the semi-discrete problem (3). Furthermore, the stability results are developed. We use the time-graded mesh having subintervals I n = [ t n , t n + 1 ] , n = 0 , , N 1 , with 0 = t 0 < < t N = T . This has the following grid points:
t n = ( n τ ) γ , 0 n N , for γ 1 , with τ = T 1 / γ N .
Let τ n = t n + 1 t n denote the stepsize of the n-th subinterval I n . The following properties (see [24,33]) hold for n 1 ,
t n + 1 2 γ t n ,
γ τ t n 1 1 / γ τ n γ τ t n + 1 1 1 / γ ,
τ n τ n 1 C γ τ 2 min { 1 , t n + 1 1 2 / γ }
τ n τ max γ T / N .
where
τ max = max 1 j n 1 τ j
We note that Equation (3) can be reformulated in the form of the Volterra integral equation
u ( t ) u 0 = 1 Γ ( α ) 0 t ( t s ) α 1 A β 2 u ( s ) + g ( u ( s ) ) d s = 0 I t α g ( u ( t ) ) A β 2 0 I t α u ( t ) ,
where
a I t α w ( t ) = 1 Γ ( α ) a t ( t s ) α 1 w ( s ) d s , a I t α g ( w ( t ) ) = 1 Γ ( α ) a t ( t s ) α 1 g ( w ( s ) ) d s .
Now, let u n : = u ( t n ) and g ( u n ) : = g ( u ( t n ) ) . Estimating t = t n and t n + 1 , we obtain the difference in successive terms as
u ( t n + 1 ) u ( t n ) = 0 I t n + 1 α g ( u n + 1 ) 0 I t n α g ( u n ) A β 2 0 I t n + 1 α u ( t n + 1 ) 0 I t n α u ( t n ) = t n I t n + 1 α g ( u n + 1 ) A β 2 t n I t n + 1 α u ( t n + 1 ) + Q n , u e + Q n , g e ,
where
Q n , u e = A β 2 0 I t n α u ( t n + 1 ) u ( t n )
Q n , g e = 0 I t n α g ( u n + 1 ) g ( u n )
That is,
u ( t n + 1 ) u ( t n ) = 1 Γ ( α ) A β 2 t n t n + 1 ( t n + 1 s ) α 1 u ( s ) d s + 1 Γ ( α ) t n t n + 1 ( t n + 1 s ) α 1 g ( u ( s ) ) d s + H n e ,
where
H n e = Q n , u e + Q n , g e .
If we replace u ( s ) and g ( u ) with linear interpolants over the interval I n , that is,
u ( s ) u n + ( s t n ) u n + 1 u n τ n , s [ t n , t n + 1 ] ,
and
g ( u ( s ) ) g ( u n ) + ( s t n ) g ( u n + 1 ) g ( u n ) τ n , s [ t n , t n + 1 ] .
We obtain
u n + 1 u n = α τ n α Γ ( α + 2 ) A β 2 u n τ n α Γ ( α + 2 ) A β 2 u n + 1 + α τ n α Γ ( α + 2 ) g n + τ n α Γ ( α + 2 ) g n + 1 + H n e .
At a glance, we see that (12) is an implicit scheme. In order to reduce the computation burden, we thus go back to (11) and approximate the nonlinear function, g ( u ) , on the interval [ t n , t n + 1 ] by a constant polynomial to obtain the following predictor–corrector scheme after some simplifications:
Γ ( α + 2 ) I + τ n α A β 2 u n + 1 p = Γ ( α + 2 ) I α τ n α A β 2 u n + τ n α ( α + 1 ) g ( u n ) + Γ ( α + 2 ) H n e Γ ( α + 2 ) I + τ n α A β 2 u n + 1 = Γ ( α + 2 ) I α τ n α A β 2 u n + τ n α α g ( u n ) + g ( u n + 1 p ) + Γ ( α + 2 ) H n e .
where I is the identity matrix.
Using linear approximations for both u ( t ) and g ( u ( t ) ) in Equations (9) and (10), the history term H n e is approximated as
H n e H n a = j = 0 n a j , n A β 2 u j + g ( u j ) ,
where
a j , n = 1 Γ ( α + 2 ) τ 0 1 τ n + 1 , 0 α ( τ n + 1 , 1 α τ 0 ) + τ n + 1 , 1 1 + α + τ n , 1 τ n , 0 α α τ 0 τ n , 0 α τ n , 1 1 + α , j = 0 , τ j 1 1 τ n + 1 , j 1 1 + α τ n , j 1 1 + α + τ j 1 τ n + 1 , j + 1 1 + α τ n , j + 1 1 + α τ n + 1 , j α τ j 1 1 τ n + 1 , j 1 + τ j 1 τ n + 1 , j + 1 1 j n 1 , + τ n , j α τ j 1 1 τ n , j 1 + τ j 1 τ n , j + 1 , τ n 1 1 τ n + 1 , n 1 1 + α τ n + 1 , n 1 τ n α τ n 1 1 + α α τ n 1 τ n α , j = n .
             τ n , j = t n t j .

3. Error and Stability Analysis

Here, we carry out the error analysis and discuss the stability property of the proposed scheme (13). For the error analysis, it is assumed that u, the solution to (1), satisfies
| | u ( t ) | | 1 M and | | u ( ) ( t ) | | 1 M t σ + α / 2 , = 1 , 2 , 3 ,
with the regularity parameter σ ( 0 , 1 ) to be determined from the error analysis. This is in line with the published results [27,30,31,33,34]. The constant M is positive, denotes the partial derivative of an appropriate order of u with respect to t, and | | · | | is the Sobolev norm on H ( Ω ) . Of course, this reduces to the L 2 norm whenever = 0 . The stability analysis of this article considers the initial data perturbations, i.e., the sensitivity of the numerical solutions to the small changes in the initial data. Function g ( u ) and its derivatives with respect to u, g ( r ) ( u ) for r = 1 , 2 , are assumed to be Lipschitz in the time domain Ω × [ 0 , T ] .

3.1. Error Analysis

Lemma 1. 
Given any positive sequence { a j } and for γ 1 , we have
j = 2 n 1 a j τ j 3 τ j 1 3 L j 1 , n L j , n + 1 C τ t n α 1 / γ max j = 2 n 1 ( a j τ j 2 ) ,
where
L j , n = 1 2 Γ ( α ) t j t j + 1 ( s t j ) ( t j + 1 s ) ( t n s ) α 1 d s .
Proof. 
Cf. [31]. □
Lemma 2. 
For 1 n N , γ > 2 σ + α / 2 + 1 and τ is sufficiently small,
1 Γ ( α ) j = 0 n 1 t j t j + 1 ( t n + 1 s ) α 1 ( t n s ) α 1 E 1 , u ( s ) d s 1 C τ 2 t n σ + 3 α / 2 2 / γ ,
where
E 1 , u ( s ) = 1 2 t j s ( w s ) 2 u ( w ) d w ( s t j ) 2 τ j t j t j + 1 ( w t j ) 2 u ( w ) d w .
Proof. 
Let
I 1 , u = 1 Γ ( α ) j = 0 n 1 t j t j + 1 ( t n + 1 s ) α 1 ( t n s ) α 1 E 1 , u ( s ) d s
| | I 1 , u | | 1 = 1 Γ ( α ) j = 0 n 1 t j t j + 1 ( t n + 1 s ) α 1 ( t n s ) α 1 E 1 , u ( s ) d s 1 2 Γ ( α ) j = 0 n 1 t j t j + 1 ( t n s ) α 1 | | E 1 , u ( s ) | | 1 d s .
For n = 1 , s ( t 0 , t 1 )
| | E 1 , u ( s ) | | 1 1 2 0 s ( w s ) 2 | | u ( w ) | | 1 d w + s t 1 0 t 1 w 2 | | u ( w ) | | 1 d w M 0 s w 2 w σ + α 2 3 + M 2 t 1 s 0 t 1 w 2 w σ + α 2 3 d w C max { s σ + α 2 , s t 1 σ + α 2 1 } .
Noting that
t n s = n γ t 1 s n γ n γ ( t 1 s ) ,
we obtain
0 t 1 ( t n s ) α 1 | | E 1 , u ( s ) | | 1 d s C n γ ( α 1 ) 0 t 1 ( t 1 s ) α 1 max { s σ + α 2 , s t 1 σ + α 2 1 } d s C t n α 1 T 1 α N γ ( 1 α ) 0 t 1 ( t 1 s ) α 1 max { s σ + α 2 , s t 1 σ + α 2 1 } d s C t n α 1 τ γ ( 1 α ) 0 t 1 ( t 1 s ) α 1 max { s σ + α 2 , s t 1 σ + α 2 1 } d s C t n α 1 τ γ ( 1 α ) t 1 σ + 3 α 2 = C t n α 1 τ 2 t 1 σ + α 2 2 γ + 1 C τ 2 t n σ + 3 α 2 2 γ , for γ 2 ( σ + α 2 + 1 ) .
Let s ( t j , t j + 1 ) , j 1 and n 2 ,
| | I 1 , u | | 1 2 Γ ( α ) 0 t 1 ( t n s ) α 1 | | E 1 , u ( s ) | | 1 d s + j = 1 n 1 t j t j + 1 ( t n s ) α 1 | | E 1 , u ( s ) | | 1 d s .
For s ( t j , t j + 1 ) and j 1
| | E 1 , u ( s ) | | 1 = 1 2 t j s ( w s ) 2 | | u ( w ) | | 1 d w + ( s t j ) 2 τ j t j t j + 1 ( w t j ) 2 | | u ( w ) | | 1 d w C τ j 2 t j t j + 1 | | u ( w ) | | 1 d w C τ 3 t j + 1 3 3 / γ s σ + α / 2 3 .
Therefore,
2 Γ ( α ) j = 1 n 1 t j t j + 1 ( t n s ) α 1 | | E 1 , u ( s ) | | 1 d s C τ 3 j = 1 n 1 t j + 1 3 3 / γ t j t j + 1 ( t n s ) α 1 s σ + α 2 3 d s C τ 3 t n 3 3 / γ t 1 t n ( t n s ) α 1 s σ + α 2 3 d s C τ 3 t n σ + 3 α 2 3 γ .
Hence, the following bound is obtained
| | I 1 , u | | 1 C τ 2 t n σ + 3 α 2 2 γ , γ > 2 σ + α 2 + 1
Lemma 3. 
Assume 1 n N and let τ be sufficiently small,
j = 0 n 1 u ( t j ) L j , n L j , n + 1 1 C τ 2 t n σ + 3 α 2 2 γ ,
holds if γ > max 2 σ + α 2 + 1 , 2 σ + 3 α 2 3 .
Proof. 
We begin by letting
I 2 , u = j = 0 n 1 u ( t j ) L j , n L j , n + 1 = j = 1 n 1 u ( t j + 1 ) u ( t j ) L j , n + u ( t 0 ) L 0 , n L 0 , n + 1 j = 1 n 1 u ( t j + 1 ) τ j + 1 3 τ j 3 1 L j , n + j = 2 n 1 u ( t j ) τ j 3 τ j 1 3 L j 1 , n L j , n + 1 u ( t 1 ) L 1 , n + 1 + u ( t n ) τ n 3 τ n 1 3 L n 1 , n = η 1 η 2 + η 3 η 4 ,
where
η 1 = j = 1 n 1 [ u ( t j + 1 ) u ( t j ) ] L j , n , η 2 = j = 1 n 1 u ( t j + 1 ) τ j + 1 3 τ j 3 1 L j , n , η 3 = j = 2 n 1 u ( t j ) τ j 3 τ j 1 3 L j 1 , n L j , n + 1 , η 4 = u ( t 0 ) [ L 0 , n + 1 L 0 , n ] + u ( t 1 ) L 1 , n + 1 u ( t n ) τ n 3 τ n 1 3 L n 1 , n .
For η 1 , we have
| | η 1 | | 1 C j = 1 n 1 τ j 2 t j t j + 1 | | u ( w ) | | 1 d w t j t j + 1 ( t n s ) α 1 d s C j = 1 n 1 τ j 2 t j + 1 σ + α / 2 2 t j t j + 1 ( t n s ) α 1 d s C τ 2 t 1 t n ( t n s ) α 1 s σ + α / 2 2 / γ d s C τ 2 t n σ + 3 α 2 2 γ , holds if γ > 2 σ + α 2 + 1 .
Noting that
τ j + 1 3 τ j 3 1 = τ j + 1 3 τ j 3 τ j 3 C τ j 1 ( τ j + 1 τ j ) C τ 2 τ j 1 t j + 2 1 2 / γ .
Therefore,
| | η 2 | | 1 C τ 2 j = 1 n 1 τ j 1 t j + 2 1 2 / γ | | u ( t j + 1 ) | | 1 τ j 2 t j t j + 1 ( t n s ) α 1 d s C τ 2 j = 1 n 1 τ j t j + 2 1 2 / γ | | u ( t j + 1 ) | | 1 t j t j + 1 ( t n s ) α 1 d s C τ 3 j = 1 n 1 t j + 1 1 1 / γ t j + 2 1 2 / γ t j + 1 σ + α / 2 2 t j t j + 1 ( t n s ) α 1 d s C τ 3 max j = 1 n 1 t j + 2 1 2 / γ t 1 t n s σ + α / 2 1 / γ 1 ( t n s ) α 1 d s C τ 3 max j = 1 n 1 t j + 2 1 2 / γ t n σ + 3 α / 2 1 / γ 1 .
Next, we bound η 3 using Lemma 1.
| | η 3 | | 1 C j = 2 n 1 t j σ + α / 2 2 τ j 3 τ j 1 3 L j 1 , n L j , n + 1 C τ t n α 1 / γ max j = 2 n 1 t j σ + α / 2 2 τ j 2 C τ 3 t n α + 2 3 / γ max j = 2 n 1 t j σ + α / 2 2 .
To estimate | | η 4 | | 1 , we begin with the bound
| | u ( t 0 ) [ L 0 , n + 1 L 0 , n ] | | 1 = 1 2 Γ ( α ) u ( t 0 ) 0 t 1 s ( t 1 s ) [ ( t n + 1 s ) α 1 ( t n s ) α 1 ] d s 1 C τ 1 2 0 t 1 ( t n s ) α 1 | | u ( t 0 ) | | 1 d s C τ 2 0 t n s σ + α / 2 2 / γ ( t n s ) α 1 d s C τ 2 t n σ + 3 α / 2 2 / γ .
In addition,
| | u ( t 1 ) L 1 , n + 1 | | 1 C τ 1 2 t 1 σ + σ / 2 2 t 1 t 2 ( t n + 1 s ) α 1 d s C τ 2 t 1 t 2 s σ + α / 2 2 / γ ( t n s ) α 1 d s C τ 2 t n σ + 3 α / 2 2 / γ .
Hence, we obtain
τ n 3 τ n 1 3 u ( t n ) L n 1 , n 1 C τ n 3 τ n 1 | | u ( t n ) | | 1 t n 1 t n ( t n s ) α 1 d s C τ 2 t n + 1 3 t n 1 1 / γ 1 t n σ + α / 2 2 3 / γ t n 1 t n ( t n s ) α 1 d s C τ 2 t n + 1 3 t n 1 t n s σ + α / 2 2 / γ 3 ( t n s ) α 1 d s C τ 2 t n σ + 3 α / 2 2 / γ 3 t n + 1 3 C τ 2 t n + 1 σ + 3 α / 2 2 / γ , for γ > 2 σ + 3 α / 2 3 .
Thus, for γ > 2 σ + 3 α / 2 3 , we obtain
| | η 4 | | 1 C τ 2 t n + 1 σ + 3 α / 2 2 / γ .
Combining all the bounds, we finally obtain
| | I 2 , u | | 1 | | η 1 | | 1 + | | η 2 | | 1 + | | η 3 | | 1 + | | η 4 | | 1 C τ 2 t n + 1 σ + 3 α / 2 2 / γ + C τ 3 t n σ + 3 α / 2 1 / γ 1 max j = 1 n 1 t j + 1 1 2 / γ + C τ 3 t n α + 2 3 / γ max j = 2 n 1 t j σ + α / 2 2 C τ 2 t n + 1 σ + 3 α / 2 2 / γ , for γ > max 2 σ + α / 2 + 1 , 2 σ + 3 α / 2 3 ,
with σ + 3 α / 2 > 3 and τ is sufficiently small. □
Lemma 4. 
For 1 n N , γ > max 2 σ + α / 2 + 1 , 2 σ + 3 α / 2 3 and for a sufficiently small τ, the following error bound arises
| | Q n , u e Q n , u a | | 1 C τ 2 t n + 1 σ + 3 α / 2 2 / γ ,
where Q n , u e is given by (9) and
Q n , u a = 1 Γ ( α ) j = 0 n t j t j + 1 ( t n + 1 s ) α 1 ( t n s ) α 1 u ˜ ( s ) d s
with
u ˜ ( s ) = u j + ( s t j ) u j + 1 u j τ j .
Proof. 
We begin the proof by observing that
Q n , u e Q n , u a = 1 Γ ( α ) j = 0 n t j t j + 1 ( t n + 1 s ) α 1 ( t n s ) α 1 ( u u ˜ ) ( s ) d s ,
Now, u ( s ) u ˜ ( s ) = E 1 , u ( s ) + E 2 , u ( s ) , where
E 1 , u ( s ) = 1 2 t j s ( w s ) 2 u ( w ) d w ( s t j ) 2 τ j t j t j + 1 ( w t j ) 2 u ( w ) d w
and
E 2 , u ( s ) = 1 2 ( s t j ) ( s t j + 1 ) u ( t j ) .
Then,
| | Q n , u e Q n , u a | | 1 | | I 1 , u | | 1 + | | I 2 , u | | 1
To complete the proof, we use the results in Lemmas 2 and 3. □
Lemma 5. 
Let g ( r ) ( u ) be Lipschitz in u for r = 0 , 1 , 2 . Then, for 1 n N ,
γ > max 2 σ + α / 2 + 1 , 2 σ + 3 α / 2 3 with σ + 3 α / 2 > 3 and for a sufficiently small τ, we have
| | Q n , g e Q n , g a | | 1 C τ 2 t n + 1 σ + 3 α / 2 2 / γ ,
where Q n , g e is given by (10) and
Q n , g a = 1 Γ ( α ) j = 0 n t j t j + 1 [ ( t n + 1 s ) α 1 ( t n s ) α 1 ] g ˜ 2 ( u ) ( s ) d s
with
g ˜ 2 ( u ( s ) ) = g ( u j ) + ( s t j ) g ( u j + 1 ) g ( u j ) τ j .
Proof. 
Q n , g e Q n , g a = 1 Γ ( α ) j = 0 n t j t j + 1 [ ( t n + 1 s ) α 1 ( t n s ) α 1 ] g ( u ) g ˜ 2 ( u ) ( s ) d s ,
Let
g ( u ( s ) ) g ˜ 2 ( u ( s ) ) = E 1 , g ( s ) + E 2 , g ( s ) + E 3 , g ( s ) ,
where
E 1 , g ( s ) = u ( s ) u ( t j ) s t j τ j u ( t j + 1 ) u ( t j ) g ( u ( t j ) ) E 2 , g ( s ) = 1 2 u ( s ) u ( t j ) 2 s t j τ j u ( t j + 1 ) u ( t j ) 2 g ( u ( t j ) ) E 3 , g ( s ) = 1 2 u ( t j ) u ( s ) u ( s ) u ( w ) 2 g ( u ( w ) ) d u ( w ) s t j 2 τ j u ( t j ) u ( t j + 1 ) ( u ( w ) u ( t j ) ) 2 g ( u ( t j ) ) d u ( w ) .
Thus, we have
| | E 1 , g ( s ) | | 1 ( u ( s ) u ˜ ( s ) ) | | g ( u ( t j ) ) | | 1 M | | u ( s ) u ˜ ( s ) | | 1
and
| | E 2 , g | | 1 1 2 u ( s ) u ( t j ) 2 s t j τ j u ( t j + 1 ) u ( t j ) 2 1 | | g ( u ( t j ) ) | | 1 max | | u ( s ) u ( t j ) | | 1 , s t j τ j | | u ( t j + 1 ) u ( t j ) | | 1 | | u ( s ) u ˜ ( s ) | | 1 | | g ( u ( t j ) ) | | 1 max t j s | | u ( w ) | | 1 d w , s t j τ j t j t j + 1 | | u ( w ) | | 1 d w | | u ( s ) u ˜ ( s ) | | 1 | | g ( u ( t j ) ) | | 1 M τ j t j + 1 σ + α / 2 1 | | u ( s ) u ˜ ( s ) | | 1 M τ t j + 1 σ + α / 2 1 / γ | | u ( s ) u ˜ ( s ) | | 1 .
Moreover,
| | E 3 , g ( s ) | | 1 M 1 u ( t j ) u ( s ) u ( s ) u ( w ) 2 d u ( w ) 1 + M 2 s t j 2 τ j u ( t j ) u ( t j + 1 ) u ( w ) u ( t j ) 2 d u ( w ) 1 M ( s t j ) 3 s 3 σ + 3 α / 2 3 + s t j τ j τ 3 t j + 1 3 σ + 3 α / 2 3 / γ C max ( s t j ) 3 s 3 σ + 3 α / 2 3 , s t j τ j τ 3 t j + 1 3 σ + 3 α / 2 3 / γ
Therefore,
| | Q n , g e Q n , g a | | 1 | | I 1 , g | | 1 + | | I 2 , g | | 1 + | | I 3 , g | | 1
where the terms in the RHS are estimated as
| | I 1 , g | | 1 1 Γ ( α ) j = 0 n 1 t j t j + 1 ( t n + 1 s ) α 1 ( t n s ) α 1 | | E 1 , g ( s ) | | 1 d s C τ 2 t n + 1 σ + 3 α / 2 2 / γ , for γ > max 2 σ + α / 2 + 1 , 2 σ + 3 α / 2 3 ,
| | I 2 , g | | 1 1 Γ ( α ) j = 0 n 1 t j t j + 1 ( t n + 1 s ) α 1 ( t n s ) α 1 | | E 2 , g | | 1 d s C τ j = 0 n 1 t j + 1 σ + α / 2 1 / γ t j t j + 1 ( t n + 1 s ) α 1 ( t n s ) α 1 | | u ( s ) u ˜ ( s ) | | 1 d s C τ t n σ + α / 2 1 / γ j = 0 n 1 t j t j + 1 ( t n + 1 s ) α 1 ( t n s ) α 1 | | u ( s ) u ˜ ( s ) | | 1 d s , for γ > 1 σ + α / 2 C τ t n σ + α / 2 1 / γ τ 2 t n σ + 3 α / 2 2 / γ , = C τ 3 t n 2 σ + 2 α 3 / γ , for γ > max 1 σ + α / 2 , 2 σ + α / 2 + 1 , 2 σ + 3 α / 2 3
and
| | I 3 , g | | 1 2 Γ ( α ) j = 0 n 1 t j t j + 1 ( t n s ) α 1 | | E 3 , g | | 1 d s C j = 0 n 1 t j t j + 1 ( t n s ) α 1 max ( s t j ) 3 s 3 σ + 3 α / 2 3 , s t j 2 τ j τ 3 t j + 1 3 σ + 3 α / 2 3 / γ C τ 3 j = 0 n 1 t j t j + 1 ( t n s ) α 1 s 3 σ + 3 α / 2 3 / γ d s C τ 3 t n 3 σ + 3 α / 2 3 / γ .
The proof is completed using the estimates for | | I 1 , g | | 1 , | | I 2 , g | | 1 and | | I 3 , g | | 1 in Equation (16), where for a sufficiently small τ , the τ 3 terms are assumed to be negligible. □
Lemma 6. 
Assume the conditions given in Lemma 5. Then, for a sufficiently small τ, the error bound
1 Γ ( α ) t n t n + 1 ( t n + 1 s ) α 1 ( g ( u ) g ˜ 1 ( u ) ) ( s ) d s 1 C τ t n + 1 σ + α 1 / γ
holds uniformly on [ t n , t n + 1 ] .
Proof. 
Noting that
g ( u ( s ) ) g ˜ 1 ( u ( t n ) ) = ( u ( s ) u ( t n ) ) g ( u ( t n ) ) + u ( t n ) u ( s ) ( u ( s ) u ( w ) ) g ( u ( w ) ) d u ( w )
and
| | u ( s ) u ( t n ) | | 1 = | | t n s u ( w ) d w | | 1 ( s t n ) s σ + α / 2 1 ,
we have
1 Γ ( α ) t n t n + 1 ( t n + 1 s ) α 1 ( g ( u ) g ˜ 1 ( u ) ) ( s ) d s 1 | | I 1 , g 1 | | 1 + | | I 2 , g 1 | | 1 ,
where
| | I 1 , g 1 | | 1 t n t n + 1 ( t n + 1 s ) α 1 | | u ( s ) u ( t n ) | | 1 | | g ( u ( t n ) ) | | 1 d s C t n t n + 1 ( t n + 1 s ) α 1 ( s t n ) s σ + α / 2 1 d s C τ t n + 1 1 1 / γ t n t n + 1 ( t n + 1 s ) α 1 s σ + α / 2 1 d s C τ t n + 1 σ + 3 α / 2 1 / γ .
Furthermore,
| | I 2 , g 1 | | 1 1 Γ ( α ) t n t n + 1 ( t n + 1 s ) α 1 | | u ( t n ) u ( s ) ( u ( s ) u ( w ) ) g ( u ( w ) ) d u ( w ) | | 1 d s C t n t n + 1 ( t n + 1 s ) α 1 ( s t n ) 2 s 2 σ + α 2 d s C τ n 2 t n t n + 1 ( t n + 1 s ) α 1 s 2 σ + α 2 d s C τ 2 t n + 1 2 σ + 2 α 2 / γ .
We complete the proof using the bounds for | | I 1 , g 1 | | 1 and | | I 2 , g 1 | | 1 . □
Lemma 7. 
Assume the conditions of Lemma 5. Then, τ is sufficiently small and we have the estimate
1 Γ ( α ) t n t n + 1 ( t n + 1 s ) α 1 ( g ( u ) g ˜ 2 ( u ) ) ( s ) d s 1 C τ 2 t n + 1 σ + 3 α / 2 2 / γ ,
for
γ > max 2 σ + α / 2 + 1 , 2 σ + 3 α / 2 3
with σ + 3 α / 2 > 3 .
Proof. 
The proof follows from Lemma 5 and is omitted. □
Theorem 1. 
Assume the conditions in Lemmas 4–7. Then, the error bounds
| | u n + 1 p u ( t n + 1 ) | | 1 C τ and | | u n + 1 u ( t n + 1 ) | | 1 C τ 2
hold uniformly on 0 t n T , for a sufficiently small τ, where u, u p are the solutions obtained from the predictor and corrector schemes in (13).
Proof. 
Noting that A β 2 is symmetric, positive definite, and using the error bounds obtained in Lemmas 4–7, we obtain
u ( t n + 1 ) u n + 1 p u ( t n ) u n + 1 Γ ( α ) A β 2 t n t n + 1 ( t n + 1 s ) α 1 ( u u ˜ ) ( s ) d s 1 + 1 Γ ( α ) t n t n + 1 ( t n + 1 s ) α 1 ( g ( u ) g 1 ˜ ( u ) ) ( s ) d s 1 + Q n , u e Q n , u a 1 + Q n , g e Q n , g a 1 u ( t n ) u n + C τ 2 t n + 1 σ + 3 α / 2 2 / γ + C τ t n + 1 σ + α 1 / γ u ( t n ) u n + C τ t n + 1 σ + α 1 / γ
Similarly,
u ( t n + 1 ) u n + 1 u ( t n ) u n + C τ 2 t n + 1 σ + 3 α / 2 2 / γ .
The proof is completed through mathematical induction. □

3.2. Stability Analysis

Definition 1. 
The scheme given by (13) is said to be stable if there is K > 0 , independent of τ and n, so that
| | u n u ^ n | | K | | u 0 u ^ 0 | | , n = 1 , 2 , , M .
where u n and u ^ n satisfy (13) with the initial data u 0 and u ^ 0 .
Lemma 8. 
If 0 < α 1 and t j = ( j τ ) γ , j = 0 , 1 , , n , γ 1 and τ = T 1 / γ N . Then, the following estimate holds
a j , n K α τ n τ 0 ( t n + 1 t 0 ) α , j = 0 , τ n τ j 1 ( t n + 1 t j 1 ) α , 1 j n 1 , τ n 1 α , j = n ,
where
K α = max α + 1 Γ ( α + 2 ) , 2 ( α + 1 ) Γ ( α + 2 ) , 2 α + 1 α + 1 .
Proof. 
For j = 0 , we have
τ 0 Γ ( α + 2 ) a 0 , n = ( t n + 1 t 1 ) α + 1 ( t n t 1 ) α + 1 + ( t n t 0 ) α [ ( t n ( α + 1 ) t 1 ] ( t n + 1 t 0 ) α [ ( t n + 1 ( α + 1 ) t 1 ] ( t n + 1 t 1 ) α + 1 ( t n t 1 ) α + 1
With ξ ( t n , t n + 1 ) and using the MVT, we have
τ 0 Γ ( α + 2 ) a 0 , n ( α + 1 ) τ n ( ξ t 1 ) α
which implies
a 0 , n ( α + 1 ) Γ ( α + 2 ) τ n τ 0 ( t n + 1 t 1 ) α .
For 1 j n 1 ,
Γ ( α + 2 ) a j , n = 1 τ j 1 [ ( t n + 1 t j 1 ) α + 1 ( t n t j 1 ) α + 1 + ( t n t j ) α ( t n t j 1 ) ( t n + 1 t j ) α ( t n + 1 t j 1 ) ] + 1 τ j [ ( t n + 1 t j + 1 ) α + 1 ( t n t j + 1 ) α + 1 + ( t n t j ) α ( t n t j + 1 ) ( t n + 1 t j ) α ( t n + 1 t j + 1 ) ] 1 τ j 1 [ ( t n + 1 t j 1 ) α + 1 ( t n t j 1 ) α + 1 + ( t n + 1 t j + 1 ) α + 1 ( t n t j + 1 ) α ] .
Again, applying the MVT, we have
a j , n 2 ( α + 1 ) Γ ( α + 2 ) τ n τ j 1 ( t n + 1 t j 1 ) α
For j = n ,
τ n 1 Γ ( α + 2 ) a n , n = ( τ n + τ n 1 ) α + 1 τ n 1 α + 1 τ n α [ τ n + ( α + 1 ) τ n 1 ] = τ n 1 α + 1 [ α ( α + 1 ) 2 ! ( τ n 1 τ n ) 1 α + ( α 1 ) α ( α + 1 ) 3 ! ( τ n 1 τ n ) 2 α + 1 ] τ n 1 α + 1 [ α ( α + 1 ) 2 ! + ( α 1 ) α ( α + 1 ) 3 ! + 1 ] = ( 2 α + 1 α 2 ) τ n 1 α + 1 ,
where we have used the generalized binomial theorem to arrive at the last inequality and the fact that τ j is non-decreasing.
Therefore, a n , n K α τ n 1 α .
Lemma 9. 
Assume that 0 < α 1 and
a j , n = K α τ n τ j 1 ( t n + 1 t j 1 ) α , ( j = 1 , 2 , , n 1 ) ,
and
a 0 , n = K α τ n τ 0 ( t n + 1 t 1 ) α
for t j = ( j τ ) γ , j = 0 , 1 , , n , n = 1 , 2 , , M . Let g 0 be a positive number and assume the sequence { ψ j } satisfies
ψ 0 g 0 , ψ n j = 0 n 1 a j , n ψ j + C 0 g 0 ,
Then,
ψ n C 0 g 0 , n = 1 , 2 , , M .
Proof. 
The proof uses a modification of that of [35] (Lemma 3.3). □
Theorem 2. 
Suppose that u j ( j = 1 , 2 , , N ) (3) due to Scheme (13) and where g ( u ) is Lipschitz in Ω × ( 0 , T ] (with respect to u). Then, (13) is stable.
Proof. 
We start by considering a history-term perturbation in the form
H n a = j = 0 n 1 a j , n A β 2 u j + g ( u j + u j ) g ( u j ) + a n , n A β 2 u n + g ( u n + u n ) g ( u n ) .
By using the positive definiteness of A β 2 , the fact that g ( u ) is Lipschitz continuous, and Lemma 8, we obtain
| | H n a | | K j = 0 n 1 a j , n | | u j | | + τ n α | | u n | | .
Here, K is assumed to be a positive constant. The perturbation of Equation (13) works out to be
u n + 1 p = Γ ( α + 2 ) I + τ n α A β 2 1 Γ ( α + 2 ) I α τ n α A β 2 u n + τ n α ( α + 1 ) ( g ( u n + u n ) g ( u n ) ) + Γ ( α + 2 ) I + τ n α A β 2 1 Γ ( α + 2 ) H n a , u n + 1 = Γ ( α + 2 ) I + τ n α A β 2 1 Γ ( α + 2 ) I α τ n α A β 2 u n + τ n α α ( g ( u n + u n ) g ( u n ) ) , + Γ ( α + 2 ) I + τ n α A β 2 1 τ n α ( g ( u n + 1 p + u n + 1 p ) g ( u n + 1 p ) ) + Γ ( α + 2 ) H n a .
By the positive definiteness of A β 2 , it follows that 0 < C < 1 , where
C = Γ ( α + 2 ) I + τ n α A β 2 1 Γ ( α + 2 ) I α τ n α A β 2 .
Therefore,
| | u n + 1 p | | C | | u n | | + K 1 τ max α | | u n | | + τ max α | | u n | | + j = 0 n 1 a j , n | | u j | | , | | u n + 1 | | C | | u n | | + K 2 τ max α ( | | u n | | + | | u n + 1 p | | ) + τ max α | | u n | | + j = 0 n 1 a j , n | | u j | | ,
where C , K 1 , K 2 are constants, with τ max = max τ j j = 0 n . We show the remaining part by employing induction. For n = 0 and a sufficiently small τ max , it follows that
| | u 1 p | | | | u 0 | | and | | u 1 | | | | u 0 | | .
Suppose that
| | u j | | | | u 0 | | , j = 1 , 2 , , n .
We consider j = n + 1 , for u n + 1 p , that is,
| | u n + 1 p | | C | | u n | | + K 1 τ max α | | u n | | + τ max α | | u n | | + j = 0 n 1 a j , n | | u j | | C 0 | | u n | | + K 1 j = 0 n 1 a j , n | | u j | | | | u 0 | | ,
where 0 < C 0 = C + 2 K 1 τ max α < 1 for a sufficiently small τ max , where Lemma 9 has been used.
We have
| | u n + 1 | | C | | u n | | + K 2 τ m a x α ( | | u n | | + | | u n + 1 p | | ) + τ m a x α | | u n | | + j = 0 n 1 a j , n | | u j | | C 1 | | u n | | + K 2 j = 0 n 1 a j , n | | u j | | | | u 0 | | ,
where 0 < C 1 = C + 3 K 2 τ max α < 1 . This completes the proof. □

4. Numerical Illustrations

Here, we corroborate the analysis through the empirical study of the convergence rate for different test problems. For the examples that we consider in this section, the convergence rate ( C R ) is given by
C R = log 2 Error M 2 / Error M ,
where
Error M = u M u M 2
and u M is the vector of the solution with M mesh points. The numerical examples are the same as those given in Biala and Khaliq [18] and the results for γ = 1 can be found there.
Example 1. 
We consider
c D 0 , t α u = ( Δ ) β 2 u + g ( u ) , t ( 0 , 1 ] , x [ 0 , 1 ]
φ ( x ) = x 2 ( 1 x ) 2 .
As the initial data φ ( x ) C ( Ω ) H 0 1 ( Ω ) , the regularity property in (15) holds true for any σ I σ = 0 , α 4 . Consider the problem g ( u ) = 0 , whose solution can be seen to be
u ( x , t ) = n = 1 4 12 + n 2 π 2 1 + ( 1 ) n n 5 π 5 E α ( ( n π ) β t α ) sin ( n π x ) ,
Here, E α is the one-parameter Mittag–Leffler function. Table 1 and Table 2 show the error and the convergence rates when g ( u ) = 0 and g ( u ) = u 2 , respectively, using the L 2 norm. We used a small step size of d x = 0.001 so that the error in time is dominant. By Theorem 1, we expect to have O ( τ 2 ) convergence for
γ > max 2 σ + α / 2 + 1 , 2 σ + 3 α / 2 3 = max 8 3 α + 4 , 8 7 α 12 .
In fact, the second term in the maximum function is not necessary since 0 < α 1 . In order not to pepper the text with so many tables with different values of γ > 8 3 α + 4 , we show the results for only two values of γ (one that is slightly greater than 8 3 α + 4 and another that is slightly lower) to validate our theoretical order of convergence. We observe that with γ = 8 3 α + 5 , the O ( τ 3 / 2 + ϵ ) for some ϵ ( 0 , 1 / 2 ) is achieved. However, for γ = 8 3 α + 3 , the O ( τ 2 ) is obtained, which corroborates our theoretical analysis. These observations are further depicted in Figure 1 and Figure 2, where we fit a linear line for the logarithm (base 10) of M 1 and the corresponding errors. The values of the slope in these figures that depict the rates of convergence for different values of α and γ further support our theoretical observations.
Example 2. 
Let us consider a two-dimensional time–space reaction–diffusion problem of fractional order
c D 0 , t α u = ( Δ ) β 2 u + g ( u ) , t ( 0 , 1 ] , ( x , y ) [ 0 , 1 ] × [ 0 , 1 ]
u ( x , y , 0 ) = x y ( 1 x ) ( 1 y )
with a Dirichlet homogeneous boundary condition. We first solve the problem when g ( u ) = 0 . The solution in this case is given in Yang et al. [36], by
u ( x , y , t ) = n = 1 m = 1 E α λ n , m β 2 t α c n , m ϕ n , m ( x , y ) . λ n , m = ( n 2 + m 2 ) π 2 , ϕ n , m ( x , y ) = 2 sin ( n π x ) sin ( n π y ) , c n , m = 0 1 0 1 u v ( 1 u ) ( 1 v ) ϕ n , m ( u , v ) d u d v .
A space-step size of d x = 0.008 (for CPU memory constraints) is used in this problem. Similar to the 1D problem, the results here (see Table 3 and Table 4) show that the O ( τ 2 ) order of convergence is achieved when γ > 8 3 α + 4 . Figure 3 shows the exact and numerical solutions with β = 1.4 and α = 0.2 .

5. Conclusions

In this work, we developed a numerical scheme over time-graded meshes for nonlinear time–space fractional reaction–diffusion equations. The analysis uses the regularity properties of the solutions of the proposed equations and an O ( τ 2 ) order of convergence is achieved. The regularity properties of the solution to this class of problem are used to improve the convergence properties of the proposed numerical scheme on time-graded meshes. The stability results are discussed and proved. Furthermore, the sharp error estimates for an optimal O ( τ 2 ) rate of convergence are proved. Some examples are provided to demonstrate the efficiency and accuracy of our proposed scheme across different values of the fractional order α .

Author Contributions

Conceptualization, Y.O.A., T.A.B., O.S.I., A.Q.M.K. and B.A.W.; Methodology, Y.O.A., T.A.B. and O.S.I.; Software, Y.O.A. and T.A.B.; Formal analysis, Y.O.A., T.A.B. and O.S.I.; Writing—original draft, Y.O.A., T.A.B. and O.S.I.; Writing—review & editing, A.Q.M.K. and B.A.W.; Supervision, O.S.I., A.Q.M.K. and B.A.W. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Di Giuseppe, E.; Moroni, M.; Caputo, M. Flux in porous media with memory: Models and experiments. Transp. Porous Media 2010, 83, 479–500. [Google Scholar] [CrossRef]
  2. Baleanu, D.; Diethelm, K.; Scalas, E.; Trujillo, J.J. Fractional Calculus: Models and Numerical Methods; World Scientific: Singapore, 2016; Volume 5. [Google Scholar]
  3. Podlubny, I. Fractional Differential Equations. In Mathematics in Science and Engineering 198; Academic Press: Cambridge, MA, USA, 1999. [Google Scholar]
  4. Caputo, M.; Mainardi, F. A new dissipation model based on memory mechanism. Pure Appl. Geophys. 1971, 91, 134–147. [Google Scholar] [CrossRef]
  5. Fomin, S.; Chugunov, V.; Hashida, T. Mathematical modeling of anomalous diffusion in porous media. Fract. Differ. Calc. 2011, 1, 1–28. [Google Scholar] [CrossRef]
  6. Metzler, R.; Klafter, J. The random walk’s guide to anomalous diffusion: A fractional dynamics approach. Phys. Rep. 2000, 339, 1–77. [Google Scholar] [CrossRef]
  7. Cloot, A.; Botha, J. A generalised groundwater flow equation using the concept of non-integer order derivatives. Water SA 2006, 32, 1–7. [Google Scholar] [CrossRef] [Green Version]
  8. Iaffaldano, G.; Caputo, M.; Martino, S. Experimental and theoretical memory diffusion of water in sand. Hydrol. Earth Syst. Sci. Discuss. 2005, 2, 1329–1357. [Google Scholar] [CrossRef] [Green Version]
  9. Podlubny, I. Fractional-order systems and fractional-order controllers. Inst. Exp. Physics, Slovak Acad. Sci. Kosice 1994, 12, 1–18. [Google Scholar]
  10. Iyiola, O.; Oduro, B.; Zabilowicz, T.; Iyiola, B.; Kenes, D. System of time fractional models for COVID-19: Modeling, analysis and solutions. Symmetry 2021, 13, 787. [Google Scholar] [CrossRef]
  11. Biala, T.; Khaliq, A. A fractional-order compartmental model for the spread of the COVID-19 pandemic. Commun. Nonlinear Sci. Numer. Simul. 2021, 98, 1–19. [Google Scholar] [CrossRef]
  12. Biala, T.; Afolabi, Y.; Khaliq, A. How Efficient is Contact Tracing in Mitigating the Spread of COVID-19? A Mathematical Modeling Approach. Appl. Math. Model. 2022, 103, 714–730. [Google Scholar] [CrossRef]
  13. Furati, K.; Sarumi, I.; Khaliq, A. Fractional model for the spread of Covid-19 subject to government intervention and public perception. Appl. Math. Model. 2021, 95, 89–105. [Google Scholar] [CrossRef] [PubMed]
  14. Lu, Z.; Yu, Y.; Chen, Y.; Ren, G.; Xu, C.; Wang, S.; Yin, Z. A fractional-order SEIHDR model for COVID-19 with inter-city networked coupling effects. Nonlinear Dyn. 2020, 101, 1717–1730. [Google Scholar] [CrossRef] [PubMed]
  15. Ilic, M.; Liu, F.; Turner, I.; Anh, V. Numerical approximation of a fractional-in-space diffusion equation (II)- with nonhomogenous boundary conditions. Fract. Calc. Appl. Anal. 2006, 9, 333–349. [Google Scholar]
  16. Iyiola, O.; Wade, B. Exponential integrator methods for systems of non-linear space-fractional models with super-diffusion processes in pattern formation. Comput. Math. Appl. 2018, 75, 3719–3736. [Google Scholar] [CrossRef]
  17. Mustapha, K.; Abdallah, B.; Furati, K.M. A Discontinuous Petrov–Galerkin Method for Time-Fractional Diffusion Equations. SIAM J. Numer. Anal. 2014, 52, 2512–2529. [Google Scholar] [CrossRef] [Green Version]
  18. Biala, T.; Khaliq, A. Parallel algorithms for nonlinear time–space fractional parabolic PDEs. J. Comput. Phys. 2018, 375, 135–154. [Google Scholar] [CrossRef]
  19. Lin, Y.; Xu, C. Finite Difference/Spectral approximations for the time-fractional diffusion equation. J. Comput. Phys. 2007, 225, 1533–1552. [Google Scholar] [CrossRef]
  20. Iyiola, O.; Asante-Asamani, E.; Furati, K.; Khaliq, A.; Wade, B. Efficient time discretization scheme for nonlinear space fractional reaction-diffusion equations. Int. J. Comput. Math. 2018, 95, 1–17. [Google Scholar] [CrossRef]
  21. Alikhanov, A.A. A new difference scheme for the time fractional diffusion equation. J. Comput. Phys. 2015, 280, 424–438. [Google Scholar] [CrossRef] [Green Version]
  22. Gao, G.; Sun, Z.; Zhang, H. A new fractional numerical differentiation formula to approximate the Caputo fractional derivative and its applications. J. Comput. Phys. 2014, 259, 33–50. [Google Scholar] [CrossRef]
  23. Jin, B.; Li, B.; Zhou, Z. An analysis of the Crank-Nicolson method for subdiffusion. IMA J. Numer. Anal. 2017, 38, 518–541. [Google Scholar] [CrossRef] [Green Version]
  24. Brunner, H. The Numerical Solution of Weakly Singular Volterra Integral Equations by Collocation on Graded Meshes. Math. Comput. 1985, 45, 417–437. [Google Scholar] [CrossRef]
  25. Zhang, Y.; Sun, Z.; Liao, H. Finite difference methods for the time fractional diffusion equation on non-uniform meshes. J. Comput. Phys. 2014, 265, 195–210. [Google Scholar] [CrossRef]
  26. Lyu, P.; Vong, S. A high-order method with a temporal nonuniform mesh for a time-fractional Benjamin-Bona-Mahony equation. J. Sci. Comput. 2019, 80, 1607–1628. [Google Scholar] [CrossRef]
  27. Stynes, M.; OŔiordan, E.; Gracia, J.L. Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation. SIAM J. Numer. Anal. 2017, 55, 1057–1079. [Google Scholar] [CrossRef] [Green Version]
  28. Liao, H.; Li, D.; Zhang, J. Sharp error estimate of the nonuniform L1 formula for linear reaction-subdiffusion equations. SIAM J. Numer. Anal. 2018, 56, 1112–1133. [Google Scholar] [CrossRef]
  29. Kopteva, N. Error analysis of the L1 method on graded and uniform meshes for a fractional- derivative problem in two and three dimensions. Math. Comp. 2019, 88, 2135–2155. [Google Scholar] [CrossRef]
  30. Wang, K.; Zhou, Z. High-Order Time-Stepping Schemes for semilinear subdiffusion equations. SIAM J. Numer. Anal. 2020, 58, 3226–3250. [Google Scholar] [CrossRef]
  31. Mustapha, K. An L1 Approximation for a Fractional Reaction-Diffusion Equation, a Second-Order Error Analysis over Time-Graded Meshes. SIAM J. Numer. Anal. 2020, 58, 1319–1338. [Google Scholar] [CrossRef]
  32. Simpson, D.P. Krylov Subspace Methods for Approximating Functions of Symmetric Positive Definite Matrices with Applications to Applied Statistics and Anomalous Diffusion. Ph.D. Thesis, Queensland University of Technology, Brisbane, Australia, 2008. [Google Scholar]
  33. McLean, W.; Mustapha, K. A second-order accurate numerical method for a fractional wave equation. Numer. Math. 2007, 105, 481–510. [Google Scholar] [CrossRef]
  34. McLean, W. Regularity of solutions to a time-fractional diffusion equation. ANZIAM J. 2010, 52, 123–138. [Google Scholar] [CrossRef] [Green Version]
  35. Li, C.; Yi, Q.; Chen, A. Finite difference methods with non-uniform meshes for nonlinear fractional differential equations. J. Comput. Phys. 2016, 316, 614–631. [Google Scholar] [CrossRef]
  36. Yang, Q.; Turner, I.; Liu, F.; Ilić, M. Novel numerical methods for solving the time-space fractional diffusion equation in two dimensions. SIAM J. Sci. Comput. 2011, 33, 1159–1180. [Google Scholar] [CrossRef]
Figure 1. Log-log error plots for Example 1 with g ( u ) = 0 , showing the rate of convergence of the scheme.
Figure 1. Log-log error plots for Example 1 with g ( u ) = 0 , showing the rate of convergence of the scheme.
Fractalfract 07 00040 g001
Figure 2. Log -log error plots for Example 1 with g ( u ) = u 2 , showing the rate of convergence of the scheme.
Figure 2. Log -log error plots for Example 1 with g ( u ) = u 2 , showing the rate of convergence of the scheme.
Fractalfract 07 00040 g002
Figure 3. Plots of exact (left) and numerical solutions (right) with β = 1.4 and α = 0.2 . .
Figure 3. Plots of exact (left) and numerical solutions (right) with β = 1.4 and α = 0.2 . .
Fractalfract 07 00040 g003
Table 1. g ( u ) = 0 with β = 1.2 .
Table 1. g ( u ) = 0 with β = 1.2 .
α = 0.4 α = 0.8
γ = 8 3 α + 5 γ = 8 3 α + 3 γ = 8 3 α + 5 γ = 8 3 α + 3
M ErrorCRErrorCRErrorCRErrorCR
10 1.324 × 10 3 4.925 × 10 4 1.094 × 10 3 9.290 × 10 4
20 3.884 × 10 4 1.7690 1.214 × 10 4 2.0201 2.896 × 10 4 1.9167 2.340 × 10 4 1.9891
40 1.143 × 10 4 1.7648 2.990 × 10 5 2.0217 7.657 × 10 5 1.9194 5.854 × 10 5 1.9989
80 3.367 × 10 5 1.7629 7.328 × 10 6 2.0287 2.021 × 10 5 1.9216 1.459 × 10 5 2.0046
160 9.889 × 10 6 1.7677 1.751 × 10 6 2.0656 5.300 × 10 6 1.9311 3.599 × 10 6 2.0193
320 2.859 × 10 6 1.7901 3.707 × 10 7 2.2394 1.354 × 10 6 1.9684 8.518 × 10 7 2.0792
Table 2. g ( u ) = u 2 with β = 1.6 .
Table 2. g ( u ) = u 2 with β = 1.6 .
α = 0.3 α = 0.7
γ = 8 3 α + 5 γ = 8 3 α + 3 γ = 8 3 α + 5 γ = 8 3 α + 3
M ErrorCRErrorCRErrorCRErrorCR
10 9.400 × 10 4 2.356 × 10 4 6.739 × 10 4 3.985 × 10 4
20 2.804 × 10 4 1.7450 5.495 × 10 5 2.1002 1.854 × 10 4 1.8621 1.009 × 10 4 1.9823
40 8.348 × 10 5 1.7482 1.259 × 10 5 2.1258 5.089 × 10 5 1.8649 2.523 × 10 5 1.9991
80 2.454 × 10 5 1.7664 2.655 × 10 6 2.2456 1.398 × 10 5 1.8638 6.285 × 10 6 2.0051
160 6.914 × 10 6 1.8276 3.816 × 10 7 2.7986 3.839 × 10 6 1.8646 1.564 × 10 6 2.0071
320 1.709 × 10 6 2.0161 9.071 × 10 8 2.0728 1.052 × 10 6 1.8676 3.885 × 10 7 2.0090
Table 3. g ( u ) = 0 with β = 1.4 .
Table 3. g ( u ) = 0 with β = 1.4 .
α = 0.2 α = 0.6
γ = 8 3 α + 5 γ = 8 3 α + 3 γ = 8 3 α + 5 γ = 8 3 α + 3
M ErrorCRErrorCRErrorCRErrorCR
10 1.369 × 10 2 3.031 × 10 3 1.316 × 10 2 6.117 × 10 3
20 4.178 × 10 3 1.7126 6.349 × 10 4 2.2554 3.441 × 10 3 1.9352 1.368 × 10 3 2.1609
40 1.270 × 10 3 1.7176 1.299 × 10 4 2.2886 9.647 × 10 4 1.8346 3.370 × 10 4 2.0209
80 3.804 × 10 4 1.7395 1.953 × 10 5 2.7343 2.666 × 10 4 1.8556 7.703 × 10 5 2.1293
160 1.076 × 10 4 1.8225 6.519 × 10 6 1.5829 6.912 × 10 5 1.9473 1.234 × 10 5 2.6424
Table 4. g ( u ) = u 3 with β = 1.8 .
Table 4. g ( u ) = u 3 with β = 1.8 .
α = 0.5 α = 0.9
γ = 8 3 α + 5 γ = 8 3 α + 3 γ = 8 3 α + 5 γ = 8 3 α + 3
M ErrorCRErrorCRErrorCRErrorCR
10 1.049 × 10 2 6.170 × 10 3 5.917 × 10 3 2.970 × 10 2
20 2.962 × 10 3 1.8250 6.315 × 10 4 3.2883 1.551 × 10 3 1.9312 1.697 × 10 3 4.1296
40 8.509 × 10 4 1.7993 1.549 × 10 4 2.0278 4.060 × 10 4 1.9341 8.376 × 10 5 4.3404
80 2.454 × 10 4 1.7937 3.790 × 10 5 2.0307 1.047 × 10 4 1.9556 1.807 × 10 5 2.2127
160 7.107 × 10 5 1.7880 9.292 × 10 6 2.0244 2.690 × 10 5 1.9600 4.539 × 10 6 1.9933
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Afolabi, Y.O.; Biala, T.A.; Iyiola, O.S.; Khaliq, A.Q.M.; Wade, B.A. A Second-Order Crank-Nicolson-Type Scheme for Nonlinear Space–Time Reaction–Diffusion Equations on Time-Graded Meshes. Fractal Fract. 2023, 7, 40. https://doi.org/10.3390/fractalfract7010040

AMA Style

Afolabi YO, Biala TA, Iyiola OS, Khaliq AQM, Wade BA. A Second-Order Crank-Nicolson-Type Scheme for Nonlinear Space–Time Reaction–Diffusion Equations on Time-Graded Meshes. Fractal and Fractional. 2023; 7(1):40. https://doi.org/10.3390/fractalfract7010040

Chicago/Turabian Style

Afolabi, Yusuf O., Toheeb A. Biala, Olaniyi S. Iyiola, Abdul Q. M. Khaliq, and Bruce A. Wade. 2023. "A Second-Order Crank-Nicolson-Type Scheme for Nonlinear Space–Time Reaction–Diffusion Equations on Time-Graded Meshes" Fractal and Fractional 7, no. 1: 40. https://doi.org/10.3390/fractalfract7010040

APA Style

Afolabi, Y. O., Biala, T. A., Iyiola, O. S., Khaliq, A. Q. M., & Wade, B. A. (2023). A Second-Order Crank-Nicolson-Type Scheme for Nonlinear Space–Time Reaction–Diffusion Equations on Time-Graded Meshes. Fractal and Fractional, 7(1), 40. https://doi.org/10.3390/fractalfract7010040

Article Metrics

Back to TopTop