In this paper, we consider the following two-dimensional nonlinear fractional reaction- subdiffusion equation
with boundary and initial conditions
where, , is sufficiently smooth function. For simplicity, we assume coefficients, and are positive constants in this paper. In fact, our method and its corresponding theoretic result are also valid for variable coefficients. is the Riemann-Liouville time fractional derivative of order defined by 
where denotes the Riemann-Liouville fractional integral operator defined as 
In addition, we assume that the nonlinear source term satisfies the Lipschitz condition with respect to, i.e., there exists a positive constant such that
Problem (1) can be considered as a model for reaction-diffusion phenomena with anomalous diffusion, which has been widely applied in various fields of science and engineering. Generally, solutions of (1) can’t be obtained by analytical approach. So, there are various numerical methods developed for solving (1). Li and Ding  proposed higher order finite difference methods for solving 1D linear reaction and anomalous- diffusion equations. Zhuang, Liu and Anh, et al.  presented an implicit finite element method for solving 1D nonlinear fractional reaction-subdiffusion process. Dehghan, Abbaszadeh and Mohebbi  analyzed a meshless Galerkin method with radial basis functions of 2D linear fractional reaction-subdiffusion process. Yu, Jiang and Xu  derived an implicit compact finite difference scheme for solving 2D nonlinear fractional reaction-subdiffusion equation.
Alternating direction implicit (ADI) method was proposed by Peaceman, Rachford and Douglas    in 1950’s for multidimensional differential equations of integer order, which could reduce original multidimensional problem into a sequences of one- dimensional problems. Since the first ADI based finite difference (FD) scheme presented for 2D space fractional diffusion equation by Meerschaert, Scheffler and Tadjeran  , there are many literatures about various multidimensional fractional differential equations numerically solved by ADI technique. The following problem is always discussed:
where is Caputo fractional derivative of order defined as 
In the case of, Zhang and Sun  presented an ADI FE scheme and analyzed its stability and convergence property by energy method. Cui  constructed a compact ADI FD scheme and discussed its stability by Fourier method. In the case of, Zhang, Sun and Zhao  proposed a Crank-Nicolson compact ADI FD scheme, where stability and convergence property are also proved by energy method. Wang and Vong  presented another compact ADI FD scheme, where the original equation was first transformed into an equivalent form and then discretized by compact FD scheme combining with ADI technique. Li, Xu and Luo  presented an ADI finite element (FE) scheme for Equation (5), where stability and L2 error estimate are analyzed. ADI orthogonal spline collocation (OSC) scheme was developed by Fairweather, Yang and Xu, et al.  for solving Equation (5); stability and error estimate in various norms are given. The equivalent form of Equation (5) as following
is also often studied. For instances, Cui  presented a compact ADI FD scheme, where the local truncation error was analyzed and the stability was discussed by the Fourier method. Furthermore, the author analyzed the cause of low time accuracy when and gave a remedy. Zhang and Sun  proposed a Crank-Nicol- son compact ADI FD scheme for Equation (6), where stability and two error estimates are proved rigorously by energy method. Yao, Sun and Wu, et al.  first transformed Equation (6) into an equivalent form of Caputo fractional derivative and then derived a compact ADI FD scheme. Their numerical experiments show better numerical performance than the scheme in  . Compact ADI FD schemes were also derived for solving 2D/3D linear time fractional convection-diffusion equations    and 2D linear time fractional diffusion equations of distributed-order   .
There are also lots of ADI based numerical methods for multidimensional space fractional differential equations. Fast iterative ADI FD schemes   are designed for 2D/3D linear space fractional diffusion equations, which are first order accuracy in both time and space and have the advantage of low computational work and low memory storage. High order accurate ADI FD schemes are proposed for 2D linear space fractional diffusion equations   and two-sided space fractional convection-dif- fusion equations  , which are based on weighted and shifted Grünwald operators or Lubich operators approximating Riemann-Liouville fractional derivatives respectively. Spectral direction splitting methods  are derived for 2D space fractional differential equations. Semi-implicit alternating direction FD scheme  and ADI FE scheme  are used for solving 2D fractional Fitz Hugh-Nagumo monodomain model, which consists of a coupled 2D space fractional nonlinear reaction-diffusion equation and an ordinary differential equation, on irregular domain and rectangle domain respectively. ADI Galerkin-Legendre spectral method  is developed for 2D Riesz space fractional nonlinear reaction-diffusion equation.
Most of the above mentioned works contribute on linear fractional differential equations and finite difference method combined with ADI technique. A few work consider ADI FEM   or nonlinear fractional differential equations    . Compared with FD method, FE method has the advantage of easily handling variable coefficients problem and boundary conditions. And many realistic problems involve nonlinear fractional differential equations. Based on these motivations, our attention in this paper will focus on developing ADI FE schemes for efficiently solving a class of nonlinear time fractional differential equations. This is the first time ADI FE scheme proposed and analyzed rigorously for nonlinear fractional differential equations. We will use problem (1.1) as a model problem to illustrate our approach.
The outline of this paper is organized as follows. In Section 2, we introduce some preliminaries and notations which will be used later. The formulation of ADI finite element method for nonlinear time fractional reaction-subdiffusion equation is presented in Section 3. The stability and error estimates of the proposed method are discussed in Section 4. In convenience of computation, we give the matrix form of ADI finite element scheme in Section 5. Some numerical experiments are displayed in Section 6. It aims to confirm our theoretical results. In the end, some concluding remarks are given in Section 5.
In the following, and denote generic positive constants independent of, and, and their value will not be the same in different equations or inequalities.
2. Preliminary and Notations
Let be a bounded and open domain in. Denote the inner product and norm on by
Recall that the Sobolev space is the closure of in the norm
Denote by the closure of in the norm; it is well known that an equivalent norm on is
If is a normed space with norm, recall that
Denote. Let be the finite dimensional subspace of and satisfy
In convenience, the following notations will be used. For a positvie integer, let, where is the time step. Denote the grid function
Define linear operator by
and its variational form
and corresponding energy norm by
Obviously, we have
Some useful lemmas are given as follows.
Lemma 1.  Let
then satisfy the following properties
Lemma 2.  If, then
where, is bounded by
We state here for convenience the discrete version of Gronwall’s inequality.
Lemma 3.  Suppose that and are nonnegative functions defined for, and that is nondecreasing. If
where is a positive constant, then
3. Formulation of ADI FEM
Integrating both sides of (1) with respect to the time variable from to, and noticing the definition (3) and (7), we can obtain
Applying Lemma 2 and the following integration formula
Equation (8) is equivalent to
where the remainder term can be bounded by
The weak form of Equation (9) is: find such that for any,
Then the finite element approximation to Equation (11) is: find such that for any,
The choice of the initial values will be discussed later.
Let, then, (12) can be transformed into
so that the alternating-direction Galerkin scheme of (1) can be defined as, for ,
Remark 1. Numerical experiments in Section 6 demonstrate that the ADI Galerkin finite element scheme (14) has bad numerical performance for. Similar phenomena have been reported in   and  , where compact ADI finite difference schemes are designed respectively for solving 2D linear time fractional sub-diffu- sion equations. To construct ADI finite element method, the third term
in left hand side of (14) is extra added. Its effect on temporal accuracy cannot be ignored when. To balance it, we can add a correction term
on the right hand side of (14). By the way, the similar remedy was also adopted by   in finite difference framework.
In conclusion, when, it is beneficial to use the following scheme for
is defined by (15). And for the approximation is still given by (14).
In the following, we will focus our attention on the ADI Galerkin finite element scheme (14).
4. Stability and Error Estimate
Firstly, we introduce some notation and lemmas. Given a smooth function, define by
The operator defined in (18) has the following approximate properties:
Lemma 4.  Let for, , there exists a constant, independent of meshsize, such that
Lemma 5.  If denote the operators and, then
Lemma 6. Let and are defined in (7) and Lemma 1, respectively. For any positive integer and, let
then it holds that
Proof. When, it is easily get. Now we consider. Note that
and, we get
By direct algebraic calculation, and noticing the properties of in Lemma 1, we obtain
The proof is completed.
Next, we consider the stability of the ADI Galerkin method (14). Define the following problem dependent norm for any,
which will be used in stability analysis.
Assume the initial value has an error, which will cause the solution of ADI Galerkin scheme (14), , has a perturbation. Then the stability property of ADI Galerkin methods can be represented as following:
Theorem 1. The ADI Galerkin method (14) is stable with respect to initial value in the problem dependent norm (22), i.e., there exist a positive constant which independent of such that
holds for any.
Proof. The equivalent form of (14) is:
Then the perturbation equation of (23) can be written as
Taking in (24), employing Schwartz inequality and the notation, yields
Further, by Lipschitz property of and Schwartz inequality, we can estimate the right hand side of (25) as
Summing for in (25) and using (26) and Lemma 6, we can write
For sufficiently small, it follows from (27) that
Using discrete Gronwall inequality, we have
Let, we conclude. The proof is completed.
Theorem 2. Let and denote the solutions of (1) and (14) respectively.
Assume that, where. Then, for sufficiently small,
provided the initial value satisfy
Proof. Denote, then. With as in (18) we have from (11) that
Subtracting (28) from (23) leads to
Taking, by Schwartz inequality and using the notation, then (29) reduces to
Now, by Young inequality and Schwartz inequality, we can write
Substitute the above three inequalities into the right hand side of (30), and sum for. Employing Lemma 6, we can write
By discrete Gronwall inequality, if sufficiently small, we can estimate
Consequently, we obtain
Using the triangle inequality and (31), we get
It remains to estimate terms on the right-hand side of (32). Firstly, by Lemma 3, we can conclude
By (10) and, we easily get
Secondly, using calculus equality
and Hölder inequality, we have
Further, combined with Lemma 3, we can write
similar as (36), we can prove
Additionally, according to Lemma 3 and Lemma 4, we have
As a result, we obtain
Similar as (37), we can write
Combining (33)-(38), and recalling gives
provided is chosen so that and are. This completes the proof.
Remark 2. Although in our theoretic analysis, we only obtain accuracy in temporal error, numerical experiments in Section 6 demonstrate that temporal error is if. Checking the above analysis process, we will find the error estimates in (37) and (38) are obstacles to obtain temporal error theoretically. This suggests we may take alternative analysis techniques to improve the error estimates in (37) and (38).
Note that we may choose the initial approximation as, where from (18),
This involves an elliptic problem to be solved. With this choice,
In practical computations, it is often sufficient to take as interpolants of in.
5. Matrix Form of ADI FEM
Equations (14) define the ADI finite element method in inner product form. To describe the algebraic problem to which these equations lead, suppose,
where and are finite-dimensional subspaces of, and let
be a tensor product basis for, where and are bases for the subspaces and, respectively. Let
For convenience, we denote
If in (14) we choose, and then
We define the matrices
with defined similarly. Then the matrix form of (39) is
where, from (14), the components of are given by
and denotes the tensor product. We may factor (40) in the form
which is equivalent to
where and denote the identity matrices of order and, respectively. Thus we determine by solving two sets of independent one-dimensional problems, first
in the -direction, where, followed by
in the -direction, where. Clearly the computation of each of the vectors and is highly parallel.
6. Numerical Experiments
In this section, two numerical examples are given to demonstrate the effectiveness and accuracy of the ADI Galerkin finite element methods. In all numerical examples, we take the linear tensor product basis
where. is constructed in a similar way. In all examples, we will take the same spatial step size in each direction, i.e..
In our numerical simulation, we present the errors in norm
and numerical convergence orders are computed by
Example 1. Consider the following problem
with initial and boundary conditions
where. The exact solution of this equation is
Table 1 shows the L2 norm errors and the temporal convergence orders for = 0.6, 0.7, 0.8, 0.9 with fixed using ADI Galerkin finite element scheme (14). Computational results in Table 1 illustrate that our ADI Galerkin finite element scheme has first order accuracy in time, which is better than our predicted order in Theorem 2.
To investigate the necessary of the correction term (16) in ADI scheme (17) if
, we have compared the two ADI Galerkin finite element scheme (14) and (17)
for = 0.1, 0.2, 0.3, 0.4, with fixed. For saving space, we only present the comparison results of in Table 2, the others cases are similar. The first three columns of Table 2 show the L2 norm errors and its corresponding temporal convergence orders using ADI scheme (14). Computational results demonstrate the temporal
Table 1. L2 errors and temporal convergence orders, fixing h = π/64.
Table 2. Comparison between ADI scheme (14) and (17) for α = 0.1, h = π/64.
convergence order is 0.1, which also suggest our theoretical order proved in Theorem 2 is optimal. The last three columns of Table 2 display the L2 norm errors and its corresponding temporal convergence orders using ADI scheme (17). As we can see, the experiment convergence order is approximately one now. It indicates that the correction term (16) in ADI scheme (17) is beneficial to keep the temporal accuracy. In a word, Table 2 shows ADI scheme (17) has better numerical performance than ADI
scheme (14) if.
Table 3 shows the L2 norm errors and the spatial convergence orders for = 0.6, 0.7, 0.8, 0.9 with fixed using ADI Galerkin finite element scheme (14). The numerical results in Table 3 suggest our ADI Galerkin finite element scheme gets second order accuracy in space, which is consistent with our theoretical result in Theorem 2.
Example 2. Consider the following problem
with initial and boundary conditions
where,. The exact solution of this equation is
The nonlinearity of Example 2 is stronger than Example 1. In the following simulation, we have to take much smaller temporal step than Example 1 to obtain good numerical performance. Table 4 shows the L2 norm errors and the temporal convergence orders for = 0.6, 0.7, 0.8, 0.9 with fixed using ADI Galerkin finite element scheme (14). Computational results in Table 4 illustrate that our ADI Galerkin finite element scheme is first order temporal accuracy again, if.
Table 5 displays the comparison between the two ADI Galerkin finite element scheme (14) and (17) for, with fixed. From Table 5, we can conclude that: 1) The first three columns show the order of temporal accuracy of ADI scheme (14) is 0.2 for, confirming our theoretical result in Theorem 2 again; 2) The last three columns illustrate the ADI scheme (17) has first order temporal accuracy, indicating the essentiality of the correction term (16) to keep temporal accuracy once again; 3) To obtain the same level of temporal truncation error, ADI scheme (17) can take much smaller time step than ADI scheme (14) due to its higher accuracy. It also suggests that high order of temporal accuracy is important to time-dependent nonlinear problem.
Table 3. L2 errors and spatial convergence orders, fixing τ = 1/5000.
Table 4. L2 errors and temporal convergence orders, fixing h = π/64.
Table 5. Comparison between ADI scheme (14) and (17) for α = 0.2, h = π/128.
Table 6. L2 errors and spatial convergence orders, fixing τ = 1/5000.
Table 6 shows the L2 norm errors and the spatial convergence orders for = 0.6, 0.7, 0.8, 0.9 with fixed using ADI Galerkin finite element scheme (14). It can be shown from Table 6 that the numerical results are in accordance with the theoretical analysis in Theorem 2.
In this work, we proposed an alternating direction Galerkin finite element method for 2D nonlinear time fractional reaction sub-diffusion equation in the Riemann-Liouville type. The stability and convergence of the method are proved for. Numerical results show that our error estimate is optimal. Like the other ADI schemes   , the temporal accuracy of our ADI Galerkin finite element method becomes very low if. A remedy is introduced by adding a correction term (15) to keep temporal accuracy.
We thank the Editor and the referee for their comments. Research of P. Zhu is funded by the Natural Science Foundation of Zhejiang province, China (Grant No. LY15A010018). This support is greatly appreciated.
 Li, C. and Ding, H. (2014) Higher Order Finite Difference Method for the Reaction and Anomalous-Diffusion Equation. Applied Mathematical Modelling, 38, 3802-3821.
 Zhuang, P., Liu, F., Anh, V., et al. (2009) Stability and Convergence of an Implicit Numerical Method for the Nonlinear Fractional Reaction-Subdiffusion Process. SIAM Journal on Mathematical Analysis, 74, 645-667.
 Dehghan, M., Abbaszadeh, M. and Mohebbi, A. (2015) Error Estimate for the Numerical Solution of Fractional Reaction-Subdiffusion Process Based on a Meshless Method. Journal of Computational and Applied Mathematics, 280, 14-36.
 Yu, B., Jiang, X. and Xu, H. (2015) A Novel Compact Numerical Method for Solving the Two-Dimensional Non-Linear Fractional Reaction-Subdiffusion Equation. Numerical Algorithms, 68, 923-950.
 Douglas, J. and Gunn, J.E. (1964) A General Formulation of Alternating Direction Method, Part I. Parabolic and Hyperbolic Problem. Numerische Mathematik, 6, 428-453.
 Peaceman, D.W. and Rachford, H.H. (1955) The Numerical Solutions of Parabolic and Elliptic Differential Equations. Journal of the Society for Industrial and Applied Mathematics, 3, 28-41.
 Meerschaert, M.M., Scheffler, H.P. and Tadjeran, C. (2006) Finite Difference Methods for Two-Dimensional Fractional Dispersion Equation. Journal of Computational Physics, 211, 249-261.
 Zhang, Y. and Sun, Z. (2011) Alternating Direction Implicit Schemes for the Two-Dimensional Fractional Sub-Diffusion Equation. Journal of Computational Physics, 230, 8713-8728.
 Cui, M. (2013) Convergence Analysis of High-Order Compact Alternating Direction Implicit Schemes for the Two-Dimensional Time Fractional Diffusion Equation. Numerical Algorithms, 62, 383-409.
 Zhang, Y., Sun, Z. and Zhao, X. (2012) Compact Alternating Direction Implicit Scheme for the Two-Dimensional Fractional Diffusion-Wave Equation. SIAM Journal on Numerical Analysis, 50, 1535-1555.
 Wang, Z. and Vong, S. (2015) A High-Order ADI Scheme for the Two-Dimensional Time Fractional Diffusion-Wave Equation. International Journal of Computer Mathematics, 92, 970-979.
 Li, L., Xu, D. and Luo, M. (2013) Alternating Direction Implicit Galerkin Finite Element Method for the Two-Dimensional Fractional Diffusion-Wave Equation. Journal of Computational Physics, 255, 471-485.
 Fairweather, G., Yang, X., Xu, D. and Zhang, H. (2015) An ADI Crank-Nicolson Orthogonal Spline Collocation Method for the Two Dimensional Fractional Diffusion-Wave Equation. Journal of Scientific Computing, 65, 1217-1239.
 Cui, M. (2012) Compact Alternating Direction Implicit Method for Two-Dimensional Time Fractional Diffusion Equation. Journal of Computational Physics, 231, 2621-2633.
 Zhang, Y. and Sun, Z. (2014) Error Analysis of a Compact ADI Scheme for the 2D Fractional Subdiffusion Equation. Journal of Scientific Computing, 59, 104-128.
 Yao, W., Sun, J., Wu, B. and Shi, S. (2015) Numerical Simulation of a Class of Fractional Subdiffusion Equations via the Alternating Direction Implicit Method. Numerical Methods for Partial Differential Equations, 32, 531-547.
 Wang, Y. and Wang, T. (2016) Error Analysis of a High-Order Compact ADI Method for Two-Dimensional Fractional Convection-Subdiffusion Equations. Calcolo, 53, 301-330.
 Wang, Z. and Vong, S. (2014) A High-Order Exponential ADI Scheme for Two Dimensional Time Fractional Convection-Diffusion Equations. Computers & Mathematics with Applications, 68, 185-196.
 Zhai, S., Feng, X. and He, Y. (2014) An Unconditionally Stable Compact ADI Method for Three-Dimensional Time-Fractional Convection-Diffusion Equation. Journal of Computational Physics, 269, 138-155.
 Gao, G. and Sun, Z. (2015) Two Alternating Direction Implicit Difference Schemes with the Extrapolation Method for the Two-Dimensional Distributed-Order Differential Equations. Computers & Mathematics with Applications, 69, 926-948.
 Gao, G. and Sun, Z. (2016) Two alternating Direction Implicit Difference Schems for Two-Dimensional Distributed-Order Fractional Diffusion Equations. Journal of Scientific Computing, 69, 506-531.
 Wang, H. and Wang, K. (2011) An O(Nlog2N) Alternating-Direction Finite Difference Method for Two-Dimensional Fractional Diffusion Equations. Journal of Computational Physics, 230, 7830-7839.
 Wang, H. and Du, N. (2014) Fast Alternating-Direction Finite Difference Methods for Three-Dimensional Space-Fractional Diffusion Equations. Journal of Computational Physics, 258, 305-318.
 Tian, W., Zhou, H. and Deng, W. (2015) A Class of Second Order Difference Approximations for Solving Space Fractional Diffusion Equations. Mathematics of Computation, 294, 1703-1727.
 Chen, M. and Deng, W. (2014) A Second-Order Numerical Method for Two-Dimensional Two-Sided Space Fractional Convection Diffusion Equation. Applied Mathematical Modelling, 38, 3244-3259.
 Song, F. and Xu, C. (2015) Spectral Direction Splitting Methods for Two-Dimensional Space Fractional Diffusion Equations. Journal of Computational Physics, 299, 196-214.
 Liu, F., Zhuang, P., Turner, I., Anh, V. and Burrage, K. (2015) A Semi-Alternating Direction Method for a 2-D Fractional Fitz Hugh-Nagumo Monodomain Model on an Approximate Irregular Domain. Journal of Computational Physics, 293, 252-263.
 Bu, W., Tang, Y., Wu, Y. and Yang, J. (2015) Crank-Nicolson ADI Galerkin Finite Element Method for Two-Dimensional Fractional FitzHugh-Nagumo Mondomain Model. Applied Mathematics and Computation, 257, 355-364.
 Zeng, F., Liu, F., Li, C., et al. (2014) A Crank-Nicolson ADI Spectral Method for a Two-Dimensional Riesz Space Fractional Nonlinear Reaction-Diffusion Equation. SIAM Journal on Numerical Analysis, 52, 2599-2622.
 Fernandes, R.I. and Fairweather, G. (1991) An Alternating Direction Galerkin Method for a Class of Second-Order Hyperbolic Equations in Two Space Variables. SIAM Journal on Numerical Analysis, 28, 1265-1281.