Rational Krylov methods for functions of matrices with applications to fractional partial differential equations

TL;DR

This paper introduces a new pole selection strategy for rational Krylov methods to efficiently approximate functions of positive definite matrices, especially fractional powers, with promising numerical results for fractional PDEs.

Contribution

It proposes a novel pole choice for rational Krylov methods tailored for fractional matrix functions, improving reliability and efficiency in solving fractional PDEs.

Findings

Numerical experiments confirm the effectiveness of the new pole selection.

The approach is promising for solving fractional partial differential equations.

The method improves approximation accuracy for fractional matrix functions.

Abstract

In this paper, we propose a new choice of poles to define reliable rational Krylov methods. These methods are used for approximating function of positive definite matrices. In particular, the fractional power and the fractional resolvent are considered because of their importance in the numerical solution of fractional partial differential equations. The results of the numerical experiments we have carried out on some fractional models confirm that the proposed approach is promising.

Equations

$$\left\{\begin{array}[]{ll}(-\Delta)^{\frac{\alpha}{2}}u=s(\mathbf{x},u),&\mathbf{x}\in\Omega,\\ u(\mathbf{x})=0,&\mathbf{x}\in\partial\Omega,\end{array}\right.$$

$$\left\{\begin{array}[]{ll}\frac{\partial u}{\partial t}=-\mu(-\Delta)^{\frac{\alpha}{2}}u+s(\mathbf{x},u,t),&(\mathbf{x},t)\in\Omega\times(t_{0},T],\\ u(\mathbf{x},t_{0})=u_{0}(\mathbf{x}),&\mathbf{x}\in\Omega,\\ u(\mathbf{x},t)=0\quad\mbox{(or}\quad\partial_{n}u(\mathbf{x},t)=0),&(\mathbf{x},t)\in\partial\Omega\times(t_{0},T],\end{array}\right.$$

$$(-\Delta {)}^{\alpha /2}u=\sum _{p=1}^{\infty }{\mu }_p^{\alpha /2}{c}_p{\phi }_p,u=\sum _{p=1}^{\infty }{c}_p{\phi }_p,$$

$$f(A_n)=\frac{1}{2\pi i}{\int }_{\Gamma }f(z)(zI-A_n{)}^{-1}dz.$$

$$\begin{array}{rl}f(A_n)\mathbf{v}=& \frac{1}{2\pi i}{\int }_{\Gamma }f(z)(zI-A_n{)}^{-1}\mathbf{v}dz\\ \approx & \frac{1}{2\pi i}{\int }_{\Gamma }f(z)V_k(zI-V_k^T{A}_n{V}_k{)}^{-1}{V}_k^T\mathbf{v}dz\\ =& V_{k}f(V_k^T{A}_n{V}_k)V_k^T\mathbf{v}.\end{array}$$

$$q_{k-1}(z)={\prod }_{j=1}^{k-1}({\sigma }_j-z).$$

$${\mathcal{Q}}_k(A_n,\mathbf{v})={\left[q_{k-1}(A_n)\right]}^{-1}{\mathcal{K}}_k(A_n,\mathbf{v}),$$

$${\mathcal{K}}_k(A_n,\mathbf{v})=\mathrm{Span}\{\mathbf{v},A_n\mathbf{v},\dots ,A_n^{k-1}\mathbf{v}\}$$

$$C_j=\left({\mu }_j{\sigma }_j{A}_n-I\right)({\sigma }_{j}I-A_n{)}^{-1},$$

$${\mathcal{Q}}_k(A_n,\mathbf{v})=\mathrm{Span}\{\mathbf{v},C_1\mathbf{v},\dots ,C_{k-1}\cdots C_2{C}_1\mathbf{v}\}.$$

$${\mathcal{W}}_{2k-1}(A_n,\mathbf{v})=\mathrm{Span}\{\mathbf{v},A_n^{-1}\mathbf{v},A_n\mathbf{v},\dots ,A_n^{-(k-1)}\mathbf{v},A_n^{k-1}\mathbf{v}\},$$

$$(\mu_{j},\sigma_{j})=\left\{\begin{array}[]{ll}(1,\infty),&\mbox{for $j$ even},\\ (0,0),&\mbox{for $j$ odd.}\end{array}\right.$$

$${\mathcal{W}}_k(A_n,\mathbf{v})=\mathrm{Span}\{\mathbf{v},(\sigma I-A_n{)}^{-1}\mathbf{v},\dots ,(\sigma I-A_n{)}^{-(k-1)}\mathbf{v}\},$$

$${\mathcal{W}}_k(A_n,\mathbf{v})=\mathrm{Span}\{\mathbf{v},({\xi }_{1}I+A_n{)}^{-1}\mathbf{v},\dots ,({\xi }_{k-1}I+A_n{)}^{-1}\cdots ({\xi }_{1}I+A_n{)}^{-1}\mathbf{v}\}.$$

$${\mathbf{w}}_j=({\xi }_{j}I+A_n{)}^{-1}{\mathbf{v}}_j$$

$${\mathbf{v}}_j$$

$${\mathbf{v}}_k$$

$$V_k^T{A}_n{V}_k=(I-H_k{D}_k)H_k^{-1}-h_{k+1,k}{V}_k^T{A}_n{\mathbf{v}}_{k+1}{\mathbf{e}}_k^T{H}_k^{-1},$$

$$h_{k+1,k}{V}_k^T{A}_n{\mathbf{v}}_{k+1}{\mathbf{e}}_k^T{H}_k^{-1}=h_{k+1,k}{V}_k^T{V}_{k+1}{H}_{k+1}{\mathbf{e}}_{k+1}{\mathbf{e}}_k^T{H}_k^{-1}=O_k,$$

$$V_k^T{A}_n{V}_k=(I-H_k{D}_k)H_k^{-1}.$$

$$f(A_n)\mathbf{v}\approx \beta V_{k}f(V_k^T{A}_n{V}_k){\mathbf{e}}_1,$$

$$\text{given }f(A_n),A_n,\mathbf{v},\text{ find }{R}_{k+\ell ,k}(\cdot )\text{ to minimize }\Vert f(A_n)\mathbf{v}-R_{k+\ell ,k}(A_n)\mathbf{v}{\Vert }_2^2,$$

$$z^{-\alpha /2}\approx \sum _{j=1}^k\frac{2\sin (\frac{\alpha }{2}\pi ){\tau }^{1-\alpha /2}}{\pi }\frac{{\omega }_j}{1+{\theta }_j}{\left(\frac{\tau (1-{\theta }_j)}{1+{\theta }_j}+z\right)}^{-1}\triangleq R_{k-1,k}\left(z\right),$$

$${ϵ}_r$$

$${\eta }_j$$

$$R_{k-1,k}(z)=\frac{p_{k-1}(z)}{q_k(z)}=\frac{\chi {\prod }_{r=1}^{k-1}(z+{ϵ}_r)}{{\prod }_{j=1}^k(z+{\eta }_j)},$$

$$\chi =\frac{{\eta }_k}{{\tau }^{\alpha /2}}\frac{\left(\genfrac{}{}{0ex}{}{k+\text{nicefrac}\alpha 2-1}{k-1}\right)}{\left(\genfrac{}{}{0ex}{}{k-\text{nicefrac}\alpha 2}{k}\right)}\prod _{j=1}^{k-1}\frac{{\eta }_j}{{ϵ}_j}.$$

$${\tilde{\tau }}_k\triangleq {\lambda }_{\text{min}}{\left(\frac{\alpha /2}{2ke}\right)}^2\exp \left(2W\left(\frac{4k^{2}e}{(\alpha /2{)}^2}\right)\right)$$

$$\tau=\tau_{k}:=\left\{\begin{array}[]{ll}{\tilde{\tau}}_{k},&\text{ if }k\leq\overline{k},\\ \displaystyle\left({\tilde{\sigma}}_{k}+\sqrt{{\tilde{\sigma}}_{k}^{2}+(\lambda_{\text{min}}\lambda_{\text{max}})^{\nicefrac{{1}}{{2}}}}\right)^{2},&\text{ if }k>\overline{k},\\ \end{array}\right.$$

$${\tilde{\sigma }}_k\triangleq -\frac{\text{nicefrac}\alpha 2}{8k}\ln \left(\frac{{\lambda }_{\text{max}}}{{\lambda }_{\text{min}}}\right){\lambda }_{\text{max}}^{\text{nicefrac}12},$$

$$\overline{k}=\frac{(\text{nicefrac}\alpha 2{)}^2}{8}{\left(\frac{{\lambda }_{\text{max}}}{{\lambda }_{\text{min}}}\right)}^{\text{nicefrac}12}\left[\ln \left(\frac{{\lambda }_{\text{max}}}{{\lambda }_{\text{min}}}\right)+2\right].$$

Full text

††journal: Journal of Computational Physics

url]lidia.aceto@unipi.it

url]bertaccini@mat.uniroma2.it

url]fabio.durastante@di.unipi.it

url]novati@units.it

1 Introduction

This paper deals with some computational issues concerning the use of Krylov-type methods for computing certain functions of matrices occurring in the solution of fractional partial differential equations. In particular, we focus on the following model problems:

[TABLE]

and

[TABLE]

where $\alpha\in(1,2]$, $\mu>0,$ and $-(-\Delta)^{\alpha/2}$ is the fractional Laplacian operator with homogeneous Dirichlet (or Neumann) boundary conditions. As usual, this operator can be defined using the spectral decomposition of the Laplacian, that is,

[TABLE]

where $\{\varphi_{p}\}_{p=1}^{\infty}$ are the eigenfunctions of $(-\Delta)$ and $\{\mu_{p}\}_{p=1}^{\infty}$ are the corresponding positive real eigenvalues. In this view, any positive definite matrix representing an appropriate discretization $A_{n}\in\mathbb{R}^{n\times n}$ of the Laplacian operator (that is, $\mathbf{x}^{T}A_{n}\mathbf{x}>0$ for any nonzero $\mathbf{x}\in\mathbb{C}^{n}$) can be used to approximate $(-\Delta)^{\alpha/2}$ by means of $A_{n}^{\alpha/2}.$ This approach is known as Matrix Transfer Technique (MTT) [1, 2, 3, 4, 5, 6]. Here we use finite differences, finite volumes and finite elements methods to generate the underlying approximation $A_{n}$.

For what concerns the steady state problem (1), the solution can be approximated by any method able to compute the action of $A_{n}^{-\alpha/2}.$ On the other hand, the solution of (2) requires also the discretization in time. Since $A_{n}$ inherits the stiffness of the Laplacian operator, the time integration of (2) needs to be implicit. Therefore, at each time step, we need to compute the action of matrix functions of the type $\left(I+\nu A_{n}^{\alpha/2}\right)^{-1},$ where $I$ denotes the identity matrix of size $n$ and $\nu>0$ is a parameter that depends on the stepsize and the integrator. Assuming that the stiffness is only due to the Laplacian operator, here we consider IMEX-type methods, which allow an explicit treatment of the forcing term in (2).

In this work, we analyze some Krylov-type methods approximating the action of the product of the above matrix functions by a given vector. We compare the behavior of the polynomial, shift-and-invert, extended, and rational Krylov methods in terms of accuracy and computational cost on various test problems and different discretizations. Our aim is to show that our rational Krylov method can outperform the other approaches provided that the poles of the underlying rational forms are suitably defined. To this purpose, we define the poles using the Gauss-Jacobi approach for the computation of $A_{n}^{-\alpha/2}$ as in [7, 8, 9], and then extend this idea for the computation of $\left(I+\nu A_{n}^{\alpha/2}\right)^{-1}\mathbf{v}$ giving also some details for a reliable implementation. We prove that our new poles are real and simple, so that they can be computed with a root-finder less prone to error amplification. We remark that the computation of $\left(I+\nu A_{n}^{\alpha/2}\right)^{-1}\mathbf{v}$ with Krylov methods has been considered in [10], but by using only the shift-and-invert approach.

The paper is organized as follows. In Section 2 we recall the basic features of Krylov methods of polynomial and rational types tailored for approximating function of matrices, while in Section 3 we give a theoretical and computational analysis of the poles of the proposed rational methods. In Section 4 we discuss some issues of the MTT based on finite differences, finite elements, and finite volumes approximation applied to the Laplacian operator. Finally, Section 5 is devoted to some numerical experiments in which we compare the performances of our approach with the other Krylov methods presented in Section 2.

2 Rational Krylov Methods

We recall that for a given matrix $A_{n}\in\mathbb{R}^{n\times n}$ and a function $f$ that is analytic on and inside a closed contour $\Gamma$ that encloses the spectrum of $A_{n}$, the matrix function $f(A_{n})$ is defined as

[TABLE]

Let $V_{k}$ be an orthogonal matrix whose columns $\mathbf{v}_{1},\dots,\mathbf{v}_{k}$ span an arbitrary Krylov subspace $\mathcal{W}_{k}(A_{n},\mathbf{v})$ of dimension $k$. We obtain an approximation of $f(A_{n})\mathbf{v}$ by

[TABLE]

Different methods for the approximation of matrix functions are obtained for different choices of the projection spaces $\mathcal{W}_{k}(A_{n},\mathbf{v})$. Given a set of scalars $\left\{\sigma_{1},\dots,\sigma_{k-1}\right\}\subset\overline{\mathbb{C}}$ (the extended complex plane), that are not eigenvalues of $A_{n}$, let

[TABLE]

The rational Krylov subspace of order $k$ associated with $A_{n}$, $\mathbf{v}$ and $q_{k-1}$ is defined by

[TABLE]

where

[TABLE]

is the standard polynomial Krylov space. By defining the matrices

[TABLE]

where $\left\{\mu_{1},\dots,\mu_{k-1}\right\}\subset\overline{\mathbb{C}}$ are such that $\sigma_{j}\neq$ $\mu_{j}^{-2}$, it is known that the rational Krylov space can also be written as follows (see [11])

[TABLE]

This general formulation allows to recast most of the classical Krylov methods in terms of a rational Krylov method with a specific choice of $\sigma_{j}$ and $\mu_{j}$. In particular, the standard (polynomial) Krylov method in which $\mathcal{W}_{k}(A_{n},\mathbf{v})=\mathcal{K}_{k}(A_{n},\mathbf{v})$ can be recovered by defining $\mu_{j}=1$ and $\sigma_{j}=\infty$ for each $j$. The extended Krylov method (see [12, 13]), in which

[TABLE]

is obtained by setting

[TABLE]

The shift-and-invert (single pole) rational Krylov, see [14, 15], where

[TABLE]

is defined by taking $\mu_{j}=0$ and $\sigma_{j}=\sigma$ for each $j$.

As for the use of Krylov methods for the computation of functions of operators involving $(-\Delta)^{\alpha/2}$, that is for our cases of interest, we quote here [5] and [10] in which the standard Krylov method and the shift-and-invert approach are investigated, respectively. On the same line, the computation of the matrix square root ($\alpha=1$) is done by using the extended Krylov method in [16], see also [17], and the rational Lanczos approximation in [18].

In this work we consider a rational Krylov approach in which $\mu_{j}=0$ for each $j$, and where $-\sigma_{j}=\xi_{j}>0$ are suitably defined (see Section 3). In this view, the method here presented is a rational Krylov method in which

[TABLE]

Starting from $\mathbf{v}_{1}=\mathbf{v/}\beta$, where $\beta=\|\mathbf{v}\|_{2}$, we determine $\mathbf{v}_{j+1}$ by orthogonalizing the vector

[TABLE]

against $\mathbf{v}_{1},\dots,\mathbf{v}_{j}$, followed by normalization. In this way, a sequence of vectors $\{\mathbf{v}_{j}\}_{j=1}^{k}$ is generated such that

[TABLE]

obtaining the following Arnoldi–like decomposition:

[TABLE]

where $D_{k}=\operatorname{diag}(\{\xi_{j}^{-1}\}_{j=1}^{k})$ and $H_{k}$ is the Hessemberg matrix $H_{k}=[h_{i,j}]$. Since

[TABLE]

we find the following expression for the projected matrix

[TABLE]

Finally, we approximate $f(A_{n})\mathbf{v}$ as

[TABLE]

where $\mathbf{e}_{1}=(1,0,\dots,0)^{T}\in\mathbb{R}^{k}.$

As is well known, rational Krylov methods are generally quite fast (in terms of iterations with respect to the polynomial counterpart) whenever $A_{n}$ represents an unbounded self-adjoint operator. In this view, the most computational demanding part of these methods is typically due to the solution of the linear systems (5). The computational cost of an iterative solver for (5), whenever $A_{n}$ is large, is highly problem-dependent and often influences very much the performance of these algorithms and therefore their competitiveness with respect to polynomial methods. Anyway, in the examples considered here, we always deal with sparse and, in particular, banded matrices $A_{n}$. In this sense, we compute the approximate solution of (5) by means of sparse direct solvers. Alternatively, one may consider iterative solvers with preconditioners able to handle multiple shifts (see e.g. [19, 20, 21, 22]).

In Section 3 we propose to select the poles $\xi_{j}$ by exploiting some insights on the functions to be approximated here, i.e. $f(z)=z^{-\alpha/2}$ and $f(z)=(1+\nu z^{\alpha/2})^{-1}$, by working with suitable rational approximations. Clearly, one could look for an automatic selection of the poles that avoids the use of such information. An example of application of this strategy is represented by the RKFIT algorithm [23]. It can be used for computing a rational approximation for a matrix function $f(A_{n})$ by solving a rational least square problem of the form

[TABLE]

where $f$ is a matrix function (or an approximation), and $R_{k+\ell,k}$ is a rational function of type $(k+\ell,k)$ and $\ell\geq-k$. Thus, the minimum is taken on the roots of the polynomials defining $R_{k+\ell,k}$. To effectively use this procedure, one needs to find some starting poles for $R_{k+\ell,k}(z)$ and the values of $k$ and $\ell.$

3 Poles Selection

Our proposal for selecting the poles for the construction of the rational Krylov subspace (4) relies on the rational approximation of $z^{-\alpha/2}$ proposed in [8, 9, 7]. In particular, following [8, eq. (5)], we get

[TABLE]

where $\omega_{j}$ and $\theta_{j}$ are, respectively, the weights and nodes of the Gauss–Jacobi quadrature formula with weight function $(1-x)^{-\frac{\alpha}{2}}(1+x)^{\frac{\alpha}{2}-1}$ and $\tau$ is a positive real parameter that should be defined suitably. In addition, denoting by $\zeta_{r}$ the $r$th zero of the Jacobi polynomial ${\mathcal{P}}_{k-1}^{({\alpha}/{2},1-{\alpha}/{2})}\left(z\right)$ and setting

[TABLE]

from [8, Proposition 1] we can express $R_{k-1,k}(z)$ as the rational function

[TABLE]

where

[TABLE]

It is worth noting that $\{\theta_{j}\}_{j=1}^{k}$ are the zeros of the Jacobi polynomial ${\mathcal{P}}_{k}^{(-{\alpha}/{2},{\alpha}/{2}-1)}\left(z\right).$ In this context, this implies that the values $\eta_{j}$ are all real and simple and therefore the roots of $q_{k}(z)$ as well.

Given a positive definite matrix $A_{n}$ with spectrum $\sigma(A_{n})\subseteq[\lambda_{\min},\lambda_{\max}],$ a specific rational approximation for the matrix function $f(A_{n})=A_{n}^{-\alpha/2}$ can be deduced from the above, provided that the parameter $\tau$ is selected. To this purpose, setting

[TABLE]

with $W(\cdot)$ denoting the Lambert-W function, in [7, Propositions 3.4, 4.1] it has been shown that for the matrix function $A_{n}^{-\alpha/2}$ a reliable value of $\tau$ is given by

[TABLE]

where

[TABLE]

and

[TABLE]

This choice is the result of an analysis on the Padé-type approximations of $z^{-\alpha/2}.$ Now, by using the above arguments, we can get an approximation also for $f(z)=(1+\nu z^{\nicefrac{{\alpha}}{{2}}})^{-1}.$ Indeed, we can write

[TABLE]

For the roots of $\tilde{q}_{k}(z)$ we observe the following useful result.

Proposition 1**.**

All the roots of the polynomial $\tilde{q}_{k}(z)$ are real and simple.

To prove this we need two auxiliary results. First, there is an interlacing of the zeros of Jacobi polynomials of different orders and weights.

Theorem 1**.**

[24, Theorem 2.3]** Let $\beta,\gamma>-1.$ For any $t,s\in[0,2],$ let

$-1<\theta_{1}<\theta_{2}<\ldots<\theta_{n}<1$* be the zeros of ${\mathcal{P}}^{(\beta,\gamma)}_{n},$ and*

• 2.

$-1<\zeta_{1}<\zeta_{2}<\ldots<\zeta_{n-1}<1$* be the zeros of ${\mathcal{P}}^{(\beta+t,\gamma+s)}_{n-1}$.*

Then

[TABLE]

Theorem 2**.**

[25, Theorem 8]** Let $p$ and $q$ be real polynomials. Then $p$ and $q$ have strictly alternating roots if and only if all polynomials in the space

[TABLE]

have real and simple roots.

Proof of Proposition 1.

To apply Theorem 2, we need to show that the roots of the polynomials $p_{k-1}(z)$ and $q_{k}(z)$ are strictly interlaced. By construction we know that the roots of $p_{k-1}(z)$ are the real values $\{-\epsilon_{r}\}_{r=1}^{k-1}$ in (6), while the roots of $q_{k}(z)$ are the real values $\{-\eta_{j}\}_{j=1}^{k}$ in (7). By Theorem 1 applied to the Jacobi polynomials ${\mathcal{P}}_{k}^{(-{\alpha}/{2},{\alpha}/{2}-1)}\left(z\right)$ and ${\mathcal{P}}_{k-1}^{({\alpha}/{2},1-{\alpha}/{2})}\left(z\right)$ we have

[TABLE]

Then, for the interlacing between the $\{-\epsilon_{r}\}_{r=1}^{k-1}$ and the $\{-\eta_{j}\}_{j=1}^{k}$, we need to prove that $-\eta_{1}<-\epsilon_{1}<-\eta_{2}<-\epsilon_{2}<\ldots<-\eta_{k-1}<-\epsilon_{k-1}<-\eta_{k}$ or, equivalently,

[TABLE]

Now, using (6) and (7) for $i=1,\dots,k-1,$ we have

[TABLE]

Since $\tau>0,$ $(1+\zeta_{i})>0$ and $(1+\theta_{i+1})>0,$ the previous inequalities are satisfied if and only if $\zeta_{i}<\theta_{i+1},$ which is true by (8). Similar arguments lead to verify that $\eta_{i}>\varepsilon_{i}.$ By using Theorem 2, the proof is complete. ∎

In addition, since all the coefficients of ${\tilde{q}}_{k}(z)$ are strictly positive by construction, according to the Descartes’ rule of signs, we are also sure that all its roots are negative.

In consideration of all the above arguments, the poles that we use for the computation of the proposed rational Krylov methods are positive, real and simple. In fact, when $f(A_{n})=A_{n}^{-\alpha/2},$ in (4) we take as poles the opposite of the roots of $q_{k}(z)$ (i.e., $\xi_{j}=\eta_{j}$), while, when $f(A_{n})=(I+\nu A_{n}^{\nicefrac{{\alpha}}{{2}}})^{-1},$ we take as poles the opposite of the roots of ${\tilde{q}}_{k}(z).$ Thus, one of the main features of our proposal is represented by the possibility to work only with real arithmetic. From now on, we refer to the rational Krylov methods based on these poles as Krylov Jacobi.

Just to provide an example of computed poles, in Figure 1 we focus on the matrix $A_{n}$ given by the centered differences discretization of the 1D Laplacian with $n=2^{12}$ and $\alpha=1.2.$ In this case, for $k=10,20,30,$ we computed the roots $\{\theta_{j}\}_{j=1}^{k}$ and $\{\zeta_{r}\}_{r=1}^{k-1}$ by the JACOBI_POLYNOMIAL Matlab code; see [26]. On the other hand, we used the package MPSolve to compute the roots of $\tilde{q}_{k}(z)$; see [27, 28].

4 Solving the model problems by MTT using Krylov methods

Computing efficiently functions of matrix-vector products by Krylov methods is beneficial for many applications. However, in order to provide reasonable comparisons, here we consider one class of mathematical tools which is quite wide and important because it forms the ground of many numerical models: (evolutionary) fractional partial differential equations.

As a case study, let us consider the $d$–dimensional Laplacian operator

[TABLE]

on the domain $\Omega$ subject to Dirichlet (or Neumann) boundary conditions. The Matrix Transfer Technique (MTT) introduced in [3, 4] considers approximations to the fractional Laplacian of the form

[TABLE]

where $A_{n}$ is a suitable matrix approximation of the standard Laplacian. In this work, we consider the following discretizations:

finite differences (FD);

• 2.

finite elements methods (FEM) with first order Lagrangian elements on a triangular tessellation;

• 3.

cell–centered finite volume (FV) methods with a piecewise linear basis, [29].

For FEM and FV we get $A_{n}=M_{n}^{-1}K_{n},$ where $K_{n}$ denotes the stiffness matrix, while $M_{n}$ is the mass matrix for the FEM or the diagonal matrix of the cell volumes for FV. In these cases, even if $A_{n}$ is generally not symmetric, it is similar to a positive definite matrix. In fact, $M_{n}^{\frac{1}{2}}A_{n}M_{n}^{-\frac{1}{2}}=M_{n}^{-\frac{1}{2}}K_{n}M_{n}^{-\frac{1}{2}}.$ Whenever the mass matrix $M_{n}$ is not trivial, i.e. $A_{n}$ is not a multiple of $K_{n},$ we should solve linear systems with coefficient matrices $(\xi_{j}I+M_{n}^{-1}K_{n}),$ by means of a direct method. This can be inappropriate because it would require assembling the matrix $A_{n}=M_{n}^{-1}K_{n}$ that can be dense. Therefore, instead of computing the basis of the rational Krylov space with $A_{n},$ we generate it working with $\hat{A}_{n}=\xi_{j}M_{n}+K_{n}$. This is completely analogous to what is usually done for generalized eigenvalues for the pencil $(K_{n},M_{n})$ with rational Krylov methods, see, e.g., [30, 31]. Then, with the attained basis, we approximate the matrix–vector product $f(A_{n})\mathbf{v}$ as in (3).

Finally, we stress that direct solvers for sparse linear systems can benefit by the use of permutations to reduce fill–in. This can be mostly useful when unstructured matrices are used for either FEM or FV discretizations.

4.1 Integration in time

If we consider the solution of the problem (2) semidiscretized with respect to the space variables, we get the initial value problem

[TABLE]

where ${\bf y}(t):\mathbb{R}\rightarrow\mathbb{R}^{n}$ and ${\bf y_{0}}\in\mathbb{R}^{n}$ (with respect to problem (2), to simplify notation, here we set $\mu=1$ and $s(\mathbf{x},u,t)=s(\mathbf{x},t)$). To integrate (10) in time, we apply a Linear Multistep Method (LMM) over an equispaced grid $t_{j}=t_{0}+j\delta_{t},j=0,1,\ldots,n_{t},$ with $\delta_{t}=(T-t_{0})/n_{t}$ the selected stepsize. Setting $\mathbf{g}(\mathbf{y},t)=-A_{n}^{\frac{\alpha}{2}}\mathbf{y}(t)+\mathbf{s}(t),$ we get

[TABLE]

with

[TABLE]

Since the fractional Laplacian operator is unbounded and give stiffness , implicit schemes are preferable to avoid possible severe stepsize restrictions in order to satisfy stability requirements; see, e.g., [32]. Using implicit schemes, we need to solve linear systems of the form

[TABLE]

that can be recast in the computation of a function of matrix times vector, that is,

[TABLE]

where

[TABLE]

and

[TABLE]

By using similar arguments, we can also deal with semi–linear problems having non–linear forcing term

[TABLE]

Assuming that the function $\mathbf{s}$ is not responsible for further stiffness, we can use implicit–explicit methods (IMEX), similarly to [8]. A generic $\ell$-step IMEX method for (13) can be written as

[TABLE]

For further details on IMEX methods see, e.g.,  [33, 34]. Therefore, we need again to compute matrix functions of the form (12).

5 Numerical Experiments

The examples in this section are collected in two groups. In Section 5.1 we deal with the discrete version of the steady state problem (1), while in Section 5.2 the time dependent problem (2).

All the numerical experiments are performed on a laptop running Linux with 8 Gb memory and CPU Intel® Core™  i7-4710HQ CPU with clock 2.50 GHz. The codes are written and executed in MATLAB R2018a. The following external codes are used in our routines:

MPSolve package for the computation of the polynomial zeros; see [27],

• 2.

\mlttfamilyrat_krylov for the construction of the basis of the rational Kyrlov methods; see [11],

• 3.

\mlttfamilyEKS class for the construction of the basis of the extended Krylov method; see [13],

• 4.

FENICS library [35] v.2018.1 to assemble the finite elements matrices.

• 5.

FiPy library [36] v.3.1.3-dev2-g11937196 to assemble the finite volume matrices.

The auxiliary linear systems are solved here by Matlab’s standard backslash. Whenever is feasible, the reference solution $\mathbf{u}^{*}$ for the various problem is computed directly by the Schur–Parlett algorithm. We denote by $\varepsilon=\|\mathbf{u}^{*}-\mathbf{u}\|_{2}/\|\mathbf{u}^{*}\|_{2}$ the relative error. All the timings are measured in seconds averaged on one hundred runs. For the definition of the poles, the interval $[\lambda_{\min},\lambda_{\max}]$ is always assumed explicitly known.

5.1 Stationary problems

We report the results for the steady state problem (1) discretized using second order centered differences for the 1D and 2D Laplacian on the domains $[0,1]$ and $[0,1]^{2}$, respectively. The right–hand side terms for the two problems are $f(x)=\sin(\pi x)$, and $f(x,y)=\sin(\pi x)\sin(\pi y)$. We observe the expected behavior: rational Krylov Jacobi algorithms outperform (in term of iterations) the polynomial Krylov; see Figure 2 for the behavior for the one dimensional problem, and Figure 3 for the two dimensional. In both cases the shift parameter $\sigma$ for the shift–and–invert method is computed as $\sigma=\sqrt{\lambda_{\min}(A_{n})\lambda_{\max}(A_{n})}$, cf. [18].

5.2 Time–dependent problems

In this section, we focus on the problem (2). First, we set $s=0,$ $\mu=1,$ $T=1,$

[TABLE]

on the unit square $\Omega=[0,1]^{2}$ and consider Dirichlet boundary conditions. We refer to this test problem as P1. By using finite differences (FD), finite volumes (FV) and finite elements methods (FEM) we discretize this problem and test on it the underlying rational Krylov methods. Then, we consider the problem (2) with $s$ a nonzero polynomial function of the solution, producing the fractional Allen–Cahn equation, that is

[TABLE]

Here $\mu>0$ is a small parameter defining the thickness of the interface separating the different phases and $u_{0}\in\mathbb{L}^{2}(\Omega).$ In our numerical experiments we choose $\Omega=[0,1]^{2},[t_{0},T]=[0,4],$ and

[TABLE]

Refer to this test problem as P2. We apply finite volumes and finite elements methods to discretize this problem. In both cases the obtained matrices $A_{n}$ are positive definite, and thus the analysis in Section 3 applies straightforwardly.

P1 by using FD

We consider the five–point stencil of the FD discretization for the Laplacian operator on the grid with $n_{x}=n_{y}=n_{t}=2^{6}$ points. The method for marching in time is the implicit Euler. Therefore, as discussed in Section 4, in order to advance in time, we need to compute the matrix function

[TABLE]

where $c$ is a constant. In Figure 4 the relative error and execution time comparison for the first iterate of the method is reported. The poles for the RKFIT algorithm are computed on the complete matrix function on a reduced size grid ($n_{x}=n_{y}=2^{4}$ and $n_{t}=2^{6}$) from the initial guess $\{\sigma_{j}=+\infty\}_{j=1}^{k}$. Note that the RKFIT can compute only $k=10$ different nodes. We recall also that for the extended Krylov algorithm the number $k$ represent the size of the Krylov subspaces, thus the number of linear auxiliary linear systems solved is $\lfloor{k}/{2}\rfloor$. This explain the behavior with respect to time: the Krylov space of the same dimension for the extended Krylov costs roughly a half of the rational subspace of the same size.

P1 by using FV

We consider the cell–centered FV discretization of the problem (2) on the same grid used for finite differences. In particular, $n_{x}=n_{y}=n_{t}=2^{6}$. Results in Figure 5 are comparable with those of the FD discretization. The RKFIT algorithm produces less than the requested number of nodes and at a greater computational effort. By fixing the size of the Krylov space, we observe that the Krylov Jacobi method is more accurate than the extended Krylov method which, on the other hand, keeps the same ratio with respect to the achieved timings.

P1 by using FEM

Consider here the FEM discretization of problem (2) on the structured mesh made by triangular cells with straight sides, i.e., each side of the squared domain $\Omega$ is divided into $n_{x}$ and $n_{y}$ rectangles, and then divided into a pair of triangles, with Lagrangian elements of order $1$. From Figure 6 we observe that the RKFIT algorithm produces again less than the requested number of nodes and at a greater computational effort. Considering Krylov spaces of the same size, the Krylov Jacobi method is more accurate than the extended Krylov. Moreover, the latter in this example is only marginally better with respect to the execution time.

P2 by using FV

Consider the cell–centered FV for (16) over a uniform quadrangular mesh with $n_{x}=n_{y}=64$ intervals on the unit square $[0,1]^{2}.$ To solve this problem we consider the application of the IMEX backward Euler method with $\delta_{t}=10^{-2}.$

The results are collected in Figure 7, from which we observe that the Krylov Jacobi method achieves the better accuracy already with $k=5$ nodes. The RKFIT Algorithm, tuned on the problem with $n_{x}=n_{y}=32$, fails to compute acceptable nodes in this case, see the example of the errors for $\alpha=1.8$ given in Figure 8.

P2 by using FEM

Consider the discretization of (16) by means of Lagrange linear elements over a uniform finite element mesh on the unit square $[0,1]^{2},$ i.e., a mesh consisting of triangular cells with straight sides dividing each side of the square into $n_{x}$ and $n_{y}$ rectangles, each divided into a pair of triangles, for $n_{x}=n_{y}=80.$ The integration method in time is the IMEX backward Euler with $\delta_{t}=10^{-2}.$ The performance of the considered rational Krylov method are shown in Figure 9. Note that the behavior is analogous to the one observed for the time–dependent problem (2) (compare the results with those in Figure 6).

6 Conclusion and perspectives

We considered certain matrix functions involving the fractional power of real positive definite matrices computed by rational Krylov methods with application to the numerical approximation of fractional-in-space partial differential equations. We proposed a new choice of the poles, real and simple by construction, that appears to be effective and faster with respect to polynomial, shift-and-invert, extended Krylov methods, and RKFIT.

Within the same framework, one can apply the proposed rational Krylov Jacobi method also to fractional linear multistep methods [37] for solving (2) in which the first order derivative in time is replaced by a fractional derivative of order $\beta,$ with $\beta\in(0,1).$

In principle, the proposed method can work on irregular domains, Robin or Dirichlet BC, see, e.g., [38, 39, 40], and can be used also in contexts of adaptivity, provided that these generate a symmetric positive definite matrix. In the latter case, studying how the selected poles vary as the mesh is changed would be of interest, and could also open new alternative approaches.