A Numerical Solution Using an Adaptively Preconditioned Lanczos Method for a Class of Linear Systems Related with the Fractional Poisson Equation

Abstract

This study considers the solution of a class of linear systems related with the fractional Poisson equation (FPE) (−∇²)^𝛼/2𝜑=𝑔(𝑥,𝑦) with nonhomogeneous boundary conditions on a bounded domain. A numerical approximation to FPE is derived using a matrix representation of the Laplacian to generate a linear system of equations with its matrix 𝐴 raised to the fractional power 𝛼/2. The solution of the linear system then requires the action of the matrix function 𝑓(𝐴)=𝐴^−𝛼/2 on a vector 𝑏. For large, sparse, and symmetric positive definite matrices, the Lanczos approximation generates 𝑓(𝐴)𝑏≈𝛽₀𝑉_𝑚𝑓(𝑇_𝑚)𝑒₁. This method works well when both the analytic grade of 𝐴 with respect to 𝑏 and the residual for the linear system are sufficiently small. Memory constraints often require restarting the Lanczos decomposition; however this is not straightforward in the context of matrix function approximation. In this paper, we use the idea of thick-restart and adaptive preconditioning for solving linear systems to improve convergence of the Lanczos approximation. We give an error bound for the new method and illustrate its role in solving FPE. Numerical results are provided to gauge the performance of the proposed method relative to exact analytic solutions.

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 𝑁log𝑁 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 𝐻∈(0,1) 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:−∇²ℎ(𝑥^𝛼/2)ℎ(=𝑥)0=i𝑔f(𝑥𝑥)∉i(f0,𝑥𝑇∈),(0,𝑇),(1.1)for a certain integrable function ℎ(𝑥) defined on [0,𝑇], where 𝑔∶[0,𝑇]→ℝ. 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:−∇²^𝛼/2𝜑=𝑔(𝑥,𝑦),0<𝑥<𝑎,0<𝑦<𝑏,(1.2)subject to−𝑘₁𝜕𝜑𝜕𝑥+ℎ¹𝑘𝜑=𝑓₁(𝑦)at𝑥=0,₂𝜕𝜑𝜕𝑥+ℎ²−𝑘₃𝜑=𝑓₂(𝑦)at𝑥=𝑎,𝜕𝜑𝜕𝑦+ℎ³𝑘𝜑=𝑓₃(𝑥)at𝑦=0,₄𝜕𝜑𝜕𝑦+ℎ⁴𝜑=𝑓₄(𝑥)at𝑦=𝑏.(1.3)We 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 (𝑘_𝑖=0,ℎ_𝑖=1,𝑖=1,…,4) and type II boundary conditions (ℎ_𝑖=0,𝑘_𝑖=1,𝑖=1,…,4). The latter case has to be analysed separately with care since 0 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 equation𝐴^𝛼/2Φ=𝑏(1.4)must 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 𝑏=𝑏₁+𝐴^𝛼/2−1𝑏₂, 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, 𝐴∈SPD.
The general solution of (1.4) can be written asΦ=𝐴^−𝛼/2𝑏=𝐴^−𝛼/2𝑏₁+𝐴⁻¹𝑏₂,(1.5)and one notes the need to determine both the action of the matrix function 𝑓(𝐴)=𝐴^−𝛼/2 on the vector 𝑏₁ and the action of the standard inverse on 𝑏₂, where the matrix 𝐴 can be large and sparse.
In the case where 𝛼=2, 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 0<𝛼<2 and indicate how the spectral information gathered from at first solving 𝐴Φ₂=𝑏₂ can be recycled to obtain the complete solution Φ=Φ₁+Φ₂ in (1.5), where Φ₁=𝑓(𝐴)𝑏₁=𝐴^−𝛼/2𝑏₁.
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 approximation𝑓𝐴𝑣≈𝑉_𝑚𝑓𝑇_𝑚𝑒₁,𝑣=‖‖𝑣‖‖𝑉_𝑚𝑒₁,(1.6)where𝐴𝑉_𝑚=𝑉_𝑚𝑇_𝑚+𝛽_𝑚𝑣_𝑚+1𝑒^𝑇𝑚(1.7)is the usual Lanczos decomposition, and the columns of 𝑉_𝑚 form an orthonormal basis for Krylov subspace 𝒦_𝑚(𝐴,𝑣)={𝑣,𝐴𝑣,…,𝐴^𝑚−1𝑣}. 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 𝑓(𝑡)=1/𝑡, 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 SPD 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), span{𝑞₁,…,𝑞_𝑘} for the set of eigenvectors associated with the singular part of 𝜎(𝐴) with respect to 𝑓(𝑡). Form 𝑄_𝑘=[𝑞₁,…,𝑞_𝑘] and set Λ_𝑘=diag{𝜆₁,𝜆₂,…,𝜆_𝑘}, where 𝜆_𝑖 are the eigenvalues associated with the eigenvectors 𝑞_𝑖,𝑖=1,…,𝑘. 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 ℓ, 𝐴𝑉_ℓ=𝑉_ℓ𝑇_ℓ+𝛽_ℓ𝑣_ℓ+1𝑒^𝑇ℓ [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., 𝛼=2 in system (1.2)) in the case of homogeneous boundary conditions to obtain the following matrix form:1ℎ²𝐴Φ=𝑔̃,(2.1)where it is assumed that (1/ℎ²)𝐴=𝑚(−∇²) 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:1ℎ^𝛼𝐴^𝛼/2Φ=𝑔̃.(2.2)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:1ℎ²𝐴Φ−𝑏=0,(2.3)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 representation𝑚((−∇²)𝜑)=1ℎ²𝐴Φ−𝑏(2.4)is 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:−∇²^𝛼/2=−∇²^𝛼/2−1(−∇²),(2.5)and its matrix representation as𝑚(−∇²^𝛼/2)=𝑚(−∇²^𝛼/2−1)𝑚((−∇²)).(2.6)Hence, the matrix representation for FPE is1ℎ^𝛼𝐴^𝛼/2Φ=𝑔̃+1ℎ^𝛼−2𝐴^𝛼/2−1𝑏.(2.7)Assuming that 𝐴 has an inverse, the solution of this equation isΦ=ℎ^𝛼𝐴^−𝛼/2𝑔̃+ℎ²𝐴⁻¹𝑏.(2.8)

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 0<𝑞=𝛼/2<1. 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 by𝐴𝑉_ℓ=𝑉_ℓ𝑇_ℓ+𝛽_ℓ𝑣_ℓ+1𝑒^𝑇ℓ=𝑉_ℓ+1𝑇_ℓ,(2.9)where 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 ‖𝑢_ℓ−𝒫_ℓ𝑢_ℓ‖/‖𝑢_ℓ‖<10^−𝑡, 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 𝑡_𝑖=𝑇_𝑖𝑡_𝑖−1, for 𝑖=1,…,ℓ, then ‖𝑢_ℓ−𝒫_ℓ𝑢_ℓ‖/‖𝑢_ℓ‖=|𝑒^𝑇𝑙+1𝑡_𝑙|/‖𝑡_𝑙‖.

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 decomposition𝐴𝑉₁=𝑉₁𝑇₁+𝛽₁𝑢₁𝑒^𝑇,(2.10)where 𝑉₁=[𝑣¹⁽¹⁾,…,𝑣^ℓ(1)], 𝑣¹⁽¹⁾=𝑟₀/𝛽₀, 𝑇₁ is tridiagonal, 𝑢₁=𝑣^(1)ℓ+1, 𝛽₀=‖𝑟₀‖, 𝛽₁=𝛽^ℓ(1), and 𝑒^𝑇=𝑒^𝑇ℓ.(ii)Obtain approximate solution Φ²⁽¹⁾=𝑉₁𝑇¹⁻¹𝑉^𝑇1𝑟₀, so thatΦ²⁽¹⁾=Φ²⁽⁰⁾+𝑉₁𝑇¹⁻¹𝑉^𝑇1𝑟₀,(2.11)and residual𝑟₁=𝑏−𝐴Φ²⁽¹⁾=𝑟₀−𝐴Φ²⁽¹⁾=−𝛽₁𝑒^𝑇𝑇¹⁻¹𝑉^𝑇1𝑟₀𝑢¹.(2.12)Test if ‖𝑟₁‖<𝜀. If yes, stop; otherwise, continue to cycle 2.

Cycle 2
(i)Find eigenvalue decomposition of 𝑇₁, that is, 𝑇₁𝑌=𝑌Λ, where Λ=diag{𝜃₁,…,𝜃_ℓ}.(ii)Select the 𝑘 orthonormal (ON) eigenvectors, 𝑌₁, of 𝑇₁ corresponding to the 𝑘 smallest in magnitude eigenvalues of 𝑇₁ and form the Ritz vectors𝑊₁=𝑉₁𝑌₁=[𝑤₁,…,𝑤_𝑘],(2.13)where 𝑤_𝑖 are ON, and let the associated Ritz values be stored in the diagonal matrix Λ₁=diag{𝜃₁,…,𝜃_𝑘}.(iii)Set 𝑉₂=[𝑊₁,𝑢₁] and generate the thick-restart Lanczos decomposition𝐴𝑉₂=𝑉₂𝑇₂+𝛽₂𝑢₂𝑒^𝑇,(2.14)where 𝑉₂=[𝑤₁,…,𝑤_𝑘,𝑣¹⁽²⁾,…,𝑣^ℓ(2)], 𝑣¹⁽²⁾=𝑢₁, 𝑢₂=𝑣^(2)ℓ+1, and𝑇₂=⎡⎢⎢⎢⎢⎢⎣Λ¹𝛽𝛽₂𝑠₁0⋯₂0𝑠^𝑇1𝛼𝛽_𝑘+1𝑘+1𝛽_𝑘+10⋯…⋱⋱⋱…⎤⎥⎥⎥⎥⎥⎦,with𝑠¹=𝑌^𝑇1𝑒.(2.15)(iv)Obtain approximate solution Φ²⁽²⁾=𝑉₂𝑇²⁻¹𝑉^𝑇2𝑟₁, so thatΦ²⁽²⁾=Φ²⁽¹⁾+𝑉₂𝑇²⁻¹𝑉^𝑇2𝑟₁,(2.16)and residual𝑟₂=𝑏−𝐴Φ²⁽²⁾=𝑟₁−𝐴Φ²⁽²⁾=−𝛽₂𝑒^𝑇𝑇²⁻¹𝑉^𝑇2𝑟₁𝑢².(2.17)Test if ‖𝑟₂‖<𝜀. If yes, stop; otherwise, continue to the next cycle.

Cycle (𝑗+1)
(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 𝑉_𝑗+1=[𝑊_𝑗,𝑢_𝑗] and generate thick-restart Lanczos decomposition𝐴𝑉_𝑗+1=𝑉_𝑗+1𝑇_𝑗+1+𝛽_𝑗+1𝑢_𝑗+1𝑒^𝑇,(2.18)where 𝑇_𝑗+1 has similar form as 𝑇₂.(iv)Obtain approximate solution Φ^2(𝑗+1)=𝑉_𝑗+1𝑇^−1𝑗+1𝑉^𝑇𝑗+1𝑟_𝑗, so thatΦ^2(𝑗+1)=Φ^2(𝑗)+Φ^2(𝑗+1)=Φ²⁽⁰⁾+^𝑗+1_𝑖=1𝑉_𝑖𝑇^𝑖−1𝑉^𝑇𝑖𝑟_𝑖−1,(2.19)and residual𝑟_𝑗+1=𝑏−𝐴Φ^2(𝑗+1)=−𝛽_𝑗+1𝑒^𝑇𝑇^−1𝑗+1𝑉^𝑇𝑗+1𝑟_𝑗𝑢^𝑗+1.(2.20)Test if ‖𝑟_𝑗+1‖<𝜀. 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 {𝜃_𝑖,𝑤_𝑖}^𝑘𝑖=1, 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 residual‖‖𝐴𝑤_𝑖−𝜃_𝑖𝑤_𝑖‖‖≤𝛽_𝑗||𝑒^𝑇𝑦_𝑖||<𝜀₂.(2.21)The 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 Λ_𝑗=diag{𝜃₁,…,𝜃_𝑝}, and form𝑀^𝑗−1=𝛾𝑄_𝑗Λ^𝑗−1𝑄^𝑇𝑗+𝐼−𝑄_𝑗𝑄^𝑇𝑗,(2.22)where 𝛾=(𝜃_min+𝜃_max)/2. 𝜃_min, 𝜃_max are the current estimates of the smallest and largest eigenvalues of 𝐴, respectively, obtained from the restart process. Then, 𝐴_𝑗=𝐴𝑀^𝑗−1 has the same eigenvectors as 𝐴; however its eigenvalues {𝜆_𝑖}^𝑝𝑖=1 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 Λ=diag(Λ_𝑖,𝑖=1,…,𝑗), we have𝑀⁻¹=𝑀^𝑗−1⋯𝑀²⁻¹𝑀¹⁻¹=𝛾𝑄Λ⁻¹𝑄^𝑇+𝐼−𝑄𝑄^𝑇.(2.23)Thus, during the cycles (say cycle 𝑗+1) the adaptively preconditioned, thick- restart Lanczos decomposition𝐴𝑀⁻¹𝑉_𝑗+1=𝑉_𝑗+1𝑇_𝑗+1+𝛽_𝑗+1𝑢_𝑗+1𝑒^𝑇(2.24)is 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 {𝜃_𝑖,𝑞_𝑖}^𝑘𝑖=1 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 𝒱=span{𝑞₁,𝑞₂,…,𝑞_𝑘} 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 𝑀⁻¹=𝛾𝑄_𝑘Λ^𝑘−1𝑄^𝑇𝑘+𝐼−𝑄_𝑘𝑄^𝑇𝑘, where 𝑄_𝑘=[𝑞₁,…,𝑞_𝑘]. The important observation at this point is the following relationship between 𝑓(𝐴) and 𝑓(𝐴𝑀⁻¹).

Proposition 2.1. Let span{𝑞₁,𝑞₂,…,𝑞_𝑘} be an eigenspace of symmetric matrix 𝐴 such that 𝐴𝑄_𝑘=𝑄_𝑘Λ_𝑘, with 𝑄_𝑘=[𝑞₁,𝑞₂,…,𝑞_𝑘] and Λ_𝑘=diag(𝜇₁,…,𝜇_𝑘). Define 𝑀=(1/𝛾)𝑄_𝑘Λ_𝑘𝑄^𝑇𝑘+𝐼−𝑄_𝑘𝑄^𝑇𝑘, then for 𝑣∈ℝ^𝑛, 𝑓(𝐴)𝑣=1𝑓(𝛾)𝑓(𝐴𝑀⁻¹)𝑓(𝛾𝑀)𝑣.(2.25)

Proof. Let 𝑊𝑊^𝑇=𝐼−𝑄_𝑘𝑄^𝑇𝑘,𝑊^𝑇𝐴𝑊=𝐵, then 𝑀𝑀⁻¹=𝑄_𝑘𝑄^𝑇𝑘+𝑊𝑊^𝑇=𝐼=𝑀⁻¹𝑀. Furthermore,
𝑀⁻¹𝐴=𝑄_𝑘𝑊𝛾0𝐼𝐵0𝑄^𝑇𝑘𝑊^𝑇=𝐴𝑀⁻¹.(2.26)Thus,𝑓𝐴𝑀⁻¹=𝑄_𝑘𝑊𝑓0𝛾𝐼𝑓0𝐵𝑄^𝑇𝑘𝑊^𝑇.(2.27)By noting that𝑓𝐴=𝑄_𝑘𝑊𝑓Λ^𝑘0𝑓0𝐵𝑄^𝑇𝑘𝑓𝛾𝑀=𝑄_𝑘𝑊^𝑇,𝑊𝑓Λ^𝑘0𝑓0𝛾𝐼𝑄^𝑇𝑘𝑊^𝑇,(2.28)we obtain the main result𝑓(𝐴)𝑓(𝛾)[𝑓(𝛾𝑀)]⁻¹=𝑄_𝑘𝑊𝑓0𝛾𝐼𝑓0𝐵𝑄^𝑇𝑘𝑊^𝑇=𝑓(𝐴𝑀⁻¹).(2.29)

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 𝑓(𝐴)𝑣=1𝑓(𝛾)𝑓(𝐴𝑀¹⁻¹𝑀²⁻¹⋯𝑀^𝑗−1)𝑓(𝛾𝑀₁𝑀₂⋯𝑀_𝑗)𝑣.(2.30)

Proof. Let 𝑄=[𝑄₁,𝑄₂,…,𝑄_𝑘] and Λ=diag(Λ_𝑖,𝑖=1,…,𝑗), then observe that 𝑀=𝑀₁𝑀₂⋯𝑀_𝑗=(1/𝛾)𝑄Λ𝑄^𝑇+𝐼−𝑄𝑄^𝑇 and 𝑓(𝐴)=𝑓(𝐴𝑀⁻¹)(1/𝑓(𝛾))𝑓(𝛾𝑀).

Corollary 2.3. Under the hypothesis of Proposition 2.1, one notes the equivalent form of (2.25) as 𝑓(𝐴)𝑣=𝑄_𝑘𝑓(Λ_𝑘)𝑄^𝑇𝑘𝑣+𝑓(𝐴𝑀⁻¹)(𝐼−𝑄_𝑘𝑄^𝑇𝑘)𝑣,(2.31)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 as𝐴^−𝑞𝑔̃=^𝑘_𝑖=1𝜃^𝑖−𝑞𝑞_𝑖𝑞^𝑇𝑖𝑔̃+(𝐴𝑀⁻¹)^−𝑞𝑔̂,(2.32)where 𝑔̂=(𝐼−𝑄_𝑘𝑄^𝑇𝑘)𝑔̃. First note that if 𝐴∈SPD, then 𝐴𝑀⁻¹∈SPD. We expand the Lanczos decomposition 𝐴𝑀⁻¹𝑉_ℓ=𝑉_ℓ𝑇_ℓ+𝛽_ℓ𝑣_ℓ+1𝑒^𝑇ℓ to the analytic grade ℓ of 𝐴𝑀⁻¹ with 𝑣₁=𝑔̂/‖𝑔̂‖. Next perform the spectral decomposition of 𝑇_ℓ=𝑌_ℓΛ_ℓ𝑌^𝑇ℓ and set 𝑄_ℓ=𝑉_ℓ𝑌_ℓ, then compute the Lanczos approximation(𝐴𝑀⁻¹)^−𝑞𝑔̂≈𝑉_ℓ𝑇^ℓ−𝑞𝑉^𝑇ℓ𝑔̂=𝑄_ℓΛ^ℓ−𝑞𝑄^𝑇ℓ𝑔̂.(2.33)Based 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, 𝒬_𝑘=span{𝑞₁,…,𝑞_𝑘}. Return the preconditioner 𝑀=(1/𝛾)𝑄_𝑘Λ_𝑘𝑄^𝑇𝑘+𝐼−𝑄_𝑘𝑄^𝑇𝑘, where 𝑄_𝑘=[𝑞₁,…,𝑞_𝑘].
Stage 2. Compute Φ₁=𝐴^−𝑞𝑔̃ using the following strategy.
(1)Set 𝑔̂=(𝐼−𝑄_𝑘𝑄^𝑇𝑘)𝑔̃.(2)Compute Lanczos decomposition 𝐴𝑀⁻¹𝑉_ℓ=𝑉_ℓ𝑇_ℓ+𝛽_ℓ𝑣_ℓ+1𝑒^𝑇ℓ, where ℓ is the analytic grade of 𝐴𝑀⁻¹ and 𝑉_ℓ=[𝑣₁,…,𝑣_ℓ], with 𝑣₁=𝑔̂/‖𝑔̂‖.(3) Perform the spectral decomposition 𝑇_ℓ=𝑌_ℓΛ_ℓ𝑌^𝑇ℓ.(4)Compute linear system residual ‖𝑟_ℓ‖=|𝛽_ℓ𝑒^𝑇ℓ𝑌_ℓΛ^ℓ−1𝑌^𝑇ℓ𝑉^𝑇ℓ𝑔̂| and estimate 𝜆_min≈𝜇_min from 𝑇_ℓ to compute bound (3.9) 𝜇^−𝑞_min‖𝑟_ℓ‖ 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:^∞0𝑥^−𝑞𝑥+1𝑑𝑥=sin𝜋(𝑞𝜋).(3.1)By changing the variable, one can deduce the following expression, for 𝜆^−𝑞,𝜆>0:𝜆^−𝑞=sin(𝑞𝜋)(1−𝑞)𝜋^∞0𝑑𝑡𝑡^1/(1−𝑞)+𝜆.(3.2)Noting that 𝐴=𝐴𝑀⁻¹∈SPD, the spectral decomposition and the usual definition of the matrix function enable the following expression for computing 𝐴^−𝑞 to be obtained:𝐴^−𝑞=sin(𝑞𝜋)(1−𝑞)𝜋^∞0𝑡^1/(1−𝑞)𝐼+𝐴⁻¹𝑑𝑡.(3.3)

Recall that the approximate solution of the linear system 𝐴𝑥=𝑣 from 𝒦_ℓ(𝐴,𝑣) using the Galerkin approach (FOM or CG) is given by 𝑥_ℓ=𝑉_ℓ𝑇^ℓ−1𝑉^𝑇ℓ𝑣, with residual 𝑟_ℓ=𝑏−𝐴𝑥_ℓ=−(𝛽_ℓ𝑒^𝑇ℓ𝑇^ℓ−1𝑉^𝑇ℓ𝑣)𝑣_ℓ+1. 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 𝑟_ℓ(𝑡)=𝑣−(𝐴+𝑡^1/(1−𝑞)𝐼)𝑉_ℓ(𝑇_ℓ+𝑡^1/(1−𝑞)𝐼)⁻¹𝑉^𝑇ℓ𝑣 be the residual to the shifted linear system (𝐴+𝑡^1/(1−𝑞)𝐼)𝑥=𝑣, then 𝐴^−𝑞𝑣−𝑉_ℓ𝑇^ℓ−𝑞𝑉^𝑇ℓ𝑣=sin(𝑞𝜋)(1−𝑞)𝜋^∞0𝑡^1/(1−𝑞)𝐼+𝐴⁻¹𝑟_ℓ(𝑡)𝑑𝑡.(3.4)

Proof. It is known that 𝐴^−𝑞𝑣−𝑉_ℓ𝑇^ℓ−𝑞𝑉^𝑇ℓ𝑣=sin(𝑞𝜋)(1−𝑞)𝜋^∞0=sin(𝑞𝜋){𝑡^1/(1−𝑞)𝐼+𝐴⁻¹−𝑉_ℓ𝑡^1/(1−𝑞)𝐼+𝑇_ℓ⁻¹𝑉^𝑇ℓ}𝑣𝑑𝑡(1−𝑞)𝜋^∞0𝑡^1/(1−𝑞)𝐼+𝐴⁻¹{𝑣−𝑡^1/(1−𝑞)𝐼+𝐴𝑉_ℓ𝑡^1/(1−𝑞)𝐼+𝑇_ℓ⁻¹𝑉^𝑇ℓ𝑣}𝑑𝑡.(3.5)

It is interesting to observe that 𝑟_ℓ(𝑡)=−(𝛽_ℓ𝑒^𝑇ℓ(𝑇_ℓ+𝑡^1/(1−𝑞)𝐼)⁻¹𝑉^𝑇ℓ𝑣)𝑣_ℓ+1 for the Lanczos approximation, so that the vectors 𝑟_ℓ≡𝑟_ℓ(0) and 𝑟_ℓ(𝑡) are aligned; however their magnitudes are different. Note further that𝑟_ℓ(𝑡)=𝑒^𝑇ℓ(𝑇_ℓ+𝑡𝐼)⁻¹𝑒₁𝑒^𝑇ℓ𝑇^ℓ−1𝑒₁𝑟_ℓ(0).(3.6)An even more important result is the following relationship between their norms.

Proposition 3.2. Let 𝑇_ℓ have eigendecomposition 𝑇_ℓ𝑌_ℓ=𝑌_ℓΛ_ℓ, where Λ_ℓ=diag(𝜇_𝑖,𝑖=1,…,ℓ) with the 𝜇_𝑖 Ritz values for the Lanczos approximation, then for 𝑡>0, ‖‖𝑟_ℓ(𝑡)‖‖=|^ℓ_𝑖=1𝜇_𝑖𝜇_𝑖+𝑡^1/(1−𝑞)|‖‖𝑟_ℓ‖‖≤‖‖𝑟_ℓ‖‖.(3.7)

Proof. The result follows from [26], which gives the following polynomial characterisations for the residuals:𝑟_ℓ=𝜋^ℓ(𝐴)𝑣𝜋_ℓ(0),𝜋^ℓ𝑟(𝜏)=det(𝜏𝐼−𝑇_ℓ)=∏^ℓ_𝑖=1(𝜏−𝜇_𝑖),_ℓ(𝑡)=𝜋^𝑡ℓ𝐴+𝑡^1/(1−𝑞)𝐼𝑣𝜋^𝑡ℓ(0),𝜋^𝑡ℓ(𝜏)=det𝜏𝐼−𝑇_ℓ+𝑡^1/(1−𝑞)𝐼=∏^ℓ_𝑖=1𝜏−(𝜇^𝑖+𝑡^1/(1−𝑞)),(3.8)so that 𝑟_ℓ(𝑡)=𝜋_ℓ(𝐴)𝑣/𝜋^𝑡ℓ(0)=(𝜋_ℓ(0)/𝜋^𝑡ℓ(0))𝑟_ℓ=∏^ℓ𝑖=1(𝜇_𝑖/(𝜇_𝑖+𝑡^1/(1−𝑞)))𝑟_ℓ. The result follows by taking the norm and noting that 𝑡>0.

We are now in a position to formulate an error bound essential for monitoring the accuracy of the Lanczos approximation (2.33).

Proposition 3.3. Let 𝜆_min be the smallest eigenvalue of 𝐴 and 𝑟_ℓ the linear system residual obtained by solving the linear system 𝐴𝑥=𝑔̂ using FOM on the Krylov subspace 𝒦_ℓ(𝐴,𝑔̂), then for 0<𝑞<1, one has ‖‖𝐴^−𝑞𝑔̂−𝑉_ℓ𝑇^ℓ−𝑞𝑉^𝑇ℓ𝑔̂‖‖≤𝜆^−𝑞min‖‖𝑟_ℓ‖‖.(3.9)

Proof. Using the orthogonal diagonalisation 𝐴=𝑄Λ𝑄^𝑇, we obtain from Proposition 3.1 that𝐴^−𝑞𝑔̂−𝑉_ℓ𝑇^ℓ−𝑞𝑉^𝑇ℓ𝑔̂=sin(𝑞𝜋)(1−𝑞)𝜋^∞0𝑄𝑡^1/(1−𝑞)𝐼+Λ⁻¹𝑄^𝑇𝑟_ℓ(𝑡)𝑑𝑡.(3.10)The result follows by taking norms and using Proposition 3.2 to obtain‖‖𝐴^−𝑞𝑔̂−𝑉_ℓ𝑇^ℓ−𝑞𝑉^𝑇ℓ𝑔̂‖‖≤sin(𝑞𝜋)(1−𝑞)𝜋^∞0‖‖‖diag1𝑡^1/(1−𝑞)𝐼+𝜆_𝑖,𝑖=1,…,𝑛‖‖‖𝑑𝑡‖‖𝑟^ℓ‖‖.(3.11)

The importance of this result is that it relates the error in the matrix function approximation to a scalar multiple of the linear system residual. This bound can be monitored during the Lanczos decomposition to deduce whether a specified tolerance has been reached in the matrix function approximation. Another key observation from Proposition 3.3 is that it motivates us to shift the small troublesome eigenvalues of 𝐴, via some form of preconditioning, so that 𝜆_min≈1. In this way, the error in the function approximation is dominated entirely by the residual error.

4. Analytic Solution

In this section, we discuss the analytic solution of the fractional Poisson equation, which can be used to verify the numerical solution strategy outlined in Section 2. The theory depends on the definition of the operator (−∇²)^𝛼/2 via spectral representation. The one-dimensional case was discussed in Ilić et al. [8], and the salient results for two dimensions are repeated here for completeness.

4.1. Homogeneous Boundary Conditions

In operator theory, functions of operators are defined using spectral decomposition. Set Ω={(𝑥,𝑦)∣0<𝑥<𝑎,0<𝑦<𝑏}, and let ℋ be the real space ℒ₂(Ω) with real inner product ⟨𝑢,𝑣⟩=∫∫_Ω𝑢𝑣𝑑𝑆. Consider the operator𝑇𝜑=−𝜕²𝜕𝑥²+𝜕²𝜕𝑦²𝜑=(−Δ)𝜑(4.1)on 𝒟={𝜑∈ℋ∣𝜑_𝑥,𝜑_𝑦 absolutely continuous; 𝜑_𝑥,𝜑_𝑦,𝜑_𝑥𝑥,𝜑_𝑥𝑦,𝜑_𝑦𝑦∈ℒ₂(Ω),ℬ(𝜑)=0}, where ℬ(𝜑) is one of the boundary conditions in the FPE problem with right-hand side equal to zero.

It is known that 𝑇 is a closed self-adjoint operator whose eigenfunctions {𝜑_𝑖𝑗}^{∞𝑖,𝑗=1} form an orthonormal basis for ℋ. Thus, 𝑇𝜑_𝑖𝑗=𝜆^2𝑖𝑗𝜑_𝑖𝑗,𝑖,𝑗=1,2,…. For any 𝜑∈𝒟,𝜑=^∞_∞𝑖=1_𝑗=1𝑇𝜑=_∞𝑐_𝑖𝑗𝜑_𝑖𝑗,𝑐_𝑖𝑗=⟨𝜑,𝜑_𝑖𝑗⟩,_∞𝑖=1_𝑗=1𝜆^2𝑖𝑗𝑐_𝑖𝑗𝜑_𝑖𝑗.(4.2)If 𝜓 is a continuous function on ℝ, then𝜓(𝑇)𝜑=^∞_∞𝑖=1_𝑗=1𝜓(𝜆^2𝑖𝑗)𝑐_𝑖𝑗𝜑_𝑖𝑗,(4.3)provided that ∑^∞𝑖=1∑^∞𝑗=1|𝜓(𝜆^2𝑖𝑗)𝑐_𝑖𝑗|²<∞. Hence, if the eigenvalue problem for 𝑇 can be solved for the region Ω, then the FPE problem with homogeneous boundary conditions can be easily solved to give𝜑(𝑥,𝑦)=^∞_∞𝑖=1_𝑗=1⟨𝑔,𝜑_𝑖𝑗⟩𝜆^𝛼𝑖𝑗𝜑_𝑖𝑗(𝑥,𝑦).(4.4)

4.2. Nonhomogeneous Boundary Conditions

Before we proceed further, we need to specify the definition of −(∇²)^𝛼/2.

Definition 4.1. Let {𝜑_𝑖𝑗} be a complete set of orthonormal eigenfunctions corresponding to eigenvalues 𝜆^2𝑖𝑗(≠0) of the Laplacian (−Δ) on a bounded region Ω with homogeneous BCs on 𝜕Ω. Let𝔉_𝛾={𝑓∈𝒟(Ω)∣^∞_∞𝑖=1_𝑗=1|||𝜆_𝛾𝑖𝑗||||||𝑐_𝑖𝑗|||²<∞,𝑐_𝑖𝑗=⟨𝑓,𝜑_𝑖𝑗⟩,𝛾=max(𝛼,0)}.(4.5)Then, for any 𝑓∈𝔉_𝛾, (−Δ)^𝛼/2𝑓 is defined by−Δ^𝛼/2𝑓=^∞_∞𝑖=1_𝑗=1𝜆^𝛼𝑖𝑗𝑐_𝑖𝑗𝜑_𝑖𝑗.(4.6)If one of 𝜆^2𝑖𝑗=0 and 𝜑₀ is the eigenfunction corresponding to this eigenvalue, then one needs ⟨𝑓,𝜑₀⟩=0.

Proposition 4.2. (1)The operator 𝑇=(−Δ)^𝛼/2 is linear and self-adjoint, that is, for 𝑓,𝑔∈𝔉_𝛾, ⟨𝑇𝑓,𝑔⟩=⟨𝑓,𝑇𝑔⟩.(2)If 𝑓∈𝔉_𝛾, where 𝛾=max(0,𝛼,𝛽,𝛼+𝛽), then(−Δ)^𝛼/2(−Δ)^𝛽/2𝑓=(−Δ)^𝛼+𝛽^/2𝑓=(−Δ)^𝛽/2(−Δ)^𝛼/2𝑓.(4.7)

For 𝛼>0, Definition 4.1 may be too restrictive, since the functions we are interested in satisfy nonhomogeneous boundary conditions, and the resulting series may not converge or not converge uniformly.

Extension of Definition 4.1
(1)For 𝛼=2𝑚,𝑚=0,1,2,…, define (−Δ)^𝛼/2𝑓=(−Δ)^𝑚𝑓 for any 𝑓∈𝐶^2𝑚(Ω) (or other possibilities).(2)For 𝑚−1<𝛼/2<𝑚,𝑚=1,2,…, define (−Δ)^𝛼/2𝑓=𝑇𝑔, where 𝑔=(−Δ)^𝑚−1𝑓∈𝐶^2(𝑚−1)(Ω), and 𝑇 is the extension of 𝑇=(−Δ)^{𝛼/2+1−𝑚} as defined by Proposition 4.3 below.
It suffices to consider 0<𝛼<2.

Proposition 4.3. Let 𝜑_𝑖𝑗(𝑥,𝑦) be an eigenfunction corresponding to the eigenvalue 𝜆^2𝑖𝑗(≠0) of the Laplacian (−Δ) on the rectangle Ω, and let 𝜑 satisfy the BCs in problem 1. Then, if 𝑇 is an extension of 𝑇=(−Δ)^𝛼/2 (in symbols 𝑇⊂𝑇) with adjoint 𝑇^∗⊂𝑇^∗, ⟨𝜑_𝑖𝑗,𝑇𝜑⟩=𝜆^𝛼𝑖𝑗⟨𝜑_𝑖𝑗,𝜑⟩−𝜆^{𝛼−2𝑖𝑗}^𝑏0𝜑_𝑖𝑗0,𝑦𝑘₁𝑓₁𝑦𝑑𝑦+^𝑏0𝜑_𝑖𝑗𝑎,𝑦𝑘₂+𝑓₂𝑦𝑑𝑦_𝑎0𝜑_𝑖𝑗𝑥,0𝑘₃𝑓₃𝑥𝑑𝑥+^𝑎0𝜑_𝑖𝑗𝑥,𝑏𝑘₄𝑓₄𝑥𝑑𝑥(4.8)if 𝑘_𝑖≠0,𝑖=1,…,4. If 𝑘_𝑖=0, the second term on the right-hand side becomes −𝜆^{𝛼−2𝑖𝑗}^𝑏0𝜕𝜑_𝑖𝑗0,𝑦𝜕𝑥𝑓¹𝑦ℎ₁𝑑𝑦−^𝑏0𝜕𝜑_𝑖𝑗𝑎,𝑦𝜕𝑥𝑓²𝑦ℎ₂+𝑑𝑦_𝑎0𝜕𝜑_𝑖𝑗𝑥,0𝜕𝑦𝑓³𝑥ℎ₃𝑑𝑥−^𝑎0𝜕𝜑_𝑖𝑗𝑥,𝑏𝜕𝑦𝑓⁴𝑦ℎ₄𝑑𝑥.(4.9)

Proof. It is known that ⟨𝜑_𝑖𝑗,𝑇𝜑⟩==𝑇𝑇^∗𝜑_𝑖𝑗,𝜑_∗=𝑇𝜑𝜑_{𝑖𝑗𝑖𝑗},𝜑=−Δ,𝜑−Δ_𝛼/2−1=−Δ_𝛼/2−1𝜑_𝑖𝑗,𝜑𝜑_𝑖𝑗,−Δ𝜑,(4.10)where (−Δ) is the extension of (−Δ) with domain 𝒟′(Ω) that is the same as 𝒟(Ω) without ℬ(𝜑)=0, which is well documented in books on partial differential equations [7]. This is done by calculating the conjunct (concomitant or boundary form) using the Green's formula_Ω(𝑣∇²𝑢−𝑢∇²𝑣)𝑑𝑆=_𝜕Ω𝑣𝜕𝑢𝜕𝑛−𝑢𝜕𝜕𝑣𝑛𝑑𝑠.(4.11)Thus,𝜑^𝑖𝑗,𝑇𝜑⟩==−𝜆𝜆^{𝛼−2𝑖𝑗}_{𝛼−2𝑖𝑗}𝜑^𝑖𝑗,−Δ𝜑⟩_Ω=−𝜆_{𝛼−2𝑖𝑗}𝜑_𝑖𝑗∇²𝜑𝑑𝑆_Ω𝜑∇²𝜑_𝑖𝑗𝑑𝑆+_𝜕Ω𝜑^𝑖𝑗𝜕𝜑𝜕𝑛−𝜑𝜕𝜕𝜑𝑛^𝑖𝑗𝑑𝑠(4.12)which gives the result on substitution.

This result can be readily used to write down the analytic solution to the FPE problem. First, we obtain the spectral representation of the operator 𝑇 by solving the eigenvalue problem:−Δ𝜑=𝜆²𝜑,ℬ(𝜑)=0.(4.13)Knowing the eigenvalues 𝜆_𝑖𝑗 and the corresponding orthonormal (ON) eigenfunctions 𝜑_𝑖𝑗, we can use the finite-transform method with respect to 𝜑_𝑖𝑗 and Proposition 4.3 to obtain𝜆^𝛼𝑖𝑗⟨𝜑_𝑖𝑗,𝜑⟩=⟨𝑔,𝜑_𝑖𝑗⟩+𝜆^{𝛼−2𝑖𝑗}𝑏_𝑖𝑗,(4.14)where 𝜆^{𝛼−2𝑖𝑗}𝑏_𝑖𝑗 is the second term on the right hand side in Proposition 4.3. Hence,𝜑(𝑥,𝑦)=^∞_∞𝑖=1_𝑗=1⟨𝑔,𝜑_𝑖𝑗⟩𝜆^𝛼𝑖𝑗𝜑_𝑖𝑗(𝑥,𝑦)+^∞_∞𝑖=1_𝑗=1𝑏_𝑖𝑗𝜆^2𝑖𝑗𝜑_𝑖𝑗(𝑥,𝑦).(4.15)

5. Results and Discussion

In this section, we exhibit the results of applying Algorithm 2.4 to solve two FPE test problems. To assess the accuracy of our approximation, we compare the numerical solutions with the exact solution in each case.

Test Problem 1: FPE with Dirichlet Boundary Conditions
Solve (−∇²)^𝛼/2𝜑=𝑔₀/𝑘 on the unit square [0,1]×[0,1] subject to the type I boundary conditions 𝜑=0 on boundary 𝜕Ω. For this problem, the ON eigenfunctions are given by𝜑_𝑖𝑗=2sin(𝑖𝜋𝑥)sin(𝑗𝜋𝑦)^{∞𝑖,𝑗=1},(5.1)and the corresponding eigenvalues 𝜆^2𝑖𝑗=𝜋²(𝑖²+𝑗²). The analytical solution is then given from Section 4 as𝜑(𝑥,𝑦)=16𝑔⁰𝑘𝜋^2∞_∞𝑖=0_𝑗=0sin[(2𝑖+1)𝜋𝑥]sin[(2𝑗+1)𝜋𝑦](2𝑖+1)(2𝑗+1)[(2𝑖+1)²𝜋²+(2𝑗+1)²𝜋²]^𝛼/2.(5.2) For the numerical solution, a standard five-point finite-difference formula with equal grid spacing ℎ=1/𝑛 in the 𝑥 and 𝑦 directions has been used to generate the block tridiagonal matrix 𝐴∈ℝ^{(𝑛−1)2×(𝑛−1)2} given in (2.1) as𝐴=⎡⎢⎢⎢⎢⎢⎢⎢⎣−𝐵𝐼−⋱𝐵𝐼−−⋱𝐼𝐼−⋱𝐵𝐼−𝐵𝐼⎤⎥⎥⎥⎥⎥⎥⎥⎦where𝐵=⎡⎢⎢⎢⎢⎢⎢⎢⎣−41−−411⋱−41⋱−−11−41⎤⎥⎥⎥⎥⎥⎥⎥⎦∈ℝ^{(𝑛−1)×(𝑛−1)}.(5.3)The parameters used to test this model are listed in Table 1.

Table 1 

Physical parameters for test problem 1.

Test Problem 2: FPE with Mixed Boundary Conditions
Solve (−∇²)^𝛼/2𝜑=𝑔₀/𝑘 on the unit square [0,1]×[0,1] subject to type III boundary conditions−𝜕𝜑𝜕𝑥+𝐻¹𝜕𝜑𝜑=𝐻₁𝜑_∞at𝑥=0,𝜕𝑥+𝐻²−𝜕𝜑𝜑=𝐻₂𝜑_∞at𝑥=1,𝜕𝑦+𝐻³𝜕𝜑𝜑=𝐻₃𝜑_∞at𝑦=0,𝜕𝑦+𝐻⁴𝜑=𝐻₄𝜑_∞at𝑦=1,(5.4)where 𝐻_𝑖=ℎ_𝑖/𝑘. The analytical solution to this problem is given by𝜑(𝑥,𝑦)=𝜑_∞+𝑔⁰𝑘^∞_∞𝑖=1_𝑗=1𝛼_𝑖,𝑗𝑋(𝜇_𝑖,𝑥)𝑌(𝜈_𝑗,𝑦)(𝜇^2𝑖+𝜈^2𝑗)^𝛼/2𝑁_𝑥(𝜇_𝑖)𝑁_𝑦(𝜈_𝑗),(5.5)where the eigenfunctions are𝑋(𝜇_𝑖,𝑥)=𝜇_𝑖cos(𝜇_𝑖𝑥)+𝐻₁sin(𝜇_𝑖𝑥),𝑌(𝜈_𝑗,𝑦)=𝜈_𝑗cos(𝜈_𝑗𝑦)+𝐻₃sin(𝜈_𝑗𝑦),(5.6)with normalisation factors𝑁^2𝑥(𝜇_𝑖)=12𝜇^2𝑖+𝐻²¹1+𝐻²𝑁_2𝑦𝜇^2𝑖+𝐻²²+𝐻¹,(𝜈_𝑗)=12𝜈^2𝑗+𝐻²³1+𝐻⁴𝜈^2𝑗+𝐻²⁴+𝐻³.(5.7)The eigenvalues 𝜇_𝑖 are determined by finding the roots of the transcendental equation:tan(𝜇)=𝜇(𝐻¹+𝐻₂)𝜇²−𝐻₁𝐻₂(5.8)with 𝜈_𝑗 determined from a similar equation for 𝜈. Finally, 𝛼_𝑖,𝑗 is given by𝛼_𝑖,𝑗=¹⁰𝑋(𝜇_𝑖,𝜉)𝑌(𝜈_𝑗,𝜂)𝑁_𝑥(𝜇_𝑖)𝑁_𝑦(𝜈_𝑗)d𝜉d𝜂.(5.9)
For the numerical solution, a standard five-point finite-difference formula with equal grid spacing ℎ=1/𝑛 was again employed in the 𝑥 and 𝑦 directions. However in this example, additional finite-difference equations are required for the boundary nodes as a result of type III boundary conditions. The block tridiagonal matrix required in (2.8) is then similar to that exhibited for example 1, however it has dimension 𝐴∈ℝ^{(𝑛+1)2×(𝑛+1)2} and boundary blocks must be modified to account for the boundary condition contributions.
The parameter values used for this problem are listed in Table 2.

Table 2 

Physical parameters for test problem 2.

5.1. Discussion of Results for Test Problem 1

A comparison of the numerical and analytical solutions for test problem 1 is exhibited in Figure 1 for different values of the fractional index 𝛼=0.5,1.0,1.5, and 2 (with the value 2 representing the solution of the classical Poisson equation). In all cases, it can be observed that good agreement is obtained between theory and simulation, with the analytical (solid contour lines) and numerical (dashed contour lines) solutions almost indistinguishable. In fact, Algorithm 2.4 consistently produced a numerical solution within approximately 2% absolute error of the analytical solution.

(a)

(b)

(c)

(d)

Figure 1 

Comparisons of numerical (dashed lines) and analytical solutions (solid lines) for test problem 1 computed using Algorithm 2.4: (a) 𝛼=0.5, (b) 𝛼=1.0, (c) 𝛼=1.5, and (d) 𝛼=2 (classical case).

The impact of decreasing the fractional index from 𝛼=2 to 0.5 is particularly evident in Figure 2 from the shape and magnitude of the computed three-dimensional symmetric profiles. Low values of 𝛼 produce a solution exhibiting a pronounced hump-like shape, with the diffusion rate low, the magnitude of the solution high at the centre, and steep gradients evident near the boundary of the solution domain. As 𝛼 increases, the magnitude of the profile diminishes and the solution is much more diffuse and representative of a Gaussian process. These observations motivate the following remark.

(a)

(b)

(c)

(d)

Figure 2 

Numerical solutions for test problem 1: (a) 𝛼=0.5, (b) 𝛼=1.0, (c) 𝛼=1.5, and (d) 𝛼=2 (classical case).

Remark 5.1. Over ℝ^𝑛, the Riesz operator (−Δ)^𝛼/2 as defined by−Δ^𝛼/2𝑓(𝑥)=𝑔(1𝛼)_ℝ^𝑛||𝑥−𝑦||^𝛼−𝑛𝑓(𝑦)𝑑𝑦,(5.10)where 𝑓 is a 𝐶^∞- function with rapid decay at infinity, and
𝑔(𝛼)=𝜋^𝑛/22^𝛼Γ(𝛼/2)Γ(𝑛/2−𝛼/2)(5.11)is known to generate 𝛼-stable processes. In fact, the Green's function of the equation𝜕𝑝𝜕𝑡=−−Δ^𝛼/2𝑝(𝑡,𝑥),𝑇>0,𝑥∈ℝ(5.12)is the probability density function of a symmetric 𝛼-stable process. When 𝛼=1, it is the density function of a Cauchy distribution, and when 𝛼=2 it is the classical Gaussian density. As 𝛼→0, the tail of the density function is heavier and heavier. These behaviours are reflected in the numerical results given in the above example; namely, when 𝛼→2, the plots exhibit the bell shape of the Gaussian density, but when 𝛼→0, the curves are flatter, indicating very heavy tails as expected.

We now report on the performance of Algorithm 2.4 for computing the solution of the FPE. The numerical solutions shown in Figures 1 and 2 were generated using a standard five-point finite difference stencil to construct the matrix representation of the two-dimensional Laplacian operator. The 𝑥- and 𝑦-dimensions were divided equally into 31 divisions to produce the symmetric positive definite matrix 𝐴∈ℝ^900×900 having its spectrum 𝜎(𝐴)⊂[0.0205,7.9795]. One notes for this problem that the homogeneous boundary conditions necessitate only the solution Φ=ℎ^𝛼𝐴^−𝛼/2𝑔̃. Algorithm 2.4 was still employed in this case; however in stage 1 we at first solve 𝐴𝑥=𝑔̃ by the adaptively preconditioned thick restart procedure. The gathered spectral information from stage 1 is then used for the efficient computation of Φ during stage 2.

Figure 3 depicts the reduction in the residual of the linear system and the error in the matrix function approximation during both stages of the solution process for test problem 1 for the case 𝛼=1. For this test using FOM(25,10), (subspace size 25 with an additional 10 approximate Ritz vectors augmented at the front of the Krylov subspace) four restarts were required to reduce the linear system residual to ≈1×10⁻¹⁵, which represented an overall total of 110 matrix-vector products. This low tolerance was enforced to ensure that as many approximate eigenpairs of 𝐴 could be computed and then locked during stage 1 for use in stage 2. An eigenpair was deemed converged when the residual in the approximate eigenpair was less than 𝜃_max×10⁻¹⁰, where 𝜃_max is the current estimate of the largest eigenvalue of 𝐴. This process saw 1 eigenpair locked after 2 restarts, 5 locked after 3 restarts, and finally 9 locked after 4 restarts. From this figure, we also see that when subspace recycling is used for stage 2 only an additional 30 matrix-vector products are required to compute the solution Φ to an acceptable accuracy. It is also worth pointing out that the Lanczos approximation for this preconditioned matrix function reduces much more rapidly than for the case where preconditioning (dotted line) is not used. Furthermore, the Lanczos approximation in this example lies almost entirely on the curve that represents the optimal approximation obtainable from the Krylov subspace [26]. Finally, we see that the bound (3.9) can be used with confidence as a means for halting stage 2 once the desired accuracy in the bound is reached.

Figure 3 

Residual reduction for test problem 1 computed using the two-stage process outlined in Algorithm 2.4.

5.2. Discussion of Results for Test Problem 2

A comparison of the numerical and analytical solutions for test problem 2 is exhibited in Figure 4, again for the values of the fractional index 𝛼=0.5,1.0,1.5, and 2. It can be seen that the agreement between theory and simulation is more than acceptable for this case, with Algorithm 2.4 producing a numerical solution within approximately 4% absolute error of the analytical solution. However, the impact of increasing the fractional index from 𝛼=0.5 to 2 is less dramatic for problem 2.

(a)

(b)

(c)

(d)

Figure 4 

Comparisons of numerical (dashed line) and analytical solutions (solid line) for test problem 2 computed using Algorithm 2.4: (a) 𝛼=0.5, (b) 𝛼=1.0, (c) 𝛼=1.5, and (d) 𝛼=2 (classical case).

The numerical solutions shown in Figure 4 were again generated using a standard five-point finite-difference stencil to construct the matrix representation of the two-dimensional Laplacian operator. The 𝑥- and 𝑦-dimensions were divided equally into 30 divisions resulting in the symmetric positive definite matrix 𝐴∈ℝ^961×961 having its spectrum 𝜎(𝐴)⊂[0.000758,7.9795]. One notes for this problem that type II boundary conditions have produced a small eigenvalue that undoubtedly will hinder the performance of restarted FOM.

Figure 5 depicts the reduction in the residual of the linear system for computing the solution Φ₁ and the error in the matrix function approximation for Φ₂ during both stages of the solution process for test problem 2 with 𝛼=1. Using FOM(25,10), a total of nine restarts were required to reduce the linear system residual to ≈1×10⁻¹⁵, which represented an overall total of 240 matrix-vector products. One notes that this is much higher than Problem 1 and primarily due to the occurrence of small eigenvalues in 𝜎(𝐴). The thick restart process saw 1 eigenpair locked after 5 restarts, 4 locked after 6 restarts, and finally 10 locked after 9 restarts. From this figure, we also see that when subspace recycling is used for stage 2 only an additional 25 matrix-vector products are required to compute the solution Φ₂ to an acceptable accuracy, which is clearly much less than the unpreconditioned (dotted line) case. The Lanczos approximation in this example again lies almost entirely on the curve that represents the optimal approximation obtainable from the Krylov subspace. Finally, we see that the bound (3.9) can be used to halt stage 2 once the desired accuracy is reached.

Figure 5 

Residual reduction for test problem 2 computed using the two-stage process outlined in Algorithm 2.4.

6. Conclusions

In this work, we have shown how the fractional Poisson equation can be approximately solved using a finite-difference discretisation of the Laplacian to produce an appropriate matrix representation of the operator. We then derived a matrix equation that involved both a linear system solution and a matrix function approximation with the matrix 𝐴 raised to the same fractional index as the Laplacian. We proposed an algorithm based on Krylov subspace methods that could be used to efficiently compute the solution of this matrix equation using a two-stage process. During stage 1, we used an adaptively preconditioned thick restarted FOM method to approximately solve the linear system and then used recycled spectral information gathered during this restart process to accelerate the convergence of the matrix function approximation in stage 2. Two test problems were then presented to assess the accuracy of our algorithm, and good agreement with the analytical solution was noted in both cases. Future research will see higher dimensional fractional diffusion equations solved using a similar approach via the finite volume method.

Acknowledgment

This work was supported financially by the Australian Research Council Grant no. LP0348653.