An ADI Method for the Numerical Solution of 3D Fractional Reaction-Diffusion Equations

A numerical method for solving fractional partial differential equations (fPDEs) of the diffusion and reaction–diffusion type, subject to Dirichlet boundary data, in three dimensions is developed. Such fPDEs may describe fluid flows through porous media better than classical diffusion equations. This is a new, fractional version of the Alternating Direction Implicit (ADI) method, where the source term is balanced, in that its effect is split in the three space directions, and it may be relevant, especially in the case of anisotropy. The method is unconditionally stable, second-order in space, and third-order in time. A strategy is devised in order to improve its speed of convergence by means of an extrapolation method that is coupled to the PageRank algorithm. Some numerical examples are given.


Introduction
Diffusive flows, possibly with reaction, play an important role in a number of fields, such as fluid dynamics in porous media and percolation, as well as seepage and groundwater hydraulics; see, e.g., [1][2][3][4]. However, the celebrated Darcy's law, along with the assumption of continuity of the seepage flow in heterogeneous media and the ensuing diffusion partial differential equation (PDE), is not valid in the general case of real flows. In 1998, J. He [5] proposed a generalized version of the Darcy's law, containing fractional Riemann-Liouville (RL) derivatives, where p ≡ p(x, y, z, t) is the pressure, q = (q x , q y , q z ) is the fluid velocity, K = (K x , K y , K z ) is the percolation tensor (here assumed to be diagonal), and α := (α 1 , α 2 , α 3 ). In this paper, we chose to follow He's idea in order to model seepage flows in porous media, adopting fractional Riemann-Liouville derivatives, even though several other kinds of fractional derivative exist. The Riemann-Liouville fractional derivative of order α of F(x) is defined by where n ∈ N is such that n − 1 < α ≤ n, see [6], e.g., the literature concerning Fractional Calculus is still rapidly increasing. We only mention [7] for some review on progresses in this field, in particular concerning RL operators, and [8] for reviewing some general basic properties of fractional operators, in particular, RL operators. Therefore, the more general equation for seepage flow with fractional Riemann-Liouville derivatives that was proposed by J. He [5] reads x, y, z, t), (x, y, z) ∈ Ω. (2) Here, 1 ν is the specific storage coefficient (that is assumed to be constant), f (p; x, y, z, t) is a "source" or "sink term", Ω is a given open bounded domain of R 3 , with piecewise smooth boundary, where the seepage or percolation takes place, t is time, and 0 < β i < 1, 0 < α i ≤ 1, with 1 < β i + α i ≤ 2, for i = 1, 2, 3, see [2]. In (2), the operator ∂ α 1 ∂x α 1 denotes the fractional Riemann-Liouville derivative of order α 1 with respect to x. Note that, when the source term, f , depends on p, the fractional partial differential equation (fPDE) in (2) is a reaction-diffusion, rather than a pure diffusion equation. The boundary and initial conditions p(x, y, z, t)| ∂Ω = Φ(x, y, z, t), p(x, y, z, 0) = ϕ(x, y, z).
are imposed, with Φ and ϕ being suitable prescribed functions.
In this paper, we consider, rather, fractional diffusions in the nonconservative form, where all of the space fractional derivatives are of the RL type, defined as in (1), the set ν = 1, γ i := β i + α i , for i = 1, 2, 3. Below, we construct a "fractional Alternating Direction Implicit" (fADI) scheme in order to numerically solve such a problem. We will consider such a fractional reaction-diffusion equation with (positive) constant percolation coefficients, K x , K y , K z , in (2), with boundary and initial conditions, as above. Unlike the case of classical (integer order) derivatives, the two model Equations (2) and (4) do not coincide, even when the coefficients K x , K y , and K z are constants. In fact, it is not true, in general, that, e.g., This is only true under some conditions, for instance, when u ∈ C k (i.e., is sufficiently smooth) and α 1 , β 1 + α 1 ∈ [ − 1, ] for some ∈ N, ≤ k, see [9]. This restriction on the fractional orders is quite serious, since we expect that 0 < α 1 , β 1 < 1 and 1 < β 1 + α 1 < 2, in order to obtain an anomalous (sub-)diffusion equation. Otherwise, one should just face the fPDE in (2), even when taking out of the derivatives the constants K, and keeping it as it is. There is a name for such fractional equations, which is "sequential fractional differential equations", see, e.g., ([6], Section 4.2).
Other definitions of fractional derivative exist. One of the most useful in practice has been proven to be Caputo's definition. The Caputo fractional derivative is given by where n − 1 < α < n, n ∈ N. Moreover, the Grünwald-Letnikov fractional derivative is particularly relevant for us, in view of the ensuing numerical treatment. It is defined as [2,6], while the so-called shifted Grünwald-Letnikov derivative is given by where [x] denotes the integer part of x, with h being the step size to be used in the numerical schemes [6], and Γ is the Gamma function. The "normalized" Grünwald weights are defined as g α,i = (−1) i ( α i ). These coefficients will be used in the numerical approximation of the Grünwald-Letnikov fractional discrete operator in (7) or (8); see Section 2.
When applied to smooth functions, for example functions with continuous (n − 1)th-order derivatives and integrable nth-order derivatives, all of these definitions are known to coincide [6]. Therefore, in this paper, we consider all fractional derivatives appearing in the fPDEs in the sense of Caputo, while the discretizations will be made according to the shifted Grünwald-Letnikov definition.
Given the interest for models that are based on fPDEs and, needless to say, analytic solutions to most fPDEs are usually not available, a number of authors have proposed, over the years, a few numerical methods to solve them, see, e.g., [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28]. However, numerical methods that are capable of handling high-dimensional fPDEs seem to be rather few in the existing literature, not to mention their performance. It should be observed that, being the fractional derivatives nonlocal operators, discretizing fPDEs can be expected to be more demanding in terms of memory storage and CPU time, as compared to classical PDEs. All of this calls for devising effective numerical methods in order to solve fPDE problems already in dimension 2 and 3.
In this paper, we consider the three-dimensional (3D) fractional seepage flow equation in (2), on the bounded domain, Ω := (0, L x ) × (0, L y ) × (0, L z ), and the time interval (0, T]. Moreover, we assume that such an initial-boundary value problem equation has a unique, sufficiently smooth solution, when the (homogeneous) Dirichlet boundary conditions and initial condition p(x, y, z, t)| ∂Ω = Φ(x, y, z), p(x, y, z, 0) = ϕ(x, y, z), for some given functions, Φ and ϕ are imposed. A new, well "balanced", fractional version of the Alternating Direction Implicit (ADI) method is introduced, in order to numerically solve 3D diffusion and reaction-diffusion fPDEs. We call our algorithm "balanced", since it distributes the source term effect in an equal amount among the three equations of the ADI scheme, instead of putting the full source on the right-hand side of the first equation only, see Section 2. This choice seems to be effective, especially when the model presents some degree of anisotropy, and it allows us to increase the space accuracy. In [2], Liu et al. derived a scheme that is accurate to the second order in space, resorting to a Richardson extrapolation, but our balanced scheme turns out to perform better, see Table 12 below.
Our method is shown to be unconditionally stable for every fractional order of space derivatives, second-order accurate in space, and third-order accurate in time. This latter property is established in this paper adopting an extrapolation technique along with an optimization strategy that is realized through Google's PageRank algorithm. The method followed here to do this for the 3D case is similar to that used for the two-dimensional (2D) case in [29].
Here is the plan of the paper. In Section 2, we formulate a typical problem, and discuss some fractional versions of the ADI (fADI) 3D schemes. In Section 3, we establish the convergence and stability of our algorithm, while, in Section 4, some details regarding optimization and extrapolation are given, and a few numerical examples are presented. Section 5 contains the conclusions of the paper, briefly highlighting its main points.

Fractional Diffusion and Seepage Flow in Homogeneous Media
In order to construct our fADI scheme on a bounded domain, we first discretize space and time, as usual, setting x i := ih x , for i = 0, 1, 2, . . . , M x , with h x := L x /M x ; y j := jh y , for j = 0, 1, 2, . . . , M y , with h y := L y /M y ; z k := kh z , for k = 0, 1, 2, . . . , M z , with h z := L z /M z ; and, t n := nτ, for n = 0, 1, 2, . . . , N, τ := T/N. Here, M x , M y , M z , and N are positive integers, while h x , h y , h z , and τ are the space and time step sizes, respectively. The numerical approximation to p(x i , y j , z k , t n ) provided by the scheme at such points will be then denoted by p n i,j,k . Similarly, we write f n i,j,k in order to denote an approximation of f (p(x i , y j , z k , t n ); x i , y j , z k , t n ) for the source term, p 0 i,j,k = ϕ i,j,k for ϕ(x i , y j , z k , 0), the initial value, and p n 0,j, y j , t n ) = 0, and p n i,j,M z = Φ 6 (x i , y j , t n ), for the boundary conditions.
We discretize the fPDE, approximating the first-order derivative ∂p ∂t in (4) by forward finite differences. Assuming that the solution, p, has a first-order continuous space derivative and that its second-order space derivatives are integrable, the operators ∂ γ 1 ∂x γ 1 , ∂ γ 2 ∂y γ 2 , and ∂ γ 3 ∂z γ 3 in (4) can be discretized while using the shifted fractional Grünwald-Letnikov derivative in (8). Thus, we obtain the implicit finite-difference scheme It is known that the fractional discrete operator provides an O(h x ) approximation to the Grünwald-Letnikov shifted fractional derivative of order γ 1 [20]. Similarly, the discrete operators provide, respectively, O(h y ) and O(h z ) approximations to the Grünwald-Letnikov shifted fractional derivatives of order γ 2 and γ 3 . Using such operators, the implicit difference scheme in (10) can be rewritten in a more compact form, as Note that, evaluating f n i,j,k on the right-hand side at time t n , makes the case simpler when f depends on p, which is when the fPDE is a reaction-diffusion rather than a diffusion equation. In Section 3, we will show that this scheme is characterized by a local truncation error of the order , and it is unconditionally stable. However, adopting suitable strategies, we have been able to improve such orders, attaining the order O(τ) . Equation (13) (or (10)) yields a linear system of equations satisfied by p n+1 i,j,k . he corresponding system's matrix, however, is neither sparse nor band structured, as it happens in the corresponding classical ADI method. This implies that the scheme in (13) requires, at each step, the solution of a large dense linear system of equations. Therefore, we are facing a computationally demanding numerical problem, which calls for constructing suitable efficient numerical schemes, possibly unconditionally stable.
To this purpose, we exploit the idea of the classical ADI method in order to design an implicit difference scheme for each direction, but now in the framework of fractional differential equations. Indeed, we split the computations into three steps, each requiring a reduced computational load. On the first step, we solve the problem in the x-direction, on the second one, we solve it in the y-direction, and in the third one, in the z-direction.
In order to devise such a splitting, we add an extra higher-order term to the left-hand side of (13), so that we are able to factor the operator into three factors, into three factors, without affecting the overall convergence rate, thus obtaining the scheme Now, splitting (15) in the three dimensions, we obtain the unbalanced version of the scheme, which provides the solution at time t n+1 , and turns out to be of the first order in space [2]. However, we can choose to distribute the effect of the source term among the three steps in a more balanced way, obtaining, from (15), the fADI scheme which also provide the solution at time t n+1 . As in all classical ADI schemes, the system in (17) is not equivalent to that in (15), but it does within an error that is of order higher that that provided by the scheme. The algorithm (17) now performs better, as it will be shown in Section 3, below. The "unbalanced" scheme is accurate to the first order in space and it requires a little shorter computing time, as compared to the better "balanced" scheme. The latter however is accurate to the second order in space, see Theorem 1, below.
The idea of using a more balanced scheme, which equally affects all (the three) directions, results in a better method, which is clearly uneffective when f ≡ 0. The algorithm seems to be new, even within the framework of 3D classical (i.e., non-fractional) ADI methods.
Summarizing, the computations are split into three (time) fractional steps, in each of which a one-dimensional equation is solved at time, as follows: Step 1. Solve the problem in the x-direction (for each fixed pair (y j , z k )), in order to obtain the intermediate value of the solution, say p n+1/3 i,j,k , from the first equation in (17).
Step 2. Solve the problem in the y-direction (for each fixed pair (x i , z k )), in order to obtain the intermediate value p n+2/3 i,j,k , from the second equation in (17), while using the results that were obtained at Step 1.
Step 3. Solve the problem in the z-direction (for each fixed pair (x i , y j )), from the third equation in (17), using the results of Step 2.
Taking the boundary values p n+1/3 0,j,k and p n+2/3 M x ,j,k into account, we can construct the coefficients of the matrix A := {a s,t } of the linear system in (17): for each fixed (j, k), we have Subsequently, using the boundary values p n+2/3 i,0,k and p n+2/3 i,M y ,k , we obtain the matrix B := {b s,t } of the linear system in the second equation of (17). For each fixed (i, k), such a matrix turns out to be similar to A (as given in (18)). Finally, while using the boundary values p n+1 i,j,0 and p n+1 i,j,M z , we obtain the matrix C := {c s,t } of the linear system in the third equation in (17), which, again, for each fixed (i, j), is similar to A.
Of course, as in the classical ADI Peaceman-Rachford algorithm [30], used to solve "integer-order" (i.e., non-fractional) PDEs, it is necessary to provide certain boundary values to be associated to the diffusion in the x-direction, which is p n+1/3 0,j,k and p n+1/3 M x ,j,k , as well as p n+2/3 i,0,k and p n+2/3 i,M y ,k , associated to the diffusion along y, and similarly for the diffusion along z, when solving the systems with the coefficient matrices A, B, and C. For instance, we need the boundary values p n+1/3 0,j,k , p n+1/3 M x ,j,k , which can be obtained from Examining the three matrices A, B, and C, above, it can be seen that, at each time step, it is merely required to solve, for each fixed pair (j, k) (at each x-level), a linear upper triangular system of size M x − 1, or, for each fixed pair (i, k) (at each y-level), a linear upper triangular system of size M y − 1, or, for each fixed pair (i, k) (at each z-level) a linear upper triangular system of size M z − 1.
The "unbalanced" implicit difference scheme (16) is characterized by a local truncation error of , and it is unconditionally stable. These facts were established in [32]. Let us now consider our better "balanced" scheme (17). We will prove the following Theorem 1. The scheme in (17) is unconditionally stable, third-order accurate in time, and second-order accurate in space, for every choice of the (space) fractional orders, γ i , i = 1, 2, 3.

Proof.
In [33], the unbalanced scheme in (16) was analyzed and its numerical stability proved. Now consider our balanced scheme in (17). Here, we use the so called "Lubich second-order backward finite difference" [34], which is defined, at time n, in the x-direction, at the point x i+1 , and for fixed y and z, as L δ γ 1 x p Here, c ∈ {−1, 0, −1}, γ 1 represents the fractional order derivative, and the ω γ 1 c 's are the first three coefficients of the Taylor expansions of the corresponding generating function, say f a (z), see [34]. Here, a denotes the order of accuracy and, for a = 2, such a function is given by The Grünwald-Letnikov definition of the discretized Riemann-Liouville fractional derivative differs from that given by Lubich by an O(h 2 x ) term, see [34]. Thus, evaluating the fractional Riemann-Liouville derivatives at the point x i+1 , in order to approximate the space fractional derivative in the sense of Lubich of order γ 1 in the x-direction, we obtain where τ = h x . In order to be more general, we could set τ := s h x , for some s > 0, but this would not entail any loss of generality because the stability analyis would not affect the location of the roots by this. This holds for each fixed i and n, j, k. Using the approximation (21) in (17), we obtain with a truncation error of order O(h 2 x ). We now approximate the x-space operator applied to p n+1/3 i+1,j,k [33] by a second-order Taylor expansion, obtaining where D x := ∂ xx and L x : (see (29) in Section 4.1 below, and (32) in Section 4.2), where * and * * represent, respectively, the first and the second intermediate sub-steps between n and n + 1/3, and p n+1/3 i+1,j,k is implicitly defined through (22). We can then insert L x in (23)p n+1/3 i+1,j,k , obtaining Subtracting L xp n+1/3 i+1,j,k from both sides of (25), we finally obtain where we set, for short, v := p n+1/3 i+1,j,k −p n+1/3 i+1,j,k , and replaced the classical second order derivative D x with his fractional version, δ γ 1 x . A von Neumann stability analysis of (26) yields an equation for the amplification factor, σ, Similar results can be obtained for the y and the z directions. It can be proved that the roots of (27) are in the range 0 ≤ |σ| ≤ 1, for every value of K x , K y , K z , θ, h x , and τ [35]. Consequently, the scheme presented in (26) turns out to be of the second order in space, and unconditionally stable. The third-order accuracy in time can be established proceeding as in the 2D case worked out in [29]. We stress that such a result follows from the optimization strategy adopted there, in particular, from using the four-steps extrapolation method (see step 3, in Section 4.1 below) coupled with the Page-Rank algorithm.

The Numerical Implementation and Some Examples
In this section, we apply our numerical method in order to solve 3D fractional diffusion as well as reaction-diffusion problems. Examples, including reaction-diffusions, are then given in order to illustrate the performance of the algorithm.

Some Theoretical Preliminary Considerations
We set up an extrapolation technique to optimize our algorithm, thus attaining an appreciable acceleration. This can be accomplished exploiting the PageRank algorithm, which is widely used by Computer Scientists, but not so much yet in Numerical Analysis.
The Google's PageRank algorithm [36] makes it possible to index material in Internet, while using the degree of popularity of a given web page, in order to define the position in searching results through the computation of certain optimal coefficients. In this way, this algorithm provides, for instance, an optimization method to make it faster a given web search. This approach will be associated to an extrapolation technique, in order to determine the coefficients that make it optimal. Not only can a considerable amount of CPU time be saved in this way, [10], but this procedure also increases the accuracy of the overall method up to the third order in time, see [29,37]. Let us describe such a technique.
Step 1. Compute first the array Ψ = (ψ 1 , ψ 2 , ψ 3 , ψ 4 ) T , solving the system of four linear algebraic equations, where d is the so-called damping factor, in general assumed to be around 0.85 (see [38]; such a numerical value for d was determined by empirical trials), v : N being the dimension of the problem, which, in our case, is equal to 4, and L = { ij } for i, j = 1, . . . , 4, where l ij are certain nonnegative coefficients, such that, for each j, ∑ N i=1 l i,j = 1. Note that (28) is an implicit, but linear, equation for Ψ; hence, immediately solvable. An extrapolated solution, depending on p n ≡ p n i,j,k , is then used to solve the problem. The quantity p n requires evaluating certain coefficients, which can be obtained by the "PageRank accelerating algorithm" [36]. The previous algebraic system yields the optimal coefficients array Ψ.

Numerical Examples
In this section, we present a few numerical results in order to support the theoretical analysis developed in the previous sections. In the following, the symbol a D α b p will be used to define the fractional Riemann-Liouville derivative of p of order α, in the interval [a, b].

The (exact) analytical solution to Equation (30) under such initial and boundary conditions is known and it is given by
see [42]. All of this can be used in order to validate our 3D algorithm and test its performance. We will then feel confident that it might perform well also when other source terms, that might be encountered in real world problems, enter the model. Below, we will indeed replace the special forcing term in (30) with some other more realistic ones.
From Figure 1, it is rather clear (qualitatively) what is the effect of replacing ordinary derivatives, (α, β, γ) = (2, 2, 2), with fractional derivatives, in a given diffusion equation. Here, (α, β, γ) = (1.4, 1.5, 1.6), thus accounting for anomalous diffusion. In general, such a modification implies new geometric patterns in the solution, and a possibly anisotropic behavior. These features have been observed, e.g., in certain porous media through which some fluid flows [43,44]. Hereafter, T denotes the final time at which the solution is computed. In Figure 2, the numerical errors q n − p n are plotted in the L 2 and in the L ∞ norms, on a log-scale, where p n ≡ p(x i , y j , z k , t n ) denotes the exact solution to the fPDE (on the grid points that are defined in Section 2, and q n ≡ q n (τ) is its approximation that is given by our scheme through Equation (29).
In order to assess the convergence rate when space and time step sizes are refined, we write our approximate solution q n as q(h, τ) in order to display its dependence on the space step sizes, h x = h y = h z = h (assumed to be equal), and the time step size, τ. Subsequently, we define, for every N ∈ N, the time and the space convergence rate, in the two-norm, for a given method of order m > 1, as In Tables 1 and 2, the absolute numerical errors q n − p n ∞ and q n − p n 2 , as well as the corresponding convergence rates that are attained using different space step sizes, are shown. Figure 2 displays the infinity norm and the L 2 norm errors for Example 1, when the optimized difference scheme is implemented in (29), at times T = 1 and T = 2, for several values of h x = h y = h z = h = τ = 2/N, N = 8, 16, 32, 64, and fixed α, β, and γ. Here, there is evidence that the numerical errors decrease when increasing N, i.e., decreasing τ. Note that, due to the logarithmic scale, decreased errors correspond to increasing values.  Table 1. L ∞ and L 2 norm errors (in log scale) and convergence rates for Example 1, when the balanced scheme is used, at time T = 1, for several values of (α, β, γ), N.
(α, β, γ) N q n − p n ∞ r space q n − p n 2 r space r time ( In Figure 3, the discrepancies between the numerical solutions to the classical and the fractional differential equations, obtained by the ADI and the f ADI method, respectively, i.e., q n ADI − q n f ADI , are plotted in the L 2 and the L ∞ norms, on a log-scale. This is done in order to appreciate the differences that are visible when one switches from one to the other model. In Tables 3 and 4, such discrepancies, q n ADI − q n f ADI ∞ and q n ADI − q n f ADI 2 , as well as the corresponding convergence rates, which are achieved using different space step sizes, are shown. Here, again, q n (τ), given by 147 2500 p n (τ) − 29, 409 100, 000 is the extrapolated solution, and p n ≡ p n i,j,k satisfies the optimized scheme in (13). Figure 3 shows the convergence rates of the maximum norm and the L 2 norm errors for Example 1, when the optimized difference scheme that is presented in (32) is implemented, at times T = 1, for several values of N, and fixed α, β, and γ.    Table 5 shows the CPU times by our algorithm on a Dual Core Pentium with 4GB RAM and with 32bit-MATLAB programs, for different values of fractional orders and advection coefficients d x , d y , d z . The CPU time is measured in terms of clocks ticks. Table 5. CPU times in clocks ticks units (CT) for the unbalanced and the balanced versions of the fADI scheme, CT bal , and CT unbal , respectively, required attaining an error of order 10 −5 , for h = 1/100 and several values of the fractional orders, α, β, γ, and the advection coefficients, d x , d y , d z .
(α, β, γ) ( Tables 5 and 6 indicate that the our balanced ADI method's performance is comparable to that proposed in [2], although slightly worse. Recall that our method and that of Liu et al. differ in the strategies that are followed to accelerate them. However, our unbalanced method performs better than that concerning the required CPU times, which are smaller. This is due to the fact that the model problem is anisotropic, behaving rather differently along the direction x than along the other directions. Indeed, all Tables 7-9 show that when the anisotropic behavior is more pronounced in another direction (others than x), an appropriately balanced method wins over the others. Table 6. CPU times in clocks ticks units (CT) for the unbalanced and the balanced versions of the fractional Alternating Direction Implicit (fADI) scheme, CT bal , and CT unbal , respectively, required attaing an error of order 10 −5 , for N = 100 and several values h = 1 γN , and of the fractional orders, α, β, γ.  Table 7. CPU times in clocks ticks units (CT) for the unbalanced and the balanced versions of the fADI scheme, CT bal , and CT unbal , respectively, required attaining an error of order 10 −5 , for h = 1/100 and several values of the fractional orders, α, β, γ, and the advection coefficients d x , d y , d z .   (30), (α, β, γ) = (1.4, 1.5, 1.6), and the (linear) source term f (p) = p 2 , instead of the linear source f (x, y, z, t), independent of p, as defined in (31) above. This is a fractional linear reaction-diffusion problem. A forcing term like this often occurs, for instance, in modeling dissolution and precipitation phenomena in porous media [45,46].
Here, the forcing term mimics the behavior of a seepage flow in a porous medium, which can be either isotropic or even anisotropic in the x, y, and z direction (depending on the coefficients K x , K y , and K z ). In view of a numerical treatment of such a problem, we replace the ideally impulsive source with a function that is likely to have a similar behavior, see [47]. Let α = β = γ = 1.8 and the domain Γ = (0, 1) 3 , and the source term, which approximately reproduces the behavior of the δ function in (33), We impose the boundary conditions p(0, y, z, t) = p(x, 0, y, t) = p(x, y, 0, t) = 0, p(1, y, z, t) = e −t y 2 z 2 , p(x, 1, y, t) = e −t x 2.2 z 2 , p(x, y, 1, t) = e −t x 2.2 y 2 , and the initial value p 0 (x, y, z) = 0.
The analytic solution to this problem is known to be see [2]. Note that, when (α, β, γ) = (2, 2, 2), we obtain the classical linear forced diffusion equation p t = ∆p + f . We now compare the previous numerical solution with α, β, γ in the range [1.2, 1.8], computed using our balanced fADI method, with that of a classical (i.e., non-fractional) case, which we can consider rather "close" to it, taking α = β = γ = 2 and the same source term that is defined in (34).
In Figure 4, the absolute numerical errors between the exact solution and numerical solution obtained by a fADI method are shown. In Figure 5, the same is done for the absolute numerical errors between the numerical classical ADI method and the numerical solution obtained by the fADI method. Both of the errors are plotted in the L 2 norm.   Table 11 illustrates the convergence rates of the algorithm, as measured in the L 2 norm, in Example 3 for several values of τ and N, and fixed α, β, and γ. Here, p n is the exact solution, while q n ADI is the solution that was obtained by the classical ADI method, and q n f ADI is that provided by the fADI method. Table 11. L 2 norm errors, discrepancies and convergence rates for Example 3, when the balanced scheme is used, with α = β = γ = 1.8, at time T = 1, for h = 1/N and several values of N and τ. The two numbers in the rates columns refer to the first and second algorithms, respectively.
N τ q n f ADI − p n 2 q n f ADI − q n As one would expect in the fADI methods, coupled with the optimized extrapolation, a faster convergence is observed as the grid is refined. In the following, we show these results in a table of computational times, see Table 12.  Table 13 shows the CPU times and memory used by our algorithm for various grid sizes. The CPU time is evaluated at T = 1. Table 13 indicates that the balanced ADI method performs slightly worse than that one proposed in [2]. However, this little price is paid in order to obtain a third-order (instead of a second-order) in time method. It was also observed that balanced and unbalanced algorithms require approximately the same computational resources. Example 4. We consider the Fisher equation, which is the semilinear (reaction-diffusion) equation [48] ∂p on the set  [49,50], which is often adopted in population biology to model the spread of invasive species. Here, r and k are parameters, with r representing the intrinsic increase rate of the fluid, and k the environmental carrying capacity, which is the maximum sustainable fluid density. We computed the numerical solution to the initial-boundary value problem above, assuming the radially symmetric initial condition p(x, y, z, 0) = min (x,y,z)∈Ω 1 {0.8, 10 e −(x 2 +y 2 +z 2 ) }, (36) and the boundary conditions p(x, y, z, t)| ∂Ω = p y (x, ±20, z, t) = p x (±20, y, z, t) = p x (10, y, z, t) = 0.
Next, we consider the solution to a fPDE, like that in (35), but with a reaction term having space variable coefficients. Let C = 0.15, D = 0, 4, E = 1, r = 0.2, and assume that k varies in space, as follows: for (x, y, z) ∈ Ω := Ω 1 ∪ Ω 2 ∪ Ω 3 ,  In many applications to meteorology, this case may model a region that can contain, at most, a certain given amount of fluid, due to certain environmental conditions. The geometry is a silt barrier through which the fluid will eventually penetrate. We first consider the case α = β = γ = 2. In Figure 5 the solution pertaining to this case is shown, in a plan view, at time T = 90. Here, the pressure can be seen to penetrate the barrier very slowly due to classical diffusion along x. In many applications to meteorology, this case may model a region which can contains at most a certain given amount of fluid, due to certain environmental conditions. The geometry is a silt barrier through which the fluid will eventually penetrate. We first consider the case α = β = γ = 2. In Figure 5 the solution pertaining to this case is shown, in a plan view, at time T = 90. Here, due to classical diffusion along x, the pressure can be seen to penetrate the barrier very slowly. We then changed parameters, choosing α = 1.7, β = γ = 2, to represent a certain anomalous diffusion. In Figure 8, it is shown that at time T = 50, even before the "snapshot time" T = 90 (which refers to Figure 7), the pressure has penetrated significantly the barrier, and it spreads at the same time in the y direction. This is a striking peculiarity of fractional reaction-diffusion model which allows to predict the effects of controlling the fluid flow, for instance describing pollution due to a pollutant expanding from a given environment to another.  (36) and parameters α = β = γ = 2, C = 0.15, D = 0.4, r = 0.2, on the cross section in the plane (x, y), this plane is divided into two regions, the pressure propagates slowly from a region (left) to the other (right) through a silt barrier that links them [48].
We then changed parameters, choosing α = 1.7, β = γ = 2, in order to represent a certain anomalous diffusion. In Figure 8, it is shown that, at time T = 50, even before the "snapshot time" T = 90 (which refers to Figure 7), the pressure has significantly penetrated the barrier, and it spreads at the same time in the y direction. This is a striking peculiarity of fractional reaction-diffusion model, which allows for predicting the effects of controlling the fluid flow, for instance, describing pollution due to a pollutant expanding from a given environment to another.

Conclusions and Future Directions
In this paper, we introduced a new, well "balanced", fractional version of the ADI method, for numerically solving a number of 3D diffusion as well as reaction-diffusion problems for fPDEs, which are important in a variety of problems in porous media modelling. We proved that such a method is unconditionally stable for every fractional order of space derivatives, second-order accurate in space, and third-order accurate in time. The speed of convergence has been improved while adopting an extrapolation technique, coupled with the optimization method that is used by the PageRank algorithm.
Future directions of research should consider: (a) using Riesz fractional space operators, in order to symmetrize the RL operators; (b) better investigating the anisotropic behavior of solutions in modeling physical phenomena, possibly designing numerical schemes that might properly take anisotropy into account (adopting differente space step-sizes in different directions, so to make the numerical solution more efficient); (c) take the issue of imposing boundary conditions (BCs) into account, while facing nonlocal operators, when the problem is posed on bounded domains. (Here, we implicitly assumed Dirichlet BCs identically zero solutions outside the domain).