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

Authors: M. Ilić; I. W. Turner; V. Anh
Journal: Journal of Applied Mathematics and Stochastic Analysis
 (2008)

---

## Abstract

This study considers the solution of a class of linear systems related with the fractional
Poisson equation (FPE)(−∇2)α/2φ=g(x,y) 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 A raised to the fractional power α/2. The solution of the linear system then requires the action of the matrix function
f(A)=A−α/2 on a vector b. For large, sparse, and symmetric positive definite matrices, the Lanczos approximation
generates f(A)b≈β0Vmf(Tm)e1. This method works well when both the analytic grade of A with respect
to b 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.


---

## Body

## 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 H-matrices and
algorithms that they claim to be of NlogN 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 indexH∈(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:(1.1)(−∇2)α/2h(x)=g(x)ifx∈(0,T),h(x)=0ifx∉(0,T),for a certain integrable
function h(x) defined on [0,T],
where g:[0,T]→ℝ.
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:(1.2)(−∇2)α/2φ=g(x,y),0<x<a,0<y<b,subject to(1.3)−k1∂φ∂x+h1φ=f1(y)atx=0,k2∂φ∂x+h2φ=f2(y)atx=a,−k3∂φ∂y+h3φ=f3(x)aty=0,k4∂φ∂y+h4φ=f4(x)aty=b.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 (ki=0,hi=1,i=1,…,4) and type II boundary conditions (hi=0,ki=1,i=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(1.4)Aα/2Φ=bmust be solved. Note that in
(1.4), A∈ℝn×n 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 b=b1+Aα/2−1b2,
with b1∈ℝn a vector containing the discrete values of the
source/sink term, and b2∈ℝn 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 A is symmetric positive definite, that is, A∈SPD.
The general solution of (1.4) can be written
as(1.5)Φ=A−α/2b=A−α/2b1+A−1b2,and one notes the need to
determine both the action of the matrix function f(A)=A−α/2 on the vector b1 and the action of the standard inverse on b2,
where the matrix A 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 AΦ2=b2 can be recycled to obtain the complete
solution Φ=Φ1+Φ2 in (1.5), where Φ1=f(A)b1=A−α/2b1.
In literature, a majority of references deal with
the extraction of an approximation tof(A)v for scalar analytic function f(t):D⊂ℂ→ℂ 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)f(A)v≈Vmf(Tm)e1,v=∥v∥Vme1,where(1.7)AVm=VmTm+βmvm+1emTis the usual Lanczos
decomposition, and the columns of Vm form an orthonormal basis for Krylov subspace 𝒦m(A,v)={v,Av,…,Am−1v}.
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 f(t)=1/t,
and the restarted Arnoldi algorithm for computing f(A)v 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 f(A)v.
For example if AM=B,
one can only deduce f(A) from f(B) in a limited number of special cases for f(t).
In the previous work by the authors [27], we proposed a spectral splitting method f(A)v=Qf(Λ)QTv+pm(A)(I−QQT)v,
where QQT is an orthogonal projector onto the invariant
subspace associated with a set of eigenvalues on the “singular part” of the
spectrum σ(A) with respect to f(t) and I−QQT 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 pm(t) was a low-degree polynomial (of degree at most
5). Thick restarting was used to construct the projector QQT on the singular part. Unfortunately, the
computational overhead associated with constructing the projector QQT,
whilst maintaining the requirement of a low-degree polynomial approximation for f(t) 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 approximatef(A)v for monotone decreasing function f(t)=t−q.
(1)
Determine an approximately invariant subspace
(AIS),span{q1,…,qk} for the set of eigenvectors associated with
the singular part of σ(A) with respect to f(t).
Form Qk=[q1,…,qk] and set Λk=diag{λ1,λ2,…,λk},
where λi are the eigenvalues associated with the
eigenvectors qi,i=1,…,k.
The thick restarted Lanczos method discussed in [27, 29] or [30] can be used for the AIS
generation.(2)
Letv˜=(I−QkQkT)v and generate orthonormal basis for 𝒦ℓ(A,v˜).(3)
Approximatef(A)(I−QkQkT)v≈Vℓf(Tℓ)VℓTv˜ using the Lanczos decomposition to analytic
grade ℓ, AVℓ=VℓTℓ+βℓvℓ+1eℓT [31].(4)
Formf(A)v≈Qkf(Λk)QkT+Vℓf(Tℓ)VℓTv˜.
To avoid
components of any eigenvectors associated with the singular part reappearing in 𝒦ℓ(A,v˜),
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 Section2, we
use MTT to formulate the matrix representation of FPE to accommodate
nonhomogeneous boundary conditions. We also consider the approximation of the
matrix function f(A)=A−qv 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:(2.1)1h2AΦ=g˜,where it is assumed that (1/h2)A=m(−∇2) is the finite difference matrix representation of the Laplacian,
and h is the grid spacing. Φ=m(φ) is the representation of φ,
and g˜=m(g) is the representation of g.
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:(2.2)1hαAα/2Φ=g˜.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:(2.3)1h2AΦ−b=0,where b represents the discretized boundary values,
and the matrix A is the same as given above. In other words, if φ does not satisfy homogeneous boundary
conditions, then the modified representation(2.4)m((−∇2¯)φ)=1h2AΦ−bis used, where −∇2¯ 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.5)(−∇2)α/2=(−∇2)α/2−1(−∇2¯),and its matrix representation
as(2.6)m((−∇2)α/2)=m((−∇2)α/2−1)m((−∇2¯)).Hence, the matrix representation
for FPE is(2.7)1hαAα/2Φ=g˜+1hα−2Aα/2−1b.Assuming that A has an inverse, the solution of this equation
is(2.8)Φ=hαA−α/2g˜+h2A−1b.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, Φ=hαΦ1+h2Φ2,
where Φ1=A−qg˜, Φ2=A−1b,
and 0<q=α/2<1.
We note further in this context that the scalar function f(t)=t−q is monotone decreasing on σ(A),
where A∈ℝn×n is symmetric positive definite.There exists a plethora of Krylov-based methods
available in the literature for approximately solving the linear systemAΦ2=b 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 A during this linear system solve can be
recycled and used to aid the approximation of Φ1.
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 M is a good preconditioner for A,
it may not be useful for f(A) since we cannot find a relation between f(A) and f(AM−1).
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 systemAΦ2=b,
where A∈ℝn×n 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 ofA is given by(2.9)AVℓ=VℓTℓ+βℓvℓ+1eℓT=Vℓ+1T¯ℓ,where the columns of Vℓ form an orthonormal basis for 𝒦ℓ(A,b),
and ℓ is the analytic grade defined in
[31]. The analytic
grade of order t of the matrix A with respect to b is defined as the lowest integer ℓ for which ∥uℓ−𝒫ℓuℓ∥/∥uℓ∥<10−t, where 𝒫ℓ is the orthogonal projector onto the lth Krylov subspace 𝒦l and ul=Aℓb.
The grade can be computed from the Lanczos algorithm using the matrices T¯1,T¯2,…,T¯l generated during the process. If t1 is the 1st column of T¯1,
and ti=T¯iti−1,
for i=1,…,ℓ,
then ∥uℓ−𝒫ℓuℓ∥/∥uℓ∥=|el+1Ttl|/∥tl∥.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 Φ2(0) be some initial approximation to the solution Φ2 and define r0=b−AΦ2(0).Cycle 1
(i)
Generate Lanczos decomposition(2.10)AV1=V1T1+β1u1eT,where V1=[v1(1),…,vℓ(1)], v1(1)=r0/β0, T1 is tridiagonal, u1=vℓ+1(1),  β0=∥r0∥,  β1=βℓ(1),
and eT=eℓT.(ii)
Obtain
approximate solutionΦ˜2(1)=V1T1−1V1Tr0,
so that(2.11)Φ2(1)=Φ2(0)+V1T1−1V1Tr0,and residual(2.12)r1=b−AΦ2(1)=r0−AΦ˜2(1)=−β1(eTT1−1V1Tr0)u1.Test if ∥r1∥<ε.
If yes, stop; otherwise, continue to cycle 2.Cycle 2
(i)
Find
eigenvalue decomposition ofT1,
that is, T1Y=YΛ,
where Λ=diag{θ1,…,θℓ}.(ii)
Select
thek orthonormal (ON) eigenvectors, Y1,
of T1 corresponding to the k smallest in magnitude eigenvalues of T1 and form the Ritz vectors(2.13)W1=V1Y1=[w1,…,wk],where wi are ON, and let the associated Ritz values be
stored in the diagonal matrix Λ1=diag{θ1,…,θk}.(iii)
SetV˜2=[W1,u1] and generate the thick-restart
Lanczos decomposition(2.14)AV2=V2T2+β2u2eT,where V2=[w1,…,wk,v1(2),…,vℓ(2)], v1(2)=u1, u2=vℓ+1(2),
and(2.15)T2=[Λ1β2s10⋯β2s1Tαk+1βk+10⋯0βk+1⋱⋱…⋱…],withs1=Y1Te.(iv)
Obtain
approximate solutionΦ˜2(2)=V2T2−1V2Tr1,
so that(2.16)Φ2(2)=Φ2(1)+V2T2−1V2Tr1,and residual(2.17)r2=b−AΦ2(2)=r1−AΦ˜2(2)=−β2(eTT2−1V2Tr1)u2.Test if ∥r2∥<ε.
If yes, stop; otherwise, continue to the next cycle.Cycle(j+1)
(i)
Find
eigenvalue decomposition ofTj,
that is, TjY=YΛ.(ii)
Selectk orthonormal (ON) eigenvectors, Yj,
of Tj corresponding to the k smallest in magnitude eigenvalues of Tj and form the Ritz vectors Wj=VjYj.(iii)
SetV˜j+1=[Wj,uj] and generate thick-restart Lanczos
decomposition(2.18)AVj+1=Vj+1Tj+1+βj+1uj+1eT,where Tj+1 has similar form as T2.(iv)
Obtain
approximate solutionΦ˜2(j+1)=Vj+1Tj+1−1Vj+1Trj,
so that(2.19)Φ2(j+1)=Φ2(j)+Φ˜2(j+1)=Φ2(0)+∑i=1j+1ViTi−1ViTri−1,and residual(2.20)rj+1=b−AΦ2(j+1)=−βj+1(eTTj+1−1Vj+1Trj)uj+1.Test if ∥rj+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 j produce the k approximate smallest Ritz pairs {θi,wi}i=1k,
where wi=Vjyi.
We then check if any of these Ritz pairs have converged to approximate
eigenpairs of A by testing the magnitude of the upper bound on
the eigenpair residual(2.21)∥Awi−θiwi∥≤βj|eTyi|<ε2.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 thep locked Ritz vectors as columns of the matrix Qj=[q1,q2,…,qp],
set Λj=diag{θ1,…,θp},
and form(2.22)Mj−1=γQjΛj−1QjT+I−QjQjT,where γ=(θmin+θmax)/2. θmin, θmax are the current estimates of the smallest and
largest eigenvalues of A,
respectively, obtained from the restart process. Then, Aj=AMj−1 has the same eigenvectors as A;
however its eigenvalues {λi}i=1p are shifted to γ [33, 34]. Furthermore, it should be noted that these
preconditioners can be nested. If M1,M2,…,Mj is a sequence of such preconditioners, then
with Q=[Q1,Q2,…,Qj] and Λ=diag(Λi,i=1,…,j),
we have(2.23)M−1=Mj−1⋯M2−1M1−1=γQΛ−1QT+I−QQT.Thus, during the cycles (say
cycle j+1) the adaptively preconditioned,
thick- restart Lanczos decomposition(2.24)AM−1Vj+1=Vj+1Tj+1+βj+1uj+1eTis employed.Note.
The preconditionerM−1 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 systemAΦ2=b for Φ2.
At the completion of this process, the residual ∥r∥=∥b−AΦ2∥<ε,
and we have the set {θi,qi}i=1k 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{q1,q2,…,qk} associated with the smallest in magnitude
eigenvalues of A.
We now show how this spectral information can be recycled to aid with the
approximation of Φ1=f(A)g˜,
where f(t)=t−q.


#### 2.2.1. Adaptive Preconditioning

Recall from
stage 1 that we have availableM−1=γQkΛk−1QkT+I−QkQkT,
where Qk=[q1,…,qk].
The important observation at this point is the following relationship between f(A) and f(AM−1).Proposition 2.1.
Letspan{q1,q2,…,qk} be an eigenspace of symmetric matrix A such that AQk=QkΛk,
with Qk=[q1,q2,…,qk] and Λk=diag(μ1,…,μk).
Define M=(1/γ)QkΛkQkT+I−QkQkT,
then for v∈ℝn,(2.25)f(A)v=1f(γ)f(AM−1)f(γM)v.Proof.
LetWWT=I−QkQkT,WTAW=B,
then MM−1=QkQkT+WWT=I=M−1M. Furthermore,
(2.26)M−1A=(QkW)(γI00B)(QkTWT)=AM−1.Thus,(2.27)f(AM−1)=(QkW)(f(γ)I00f(B))(QkTWT).By noting that(2.28)f(A)=(QkW)(f(Λk)00f(B))(QkTWT),f(γM)=(QkW)(f(Λk)00f(γ)I)(QkTWT),we obtain the main result(2.29)f(A)f(γ)[f(γM)]−1=(QkW)(f(γ)I00f(B))(QkTWT)=f(AM−1).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.
LetM1,M2,…,Mj be a sequence of preconditioners as defined in
Proposition 2.1, then(2.30)f(A)v=1f(γ)f(AM1−1M2−1⋯Mj−1)f(γM1M2⋯Mj)v.Proof.
LetQ=[Q1,Q2,…,Qk] and Λ=diag(Λi,i=1,…,j),
then observe that M=M1M2⋯Mj=(1/γ)QΛQT+I−QQT and f(A)=f(AM−1)(1/f(γ))f(γM).Corollary 2.3.
Under the hypothesis of
Proposition2.1, one notes the equivalent form of
(2.25) as(2.31)f(A)v=Qkf(Λk)QkTv+f(AM−1)(I−QkQkT)v,which appears similar to the
idea of spectral splitting proposed in [27].We now turn our attention to the approximation ofΦ1=A−qg˜,
which by using Corollary 2.3 can be expressed as(2.32)A−qg˜=∑i=1kθi−qqiqiTg˜+(AM−1)−qg^,where g^=(I−QkQkT)g˜.
First note that if A∈SPD,
then AM−1∈SPD.
We expand the Lanczos decomposition AM−1Vℓ=VℓTℓ+βℓvℓ+1eℓT to the analytic grade ℓ of AM−1 with v1=g^/∥g^∥.
Next perform the spectral decomposition of Tℓ=YℓΛℓYℓT and set Q˜ℓ=VℓYℓ,
then compute the Lanczos approximation(2.33)(AM−1)−qg^≈VℓTℓ−qVℓTg^=Q˜ℓΛℓ−qQ˜ℓTg^.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.
SolveAΦ2=b  using the thick restarted adaptively
preconditioned Lanczos method  and generate the AIS, 𝒬k=span{q1,…,qk}.
Return the preconditioner M=(1/γ)QkΛkQkT+I−QkQkT,  where Qk=[q1,…,qk].
Stage 2.
ComputeΦ1=A−qg˜ using the following strategy.
(1)
Setg^=(I−QkQkT)g˜.(2)
Compute Lanczos decompositionAM−1Vℓ=VℓTℓ+βℓvℓ+1eℓT,
where ℓ  is the
analytic grade of AM−1 and Vℓ=[v1,…,vℓ],
with v1=g^/∥g^∥.(3)
Perform the spectral decompositionTℓ=YℓΛℓYℓT.(4)
Compute linear system residual∥rℓ∥=|βℓeℓTYℓΛℓ−1YℓTVℓTg^| and estimate λmin≈μmin from Tℓ  to compute bound (3.9) μmin−q∥rℓ∥  derived in
Section 3.(5)
If bound is small, then approximatef(AM−1)g^≈VℓTℓ−qVℓTg^ and exit to step (6),
otherwise continue the
Lanczos expansion until bound is satisfied.(6)
FormΦ1=f(A)g˜≈QkΛk−qQkTg˜+Q˜ℓΛℓ−qQ˜ℓTg^,
where Q˜ℓ=VℓYℓ.
Finally,
compose the approximate solution of FPE asΦ=hαΦ1+h2Φ2.Remarks
At stage 2, we monitor the upper bound given in Proposition3.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 k 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 (AM−1)−qg^ at the outset makes it impossible to answer
this question. Instead, we opt to provide an upper bound for the error ∥(AM−1)−qg^−VℓTℓ−qVℓTg^∥,
which is the topic of the following section.


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

Suppose that
the Lanczos decomposition ofA is given by(2.9)AVℓ=VℓTℓ+βℓvℓ+1eℓT=Vℓ+1T¯ℓ,where the columns of Vℓ form an orthonormal basis for 𝒦ℓ(A,b),
and ℓ is the analytic grade defined in
[31]. The analytic
grade of order t of the matrix A with respect to b is defined as the lowest integer ℓ for which ∥uℓ−𝒫ℓuℓ∥/∥uℓ∥<10−t, where 𝒫ℓ is the orthogonal projector onto the lth Krylov subspace 𝒦l and ul=Aℓb.
The grade can be computed from the Lanczos algorithm using the matrices T¯1,T¯2,…,T¯l generated during the process. If t1 is the 1st column of T¯1,
and ti=T¯iti−1,
for i=1,…,ℓ,
then ∥uℓ−𝒫ℓuℓ∥/∥uℓ∥=|el+1Ttl|/∥tl∥.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 Φ2(0) be some initial approximation to the solution Φ2 and define r0=b−AΦ2(0).Cycle 1
(i)
Generate Lanczos decomposition(2.10)AV1=V1T1+β1u1eT,where V1=[v1(1),…,vℓ(1)], v1(1)=r0/β0, T1 is tridiagonal, u1=vℓ+1(1),  β0=∥r0∥,  β1=βℓ(1),
and eT=eℓT.(ii)
Obtain
approximate solutionΦ˜2(1)=V1T1−1V1Tr0,
so that(2.11)Φ2(1)=Φ2(0)+V1T1−1V1Tr0,and residual(2.12)r1=b−AΦ2(1)=r0−AΦ˜2(1)=−β1(eTT1−1V1Tr0)u1.Test if ∥r1∥<ε.
If yes, stop; otherwise, continue to cycle 2.Cycle 2
(i)
Find
eigenvalue decomposition ofT1,
that is, T1Y=YΛ,
where Λ=diag{θ1,…,θℓ}.(ii)
Select
thek orthonormal (ON) eigenvectors, Y1,
of T1 corresponding to the k smallest in magnitude eigenvalues of T1 and form the Ritz vectors(2.13)W1=V1Y1=[w1,…,wk],where wi are ON, and let the associated Ritz values be
stored in the diagonal matrix Λ1=diag{θ1,…,θk}.(iii)
SetV˜2=[W1,u1] and generate the thick-restart
Lanczos decomposition(2.14)AV2=V2T2+β2u2eT,where V2=[w1,…,wk,v1(2),…,vℓ(2)], v1(2)=u1, u2=vℓ+1(2),
and(2.15)T2=[Λ1β2s10⋯β2s1Tαk+1βk+10⋯0βk+1⋱⋱…⋱…],withs1=Y1Te.(iv)
Obtain
approximate solutionΦ˜2(2)=V2T2−1V2Tr1,
so that(2.16)Φ2(2)=Φ2(1)+V2T2−1V2Tr1,and residual(2.17)r2=b−AΦ2(2)=r1−AΦ˜2(2)=−β2(eTT2−1V2Tr1)u2.Test if ∥r2∥<ε.
If yes, stop; otherwise, continue to the next cycle.Cycle(j+1)
(i)
Find
eigenvalue decomposition ofTj,
that is, TjY=YΛ.(ii)
Selectk orthonormal (ON) eigenvectors, Yj,
of Tj corresponding to the k smallest in magnitude eigenvalues of Tj and form the Ritz vectors Wj=VjYj.(iii)
SetV˜j+1=[Wj,uj] and generate thick-restart Lanczos
decomposition(2.18)AVj+1=Vj+1Tj+1+βj+1uj+1eT,where Tj+1 has similar form as T2.(iv)
Obtain
approximate solutionΦ˜2(j+1)=Vj+1Tj+1−1Vj+1Trj,
so that(2.19)Φ2(j+1)=Φ2(j)+Φ˜2(j+1)=Φ2(0)+∑i=1j+1ViTi−1ViTri−1,and residual(2.20)rj+1=b−AΦ2(j+1)=−βj+1(eTTj+1−1Vj+1Trj)uj+1.Test if ∥rj+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 j produce the k approximate smallest Ritz pairs {θi,wi}i=1k,
where wi=Vjyi.
We then check if any of these Ritz pairs have converged to approximate
eigenpairs of A by testing the magnitude of the upper bound on
the eigenpair residual(2.21)∥Awi−θiwi∥≤βj|eTyi|<ε2.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 thep locked Ritz vectors as columns of the matrix Qj=[q1,q2,…,qp],
set Λj=diag{θ1,…,θp},
and form(2.22)Mj−1=γQjΛj−1QjT+I−QjQjT,where γ=(θmin+θmax)/2. θmin, θmax are the current estimates of the smallest and
largest eigenvalues of A,
respectively, obtained from the restart process. Then, Aj=AMj−1 has the same eigenvectors as A;
however its eigenvalues {λi}i=1p are shifted to γ [33, 34]. Furthermore, it should be noted that these
preconditioners can be nested. If M1,M2,…,Mj is a sequence of such preconditioners, then
with Q=[Q1,Q2,…,Qj] and Λ=diag(Λi,i=1,…,j),
we have(2.23)M−1=Mj−1⋯M2−1M1−1=γQΛ−1QT+I−QQT.Thus, during the cycles (say
cycle j+1) the adaptively preconditioned,
thick- restart Lanczos decomposition(2.24)AM−1Vj+1=Vj+1Tj+1+βj+1uj+1eTis employed.Note.
The preconditionerM−1 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 systemAΦ2=b for Φ2.
At the completion of this process, the residual ∥r∥=∥b−AΦ2∥<ε,
and we have the set {θi,qi}i=1k of locked Ritz pairs. This spectral
information is then passed to accelerate the performance of stage 2 of the
solution process.


## 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 j produce the k approximate smallest Ritz pairs {θi,wi}i=1k,
where wi=Vjyi.
We then check if any of these Ritz pairs have converged to approximate
eigenpairs of A by testing the magnitude of the upper bound on
the eigenpair residual(2.21)∥Awi−θiwi∥≤βj|eTyi|<ε2.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 thep locked Ritz vectors as columns of the matrix Qj=[q1,q2,…,qp],
set Λj=diag{θ1,…,θp},
and form(2.22)Mj−1=γQjΛj−1QjT+I−QjQjT,where γ=(θmin+θmax)/2. θmin, θmax are the current estimates of the smallest and
largest eigenvalues of A,
respectively, obtained from the restart process. Then, Aj=AMj−1 has the same eigenvectors as A;
however its eigenvalues {λi}i=1p are shifted to γ [33, 34]. Furthermore, it should be noted that these
preconditioners can be nested. If M1,M2,…,Mj is a sequence of such preconditioners, then
with Q=[Q1,Q2,…,Qj] and Λ=diag(Λi,i=1,…,j),
we have(2.23)M−1=Mj−1⋯M2−1M1−1=γQΛ−1QT+I−QQT.Thus, during the cycles (say
cycle j+1) the adaptively preconditioned,
thick- restart Lanczos decomposition(2.24)AM−1Vj+1=Vj+1Tj+1+βj+1uj+1eTis employed.Note.
The preconditionerM−1 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 systemAΦ2=b for Φ2.
At the completion of this process, the residual ∥r∥=∥b−AΦ2∥<ε,
and we have the set {θi,qi}i=1k 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{q1,q2,…,qk} associated with the smallest in magnitude
eigenvalues of A.
We now show how this spectral information can be recycled to aid with the
approximation of Φ1=f(A)g˜,
where f(t)=t−q.


### 2.2.1. Adaptive Preconditioning

Recall from
stage 1 that we have availableM−1=γQkΛk−1QkT+I−QkQkT,
where Qk=[q1,…,qk].
The important observation at this point is the following relationship between f(A) and f(AM−1).Proposition 2.1.
Letspan{q1,q2,…,qk} be an eigenspace of symmetric matrix A such that AQk=QkΛk,
with Qk=[q1,q2,…,qk] and Λk=diag(μ1,…,μk).
Define M=(1/γ)QkΛkQkT+I−QkQkT,
then for v∈ℝn,(2.25)f(A)v=1f(γ)f(AM−1)f(γM)v.Proof.
LetWWT=I−QkQkT,WTAW=B,
then MM−1=QkQkT+WWT=I=M−1M. Furthermore,
(2.26)M−1A=(QkW)(γI00B)(QkTWT)=AM−1.Thus,(2.27)f(AM−1)=(QkW)(f(γ)I00f(B))(QkTWT).By noting that(2.28)f(A)=(QkW)(f(Λk)00f(B))(QkTWT),f(γM)=(QkW)(f(Λk)00f(γ)I)(QkTWT),we obtain the main result(2.29)f(A)f(γ)[f(γM)]−1=(QkW)(f(γ)I00f(B))(QkTWT)=f(AM−1).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.
LetM1,M2,…,Mj be a sequence of preconditioners as defined in
Proposition 2.1, then(2.30)f(A)v=1f(γ)f(AM1−1M2−1⋯Mj−1)f(γM1M2⋯Mj)v.Proof.
LetQ=[Q1,Q2,…,Qk] and Λ=diag(Λi,i=1,…,j),
then observe that M=M1M2⋯Mj=(1/γ)QΛQT+I−QQT and f(A)=f(AM−1)(1/f(γ))f(γM).Corollary 2.3.
Under the hypothesis of
Proposition2.1, one notes the equivalent form of
(2.25) as(2.31)f(A)v=Qkf(Λk)QkTv+f(AM−1)(I−QkQkT)v,which appears similar to the
idea of spectral splitting proposed in [27].We now turn our attention to the approximation ofΦ1=A−qg˜,
which by using Corollary 2.3 can be expressed as(2.32)A−qg˜=∑i=1kθi−qqiqiTg˜+(AM−1)−qg^,where g^=(I−QkQkT)g˜.
First note that if A∈SPD,
then AM−1∈SPD.
We expand the Lanczos decomposition AM−1Vℓ=VℓTℓ+βℓvℓ+1eℓT to the analytic grade ℓ of AM−1 with v1=g^/∥g^∥.
Next perform the spectral decomposition of Tℓ=YℓΛℓYℓT and set Q˜ℓ=VℓYℓ,
then compute the Lanczos approximation(2.33)(AM−1)−qg^≈VℓTℓ−qVℓTg^=Q˜ℓΛℓ−qQ˜ℓTg^.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.
SolveAΦ2=b  using the thick restarted adaptively
preconditioned Lanczos method  and generate the AIS, 𝒬k=span{q1,…,qk}.
Return the preconditioner M=(1/γ)QkΛkQkT+I−QkQkT,  where Qk=[q1,…,qk].
Stage 2.
ComputeΦ1=A−qg˜ using the following strategy.
(1)
Setg^=(I−QkQkT)g˜.(2)
Compute Lanczos decompositionAM−1Vℓ=VℓTℓ+βℓvℓ+1eℓT,
where ℓ  is the
analytic grade of AM−1 and Vℓ=[v1,…,vℓ],
with v1=g^/∥g^∥.(3)
Perform the spectral decompositionTℓ=YℓΛℓYℓT.(4)
Compute linear system residual∥rℓ∥=|βℓeℓTYℓΛℓ−1YℓTVℓTg^| and estimate λmin≈μmin from Tℓ  to compute bound (3.9) μmin−q∥rℓ∥  derived in
Section 3.(5)
If bound is small, then approximatef(AM−1)g^≈VℓTℓ−qVℓTg^ and exit to step (6),
otherwise continue the
Lanczos expansion until bound is satisfied.(6)
FormΦ1=f(A)g˜≈QkΛk−qQkTg˜+Q˜ℓΛℓ−qQ˜ℓTg^,
where Q˜ℓ=VℓYℓ.
Finally,
compose the approximate solution of FPE asΦ=hαΦ1+h2Φ2.Remarks
At stage 2, we monitor the upper bound given in Proposition3.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 k 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 (AM−1)−qg^ at the outset makes it impossible to answer
this question. Instead, we opt to provide an upper bound for the error ∥(AM−1)−qg^−VℓTℓ−qVℓTg^∥,
which is the topic of the following section.


## 2.2.1. Adaptive Preconditioning

Recall from
stage 1 that we have availableM−1=γQkΛk−1QkT+I−QkQkT,
where Qk=[q1,…,qk].
The important observation at this point is the following relationship between f(A) and f(AM−1).Proposition 2.1.
Letspan{q1,q2,…,qk} be an eigenspace of symmetric matrix A such that AQk=QkΛk,
with Qk=[q1,q2,…,qk] and Λk=diag(μ1,…,μk).
Define M=(1/γ)QkΛkQkT+I−QkQkT,
then for v∈ℝn,(2.25)f(A)v=1f(γ)f(AM−1)f(γM)v.Proof.
LetWWT=I−QkQkT,WTAW=B,
then MM−1=QkQkT+WWT=I=M−1M. Furthermore,
(2.26)M−1A=(QkW)(γI00B)(QkTWT)=AM−1.Thus,(2.27)f(AM−1)=(QkW)(f(γ)I00f(B))(QkTWT).By noting that(2.28)f(A)=(QkW)(f(Λk)00f(B))(QkTWT),f(γM)=(QkW)(f(Λk)00f(γ)I)(QkTWT),we obtain the main result(2.29)f(A)f(γ)[f(γM)]−1=(QkW)(f(γ)I00f(B))(QkTWT)=f(AM−1).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.
LetM1,M2,…,Mj be a sequence of preconditioners as defined in
Proposition 2.1, then(2.30)f(A)v=1f(γ)f(AM1−1M2−1⋯Mj−1)f(γM1M2⋯Mj)v.Proof.
LetQ=[Q1,Q2,…,Qk] and Λ=diag(Λi,i=1,…,j),
then observe that M=M1M2⋯Mj=(1/γ)QΛQT+I−QQT and f(A)=f(AM−1)(1/f(γ))f(γM).Corollary 2.3.
Under the hypothesis of
Proposition2.1, one notes the equivalent form of
(2.25) as(2.31)f(A)v=Qkf(Λk)QkTv+f(AM−1)(I−QkQkT)v,which appears similar to the
idea of spectral splitting proposed in [27].We now turn our attention to the approximation ofΦ1=A−qg˜,
which by using Corollary 2.3 can be expressed as(2.32)A−qg˜=∑i=1kθi−qqiqiTg˜+(AM−1)−qg^,where g^=(I−QkQkT)g˜.
First note that if A∈SPD,
then AM−1∈SPD.
We expand the Lanczos decomposition AM−1Vℓ=VℓTℓ+βℓvℓ+1eℓT to the analytic grade ℓ of AM−1 with v1=g^/∥g^∥.
Next perform the spectral decomposition of Tℓ=YℓΛℓYℓT and set Q˜ℓ=VℓYℓ,
then compute the Lanczos approximation(2.33)(AM−1)−qg^≈VℓTℓ−qVℓTg^=Q˜ℓΛℓ−qQ˜ℓTg^.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.
SolveAΦ2=b  using the thick restarted adaptively
preconditioned Lanczos method  and generate the AIS, 𝒬k=span{q1,…,qk}.
Return the preconditioner M=(1/γ)QkΛkQkT+I−QkQkT,  where Qk=[q1,…,qk].
Stage 2.
ComputeΦ1=A−qg˜ using the following strategy.
(1)
Setg^=(I−QkQkT)g˜.(2)
Compute Lanczos decompositionAM−1Vℓ=VℓTℓ+βℓvℓ+1eℓT,
where ℓ  is the
analytic grade of AM−1 and Vℓ=[v1,…,vℓ],
with v1=g^/∥g^∥.(3)
Perform the spectral decompositionTℓ=YℓΛℓYℓT.(4)
Compute linear system residual∥rℓ∥=|βℓeℓTYℓΛℓ−1YℓTVℓTg^| and estimate λmin≈μmin from Tℓ  to compute bound (3.9) μmin−q∥rℓ∥  derived in
Section 3.(5)
If bound is small, then approximatef(AM−1)g^≈VℓTℓ−qVℓTg^ and exit to step (6),
otherwise continue the
Lanczos expansion until bound is satisfied.(6)
FormΦ1=f(A)g˜≈QkΛk−qQkTg˜+Q˜ℓΛℓ−qQ˜ℓTg^,
where Q˜ℓ=VℓYℓ.
Finally,
compose the approximate solution of FPE asΦ=hαΦ1+h2Φ2.Remarks
At stage 2, we monitor the upper bound given in Proposition3.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 k 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 (AM−1)−qg^ at the outset makes it impossible to answer
this question. Instead, we opt to provide an upper bound for the error ∥(AM−1)−qg^−VℓTℓ−qVℓTg^∥,
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:(3.1)∫0∞x−qx+1dx=πsin(qπ).By changing the variable, one
can deduce the following expression, for λ−q,λ>0:(3.2)λ−q=sin(qπ)(1−q)π∫0∞dtt1/(1−q)+λ.Noting that A¯=AM−1∈SPD,
the spectral decomposition and the usual definition of the matrix function
enable the following expression for computing A¯−q to be obtained:(3.3)A¯−q=sin(qπ)(1−q)π∫0∞(t1/(1−q)I+A¯)−1dt.Recall that the approximate solution of the linear
systemA¯x=v from 𝒦ℓ(A¯,v) using the Galerkin approach (FOM or CG) is
given by xℓ=VℓTℓ−1VℓTv,
with residual rℓ=b−A¯xℓ=−(βℓeℓTTℓ−1VℓTv)vℓ+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.
Letrℓ(t)=v−(A¯+t1/(1−q)I)Vℓ(Tℓ+t1/(1−q)I)−1VℓTv be the residual to the shifted linear system (A¯+t1/(1−q)I)x=v,
then(3.4)A¯−qv−VℓTℓ−qVℓTv=sin(qπ)(1−q)π∫0∞(t1/(1−q)I+A¯)−1rℓ(t)dt.Proof.
It is known that(3.5)A¯−qv−VℓTℓ−qVℓTv=sin(qπ)(1−q)π∫0∞{(t1/(1−q)I+A¯)−1−Vℓ(t1/(1−q)I+Tℓ)−1VℓT}vdt=sin(qπ)(1−q)π∫0∞(t1/(1−q)I+A¯)−1{v−(t1/(1−q)I+A¯)Vℓ(t1/(1−q)I+Tℓ)−1VℓTv}dt.It is interesting to observe thatrℓ(t)=−(βℓeℓT(Tℓ+t1/(1−q)I)−1VℓTv)vℓ+1 for the Lanczos approximation, so that the
vectors rℓ≡rℓ(0) and rℓ(t) are aligned; however their magnitudes are
different. Note further that(3.6)rℓ(t)=eℓT(Tℓ+tI)−1e1eℓTTℓ−1e1rℓ(0).An even more important result is
the following relationship between their norms.Proposition 3.2.
LetTℓ have eigendecomposition TℓYℓ=YℓΛℓ,
where Λℓ=diag(μi,i=1,…,ℓ) with the μi Ritz values for the Lanczos approximation,
then for t>0,(3.7)∥rℓ(t)∥=|∏i=1ℓμiμi+t1/(1−q)|∥rℓ∥≤∥rℓ∥.Proof.
The result follows from [26], which gives the
following polynomial characterisations for the residuals:(3.8)rℓ=πℓ(A)vπℓ(0),πℓ(τ)=det(τI−Tℓ)=∏i=1ℓ(τ−μi),rℓ(t)=πℓt(A+t1/(1−q)I)vπℓt(0),πℓt(τ)=det(τI−(Tℓ+t1/(1−q)I))=∏i=1ℓ(τ−(μi+t1/(1−q))),so that rℓ(t)=πℓ(A)v/πℓt(0)=(πℓ(0)/πℓt(0))rℓ=∏i=1ℓ(μi/(μi+t1/(1−q)))rℓ.
The result follows by taking the norm and noting that t>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 A¯ and rℓ the linear system residual obtained by solving
the linear system A¯x=g^ using FOM on the Krylov subspace 𝒦ℓ(A¯,g^),
then for 0<q<1,
one has(3.9)∥A¯−qg^−VℓTℓ−qVℓTg^∥≤λmin−q∥rℓ∥.Proof.
Using the orthogonal diagonalisationA¯=QΛQT,
we obtain from Proposition 3.1 that(3.10)A¯−qg^−VℓTℓ−qVℓTg^=sin(qπ)(1−q)π∫0∞Q(t1/(1−q)I+Λ)−1QTrℓ(t)dt.The result follows by taking
norms and using Proposition 3.2 to obtain(3.11)∥A¯−qg^−VℓTℓ−qVℓTg^∥≤sin(qπ)(1−q)π∫0∞∥diag{1(t1/(1−q)I+λi),i=1,…,n}∥dt∥rℓ∥.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 Proposition3.3 is that it
motivates us to shift the small troublesome eigenvalues of A,
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 Section2. The theory depends on the definition of the operator (−∇2)α/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Ω={(x,y)∣0<x<a,0<y<b},
and let ℋ be the real space ℒ2(Ω) with real inner product 〈u,v〉=∫∫ΩuvdS.
Consider the operator(4.1)Tφ=−(∂2∂x2+∂2∂y2)φ=(−Δ)φon 𝒟={φ∈ℋ∣φx,φy absolutely continuous; φx,φy,φxx,φxy,φyy∈ℒ2(Ω),ℬ(φ)=0},
where ℬ(φ) is one of the boundary conditions in the FPE
problem with right-hand side equal to zero.It is known thatT is a closed self-adjoint operator whose
eigenfunctions {φij}i,j=1∞ form an orthonormal basis for ℋ.
Thus, Tφij=λij2φij,i,j=1,2,….
For any φ∈𝒟,(4.2)φ=∑i=1∞∑j=1∞cijφij,cij=〈φ,φij〉,Tφ=∑i=1∞∑j=1∞λij2cijφij.If ψ is a continuous function on ℝ,
then(4.3)ψ(T)φ=∑i=1∞∑j=1∞ψ(λij2)cijφij,provided that ∑i=1∞∑j=1∞|ψ(λij2)cij|2<∞.
Hence, if the eigenvalue problem for T can be solved for the region Ω,
then the FPE problem with homogeneous boundary conditions can be easily solved
to give(4.4)φ(x,y)=∑i=1∞∑j=1∞〈g,φij〉λijαφij(x,y).


### 4.2. Nonhomogeneous Boundary Conditions

Before we
proceed further, we need to specify the definition of−(∇2)α/2.Definition 4.1.
Let{φij} be a complete set of orthonormal
eigenfunctions corresponding to eigenvalues λij2(≠0) of the Laplacian (−Δ) on a bounded region Ω with homogeneous BCs on ∂Ω.
Let(4.5)𝔉γ={f∈𝒟(Ω)∣∑i=1∞∑j=1∞|λijγ||cij|2<∞,cij=〈f,φij〉,γ=max(α,0)}.Then, for any f∈𝔉γ, (−Δ)α/2f is defined by(4.6)(−Δ)α/2f=∑i=1∞∑j=1∞λijαcijφij.If one of λij2=0 and φ0 is the eigenfunction corresponding to this
eigenvalue, then one needs 〈f,φ0〉=0.Proposition 4.2.
(1)
The operatorT=(−Δ)α/2 is linear and self-adjoint, that is, for f,g∈𝔉γ, 〈Tf,g〉=〈f,Tg〉.(2)
Iff∈𝔉γ,
where γ=max(0,α,β,α+β),
then(4.7)(−Δ)α/2(−Δ)β/2f=(−Δ)(α+β)/2f=(−Δ)β/2(−Δ)α/2f.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 Definition4.1
(1)
Forα=2m,m=0,1,2,…,
define (−Δ)α/2f=(−Δ)mf for any f∈C2m(Ω) (or other possibilities).(2)
Form−1<α/2<m,m=1,2,…,
define (−Δ)α/2f=T¯g,
where g=(−Δ)m−1f∈C2(m−1)(Ω),
and T¯ is the extension of T=(−Δ)α/2+1−m as defined by Proposition 4.3 below.
It suffices to
consider0<α<2.Proposition 4.3.
Letφij(x,y) be an eigenfunction corresponding to the
eigenvalue λij2(≠0) of the Laplacian (−Δ) on the rectangle Ω,
and let φ satisfy the BCs in problem 1. Then, if T¯ is an extension of T=(−Δ)α/2 (in symbols T⊂T¯) with adjoint T¯*⊂T*,(4.8)〈φij,T¯φ〉=λijα〈φij,φ〉−λijα−2{∫0bφij(0,y)k1f1(y)dy+∫0bφij(a,y)k2f2(y)dy+∫0aφij(x,0)k3f3(x)dx+∫0aφij(x,b)k4f4(x)dx}if ki≠0,i=1,…,4.
If ki=0,
the second term on the right-hand side becomes(4.9)−λijα−2{∫0b∂φij(0,y)∂xf1(y)h1dy−∫0b∂φij(a,y)∂xf2(y)h2dy+∫0a∂φij(x,0)∂yf3(x)h3dx−∫0a∂φij(x,b)∂yf4(y)h4dx}.Proof.
It is known that(4.10)〈φij,T¯φ〉=〈T¯*φij,φ〉=〈T*φij,φ〉=〈Tφij,φ〉=〈(−Δ)(−Δ)α/2−1φij,φ〉=〈(−Δ)α/2−1φij,(−Δ)¯φ〉,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)∫∫Ω(v∇2u−u∇2v)dS=∫∂Ω(v∂u∂n−u∂v∂n)ds.Thus,(4.12)〈φij,T¯φ〉=λijα−2〈φij,(−Δ¯)φ〉=−λijα−2∫∫Ωφij∇2φdS=−λijα−2{∫∫Ωφ∇2φijdS+∫∂Ω(φij∂φ∂n−φ∂φij∂n)ds}which gives the result on substitution.This result can be readily used to write down the
analytic solution to theFPE problem. First, we obtain the spectral
representation of the operator T by solving the eigenvalue
problem:(4.13)−Δφ=λ2φ,ℬ(φ)=0.Knowing the eigenvalues λij and the corresponding orthonormal (ON)
eigenfunctions φij,
we can use the finite-transform method with respect to φij and Proposition 4.3 to obtain(4.14)λijα〈φij,φ〉=〈g,φij〉+λijα−2bij,where λijα−2bij is the second term on the right hand side in
Proposition 4.3. Hence,(4.15)φ(x,y)=∑i=1∞∑j=1∞〈g,φij〉λijαφij(x,y)+∑i=1∞∑j=1∞bijλij2φij(x,y).


## 4.1. Homogeneous Boundary Conditions

In operator
theory, functions of operators are defined using spectral decomposition. SetΩ={(x,y)∣0<x<a,0<y<b},
and let ℋ be the real space ℒ2(Ω) with real inner product 〈u,v〉=∫∫ΩuvdS.
Consider the operator(4.1)Tφ=−(∂2∂x2+∂2∂y2)φ=(−Δ)φon 𝒟={φ∈ℋ∣φx,φy absolutely continuous; φx,φy,φxx,φxy,φyy∈ℒ2(Ω),ℬ(φ)=0},
where ℬ(φ) is one of the boundary conditions in the FPE
problem with right-hand side equal to zero.It is known thatT is a closed self-adjoint operator whose
eigenfunctions {φij}i,j=1∞ form an orthonormal basis for ℋ.
Thus, Tφij=λij2φij,i,j=1,2,….
For any φ∈𝒟,(4.2)φ=∑i=1∞∑j=1∞cijφij,cij=〈φ,φij〉,Tφ=∑i=1∞∑j=1∞λij2cijφij.If ψ is a continuous function on ℝ,
then(4.3)ψ(T)φ=∑i=1∞∑j=1∞ψ(λij2)cijφij,provided that ∑i=1∞∑j=1∞|ψ(λij2)cij|2<∞.
Hence, if the eigenvalue problem for T can be solved for the region Ω,
then the FPE problem with homogeneous boundary conditions can be easily solved
to give(4.4)φ(x,y)=∑i=1∞∑j=1∞〈g,φij〉λijαφij(x,y).


## 4.2. Nonhomogeneous Boundary Conditions

Before we
proceed further, we need to specify the definition of−(∇2)α/2.Definition 4.1.
Let{φij} be a complete set of orthonormal
eigenfunctions corresponding to eigenvalues λij2(≠0) of the Laplacian (−Δ) on a bounded region Ω with homogeneous BCs on ∂Ω.
Let(4.5)𝔉γ={f∈𝒟(Ω)∣∑i=1∞∑j=1∞|λijγ||cij|2<∞,cij=〈f,φij〉,γ=max(α,0)}.Then, for any f∈𝔉γ, (−Δ)α/2f is defined by(4.6)(−Δ)α/2f=∑i=1∞∑j=1∞λijαcijφij.If one of λij2=0 and φ0 is the eigenfunction corresponding to this
eigenvalue, then one needs 〈f,φ0〉=0.Proposition 4.2.
(1)
The operatorT=(−Δ)α/2 is linear and self-adjoint, that is, for f,g∈𝔉γ, 〈Tf,g〉=〈f,Tg〉.(2)
Iff∈𝔉γ,
where γ=max(0,α,β,α+β),
then(4.7)(−Δ)α/2(−Δ)β/2f=(−Δ)(α+β)/2f=(−Δ)β/2(−Δ)α/2f.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 Definition4.1
(1)
Forα=2m,m=0,1,2,…,
define (−Δ)α/2f=(−Δ)mf for any f∈C2m(Ω) (or other possibilities).(2)
Form−1<α/2<m,m=1,2,…,
define (−Δ)α/2f=T¯g,
where g=(−Δ)m−1f∈C2(m−1)(Ω),
and T¯ is the extension of T=(−Δ)α/2+1−m as defined by Proposition 4.3 below.
It suffices to
consider0<α<2.Proposition 4.3.
Letφij(x,y) be an eigenfunction corresponding to the
eigenvalue λij2(≠0) of the Laplacian (−Δ) on the rectangle Ω,
and let φ satisfy the BCs in problem 1. Then, if T¯ is an extension of T=(−Δ)α/2 (in symbols T⊂T¯) with adjoint T¯*⊂T*,(4.8)〈φij,T¯φ〉=λijα〈φij,φ〉−λijα−2{∫0bφij(0,y)k1f1(y)dy+∫0bφij(a,y)k2f2(y)dy+∫0aφij(x,0)k3f3(x)dx+∫0aφij(x,b)k4f4(x)dx}if ki≠0,i=1,…,4.
If ki=0,
the second term on the right-hand side becomes(4.9)−λijα−2{∫0b∂φij(0,y)∂xf1(y)h1dy−∫0b∂φij(a,y)∂xf2(y)h2dy+∫0a∂φij(x,0)∂yf3(x)h3dx−∫0a∂φij(x,b)∂yf4(y)h4dx}.Proof.
It is known that(4.10)〈φij,T¯φ〉=〈T¯*φij,φ〉=〈T*φij,φ〉=〈Tφij,φ〉=〈(−Δ)(−Δ)α/2−1φij,φ〉=〈(−Δ)α/2−1φij,(−Δ)¯φ〉,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)∫∫Ω(v∇2u−u∇2v)dS=∫∂Ω(v∂u∂n−u∂v∂n)ds.Thus,(4.12)〈φij,T¯φ〉=λijα−2〈φij,(−Δ¯)φ〉=−λijα−2∫∫Ωφij∇2φdS=−λijα−2{∫∫Ωφ∇2φijdS+∫∂Ω(φij∂φ∂n−φ∂φij∂n)ds}which gives the result on substitution.This result can be readily used to write down the
analytic solution to theFPE problem. First, we obtain the spectral
representation of the operator T by solving the eigenvalue
problem:(4.13)−Δφ=λ2φ,ℬ(φ)=0.Knowing the eigenvalues λij and the corresponding orthonormal (ON)
eigenfunctions φij,
we can use the finite-transform method with respect to φij and Proposition 4.3 to obtain(4.14)λijα〈φij,φ〉=〈g,φij〉+λijα−2bij,where λijα−2bij is the second term on the right hand side in
Proposition 4.3. Hence,(4.15)φ(x,y)=∑i=1∞∑j=1∞〈g,φij〉λijαφij(x,y)+∑i=1∞∑j=1∞bijλij2φij(x,y).


## 5. Results and Discussion

In this
section, we exhibit the results of applying Algorithm2.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)α/2φ=g0/k 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(5.1)φij={2sin(iπx)sin(jπy)}i,j=1∞,and the corresponding
eigenvalues λij2=π2(i2+j2).
The analytical solution is then given from Section 4 as(5.2)φ(x,y)=16g0kπ2∑i=0∞∑j=0∞sin[(2i+1)πx]sin[(2j+1)πy](2i+1)(2j+1)[(2i+1)2π2+(2j+1)2π2]α/2.
For the
numerical solution, a standard five-point finite-difference formula with equal
grid spacing h=1/n in the x and y directions has been used to generate the block
tridiagonal matrix A∈ℝ(n−1)2×(n−1)2 given in (2.1) as(5.3)A=[B−I−IB−I⋱⋱⋱−IB−I−IB]whereB=[4−1−14−1−14−1⋱⋱−1−14]∈ℝ(n−1)×(n−1).The parameters used to test this model are listed in
Table 1.Table 1
Physical parameters for test
problem 1.
ParameterDescriptionValuekThermal conductivity1Wm−1K−1g0Source10Wm−3αFractional
index0.5,1,1.5Test Problem 2: FPE with Mixed Boundary Conditions
Solve(−∇2)α/2φ=g0/k on the unit square [0,1]×[0,1] subject to type III boundary
conditions(5.4)−∂φ∂x+H1φ=H1φ∞atx=0,∂φ∂x+H2φ=H2φ∞atx=1,−∂φ∂y+H3φ=H3φ∞aty=0,∂φ∂y+H4φ=H4φ∞aty=1,where Hi=hi/k.
The analytical solution to this problem is given by(5.5)φ(x,y)=φ∞+g0k∑i=1∞∑j=1∞αi,jX(μi,x)Y(νj,y)(μi2+νj2)α/2Nx(μi)Ny(νj),where the eigenfunctions
are(5.6)X(μi,x)=μicos(μix)+H1sin(μix),Y(νj,y)=νjcos(νjy)+H3sin(νjy),with normalisation factors(5.7)Nx2(μi)=12[(μi2+H12)(1+H2μi2+H22)+H1],Ny2(νj)=12[(νj2+H32)(1+H4νj2+H42)+H3].The eigenvalues μi are determined by finding the roots of the
transcendental equation:(5.8)tan(μ)=μ(H1+H2)μ2−H1H2with νj determined from a similar equation for ν.
Finally, αi,j is given by(5.9)αi,j=∬01X(μi,ξ)Y(νj,η)Nx(μi)Ny(νj)dξdη.
For the numerical solution, a standard five-point
finite-difference formula with equal grid spacingh=1/n was again employed in the x and y 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 A∈ℝ(n+1)2×(n+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 Table2.Table 2
Physical parameters for test problem 2.
ParameterDescriptionValuekThermal conductivity5Wm−1K−1g0Source−10Wm−3h1,h3Heat
transfer coefficient0Wm−2K−1h2,h4Heat
transfer coefficient2Wm−2K−1φ∞External temperature20KαFractional
index0.5,1,1.5


### 5.1. Discussion of Results for Test Problem 1

A comparison of
the numerical and analytical solutions for test problem 1 is exhibited in
Figure1 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.Comparisons of
numerical (dashed lines) and analytical solutions (solid lines) for test
problem 1 computed using Algorithm2.4: (a) α=0.5,
(b) α=1.0,
(c) α=1.5,
and (d) α=2 (classical case).
(a)(b)(c)(d)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.Numerical
solutions for test problem 1: (a)α=0.5,
(b) α=1.0,
(c) α=1.5,
and (d) α=2 (classical case).
(a)(b)(c)(d)Remark 5.1.
Overℝn,
the Riesz operator (−Δ)α/2 as defined by(5.10)(−Δ)α/2f(x)=1g(α)∫ℝn|x−y|α−nf(y)dy,where f is a C∞- function with rapid decay at infinity, and
(5.11)g(α)=πn/22αΓ(α/2)Γ(n/2−α/2)is known to generate α-stable processes. In fact, the Green's
function of the equation(5.12)∂p∂t=−(−Δ)α/2p(t,x),T>0,x∈ℝ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 Algorithm2.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 x- and y-dimensions were divided equally into 31 divisions to produce the symmetric positive
definite matrix A∈ℝ900×900 having its spectrum σ(A)⊂[0.0205,7.9795].
One notes for this problem that the homogeneous boundary conditions necessitate
only the solution Φ=hαA−α/2g˜.
Algorithm 2.4 was still employed in this case; however in stage 1 we at first
solve Ax=g˜ 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.Figure3 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−15,
which represented an overall total of 110 matrix-vector products. This low tolerance was
enforced to ensure that as many approximate eigenpairs of A 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−10,
where θmax is the current estimate of the largest
eigenvalue of A.
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
Algorithm2.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
Figure4, 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.Comparisons of
numerical (dashed line) and analytical solutions (solid line) for test problem
2 computed using Algorithm2.4: (a) α=0.5,
(b) α=1.0,
(c) α=1.5,
and (d) α=2 (classical case).
(a)(b)(c)(d)The numerical solutions shown in Figure4 were again
generated using a standard five-point finite-difference stencil to construct
the matrix representation of the two-dimensional Laplacian operator. The x- and y-dimensions were divided equally into 30 divisions resulting in the symmetric positive
definite matrix A∈ℝ961×961 having its spectrum σ(A)⊂[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.Figure5 depicts the reduction in the residual of
the linear system for computing the solution Φ1 and the error in the matrix function
approximation for Φ2 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−15,
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 σ(A).
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 Φ2 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 Algorithm2.4.


## 5.1. Discussion of Results for Test Problem 1

A comparison of
the numerical and analytical solutions for test problem 1 is exhibited in
Figure1 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.Comparisons of
numerical (dashed lines) and analytical solutions (solid lines) for test
problem 1 computed using Algorithm2.4: (a) α=0.5,
(b) α=1.0,
(c) α=1.5,
and (d) α=2 (classical case).
(a)(b)(c)(d)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.Numerical
solutions for test problem 1: (a)α=0.5,
(b) α=1.0,
(c) α=1.5,
and (d) α=2 (classical case).
(a)(b)(c)(d)Remark 5.1.
Overℝn,
the Riesz operator (−Δ)α/2 as defined by(5.10)(−Δ)α/2f(x)=1g(α)∫ℝn|x−y|α−nf(y)dy,where f is a C∞- function with rapid decay at infinity, and
(5.11)g(α)=πn/22αΓ(α/2)Γ(n/2−α/2)is known to generate α-stable processes. In fact, the Green's
function of the equation(5.12)∂p∂t=−(−Δ)α/2p(t,x),T>0,x∈ℝ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 Algorithm2.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 x- and y-dimensions were divided equally into 31 divisions to produce the symmetric positive
definite matrix A∈ℝ900×900 having its spectrum σ(A)⊂[0.0205,7.9795].
One notes for this problem that the homogeneous boundary conditions necessitate
only the solution Φ=hαA−α/2g˜.
Algorithm 2.4 was still employed in this case; however in stage 1 we at first
solve Ax=g˜ 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.Figure3 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−15,
which represented an overall total of 110 matrix-vector products. This low tolerance was
enforced to ensure that as many approximate eigenpairs of A 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−10,
where θmax is the current estimate of the largest
eigenvalue of A.
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
Algorithm2.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
Figure4, 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.Comparisons of
numerical (dashed line) and analytical solutions (solid line) for test problem
2 computed using Algorithm2.4: (a) α=0.5,
(b) α=1.0,
(c) α=1.5,
and (d) α=2 (classical case).
(a)(b)(c)(d)The numerical solutions shown in Figure4 were again
generated using a standard five-point finite-difference stencil to construct
the matrix representation of the two-dimensional Laplacian operator. The x- and y-dimensions were divided equally into 30 divisions resulting in the symmetric positive
definite matrix A∈ℝ961×961 having its spectrum σ(A)⊂[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.Figure5 depicts the reduction in the residual of
the linear system for computing the solution Φ1 and the error in the matrix function
approximation for Φ2 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−15,
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 σ(A).
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 Φ2 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 Algorithm2.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 matrixA 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.



---

*Source: 104525-2009-03-11.xml*