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

: 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-ﬁnite difference method, which entails using the L1 scheme on graded meshes for the temporal discretisation of each Caputo fractional derivative and using the ﬁnite 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 conﬁrm the sharpness of the error estimates are presented.


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 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.
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 α l t u(x, y, t) is defined as the temporal Caputo fractional derivative of order α l of u by For ([10] Lemma 2.2) and ( [11] Section 6), (1) has unique solution which satisfies the following regularity with weak initial singularity 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 timefractional 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. (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.
At each mesh point (x i , y j , t m ), the computed approximation to the analytical solution u will be denoted by u m i,j . Define the grid functions where (i, j) ∈Ω h , m = 1, 2, . . . , M.
The Caputo fractional derivative D α t u can be expressed as The L1 scheme, which is used to approximate the Caputo fractional derivative to obtain the discretisation of each time-fractional term q l D α l t u(x i , y j , t m ).
In ( [20] Lemma 4), the truncation error has the following estimate We discrete the initial condition (1b) in a standard way: for i, j ∈Ω h , set u 0 i,j = 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.

Inner Points
Then, the truncation error has the following estimate In summary, we can approximate (1)

Boundary Points
For brevity, we set Then, define grid functions Lemma 1. Assume u(·, ·, t m ) ∈ C 2 (Ω) for every t m ; then, there exists a constant C such that By using (10), we have Differentiating (9) with respect to x, (∂ 3 /∂x 3 )u can be expressed as Furthermore, in view of (10), we obtain Adding (σh/3)(∑ J l=1 q l D α l M u(0, y j , t m )) to the right-hand side and subtract it. Then, by truncation error (4), we have The other corner points can be treated similarly.

Corner Points
For convenience, we introduce the following functions and grid functions Lemma 2. Assume u(·, ·, t m ) ∈ C 2 (Ω) for fixed t m ; then, there exists a constant C such that Proof. For corner point (0, 0, t m ), where m = 1, . . . , M (1a) and (1c) at point (0, 0, t m ) are 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 ) Combining the above two equations and boundary conditions (18), we have 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 where we can apply (20) into the right-hand side of (19); thus, we have That means Add (2σh/(3 − σh))(∑ J l=1 q l D α l M u(0, y j , t m )) to the right-hand side and subtract it. Then, by truncation error (4), we have We can approximate (1) on corner point (0, 0, t m ) with the discrete problem The other corner points can be treated similarly.

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 From ( [20] Theorem 3), we obtain the next result, which will be used.
for some constant C independent of the mesh.

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).
The exact solution is u(x, y, t) = ( . The right-hand-side function f (x, y, t) can be computed from (35).
In Table 1, the table contains the 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]. 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 m i,j 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, . . . , 2N, m = 0, . . . , 2M, 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 m i,j using the same scheme as u m i,j . Now the maximum two-mesh differences are defined by In each Tables 1 and 2, let r = (2 − α 1 )/0.9 and α 2 = 0.1. The global convergence rate is bigger than α 1 . The convergence rate e m i,j ≤ M −min{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.

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.