1. Introduction

In recent times, the study of the fractional calculus and its applications in science and engineering has escalated [1–3]. The majority of papers dedicated to this topic discuss fractional kinetic equations of diffusion, diffusion-advection, and Fokker-Planck type to describe transport dynamics in complex systems that are governed by anomalous diffusion and nonexponential relaxation patterns [2, 3]. These papers provide comprehensive reviews of fractional/anomalous diffusion and an extensive collection of examples from a variety of application areas. A particular case of interest is the motion of solutes through aquifers discussed by Benson et al. [4, 5].

The generally accepted definition for the fractional Laplacian involves an integral representation (see [6] and the references therein) since the spectral resolution of the Laplacian operator over infinite domains is continuous; for the whole space, we use the Fourier transform and for initial value problems we use the Laplace transform in time [7]. However, when dealing with finite domains the fractional Laplacian subject to homogeneous boundary conditions is usually defined in terms of a summation involving the discrete spectrum. It is nontrivial to extend the latter definition to accommodate nonhomogeneous boundary conditions. To the best of our knowledge, there is no evidence in the literature that suggests this has been done apart from Ilić et al. [8] where the one-dimensional case was discussed. In this paper, we propose the extension to higher dimensions and illustrate the idea in the context of solving the fractional Poisson equation subjected to nonhomogeneous boundary conditions on a bounded domain.

Space fractional diffusion equations have been investigated by West and Seshadri [9] and more recently by Gorenflo and Mainardi [10, 11]. Numerical methods for these fractional equations are still under development. Hackbusch and his group [12–14] have developed the theory of -matrices and algorithms that they claim to be of complexity for computing functions of operators that are approximated by a finite difference (or other Galerkin schemes) discretisation matrix. However, the underlying theory is developed using integral representations of the matrix for separable coordinate systems and does not include a discussion of nonhomogeneous boundary conditions, which is essential for the fractional Poisson equation under investigation in this paper. Recently Ilić et al. [8, 15] proposed a matrix representation of the fractional-in-space operator to produce a system of linear ordinary differential equations (ODEs) with a matrix representation of the Laplacian operator raised to the same fractional power. This approach, which was coined the matrix transfer technique (MTT), enabled either the standard finite element, finite volume, or finite difference methods to be exploited for the spatial discretisation of the operator.

In recent years, fractional Brownian motion (FBM) with Hurst index has been used to introduce memory into the dynamics of diffusion processes. A prediction theory and other analytical results on FBM can be found in [16]. As shown in [17], a Girsanov-type formula for the Radon-Nikodym derivative of an FBM with drift with respect to the same FBM is determined by differential equations of fractional order with Dirichlet boundary conditions:for a certain integrable function defined on , where . In this study, we extend problem (1.1) and investigate the solution of a steady-state space fractional diffusion equation with sources, hereafter referred to as the fractional Poisson equation (FPE), on some bounded domain in two dimensions subject to either one or a combination of the usual (nonhomogeneous) boundary conditions of types I, II, or III imposed on the boundary . Although the method we present for solving the FPE is equally applicable to two- and three-dimensional problems and the various coordinate systems used in the solution by separation of variables, we consider only the following problem here.

FPE Problem
Solve the fractional Poisson equation in a finite rectangle:subject toWe choose such a simple region so that an analytic solution can be found, which can be used subsequently to verify our numerical approach. Note also that this system captures type I boundary conditions () and type II boundary conditions (). The latter case has to be analysed separately with care since is an eigenvalue that introduces singularities.
The use of our matrix transfer technique leads to the matrix representation of the FPE (1.2), which requires that the matrix function equationmust be solved. Note that in (1.4), denotes the matrix representation of the Laplacian operator obtained using any of the well-documented methods: finite difference, the finite volume method, or variational methods such as the Galerkin method using finite element or wavelets and , with a vector containing the discrete values of the source/sink term, and a vector that contains all of the discrete boundary condition information. We assume further that both the discretisation process and the implementation of the boundary conditions have been carried out to ensure that is symmetric positive definite, that is, .
The general solution of (1.4) can be written asand one notes the need to determine both the action of the matrix function on the vector and the action of the standard inverse on , where the matrix can be large and sparse.
In the case where , numerous authors have proposed efficient methods to deal directly with (1.5) using Krylov subspace methods and in particular, the preconditioned generalised minimum residual (GMRES) iterative method (see, e.g., the texts by Golub and Van Loan [18], Saad [19], and van der Vorst [20]). In this paper, we investigate the use of Krylov subspace methods for computing an approximate solution for a range of values and indicate how the spectral information gathered from at first solving can be recycled to obtain the complete solution in (1.5), where .
In literature, a majority of references deal with the extraction of an approximation to for scalar analytic function using Krylov subspace methods (see [21, Chapter 13] and the references therein). Druskin and Knizhnerman [22], Hochbruck and Lubich [23], Eiermann and Ernst [24], Lopez and Simoncini [25], van den Eshof et al. [26], as well as many other researchers use the Lanczos approximationwhereis the usual Lanczos decomposition, and the columns of form an orthonormal basis for Krylov subspace . However, as noted by Eiermann and Ernst [24], all basis vectors must be stored to form this approximation, which may prove costly for large matrices. Restarting the process is by no means as straightforward as for the case , and the restarted Arnoldi algorithm for computing given in [24] addresses this issue. Another issue worth pointing out is that although preconditioning linear systems is now well understood and numerous preconditioning strategies exist to accelerate the convergence of many iterative solvers based on Krylov subspace methods [19], preconditioning in many cases cannot be applied to . For example if , one can only deduce from in a limited number of special cases for .
In the previous work by the authors [27], we proposed a spectral splitting method , where is an orthogonal projector onto the invariant subspace associated with a set of eigenvalues on the “singular part” of the spectrum with respect to and an orthogonal projector onto the “regular part” of the spectrum. We refer to that part of the spectral interval where the function to be evaluated has rapid change with large values of the derivatives as the singular part (see [27] for more details). The splitting was chosen in such a way that was a low-degree polynomial (of degree at most 5). Thick restarting was used to construct the projector on the singular part. Unfortunately, the computational overhead associated with constructing the projector , whilst maintaining the requirement of a low-degree polynomial approximation for over the regular part, limits the application of the splitting method to a class of matrices that had fairly compact spectra. The method appeared to work well for applications in statistics [27, 28].
In this paper, we build upon the splitting method idea in the manner outlined as follows to approximate for monotone decreasing function .
(1) Determine an approximately invariant subspace (AIS), for the set of eigenvectors associated with the singular part of with respect to . Form and set , where are the eigenvalues associated with the eigenvectors . The thick restarted Lanczos method discussed in [27, 29] or [30] can be used for the AIS generation.(2) Let and generate orthonormal basis for .(3) Approximate using the Lanczos decomposition to analytic grade , [31].(4) Form . To avoid components of any eigenvectors associated with the singular part reappearing in , we show how this splitting strategy can be embedded in an adaptively constructed preconditioning of the matrix function.
The paper is organised as follows. In Section 2, we use MTT to formulate the matrix representation of FPE to accommodate nonhomogeneous boundary conditions. We also consider the approximation of the matrix function using the Lanczos method with thick restart and adaptive preconditioning. In Section 3, we give an upper bound on the error cast in terms of the linear system residual. In Section 4, we derive an analytic solution to the fractional Poisson equation using the spectral representation of the Laplacian, and in Section 5, we give the results of our algorithm when applied to two different problems, which highlight the importance of using our adaptively preconditioned Lanczos method. In Section 6, we give the conclusions of our work and hint at future research directions.

2. Matrix Function Approximation and Solution Strategy

The general numerical solution procedure MTT is implemented as follows. First apply a standard spatial discretisation process such as the finite volume, finite element, or finite difference method to the standard Poisson equation (i.e., in system (1.2)) in the case of homogeneous boundary conditions to obtain the following matrix form:where it is assumed that is the finite difference matrix representation of the Laplacian, and is the grid spacing. is the representation of , and is the representation of . Then, as was discussed in [15], the solution of FPE subject to homogeneous boundary conditions is approximated by the solution of the following matrix function equation:Next, we apply the same finite difference method to the homogeneous Poisson equation (i.e., Laplace's equation) with nonhomogeneous boundary conditions. The resulting equations can be written in the following matrix form:where represents the discretized boundary values, and the matrix is the same as given above. In other words, if does not satisfy homogeneous boundary conditions, then the modified representationis used, where denotes the extended definition of the Laplacian (see [8] and also refer to Section 4 for further details). Thirdly, we follow [8] to write the fractional Laplacian in the following form:and its matrix representation asHence, the matrix representation for FPE isAssuming that has an inverse, the solution of this equation is

Our aim is to devise an efficient algorithm to approximate the solution in (2.8) using Krylov subspace methods. One notes from (2.8) that the solution comprises two distinct components, , where , , and . We note further in this context that the scalar function is monotone decreasing on , where is symmetric positive definite.

There exists a plethora of Krylov-based methods available in the literature for approximately solving the linear system using, for example, conjugate gradient, FOM, or MINRES (see [19, 20]). Although preconditioning strategies are often employed to accelerate the convergence of many of these methods, we prefer not to adopt preconditioning here so that spectral information gathered about during this linear system solve can be recycled and used to aid the approximation of . As we will see, this recycling is affected through the use of thick restart [30, 32] and adaptive preconditioning [33, 34]. We emphasise that even if is a good preconditioner for , it may not be useful for since we cannot find a relation between and . Thus, many efficient solvers used for the ordinary Poisson equation cannot be employed for the FPE. The adaptive preconditioner, however, can.

We begin our presentation of the numerical algorithm by briefly reviewing the solution of the linear system , where is a symmetric positive definite using the full orthogonal method (FOM) [19] together with thick restart [27, 30, 32].

2.1. Stage 1—Thick Restarted, Adaptively Preconditioned, Lanczos Procedure

Suppose that the Lanczos decomposition of is given bywhere the columns of form an orthonormal basis for , and is the analytic grade defined in [31]. The analytic grade of order of the matrix with respect to is defined as the lowest integer for which , where is the orthogonal projector onto the th Krylov subspace and . The grade can be computed from the Lanczos algorithm using the matrices generated during the process. If is the 1st column of , and , for , then .

In each restart, or cycle, that follows, the Lanczos decomposition is carried up to the analytic grade , which could be different for different cycles. Consequently, for ease of exposition, the subscript will be suppressed so that the only subscript that appears throughout the description below refers to the cycle. Let be some initial approximation to the solution and define .

Cycle 1
(i)Generate Lanczos decompositionwhere , , is tridiagonal, , , , and .(ii)Obtain approximate solution , so thatand residualTest if . If yes, stop; otherwise, continue to cycle 2.

Cycle 2
(i)Find eigenvalue decomposition of , that is, , where .(ii)Select the orthonormal (ON) eigenvectors, , of corresponding to the smallest in magnitude eigenvalues of and form the Ritz vectorswhere are ON, and let the associated Ritz values be stored in the diagonal matrix .(iii)Set and generate the thick-restart Lanczos decompositionwhere , , , and(iv)Obtain approximate solution , so thatand residualTest if . If yes, stop; otherwise, continue to the next cycle.

Cycle
(i)Find eigenvalue decomposition of , that is, .(ii)Select orthonormal (ON) eigenvectors, , of corresponding to the smallest in magnitude eigenvalues of and form the Ritz vectors .(iii)Set and generate thick-restart Lanczos decompositionwhere has similar form as .(iv)Obtain approximate solution , so thatand residualTest if . If yes, stop; otherwise, continue cycling.

2.1.1. Construction of an Adaptive Preconditioner

Another important ingredient in the algorithm described above is the construction of an adaptive preconditioner [33, 34]. Let the thick-restart procedure at cycle produce the approximate smallest Ritz pairs , where . We then check if any of these Ritz pairs have converged to approximate eigenpairs of by testing the magnitude of the upper bound on the eigenpair residualThe eigenpairs deemed to have converged are then locked and used to construct an adaptive preconditioner that can be employed during the next cycle to ensure that difficulties such as spuriousness can be avoided.

Suppose we collect the locked Ritz vectors as columns of the matrix , set , and formwhere . , are the current estimates of the smallest and largest eigenvalues of , respectively, obtained from the restart process. Then, has the same eigenvectors as ; however its eigenvalues are shifted to [33, 34]. Furthermore, it should be noted that these preconditioners can be nested. If is a sequence of such preconditioners, then with and , we haveThus, during the cycles (say cycle ) the adaptively preconditioned, thick- restart Lanczos decompositionis employed.

Note. The preconditioner does not need to be explicitly formed; it can be applied in a straightforward manner from the stored locked Ritz pairs.

In summary, stage 1 consists of employing the adaptively preconditioned Lanczos procedure outlined above to approximately solve the linear system for . At the completion of this process, the residual , and we have the set of locked Ritz pairs. This spectral information is then passed to accelerate the performance of stage 2 of the solution process.

2.2. Stage 2—Matrix Function Approximation Using an Adaptively Preconditioned Lanczos Procedure

At the completion of stage 1, we have generated an approximately invariant eigenspace associated with the smallest in magnitude eigenvalues of . We now show how this spectral information can be recycled to aid with the approximation of , where .

2.2.1. Adaptive Preconditioning

Recall from stage 1 that we have available , where . The important observation at this point is the following relationship between and .

Proposition 2.1. Let be an eigenspace of symmetric matrix such that , with and . Define , then for ,

Proof. Let , then Furthermore,
Thus,By noting thatwe obtain the main result

The following proposition shows that, as was the case for the solution of the linear system in stage 1, these preconditioners can be nested in the case of the matrix function approximation.

Proposition 2.2. Let be a sequence of preconditioners as defined in Proposition 2.1, then

Proof. Let and , then observe that and

Corollary 2.3. Under the hypothesis of Proposition 2.1, one notes the equivalent form of (2.25) as which appears similar to the idea of spectral splitting proposed in [27].

We now turn our attention to the approximation of , which by using Corollary 2.3 can be expressed aswhere . First note that if , then . We expand the Lanczos decomposition to the analytic grade of with . Next perform the spectral decomposition of and set , then compute the Lanczos approximationBased on the theory presented to this point, we propose the following algorithm to approximate the solution of the fractional Poisson equation.

Algorithm 2.4 (Computing the Solution of the FPE Problem). Stage 1. Solve using the thick restarted adaptively preconditioned Lanczos method and generate the AIS, . Return the preconditioner , where .
Stage 2. Compute using the following strategy.
(1)Set .(2)Compute Lanczos decomposition , where is the analytic grade of and , with .(3) Perform the spectral decomposition .(4)Compute linear system residual and estimate from to compute bound (3.9) derived in Section 3.(5)If bound is small, then approximate and exit to step (6), otherwise continue the Lanczos expansion until bound is satisfied.(6) Form , where .
Finally, compose the approximate solution of FPE as .

Remarks
At stage 2, we monitor the upper bound given in Proposition 3.3 to check if the desired accuracy is achieved in the matrix function approximation. If the desired level is not attained, then it may be necessary to repeat the thick-restart procedure to determine the next smallest eigenvalues and their corresponding ON eigenvectors. In fact, this process may need to be repeated until there are no eigenvalues remaining in the “singular” part so that the accuracy of the approximation is dictated entirely by that of the linear system residual. We leave the design of this more sophisticated and generic algorithm for future research.

It is natural at this point to ask what is the accuracy of the approximation (2.33) for a given ? Not knowing at the outset makes it impossible to answer this question. Instead, we opt to provide an upper bound for the error , which is the topic of the following section.

3. Error Bounds for the Numerical Solution

At first, we note that Churchill [35] uses complex integration around a branch point to derive the following:By changing the variable, one can deduce the following expression, for :Noting that , the spectral decomposition and the usual definition of the matrix function enable the following expression for computing to be obtained:

Recall that the approximate solution of the linear system from using the Galerkin approach (FOM or CG) is given by , with residual . We note the similarity to (2.33); however a key observation is that the error in the matrix function approximation cannot be determined in such a straightforward manner as for the linear system [24]. The following proposition enables the error in the matrix function approximation to be expressed in terms of the integral expression given above in (3.3) and the residual of what is called a shifted linear system.

Proposition 3.1. Let be the residual to the shifted linear system , then

Proof. It is known that

It is interesting to observe that for the Lanczos approximation, so that the vectors