Next Article in Journal
Lyapunov-Type Inequalities for Systems of Riemann-Liouville Fractional Differential Equations with Multi-Point Coupled Boundary Conditions
Next Article in Special Issue
Probing Families of Optical Soliton Solutions in Fractional Perturbed Radhakrishnan–Kundu–Lakshmanan Model with Improved Versions of Extended Direct Algebraic Method
Previous Article in Journal
Synchronization of Discrete-Time Fractional-Order Complex-Valued Neural Networks with Distributed Delays
Previous Article in Special Issue
Doubling Smith Method for a Class of Large-Scale Generalized Fractional Diffusion Equations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Local Error Estimate of an L1-Finite Difference Scheme for the Multiterm Two-Dimensional Time-Fractional Reaction–Diffusion Equation with Robin Boundary Conditions

1
Department of Mathematics, Beijing Jiaotong University, Beijing 100044, China
2
Institute of Process Engineering, Chinese Academy of Sciences, Beijing 100190, China
*
Author to whom correspondence should be addressed.
Submission received: 27 April 2023 / Revised: 29 May 2023 / Accepted: 29 May 2023 / Published: 1 June 2023

Abstract

:
In this paper, the numerical method for a multiterm time-fractional reaction–diffusion equation with classical Robin boundary conditions is considered. The full discrete scheme is constructed with the L1-finite difference method, which entails using the L1 scheme on graded meshes for the temporal discretisation of each Caputo fractional derivative and using the finite difference method on uniform meshes for spatial discretisation. By dealing with the discretisation of Robin boundary conditions carefully, sharp error analysis at each time level is proven. Additionally, numerical results that can confirm the sharpness of the error estimates are presented.

1. Introduction

In recent years, fractional calculus, which is considered to be a generalisation of classical derivatives and integrals to non-integer order, has become a powerful modelling tool that is more flexible and precise for describing physical problems than integer calculus. The fractional system has been widely used in engineering, physical science, chemical science, biology, and a variety of other subjects, for which it has gradually become an essential component. For more details on fractional calculus, see [1,2,3,4,5].
At present, it is not generalised enough to consider a numerical solution of the initial boundary value problem with only the time fractional derivative term with the order α ( 0 , 1 ) , such as in [6]. On this basis, more attention is being paid to the summation form of the time fractional derivative with the order
0 < α L < < α 2 < α 1 < 1 .
where L is a positive integer. At the initial time, the typical solutions of such problems have a key factor that must be considered (as in [6]); this factor is weak singularity which significantly complicates analysis. Now, many time-fractional initial-boundary value problems with Robin boundary conditions are widely used in the research fields of heat equation, biomathematics, and so on [7,8,9]. That is the main reason why this type of boundary condition is considered in this paper.
The problem that we study in the spatial domain is Ω : = ( 0 , 1 ) 2 with closure Ω ¯ = [ 0 , 1 ] 2 . Define the boundary as Ω = Ω ¯ Ω . Set Q = Ω × ( 0 , T ] and Q ¯ = Ω ¯ × [ 0 , T ] where T > 0 be fixed.
Based on the above description, the purpose of this paper is to propose the following multiterm time-fractional reaction–diffusion problem numerically.
l = 1 L q l D t α l u ( x , y , t ) Δ u ( x , y , t ) + c ( x , y ) u ( x , y , t ) = F ( x , y , t ) for ( x , y , t ) Q ,
with the initial condition
u ( x , y , 0 ) = u 0 ( x , y ) for ( x , y ) Ω ,
and Robin boundary conditions
σ u ( x , y , t ) + n u ( x , y , t ) = g ( x , y , t ) for ( x , y , t ) Ω , 0 < t T .
where q l and σ are given positive constants, g and u 0 are sufficiently smooth in their respective domains, F C 1 ( Q ¯ ) , and c C 1 ( Q ¯ ) with c c 0 > 0 . u 0 C 1 ( Ω ¯ ) with σ u 0 + ( / n ) u 0 = 0 on Ω . D t α l u ( x , y , t ) is defined as the temporal Caputo fractional derivative of order α l of u by
D t α l u ( x , y , t ) : = 1 Γ ( 1 α l ) 0 t ( t s ) α l s u ( x , y , s ) d s .
For ([10] Lemma 2.2) and ([11] Section 6), (1) has unique solution which satisfies the following regularity with weak initial singularity
η x η u ( x , y , t ) + η y η u ( x , y , t ) C for η = 0 , 1 , 2 , 3 , 4 ,
ν t ν u ( x , y , t ) C ( 1 + t α 1 ν ) for ν = 0 , 1 , 2 ,
where C is some fixed constant.
In recent years, the introduction of the classic L1 scheme to the discrete Caputo derivative has received widespread attention [12,13]. To recover the convergence rate, researchers have used the L1 scheme on graded meshes [6,11,14]. Analysis of the local convergence rate is mathematically interesting [15,16], as the local convergence rate on every time node is sharper than the global one. This method has wide applicability. When considering practical problems such as [17,18,19], it can be combined with the finite difference and finite element methods in the space direction [6,20].
To avoid discrete errors at the boundary, Dirichlet boundary conditions are usually considered in the local error analysis because the boundary values are known. For Robin boundary conditions, to ensure global accuracy, we need to find a suitable boundary discretisation method. One of the novelties of this paper is mitigation of the difficulty caused by Robin boundary conditions in the discretisation process. At present, there are many papers consider global convergence of the time fractional problem with Robin boundary conditions [10,21,22]; But no local in time error analysis for multi-term time-fractional problems with Robin boundary conditions has been considered. This is our motivation for completing this paper. The highlights of this paper can be summarized as the following:
  • By using the L1-finite difference method to solve (1a), we have propose a discrete scheme at the boundary which can match the second-order central difference scheme at interior points.
  • When considering local errors, Dirichlet boundary conditions are usually considered [20,23,24]. In this paper, the boundary conditions are extended to Robin boundary conditions.
The outline of the paper is as follows. In Section 2, the L1-finite difference method that will be used to solve (1) is described. In Section 3, the local error of the L1-finite difference method is analysed. Then, in Section 4, numerical examples are given to verify the local error results.
Notation: throughout the paper, C is used as a generic constant to solve (1) numerically, and may take a different value each time it appears. Meanwhile, it is related to the information of the problem (1) but is independent of ( x , y , t ) and of any mesh.

2. L1-Finite Difference Method for (1)

We shall consider the L1-finite difference method for construction of the fully discrete scheme for the problem (1). At the boundary and inner points, the discrete scheme which can match each other’s accuracy is constructed.
Let M and N be positive integers. Set t m = T ( m / M ) r for m = 0 , 1 , , M , and denote the time step τ m = t m t m 1 for m = 1 , 2 , , M . The mesh grading r 1 will be chosen later. Set the spatial step as h = 1 / N . We divide Ω ¯ into ( N 1 ) × ( N 1 ) intervals; the mesh point is ( x i , y j ) with x i = i h and y j = j h , where ( 0 i , j N ) . Let
Ω ¯ h = { ( x i , y j ) | 0 i , j N } , Ω h = Ω ¯ h Ω , Ω h = Ω ¯ h Ω .
Let ( i , j ) Ω ¯ h represent ( x i , y j ) Ω ¯ h to simplify the notation. Similarly, set ( i , j ) Ω h and let ( i , j ) Ω h represent ( x i , y j ) Ω h and ( x i , y j ) Ω h , respectively. Thus, our mesh is
{ ( x i , y j , t m ) : i , j = 0 , 1 , , N and m = 1 , 2 , , M } .
At each mesh point ( x i , y j , t m ) , the computed approximation to the analytical solution u will be denoted by u i , j m . Define the grid functions
c i , j = c ( x i , y j ) , ( c 1 ) i , j = x c ( x i , y j ) , ( c 2 ) i , j = y c ( x i , y j ) , F i , j m = F ( x i , y j , t m ) , ( F 1 ) i , j m = x F ( x i , y j , t m ) , ( F 2 ) i , j m = y F ( x i , y j , t m ) .
where ( i , j ) Ω ¯ h , m = 1 , 2 , , M .
The Caputo fractional derivative D t α u can be expressed as
D t α u ( x , y , t ) : = 1 Γ ( 1 α ) k = 0 m 1 s = t k t k + 1 ( t s ) α s u ( x , y , s ) d s .
The L1 scheme, which is used to approximate the Caputo fractional derivative to obtain the discretisation of each time-fractional term q l D t α l u ( x i , y j , t m ) .
q l D M α l u i , j m : = q l Γ ( 1 α l ) k = 0 m 1 u i , j k + 1 u i , j k τ k + 1 s = t k t k + 1 ( t m s ) α l d s = q l Γ ( 2 α l ) k = 0 m 1 u i , j k + 1 u i , j k τ k + 1 ( t m t k ) 1 α l ( t m t k + 1 ) 1 α l
for l = 1 , 2 , , L .
In ([20] Lemma 4), the truncation error has the following estimate
l = 1 L q l D M α l u i , j m l = 1 L q l D t α l u ( x i , y j , t m ) C m min { 2 α 1 , r ( α 1 + 1 ) } .
For any grid function v = { v i , j | 0 i , j N } , the spatial difference operators are defined as follows:
δ x v i , j = 1 h ( v i , j v i 1 , j ) , δ y v i , j = 1 h ( v i , j v i , j 1 ) , 1 i , j N ,
and
δ x 2 v i , j = 1 h δ x v i + 1 , j δ x v i , j , 1 i N 1 , 0 j N , δ y 2 v i , j = 1 h δ y v i , j + 1 δ y v i , j , 1 i N , 0 j N 1 .
We discrete the initial condition (1b) in a standard way: for i , j Ω ¯ h , set u i , j 0 = u 0 ( x i , y j ) . In the following three subsections, we will discretise (1a) and (1c) on inner points, boundary points, and corner points separately.

2.1. Inner Points

In ( i , j ) Ω h , the diffusion term Δ u ( / x ) 2 u + ( / y ) 2 u in (1a) is approximated by a standard second-order discretisation
Δ u ( x i , y j , t m ) δ x 2 u i , j m + δ y 2 u i , j m .
Then, the truncation error has the following estimate
δ x 2 u ( x i , y j , t m ) 2 x 2 u ( x i , y j , t m ) + δ y 2 u ( x i , y j , t m ) 2 y 2 u ( x i , y j , t m ) C h 2 .
In summary, we can approximate (1) on ( i , j ) Ω h with the discrete problem
l = 1 L q l D M α l u i , j m δ x 2 u i , j m δ y 2 u i , j m + c i , j u i , j m = F i , j m , for 1 m M .

2.2. Boundary Points

For brevity, we set
g 1 ( x , y , t ) : = l = 1 J q l D t α l g ( x , y , t ) 2 y 2 g ( x , y , t ) + c ( x , y ) g ( x , y , t ) + x F ( x , y , t ) , g 2 ( x , y , t ) : = l = 1 J q l D t α l g ( x , y , t ) 2 x 2 g ( x , y , t ) + c ( x , y ) g ( x , y , t ) + y F ( x , y , t ) , p 1 ( x , y , t ) : = u ( x , y , t ) x c ( x , y ) + σ c ( x , y ) u ( x , y , t ) , p 2 ( x , y , t ) : = u ( x , y , t ) y c ( x , y ) + σ c ( x , y ) u ( x , y , t ) , q 1 ( x , y , t ) : = 2 h ( δ x u ( x , y , t ) σ u ( x , y , t ) + g ( x , y , t ) ) , q 2 ( x , y , t ) : = 2 h ( δ y u ( x , y , t ) σ u ( x , y , t ) + g ( x , y , t ) ) ,
and
δ x b u ( x , y , t ) : = q 1 ( x , y , t ) h 3 ( σ l = 1 J q l D M α l u ( x , y , t ) + p 1 ( x , y , t ) σ δ y 2 u ( x , y , t ) g 1 ( x , y , t ) ) , δ y b u ( x , y , t ) : = q 2 ( x , y , t ) h 3 ( σ l = 1 J q l D M α l u ( x , y , t ) + p 2 ( x , y , t ) σ δ x 2 u ( x , y , t ) g 2 ( x , y , t ) ) .
Then, define grid functions
( g 1 ) i , j m = g 1 ( x i , y j , t m ) , δ x b u i , j m : = δ x b u ( x i , y j , t m ) , ( p 1 ) i , j m = p 1 ( x i , y j , t m ) , ( g 2 ) i , j m = g 2 ( x i , y j , t m ) , δ y b u i , j m : = δ y b u ( x i , y j , t m ) , ( p 2 ) i , j m = p 2 ( x i , y j , t m ) , ( q 1 ) i , j m = p 1 ( x i , y j , t m ) , ( q 2 ) i , j m = q 2 ( x i , y j , t m ) .
Lemma 1. 
Assume u ( · , · , t m ) C 2 ( Ω ¯ ) for every t m ; then, there exists a constant C such that
δ x b u ( 0 , y j , t m ) 2 x 2 u ( 0 , y j , t m ) C h 2 + C h m m i n { 2 α 1 , r ( α 1 + 1 ) } .
Proof. 
For boundary points ( 0 , y j , t m ) , where j = 1 , , N 1 and m = 1 , , M . (1a) and (1c) at point ( 0 , y j , t m ) is
l = 1 L q l D t α l u ( 0 , y j , t m ) Δ u ( 0 , y j , t m )   +   c ( 0 , y j ) u ( 0 , y j , t m ) = F ( 0 , y j , t m ) ,
σ u ( 0 , y j , t m ) x u ( 0 , y j , t m ) = g ( 0 , y j , t m ) .
By Taylor expansion of u ( h , y j , t m ) at point ( 0 , y j , t m )
u ( h , y j , t m ) = u ( 0 , y j , t m ) + h x u ( 0 , y j , t m ) + h 2 2 2 x 2 u ( 0 , y j , t m ) + h 3 6 3 x 3 u ( 0 , y j , t m ) + C h 4 ,
using (10), we have
2 x 2 u ( 0 , y j , t m ) = q 1 ( 0 , y j , t m ) h 3 3 x 3 u ( 0 , y j , t m ) C h 2 .
Differentiating (9) with respect to x, ( 3 / x 3 ) u can be expressed as
3 x 3 u ( 0 , y j , t m ) = l = 1 L q l D t α l x u ( 0 , y j , t m ) x 2 y 2 u ( 0 , y j , t m ) + u ( 0 , y j , t m ) x c ( 0 , y j ) + c ( 0 , y j ) x u ( 0 , y j , t m ) x F ( 0 , y j , t m ) .
Furthermore, in view of (10), we obtain
3 x 3 u ( 0 , y j , t m ) = σ l = 1 J q l D t α l u ( 0 , y j , t m ) σ 2 y 2 u ( 0 , y j , t m ) + p 1 ( 0 , y j ) g 1 ( 0 , y j , t m ) .
So, substituting (13) into (11) to replace ( 3 / x 3 ) u yields
2 x 2 u ( 0 , y j , t m ) = q 1 ( 0 , y j , t m ) h 3 ( σ l = 1 J q l D t α l u ( 0 , y j , t m ) σ 2 y 2 u ( 0 , y j , t m ) + p 1 ( 0 , y j , t m ) g 1 ( 0 , y j , t m ) ) C h 2 .
Adding ( σ h / 3 ) ( l = 1 J q l D M α l u ( 0 , y j , t m ) ) to the right-hand side and subtract it. Then, by truncation error (4), we have
2 x 2 u ( 0 , y j , t m ) = q 1 ( 0 , y j , t m ) h 3 ( σ l = 1 J q l D M α l u ( 0 , y j , t m ) σ 2 y 2 u ( 0 , y j , t m ) + p 1 ( 0 , y j ) g 1 ( 0 , y j , t m ) ) C h 2 C h m m i n { 2 α 1 , r ( α 1 + 1 ) } .
For 1 j N 1 , 1 m M , we can approximate (1) on boundary point ( 0 , y j , t m ) by the discrete problem
l = 1 L q l D M α l u 0 , j m δ x b u 0 , j m δ y 2 u 0 , j m + c 0 , j u 0 , j m = F 0 , j m
The other corner points can be treated similarly.

2.3. Corner Points

For convenience, we introduce the following functions
δ x c u ( x , y , t ) : = 1 1 σ h 3 q 1 ( x , y , t ) h 3 σ l = 1 J q l D M α l u ( x , y , t ) + p 1 ( x , y , t ) g 1 ( x , y , t ) δ y c u ( x , y , t ) : = 1 1 σ h 3 q 2 ( x , y , t ) h 3 σ l = 1 J q l D M α l u ( x , y , t ) + p 2 ( x , y , t ) g 2 ( x , y , t )
and grid functions
δ x c u i , j m = δ x c u ( x i , y j , t m ) , δ y c u i , j m = δ y c u ( x i , y j , t m )
Lemma 2. 
Assume u ( · , · , t m ) C 2 ( Ω ¯ ) for fixed t m ; then, there exists a constant C such that
| δ x c u ( 0 , 0 , t m ) 2 x 2 u ( 0 , 0 , t m ) | + | δ y c u ( 0 , 0 , t m ) 2 y 2 u ( 0 , 0 , t m ) | C h 2 + C h m m i n { 2 α 1 , r ( α 1 + 1 ) } .
Proof. 
For corner point ( 0 , 0 , t m ) , where m = 1 , , M (1a) and (1c) at point ( 0 , 0 , t m ) are
l = 1 L q l D t α l u ( 0 , 0 , t m ) Δ u ( 0 , 0 , t m ) + c ( 0 , 0 ) u ( 0 , 0 , t m ) = F ( 0 , 0 , t m ) ,
σ u ( 0 , 0 , t m ) x u ( 0 , 0 , t m ) = g ( 0 , 0 , t m ) , σ u ( 0 , 0 , t m ) y u ( 0 , 0 , t m ) = g ( 0 , 0 , t m ) .
Similarly, by the Taylor expansion of u ( h , 0 , t m ) at point ( 0 , 0 , t m ) and u ( 0 , h , t m ) at point ( 0 , 0 , t m )
u ( h , 0 , t m ) = u ( 0 , 0 , t m ) + h x u ( 0 , 0 , t m ) + h 2 2 2 x 2 u ( 0 , 0 , t m ) + h 3 6 3 x 3 u ( 0 , 0 , t m ) + C h 4 , u ( 0 , h , t m ) = u ( 0 , 0 , t m ) + h y u ( 0 , 0 , t m ) + h 2 2 2 y 2 u ( 0 , 0 , t m ) + h 3 6 3 y 3 u ( 0 , 0 , t m ) + C h 4 .
Combining the above two equations and boundary conditions (18), we have
2 x 2 u ( 0 , 0 , t m ) + 2 y 2 u ( 0 , 0 , t m ) = q 1 ( 0 , 0 , t m ) + q 2 ( 0 , 0 , t m ) h 3 3 x 3 u ( 0 , 0 , t m ) h 3 3 y 3 u ( 0 , 0 , t m ) C h 2 .
Differentiating (17) with respect to x and y, respectively, we can express ( 3 / x 3 ) u + ( 3 / y 3 ) u at ( 0 , 0 , t m ) as
3 x 3 u ( 0 , 0 , t m ) + 3 y 3 u ( 0 , 0 , t m ) = 2 σ l = 1 J q l D t α l u ( 0 , 0 , t m ) + p 1 ( 0 , 0 , t m ) + p 2 ( 0 , 0 , t m ) σ 2 y 2 u ( 0 , 0 , t m ) σ 2 x 2 u ( 0 , 0 , t m ) g 1 ( 0 , 0 , t m ) g 2 ( 0 , 0 , t m ) ,
where we can apply (20) into the right-hand side of (19); thus, we have
2 x 2 u ( 0 , 0 , t m ) + 2 y 2 u ( 0 , 0 , t m ) = q 1 ( 0 , 0 , t m ) + q 2 ( 0 , 0 , t m ) h 3 ( 2 σ l = 1 J q l D t α l u ( 0 , 0 , t m ) σ 2 y 2 u ( 0 , 0 , t m ) σ 2 x 2 u ( 0 , 0 , t m ) + p 1 ( 0 , 0 , t m ) + p 2 ( 0 , 0 , t m ) g 1 ( 0 , 0 , t m ) g 2 ( 0 , 0 , t m ) ) C h 2 .
That means
2 x 2 u ( 0 , 0 , t m ) + 2 y 2 u ( 0 , 0 , t m ) = 1 1 σ h 3 [ q 1 ( 0 , 0 , t m ) + q 2 ( 0 , 0 , t m ) h 3 ( 2 σ l = 1 J q l D t α l u ( 0 , 0 , t m ) + p 1 ( 0 , 0 , t m ) + p 2 ( 0 , 0 , t m ) g 1 ( 0 , 0 , t m ) g 2 ( 0 , 0 , t m ) ) ] C h 2 .
Add 2 σ h / ( 3 σ h ) ( l = 1 J q l D M α l u ( 0 , y j , t m ) ) to the right-hand side and subtract it. Then, by truncation error (4), we have
2 x 2 u ( 0 , 0 , t m ) + 2 y 2 u ( 0 , 0 , t m ) = 1 1 σ h 3 [ q 1 ( 0 , 0 , t m ) + q 2 ( 0 , 0 , t m ) h 3 ( 2 σ l = 1 J q l D M α l u ( 0 , 0 , t m ) + p 1 ( 0 , 0 , t m ) + p 2 ( 0 , 0 , t m ) g 1 ( 0 , 0 , t m ) g 2 ( 0 , 0 , t m ) ) ] C h 2 C h m m i n { 2 α 1 , r ( α 1 + 1 ) } .
We can approximate (1) on corner point ( 0 , 0 , t m ) with the discrete problem
l = 1 L q l D M α l u 0 , 0 m δ x c u 0 , 0 m δ y c u 0 , 0 m + c 0 , 0 u 0 , 0 m = F 0 , 0 m for 1 j N 1 , 1 m M .
The other corner points can be treated similarly.

3. Error Analysis

The local error analysis of problem (1) is studied in this section. The discrete scheme is the same as in Section 2. Let
E m : = M r t m α 1 1 i f 1 r < 2 α 1 , M α 1 2 t m α 1 1 [ 1 + ln ( t m / t 1 ) ] i f r = 2 α 1 , M α 1 2 t m α 1 ( 2 α 1 ) / r i f r > 2 α 1 .
From ([20] Theorem 3), we obtain the next result, which will be used.
Lemma 3. 
(i) 
If the mesh function { V m } m = 0 M satisfies V 0 = 0 and
l = 1 L q l D M α l V m C m min { 2 α 1 , r ( α 1 + 1 ) } for m = 1 , 2 , , M ,
for some C > 0 , then | V m | C E m for m = 1 , 2 , , M .
(ii) 
If the mesh function { V m } m = 0 M satisfies V 0 = 0 and
l = 1 L q l D M α l V m C for m = 1 , 2 , , M ,
for some C > 0 , then | V m | C for m = 1 , 2 , , M .
Now, we provide the main result of this paper. For grid function { v i , j m } , set v m , Ω ¯ = max ( i , j ) Ω ¯ h | v i , j m | .
Theorem 1. 
The solution { u i , j m } of the L1-finite difference scheme satisfies
max ( x i , y j , t m ) Q ¯ | u ( x i , y j , t m ) u i , j m | C h 2 + E m
for some constant C independent of the mesh.
Proof. 
Set e i , j m : = u ( x i , y j , t m ) u i , j m , where ( i , j , m ) Q ¯ . Set ( i * , j * ) Ω ¯ h such that | e i * , j * m | = max ( i , j ) Ω ¯ | e i , j m | . Suppose that e i * , j * m 0 (the case e i * , j * m 0 can be proved similarly).
If ( i * , j * ) Ω h , by (1a) and (7) we obtain the error equation
l = 1 L q l D M α l e i , j m δ x 2 e i , j m δ y 2 e i , j m + c ( x i , y j ) e i , j m = ( l = 1 L q l D M α l l = 1 L q l D t α l ) u ( x i , y j , t m ) + ( δ x 2 2 x 2 ) u ( x i , y j , t m ) + ( δ y 2 2 y 2 ) u ( x i , y j , t m ) = R t i , j + R x i , j + R y i , j .
As a result of | e i * , j * m | = max ( i , j ) Ω ¯ | e i , j m | and e i * , j * m 0 , we have δ x 2 e i * , j * m δ y 2 e i * , j * m > 0 . Combine this with c c 0 > 0 , we have
l = 1 L q l D M α l e i * , j * m R t i , j + R x i , j + R y i , j C m min { 2 α 1 , r ( α 1 + 1 ) } + h 2 .
If i * = 0 and j * = 1 , , N 1 . Applying (1a) into (15) leads to the error equation
( 1 + σ h 3 ) l = 1 L q l D M α l e 0 , j m 2 h δ x e 1 , j m ( 1 + σ h 3 ) δ y 2 e 0 , j m + c 0 , j + h 3 6 σ h 2 + ( c 1 ) 0 , j + σ c 0 , j e 0 , j m = ( l = 1 L q l D M α l l = 1 L q l D t α l ) u ( 0 , y j , t m ) + ( δ x b 2 x 2 ) u ( 0 , y j , t m ) + ( δ y 2 2 y 2 ) u ( 0 , y j , t m ) = R t 0 , j + R x 0 , j + R y 0 , j .
It is easy to obtain 2 h δ x e 1 , j * m ( 1 + σ h 3 ) δ y 2 e 0 , j * m > 0 . When h is small enough, ( c 0 , j + h 3 ( 6 σ h 2 + ( c 1 ) 0 , j + σ c 0 , j ) ) > 0 . Compared to the time direction truncation error C m min { 2 α 1 , r ( α 1 + 1 ) } , we can omit the higher order truncation error C h m min { 2 α 1 , r ( α 1 + 1 ) } caused by boundary discretisation. Then, we have
l = 1 L q l D M α l e i * , j * m ( 1 + σ h 3 ) l = 1 L q l D M α l e i * , j * m R t 0 , j + R x 0 , j + R y 0 , j C m min { 2 α 1 , r ( α 1 + 1 ) } + h 2 .
If i * = j * = 0 . By (1a) into (23), the error equation is
( 1 + 2 h σ 3 σ h ) l = 1 L q l D M α l e 0 , 0 n 1 1 σ h 3 ( 2 h δ x e 1 , 0 m + 2 h δ y e 0 , 1 m ) + 1 1 σ h 3 h 3 ( c 1 ) 0 , 0 + ( c 2 ) 0 , 0 2 σ h 3 c ( x 0 , y 0 ) + 2 σ e 0 , 0 m + c 0 , 0 e 0 , 0 m = ( l = 1 L q l D M α l l = 1 L q l D t α l ) u ( 0 , 0 , t m ) + ( δ x c 2 x 2 ) u ( 0 , 0 , t m ) + ( δ y c 2 y 2 ) u ( 0 , 0 , t m ) = R t 0 , 0 + R x 0 , 0 + R y 0 , 0 .
Then, we have 2 h δ x e 1 , 0 m 2 h δ x e 0 , 1 m > 0 . When h is small enough, we have c ( 0 , 0 ) 1 1 σ h 3 ( c 1 ) 0 , 0 + ( c 2 ) 0 , 0 + 2 σ σ h 3 c ( x 0 , y 0 ) . Similarly, omitting higher order error caused by C h m min { 2 α 1 , r ( α 1 + 1 ) } , we have
l = 1 L q l D M α l e i * , j * m R t 0 , 0 + R x 0 , 0 + R y 0 , 0 C m min { 2 α 1 , r ( α 1 + 1 ) } + h 2 .
For other corner points, we shall obtain similar results.
For l = 1 , , L , rewrite the discretization of each Caputo derivative as
D M α l u i , j m = d m , 1 l Γ ( 2 α l ) u i , j m d m , m l Γ ( 2 α l ) u i , j 0 + 1 Γ ( 2 α l ) k = 1 m 1 u i , j m k d m , k + 1 l d m , k l
where
d m , k l : = ( t m t m k ) 1 α l ( t m t m k + 1 ) 1 α l τ m k + 1 .
The mean value theorem gives ( 1 α l ) ( t m t m k ) α l d m , k l ( 1 α l ) ( t m t m k + 1 ) α l and hence d m , k l d m , k + 1 l 0 . Then,
D M α l e i * , j * m = d m , m α l e i * , j * m k = 0 m 1 ( d m , k α l d m , k + 1 α l ) e i * , j * k . d m , m α l e m , Ω ¯ h k = 0 m 1 ( d m , k α l d m , k + 1 α l ) e k , Ω ¯ h = D M α l e m , Ω ¯ h .
By (28), (30) and (32), we have
l = 1 L q l D M α l e m , Ω ¯ h C m min { 2 α 1 , r ( α 1 + 1 ) } + h 2 .
Note that l = 1 J q l D M α l is associated with an M-matrix, so we can deal separately with the terms m min { 2 α 1 , r ( α 1 + 1 ) } and h 2 . This means for m = 1 , , M ,
u ( x i , y j , t m ) u i , j m , Ω ¯ h C h 2 + E m .
Then, we finish the proof. □
Remark 1. 
From (26), we can obtain the following global error results
max m = 1 , M u ( x i , y j , t m ) u i , j m , Ω ¯ h C max m = 1 , M h 2 + E m C h 2 + M min { 2 α 1 , r α 1 } .

4. Numerical Results

In order to prove the validity of the numerical scheme, two numerical examples are introduced. One example has a known solution, the other is unknown.
We use the full discrete scheme in Section 2 to discretize 1. In the following examples, we set mesh parameters r = ( 2 α 1 ) / 0.9 , L = 2 and 0 < α 2 < α 1 < 1 . Let the space interval N equals to the time interval M such that the error in the time direction dominates the space error. On this basis, we shall check the sharpness of Theorem (1).
Example 1. 
D t α 1 u + D t α 2 u 2 u x 2 2 u y 2 + ( 1 + x + y ) u = f ( x , y , t ) for ( x , y , t ) [ 0 , 2 ] × [ 0 , 2 ] × ( 0 , 1 ] , u ( x , y , 0 ) = ( 1 3 x 3 x 2 + 1 3 x + 1 3 ) ( 1 3 y 3 y 2 + 1 3 y + 1 3 ) for   ( x , y ) [ 0 , 2 ] × [ 0 , 2 ] , u ( 0 , y , t ) u x ( 0 , y , t ) = 0 , u ( 2 , y , t ) + u x ( 2 , y , t ) = 0 for   y [ 0 , 2 ] t ( 0 , 1 ] . u ( x , 0 , t ) u y ( x , 0 , t ) = 0 , u ( x , 2 , t ) + u y ( x , 2 , t ) = 0 for   x [ 0 , 2 ] t ( 0 , 1 ] .
The exact solution is u ( x , y , t ) = ( 1 + t α 1 + t 3 ) ( 1 3 x 3 x 2 + 1 3 x + 1 3 ) ( 1 3 y 3 y 2 + 1 3 y + 1 3 ) . The right-hand-side function f ( x , y , t ) can be computed from (35).
In Table 1, the table contains the global error, and local error is defined as
e r r o r G M , N : = max i , j Ω ¯ 1 m M | u i , j m u ( x i , y j , t m ) | , e r r o r L M , N : = max i , j Ω ¯ M / 10 m N | u i , j m u ( x i , y j , t m ) |
Then, we can compute the rate of convergence
r a t e G M , N = log 2 e r r o r G M , N e r r o r G 2 M , 2 N , r a t e L M , N = log 2 e r r o r L M , N e r r o r L 2 M , 2 N .
Example 2. 
D t α 1 u + D t α 2 u 2 u x 2 2 u y 2 + ( 1 + x + y ) u = 0 for ( x , y , t ) [ 0 , 2 ] × [ 0 , 2 ] × ( 0 , 1 ] , u ( x , y , 0 ) = ( 1 3 x 3 x 2 + 1 3 x + 1 3 ) ( 1 3 y 3 y 2 + 1 3 y + 1 3 ) for ( x , y ) [ 0 , 2 ] × [ 0 , 2 ] , u ( 0 , y , t ) u x ( 0 , y , t ) = 0 , u ( 2 , y , t ) + u x ( 2 , y , t ) = 0 for y [ 0 , 2 ] t ( 0 , 1 ] . u ( x , 0 , t ) u y ( x , 0 , t ) = 0 , u ( x , 2 , t ) + u y ( x , 2 , t ) = 0 for x [ 0 , 2 ] t ( 0 , 1 ] .
In this example, the exact solution is unknown, and we can use the two-mesh principle in [25] to check the convergence rate. Let u i , j m be the numerical solution computed on the mesh { ( x i , y j , t m ) } for i , j = 0 , , N , m = 0 , , M . The second mesh is defined as { ( x i / 2 , y j / 2 , t m / 2 ) } for i , j = 0 , , 2 N , m = 0 , , 2 M , where x i + 1 / 2 : = 1 2 ( x i + 1 + x i ) , y y + 1 / 2 : = 1 2 ( y j + 1 + y j ) and t m + 1 / 2 : = 1 2 ( t m + 1 + t m ) . Then, compute a new approximation u ^ i , j m using the same scheme as u i , j m .
Now the maximum two-mesh differences are defined by
e r r o r G M , N : = max i , j Ω ¯ 1 m M | u i , j m u ^ i , j m | , e r r o r L M , N : = max i , j Ω ¯ 9 M / 10 m M | u i , j m u ^ i , j m |
and they are used to compute the global and local rate of convergence rates
r a t e G M , N = log 2 e r r o r G M , N e r r o r G 2 M , 2 N , r a t e L M , N = log 2 e r r o r L M , N e r r o r L 2 M , 2 N .
In each Table 1 and Table 2, let r = ( 2 α 1 ) / 0.9 and α 2 = 0.1 . The global convergence rate is bigger than α 1 . The convergence rate e i , j m M m i n { 2 α 1 , r α 1 } can be found in other papers that only focus on global errors [26,27]. When the parameters r and a 1 are selected the same as in this paper, the convergence rate can be seen to be the same. The most important conclusion of this paper is the convergence rate of local errors. We can see that the rate of convergence is ( 2 α 1 ) . It is obvious that the local convergence rate in every time step is sharper than the global convergence rate. All these experimental results demonstrate the sharpness of our theoretical analysis.

5. Discussion

In this paper, we have presented a fully discrete scheme for multiple time-fractional reaction diffusion equations by using the L1 scheme in time and finite difference method in space. To the best of our knowledge, the Robin boundary conditions have not been explored much in this regard. For this type of boundary conditions, we have constructed a discrete scheme of (1) at the boundary points which can match the convergence rate of the inner points. Based on the fully discrete scheme, a detailed local error analysis for (1) is presented. The convergence rate of each time node is proven, and two numerical examples are used to verify the theoretical results. In our future work, we will consider some methods with higher convergence rates in time and consider methods such as the mixed finite element method in space.

Author Contributions

Conceptualization and writing—original draft, J.H.; writing—review and editing, J.H., X.M., J.W. and Y.H.; supervision, Y.Y. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China grant number 62173027 and 12101039.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Derakhshan, M.H. A numerical technique for solving variable order time fractional differential-integro equations. Commun. Math. 2024, 32, 129–147. [Google Scholar] [CrossRef]
  2. Zhou, Y.; Wang, J.; Zhang, L. Basic Theory of Fractional Differential Equations; World Scientific Publishing Co. Pte. Ltd.: Hackensack, NJ, USA, 2017. [Google Scholar]
  3. Hilfer, R. (Ed.) Applications of Fractional Calculus in Physics; World Scientific Publishing Co., Inc.: River Edge, NJ, USA, 2000; p. viii+463. [Google Scholar] [CrossRef]
  4. Karthikeyan, K.; Karthikeyan, P.; Chalishajar, D.N.; Raja, D.S.; Sundararajan, P. Analysis on ψ-Hilfer Fractional Impulsive Differential Equations. Symmetry 2021, 13, 1895. [Google Scholar] [CrossRef]
  5. Agarwal, R.P.; dos Santos, J.P.C.; Cuevas, C. Analytic resolvent operator and existence results for fractional integro-differential equations. J. Abstr. Differ. Equ. Appl. 2012, 2, 26–47. [Google Scholar]
  6. Stynes, M.; O’Riordan, 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]
  7. Povstenko, Y. Fundamental solutions to the fractional heat conduction equation in a ball under Robin boundary condition. Cent. Eur. J. Math. 2014, 12, 611–622. [Google Scholar] [CrossRef]
  8. Brociek, R.; Słota, D. Reconstruction of the Robin boundary condition and order of derivative in time fractional heat conduction equation. Math. Model. Nat. Phenom. 2018, 13, 5. [Google Scholar] [CrossRef]
  9. Jílek, M. Straight quantum waveguide with Robin boundary conditions. SIGMA Symmetry Integr. Geom. Methods Appl. 2007, 3, 12. [Google Scholar] [CrossRef]
  10. Huang, C.; Stynes, M. A direct discontinuous Galerkin method for a time-fractional diffusion equation with a Robin boundary condition. Appl. Numer. Math. 2019, 135, 15–29. [Google Scholar] [CrossRef]
  11. Kopteva, N. Error analysis of the L1 method on graded and uniform meshes for a fractional-derivative problem in two and three dimensions. Math. Comput. 2019, 88, 2135–2155. [Google Scholar] [CrossRef]
  12. Sun, Z.; Wu, X. A fully discrete difference scheme for a diffusion-wave system. Appl. Numer. Math. 2006, 56, 193–209. [Google Scholar] [CrossRef]
  13. Li, B.; Wang, T.; Xie, X. Analysis of the L1 scheme for fractional wave equations with nonsmooth data. Comput. Math. Appl. 2021, 90, 1–12. [Google Scholar] [CrossRef]
  14. Chaudhary, S.; Kundaliya, P.J. L1 scheme on graded mesh for subdiffusion equation with nonlocal diffusion term. Math. Comput. Simul. 2022, 195, 119–137. [Google Scholar] [CrossRef]
  15. Chen, H.; Chen, M.; Sun, T.; Tang, Y. Local error estimate of L1 scheme for linearized time fractional KdV equation with weakly singular solutions. Appl. Numer. Math. 2022, 179, 183–190. [Google Scholar] [CrossRef]
  16. Chen, H.; Stynes, M. Using complete monotonicity to deduce local error estimates for discretisations of a multi-term time-fractional diffusion equation. Comput. Methods Appl. Math. 2022, 22, 15–29. [Google Scholar] [CrossRef]
  17. Akil, M.; Issa, I.; Wehbe, A. Energy decay of some boundary coupled systems involving wave Euler-Bernoulli beam with one locally singular fractional Kelvin-Voigt damping. Math. Control Relat. Fields 2023, 13, 330–381. [Google Scholar] [CrossRef]
  18. Vijayakumar, V.; Nisar, K.S.; Chalishajar, D.; Shukla, A.; Malik, M.; Alsaadi, A.; Aldosary, S.F. A Note on Approximate Controllability of Fractional Semilinear Integrodifferential Control Systems via Resolvent Operators. Fractal Fract. 2022, 6, 73. [Google Scholar] [CrossRef]
  19. Hakkar, N.; Dhayal, R.; Debbouche, A.; Torres, D.F.M. Approximate Controllability of Delayed Fractional Stochastic Differential Systems with Mixed Noise and Impulsive Effects. Fractal Fract. 2023, 7, 104. [Google Scholar] [CrossRef]
  20. Meng, X.; Stynes, M. Barrier function local and global analysis of an L1 finite element method for a multiterm time-fractional initial-boundary value problem. J. Sci. Comput. 2020, 84, 16. [Google Scholar] [CrossRef]
  21. Oloniiju, S.D.; Goqo, S.P.; Sibanda, P. A Chebyshev pseudo-spectral method for the multi-dimensional fractional Rayleigh problem for a generalized Maxwell fluid with Robin boundary conditions. Appl. Numer. Math. 2020, 152, 253–266. [Google Scholar] [CrossRef]
  22. Wang, Y.; Ren, L. Analysis of a high-order compact finite difference method for Robin problems of time-fractional sub-diffusion equations with variable coefficients. Appl. Numer. Math. 2020, 156, 467–492. [Google Scholar] [CrossRef]
  23. Huang, C.; Liu, X.; Meng, X.; Stynes, M. Error analysis of a finite difference method on graded meshes for a multiterm time-fractional initial-boundary value problem. Comput. Methods Appl. Math. 2020, 20, 815–825. [Google Scholar] [CrossRef]
  24. Hemati, F.; Ghasemi, M.; Khoshsiar Ghaziani, R. Numerical solution of the multiterm time-fractional diffusion equation based on reproducing kernel theory. Numer. Methods Partial Differ. Equ. 2021, 37, 44–68. [Google Scholar] [CrossRef]
  25. Farrell, P.A.; Hegarty, A.F.; Miller, J.J.H.; O’Riordan, E.; Shishkin, G.I. Robust Computational Techniques for Boundary Layers; Applied Mathematics (Boca Raton); Chapman & Hall/CRC: Boca Raton, FL, USA, 2000; Volume 16. [Google Scholar]
  26. Ren, J.; Chen, H. A numerical method for distributed order time fractional diffusion equation with weakly singular solutions. Appl. Math. Lett. 2019, 96, 159–165. [Google Scholar] [CrossRef]
  27. Zhang, J.; Chen, H.; Sun, T.; Wang, J. Error analysis of a fully discrete scheme for time fractional Schrödinger equation with initial singularity. Int. J. Comput. Math. 2020, 97, 1636–1647. [Google Scholar] [CrossRef]
Table 1. Example 1 with α 2 = 0.1 and r = ( 2 α 1 ) / α 1 . .
Table 1. Example 1 with α 2 = 0.1 and r = ( 2 α 1 ) / α 1 . .
Global ErrorRateLocal ErrorRate
α 1 = 0.4 M = 649.1820 × 10 4 0.5093.4584 × 10 5 1.551
M = 1286.4522 × 10 4 0.5671.1802 × 10 5 1.568
M = 2564.3544 × 10 4 0.6123.9795 × 10 6 1.579
M = 5122.8472 × 10 4 1.3317 × 10 6
α 1 = 0.6 M = 644.6098 × 10 4 0.8068.7567 × 10 5 1.377
M = 1282.6355 × 10 4 0.8603.3700 × 10 5 1.387
M = 2561.4514 × 10 4 0.8871.2876 × 10 5 1.393
M = 5127.8446 × 10 5 4.9017 × 10 6
α 1 = 0.8 M = 642.2223 × 10 4 0.8522.1862 × 10 4 1.190
M = 1281.2316 × 10 4 0.9039.5808 × 10 5 1.195
M = 2566.5823 × 10 5 0.9434.1826 × 10 5 1.198
M = 5123.4215 × 10 5 1.8230 × 10 5
Table 2. Example 2 with α 2 = 0.1 and r = ( 2 α 1 ) / 0.9 . .
Table 2. Example 2 with α 2 = 0.1 and r = ( 2 α 1 ) / 0.9 . .
Global ErrorRateLocal ErrorRate
α 1 = 0.4 M = 326.6148 × 10 3 0.2201.1175 × 10 4 1.683
M = 645.6762 × 10 3 0.3743.4786 × 10 5 1.637
M = 1284.4626 × 10 3 0.4531.1178 × 10 5 1.568
M = 2563.2585 × 10 3 3.7696 × 10 6
α 1 = 0.6 M = 325.3413 × 10 3 0.6624.0714 × 10 4 1.248
M = 643.3751 × 10 3 0.7621.7139 × 10 4 1.391
M = 1281.9889 × 10 3 0.7516.5340 × 10 5 1.310
M = 2561.1813 × 10 3 2.6352 × 10 5
α 1 = 0.8 M = 323.7397 × 10 3 0.8021.6129 × 10 3 1.175
M = 642.1438 × 10 3 0.8777.1399 × 10 4 1.073
M = 1281.1668 × 10 3 0.9423.3927 × 10 4 1.163
M = 2566.0732 × 10 4 1.5143 × 10 4
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

Hou, J.; Meng, X.; Wang, J.; Han, Y.; Yu, Y. Local Error Estimate of an L1-Finite Difference Scheme for the Multiterm Two-Dimensional Time-Fractional Reaction–Diffusion Equation with Robin Boundary Conditions. Fractal Fract. 2023, 7, 453. https://0-doi-org.brum.beds.ac.uk/10.3390/fractalfract7060453

AMA Style

Hou J, Meng X, Wang J, Han Y, Yu Y. Local Error Estimate of an L1-Finite Difference Scheme for the Multiterm Two-Dimensional Time-Fractional Reaction–Diffusion Equation with Robin Boundary Conditions. Fractal and Fractional. 2023; 7(6):453. https://0-doi-org.brum.beds.ac.uk/10.3390/fractalfract7060453

Chicago/Turabian Style

Hou, Jian, Xiangyun Meng, Jingjia Wang, Yongsheng Han, and Yongguang Yu. 2023. "Local Error Estimate of an L1-Finite Difference Scheme for the Multiterm Two-Dimensional Time-Fractional Reaction–Diffusion Equation with Robin Boundary Conditions" Fractal and Fractional 7, no. 6: 453. https://0-doi-org.brum.beds.ac.uk/10.3390/fractalfract7060453

Article Metrics

Back to TopTop