Exact and numerical solutions of the fractional Sturm–Liouville problem

1 Introduction

In the paper, we study (from a theoretical and numerical point of view) the fractional Sturm-Liouville problem (FSLP) with the homogeneous mixed boundary conditions. The Sturm-Liouville problem in a fractional version can be derived by using different approaches. The first one consists of replacing the integer order derivative in the classical Sturm-Liouville problem by a fractional order derivative [1, 14]. However, this approach does not lead to orthogonal systems of eigenfunctions. The second approach is connected with the application of the calculus of variations [16, [17]. In this case, the obtained fractional differential equations can be interpreted as fractional Euler-Lagrange equations [2, 16, 18, 21, 22, 27]. They contain the differential operator, which is a composition of the left and the right fractional derivative [16, 17]. This feature leads to fundamental difficulties in calculating eigenvalues and deriving the exact solutions in a closed form, even in a simple case of a fractional oscillator problem in a bounded domain. Explicit solutions and eigenvalues are known, so far, only for a few FSLP, like fractional oscillator problem on unbounded domain [28], and for some singular cases like fractional Legendre and Jacobi problems [17, 31, 32], and fractional Bessel equation [23].

The FSLP, along with its eigenfunction’s system and eigenvalues, is connected to a fractional diffusion [19, 20] in a bounded domain. The term ‘fractional diffusion’ is to be understood as the application of fractional derivatives in description of processes of anomalous diffusion. Such diffusive processes appear in many fields of science and engineering, e.g., heat conduction in materials, fluid pressure in porous media, human migration, movement of proteins in cells, transport of lipids on cell membranes, transport on social networks, bacterial motility, and others. Classical results show that in order to solve diffusion equations in a bounded domain, one needs to apply the suitable orthogonal systems of functions, usually connected to a respective Sturm-Liouville problem. Therefore FSLPs, determined in bounded domains, are an emerging meaningful area of the fractional differential equations theory. The orthogonal eigenfunctions’ systems of FSLPs are and will be a useful tool in solving partial fractional differential equations connected to anomalous diffusion processes. Preliminary results are given in papers [17, 18, 19, 20] and show that by applying the eigenfunctions’ systems, we can study fractional diffusion problems with variable diffusivity and calculate the explicit solutions or establish the existence- uniqueness results and analyze the properties of solutions. Similar to the classical Sturm-Liouville theory, it appears that the existence of the purely discrete spectrum and the associated orthogonal eigenfunctions’ system is strongly connected to the singularity of FSLP, or in case of a regular FSLP, to the choice of boundary conditions. In paper [20], the proof was given for the regular FSLP subjected to the homogeneous Dirichlet conditions, now we shall present the result on the regular FSLP subjected to the homogeneous mixed boundary conditions. This new result is developed by converting the differential FSLP to the equivalent integral one and by applying the results of the Hilbert-Schmidt operators theory.

As we have mentioned, the problem of finding analytical solutions of the FSLP, containing both the left and right fractional derivative, is still a challenge for scientists. Therefore, numerical methods of solving the FSLP are being simultaneously developed. During the last few years, several numerical algorithms (based on direct or indirect methods) have been proposed to obtain approximate solutions of the fractional Euler-Lagrange equations [5, 6, 7, 8, 9, 10, 11, 30].

The same problem appears in calculating the eigenvalues of the FSLP. The most common approach to determine eigenvalues and eigenfunctions for Sturm-Liouville operators of integer and fractional order is to use a numerical method. The numerical solution of the Sturm-Liouville problem of integer order can be found in literature i.e. the Pruess, shooting and finite difference methods [4, 26]. In paper [29], the control volume method is used to determine the eigenvalues of the classical Sturm-Liouville problem. However, for FSLPs involving both the left and the right derivative, the adequate set of numerical tools still requires further and extensive work.

In our previous paper [12], we developed the numerical method for solving a fractional eigenvalue problem - the version of the FSLP with the homogeneous mixed boundary conditions. The proposed numerical scheme was based on the discretization of Caputo derivatives involving the boundary conditions. This approach allowed us to approximate the eigenfunctions keeping their orthogonality at each step of approximation. Moreover, the convergence was controlled by using the experimental rates of convergence formulas both for the eigenvalues and for the eigenvectors. The obtained rate of convergence was close to 1.

In paper [3], the FSLP with Dirichlet boundary conditions was considered. The authors analysed two approaches to the FSLP: discrete and continuous. They investigated the numerical solution of the FSLP by using the truncated Grunwald-Letnikov fractional derivative.

Now, we will construct a numerical scheme to calculate the approximate eigenvalues and eigenfunctions, by applying the approach presented in papers [7, 9, 11]. First, we transform the FSLP into an intermediate integral equation and then we discretize the obtained equation by using the numerical quadrature rule.. This method allows us to obtain the numerical scheme for which the experimental rate of convergence in all the considered examples tends to 2α. As we study the FSLP with order α > 1/2, we clearly see that the new numerical method gives better convergence than the one introduced in [12], while the orthogonality of the approximate eigenfunctions is also kept at each step of the procedure.

The paper is organized as follows. Section 2 presents the analyzed problem and recalls basic definitions and main properties of fractional differential and integral operators. In Section 3 the exact solution of the FSLP is depicted, while in Section 4 the numerical solution is given. Finally, in Section 5, we show numerical results for two examples of the FSLP, and we conclude the paper with a section containing brief conclusions.

2 Preliminaries

We recall the left and right Caputo fractional derivatives of order α ∈ (0, 1) (see e.g. [13, 15, 24])

and the left and right Riemann-Liouville fractional derivatives of order α ∈ (0, 1) ([13, 15, 24])

where the operators Ia+αandIb−α are respectively the left and the right fractional Riemann-Liouville integrals of order α > 0 defined by

We also recall the composition rules of fractional operators for the case of order α ∈ (0, 1]

and for the Riemann-Liouville derivatives

All the above rules are fulfilled for all x ∈ [a, b] when function y is a continuous one.

Now, we shall quote the general formulation of the fractional eigenvalue problem, introduced and investigated in papers [16, 17].

Let us observe that in the case α = 1 we have: cDa+1y=y′ and cDb−1y=−y′, hence Eq. (2.12) takes the classical form

with the boundary conditions:

The aim of this paper is to study FSLP subjected to a such a set of boundary conditions that leads to a purely discrete countable spectrum and to the orthogonal eigenfunctions’ system constituting the basis in the respective weighted Hilbert space.

3 Exact solutions

In this section, we shall formulate the FSLP with an equation containing the fractional differential operator (2.11). We investigate the eigenvalues and eigenfunctions’ system connected to the FSLE in the case of order α fulfilling condition: 1 ≥ α > 1/2

subject to the mixed boundary conditions in the fractional version and on the space of continuous functions

Let us observe that the above regular FSLP is a special case of the general FSLP given in Definition 2.1 when constants c₂ = d₁ = 0.

We propose a transformation of the introduced differential Sturm-Liouville problem to the integral one on the subspace of continuous functions defined below:

We start by introducing the following integral operator

and we note that on the C_B[a, b], subspace of continuous functions, the following relation is valid

In addition, it is easy to check that

therefore the intermediate integral form of the equation (2.12) looks as follows

In order to invert the integral operator on the left-hand side, we estimate the norm of the T_q operator in the C[a, b] space with the supremum norm ‖ ⋅ ‖ and obtain

Let us denote the parameter ξ as follows:

then the above calculations lead to the following lemma.

When the norm of the T_q-operator is smaller than 1, we can invert operator 1 + T_q.

For functions f ∈ C_B [a, b] we can prove the equivalence of the differential and the integral form of the fractional Sturm-Liouville problem. Namely, the following lemma is valid.

Next, we extend the T-operator to the Lw2 (a, b) - space and note that if order α fulfills 1 ≥ α > 1/2, then

The following lemma is a Straightforward corollary from Lemma 3.2.

In order to consider the spectral properties of the integral operator T, we shall convert it to the integral Hilbert-Schmidt operator. To this aim we explicitly calculate its kernel. First, we express operators T_q and T_w as integral operators with the corresponding kernels. We obtain the following formula for the kernel K_q of the T_q operator defined by Eq. (3.6)

and

with symmetric part

For the T_w operator defined by Eq. (3.4):

we calculate kernel K_w:

We now estimate the norms and values of the above kernels.

From the above lemma we obtain |K_w (x, s)| ≤ ξ_w /(b − a) and K_w ∈ L1⊗12 ([a, b] × [a, b]) by assuming q = w and denoting:

Now, we reformulate the T operator given by Eq. (3.11) in the form of Hilbert-Schmidt operator with kernel K:

The next step is the explicit calculation of the kernel K and examine its properties. It appears that it can be expressed as a series of convolutions described in the lemma below.

Now, we are ready to formulate the main theorem on the integral operator T.

The above theorem on the Hilbert-Schmidt operator T implies the following result on its spectrum and eigenfunctions.

As the eigenfunctions of the operator T belong to the C_B[a, b]-space, we can apply Lemma 3.3 on the equivalence of the differential and integral eigenvalue problem to obtain a principal result on the discrete spectrum of the studied FSLP.

The assumption 1p ∈ C[a, b] leads to two possible cases, namely p can be a positive or negative function. When function p is positive the spectrum is bounded from below which means that the fractional Sturm-Liouville operator has at least a finite number of negative and an infinite number of positive eigenvalues.

4 Numerical solution

We start the construction of numerical scheme by dividing the considered interval [a, b] into N equidistant subintervals of length Δx = (b − a)/N with the central points x_i = a + (i − 0.5)Δx for i = 1,…, N. A value of function y at node x_i we denote as y_i = y(x_i). In addition, we introduce the notation x_{i ± 0.5} = x_i ± 0.5Δx.

The intermediate integral equation Eq. (3.7) can be written as follows

with kernel K₁ given by Eq. (3.27). Next, we approximate the above integrals by the quadrature rule

where u_k are the weights of the quadrature rule (for the midpoint rectangular rule: u_k = Δx).

If we evaluate the equation at every node x_i, i = 1,…, N, then we obtain the following system of N linear algebraic equations

which can be written in the short notation for node values of eigenfunctions looks as follows:

where (K₁)_{i, k} = K₁ (x_i, x_k).

Now, the above system of equations can be rewritten in the matrix form [25]. We introduce:

1. two diagonal matrices Q and W:

   Q=diagq1,q2,...,qN,W=diagw1,w2,...,wN,(4.5)

   where

   qi=1Δx∫xi−0.5xi+0.5qxdx,wi=1Δx∫xi−0.5xi+0.5wxdx,(4.6)

2. the matrix M:

   M=Mi,k=ukK1i,k,(4.7)

   for i = 1,…., N and k = 1,…., N,

3. the vector Y:

   Y=y1,y2,...,yNT.(4.8)

By using these notations, the system of equations (4.4) takes the following matrix form

and after transformations it is a matrix eigenvalue problem:

where