Super-Exponentially Convergent Parallel Algorithm for Eigenvalue Problems with Fractional Derivatives

1 Introduction

It was recognized over the last decades that realistic models of various physical phenomena can be better described with fractional calculus; see, e.g., [44, 32, 31, 50, 25, 8, 33, 47, 20], just to mention a few. A very good overview of applications of fractional calculus is given in [49, 45].

There are various definitions of fractional derivatives; see, e.g., [37, 36, 34, 38, 32, 31, 1, 49]. In [46], a fractional derivative was defined in various ways, especially through integer derivatives using the Taylor and Fourier series. For example, one can introduce fractional derivatives for periodic functions using the Fourier series and the elementary relations

and replacing here the integer n by a real number. Alternatively one can use the Fourier or Laplace transform with the relations

and declare them valid for all real values of α.

The very popular definitions of the fractional calculus are the Riemann–Liouville and Caputo fractional derivatives and integrals; see, e.g., [10].

The Riemann–Liouville integral is defined by

where Γ⁢(α) is the gamma function and a is a fixed base point. Another notation, which emphasizes the base point, is

The following fundamental relations hold:

the latter of which is the semigroup property. These properties make possible not only the definition of fractional integration, but also the definition of fractional differentiation, by taking enough derivatives of Iα⁢f⁢(x). Computing the n-th order derivative over the integral of order (n-α), the α order derivative is obtained:

where n is the nearest integer bigger than α, i.e. n-1≤α<n∈ℤ+ or n=⌈α+1⌉.

The Caputo fractional derivative is another option for computing fractional derivatives. It was introduced by M. Caputo in 1967 in the following way:

where n-1≤α<n∈ℤ+. From these definitions, one can see that DtαaC⁢f(n)⁢(t)=Dtn+αaC⁢f⁢(t) and, for sufficient smooth f,

i.e., both derivatives coincide for functions f⁢(t) with f(k)⁢(a)=0, k=0,1,…,n-1. One can also rewrite this formula as

Both definitions differ only in the order of evaluation: whereas in the Caputo definition we first compute an ordinary derivative, then a fractional integral, in the Riemann–Liouville definition the operators are reversed. Comparing these definitions, we can see that functions which are derivable in the Caputo sense are much “fewer” than those which are derivable in the Riemann–Liouville sense.

Note that the fractional derivatives introduced above do not satisfy many properties of the classical differential calculus. But there are other definitions of fractional derivative in the literature which obey classical properties including: linearity, product rule, quotient rule, power rule, chain rule, vanishing derivatives for constant functions, Rolle’s theorem and the Mean Value Theorem; see, e.g., [22].

Using the definitions above, one can consider differential operators of fractional order, boundary and eigenvalue problems for these operators, etc., as well as various approximation methods for them; see, e.g., [31, 14, 11, 21].

The eigenvalue problem (EVP) is the problem of finding eigenpairs (eigenvalues and eigenfunctions or, in the language of mechanicians, frequencies and vibration shapes). It plays an important role in various applications concerned with vibrations and wave processes [19, 4, 41]. Popular methods such as the finite-difference method (FD), the finite element (FEM) and other variational methods, as well as spectral methods allow one to compute efficiently some lower eigenvalues only. At the same time there are applied problems requiring the computation of a great number (hundreds of thousands) of eigenvalues and eigenfunctions including eigenpairs with great indexes; see, e.g., [41, p. 273].

Over the last decade, it has been demonstrated that also eigen-oscillations of many systems in science and engineering can be modeled more accurately by employing fractional-order rather than integer-order derivatives [3, 13, 35]. In most of the fractional Sturm–Liouville formulations presented recently, the ordinary derivatives in a traditional Sturm–Liouville problem are replaced with fractional derivatives, and the resulting problems are solved using some numerical schemes such as the Adomian decomposition method [3], the fractional differential transform method [13], or using the method of the Haar wavelet operational matrix [35]. Some of the proposed algorithms are given in the literature without sufficient theoretical justification or possess the same drawbacks as the corresponding algorithms for the classical Sturm–Liouville problem. It turns out that the Sturm–Liouville problem with fractional derivatives can possess similar qualitative properties as a traditional Sturm–Liouville problem [51, 12], but some qualitative properties differ from the traditional ones.

In [9], Chen, Shen and Wang consider a spectral approximation of fractional differential equations (FDEs). A new class of generalized Jacobi functions (GJFs) is defined, which are the eigenfunctions of some fractional Jacobi-type differential operator and can serve as natural basis functions for properly designed spectral methods for FDEs. The efficient GJF Petrov–Galerkin methods for a class of prototypical fractional initial value problems (FIVPs) and fractional boundary value problems (FBVPs) of general order are constructed and analyzed.

In order to solve numerically an eigenvalue problem for fractional differential operators, we propose a new approach described below which we will refer to as the FD-method (from “functional-discrete method”, following [26, 27, 28, 5, 6]). Note that this approach for eigenvalue problems for nonlinear differential equations was applied for the first time in [28] and then continued in [15].

The FD-method is based on the perturbation and homotopy ideas. The perturbation in the case of ODE operators can be similar to that of the butt method (metodo dei tronconi; see, e.g., [7]), the Pruess method for BVPs [41, 39, 40] for the second-order ODEs, or of some methods for EVP from [2], where the coefficients of the differential equation are replaced by their piecewise constant approximations. This approach has been applied also to EVPs with multiple eigenvalues in [17].

The article is organized as follows. In Section 2, we describe the algorithm of the FD-method for EVPs in an abstract setting. Section 3 is devoted to the Fourier fractional derivative and some of its properties. In Section 4, we apply the FD-method for a differential equation with an integer highest derivative and a subordinated fractional derivative. Here we prove the superexponential convergence rate of our method independent of the definition of fractional derivatives in use and discuss two various algorithmic realizations: 1) solving the differential problems for the corrections of the FD-method, and 2) using a recursive procedure for the expansion coefficients of these corrections. Section 5 deals with the FD-method for a differential equation with a highest Riemann–Liouville fractional derivative, which is similar to the differential equation defining the generalized Jacobi functions (GJF). We prove the super-exponential convergence rate of our method. For the practical implementation we test the FD-algorithm directly and, besides, propose some new recursive procedure. Numerical examples are given to support the theoretical results.

2 The Homotopy-Based Method for EVPs in an Abstract Setting

Let us briefly explain the ideas of perturbation and homotopy for the eigenvalue problem

in a Hilbert space X with the scalar product (⋅,⋅) and with the null-element θ under the assumption that the spectrum of the operator A+B is discrete. We are looking for the eigenpair with a given fixed index n.

Let B¯ be an approximating operator for B in the sense that the eigenvalue problem

is “simpler” than problem (2.1).

Formally, a homotopy between two problems P1 and P2 with solutions u1 and u2 from some topological space X is defined to be a parametric problem PH⁢(t) with a solution u⁢(t) continuously and smoothly depending on the parameter t∈[0,1] and such that u⁢(0)=u1 and u⁢(1)=u2 (compare http://en.wikipedia.org/wiki/homotopy).

Following the homotopy idea for a given eigenpair number n, we embed our problem into the parametric family of problems

with W⁢(t)=B¯+t⁢φ⁢(B), φ⁢(B)=B-B¯ containing both problems (2.1) and (2.2), so that we obviously have

This suggests the idea to look for the solution of (2.3) in the form

where

Setting t=1 in (2.4), we obtain

provided that the series in (2.4) converge for all t∈[0,1].

The identities given in (2.5) are not suitable for a numerical algorithm, therefore we need another way to compute the corrections λn(j),un(j) which we describe below.

Substituting (2.4) into (2.3) and matching the coefficients in front of the same powers of t, we arrive at the following recurrence sequence:

with Fn(0)=0 and

For the pair λn(0),un(0) corresponding to the index j=-1 we get the so-called base eigenvalue problem

which for simplicity is assumed to have no multiple eigenvalues, to be “simpler” than the original one, and to produce the initial data for problems (2.6), (2.7). The case of multiple eigenvalues of the base problem was studied in [17].

Problems (2.6) for higher indices j≥0 are solvable provided that

from where we obtain

Under this condition the general solution of the inhomogeneous equation (2.6) with the singular operator can be represented by

with an arbitrary constant C. We choose the particular solution

satisfying the condition

The start values λn(0),un(0) for the recursion (2.6), (2.8) are the solutions of the base problem.

The truncated series

represent an algorithm to find the approximate solution λ𝑚n,u𝑚n (of rank m) to the solution of problem (2.1).

Below we give the error estimates of this method in the cases of a “dominated” fractional derivative and a “subordinated” fractional derivative.

3 The Fourier Fractional Derivative

The following differentiation and integration formulas can be easily proved for n∈ℕ, a∈ℝ:

We can generalize these formulas in a natural way for real n and introduce the differentiation and integration operators of fractional order α∈ℝ (the Fourier fractional derivative) by

cf. [24, 23]. One can see, especially for a=n⁢π, that

The functions 2⁢sin⁡(n⁢π⁢x), n=1,2,…, build an orthonormal basis in the space L2,0⁢(0,1) of functions vanishing at the ends of the interval with the norm ∥f∥L2=(∫01f2⁢(x)⁢𝑑x)1/2. An arbitrary function f⁢(x) from this space can be represented by the Fourier series

with an=2⁢∫01f⁢(x)⁢sin⁡(n⁢π⁢x)⁢𝑑x. The fractional derivative for such functions can be defined by

Given the Fourier representation (3.4), the inverse operator D-αF is defined by

One can also define the fractional derivative of a 2⁢π-periodic function using the exponential form of the Fourier series (see, e.g., www.xuru.org/fc/exponentials.asp):

Let us consider the asymptotic behavior with respect to n of the Caputo and Riemann–Liouville fractional derivatives of the function sin⁡(n⁢π⁢x). We have

where

The asymptotic behavior of φn(α) with respect to n gives the next simple lemma.

The lemma yields the estimate

which is of the same order in n as the Fourier fractional derivative (3.3). Due to (1.2) we have the same estimate for the Riemann–Liouville derivative too.

4 The Sturm–Liouville Problem with a Subordinated Fractional Derivative

Let us consider the following Sturm–Liouville problem:

where Dα denotes the Fourier, Caputo or Riemann–Liouville fractional derivatives. For shortness, we will use in this section the Fourier derivative DαF⁢u⁢(x), since the main property in use is the asymptotic (3.3) and (3.5) which are of the same order in n for all three derivatives.

If we approximate the coefficients k⁢(x),q⁢(x) by the constant 0 on the whole interval (the simplest variant), then the FD-method for (4.1) consists of the following sequence of recursive problems:

where

The solvability condition

implies

The particular solution of problem (4.2) satisfying the orthogonality condition

can be represented by

and we have

Using the orthonormality of the system 2⁢sin⁡(p⁢π⁢x) for the corrections of the eigenfunction, we obtain the estimates

and

with

The corrections to the eigenvalues are estimated by

Introducing the majorants Uj, Vj and Λj by

and replacing the inequality signs in (4.6)–(4.9) by equal signs, we obtain the following majorant system of equations:

A consequence of (4.10) is

where

The last equation in (4.11) is a recurrence equation of convolution type which can be solved by the method of generating functions (see, e.g., [42, 6]). We successively obtain

where

and rn∈(0,1) for n large enough.

From the definition of the majorant sequences we obtain the accuracy estimates

It is easy to see that

Therefore,

Now, the well-known Stirling’s asymptotic formula Γ⁢(t+1)≍2⁢π⁢t⁢(t/e)t implies

Analogously we obtain the corresponding estimate for the eigenvalues.

Thus, we come to the following assertion.