Exact and Numerical Solution of the Fractional Sturm–Liouville Problem with Neumann Boundary Conditions

Abstract

In this paper, we study the fractional Sturm–Liouville problem with homogeneous Neumann boundary conditions. We transform the differential problem to an equivalent integral one on a suitable function space. Next, we discretize the integral fractional Sturm–Liouville problem and discuss the orthogonality of eigenvectors. Finally, we present the numerical results for the considered problem obtained by utilizing the midpoint rectangular rule.

1. Introduction

The Sturm–Liouville boundary value problem (eigenvalue problem) is an important issue in the field of differential equations. These problems can be regular or singular at each endpoint of the considered interval. The classical Sturm–Liouville theory is still a branch of mathematics, where recent investigations yield meaningful results, both theoretical and applicable (compare [1,2] and references therein). Eigenvalue problems arise in a large number of disciplines of sciences and engineering (such as applied mathematics, classical physics or quantum mechanics).

Therefore, one of the most important problems, in the extended Sturm–Liouville theory, including fractional differential operators, is to understand how to to construct fractional analogues of a classical Sturm–Liouville operator and how its spectrum and eigenfunctions behave for various types of operator and boundary conditions (for example, for fractional Dirichlet, Neumann, Robin or mixed conditions).

It is well known from the classical spectral theory that the Dirichlet Laplacian on a bounded domain always has a purely discrete spectrum, while the Neumann Laplacian on a bounded domain may have an essential spectrum if the boundary is not smooth. For this reason, one can say that the Neumann eigenvalue problem is more subtle than the Dirichlet one [3].

This paper is a part of the project where we construct the fractional Sturm–Liouville Equation as an Euler-Lagrange Equation for a suitable action functional, including fractional derivatives. Here, we focus on the examination of the fractional Sturm–Liouville problem (FSLP) in a bounded interval $[a,b]$ with homogeneous fractional Neumann boundary conditions. The problem is analyzed under assumption: $1/2<\alpha \le 1$, where $\alpha$ is the order of fractional derivatives. Let us point out that by applying the fractional calculus and standard fractional derivatives, one can construct different types of the Sturm–Liouville operators and various types of boundary conditions (see, for example, results in papers [4,5,6,7,8,9,10,11,12]).

Fractional eigenvalue problems have also been considered within the framework of tempered fractional calculus [12,13] and conformable fractional calculus [14,15,16,17]. Recently, a fractional Sturm–Liouville operator containing composite fractional derivatives has been proposed in paper [18], and Prabhakar derivatives were applied in the construction of fractional eigenvalue problems subjected to homogeneous Dirichlet or mixed boundary conditions in papers [19,20].

However, it is important to be aware that not every construction of fractional eigenvalue problem leads to a purely discrete real spectrum of the FSLP or orthogonal eigenfunctions’ system, which is complete in the corresponding Hilbert space. For instance, fractional Sturm–Liouville operators proposed in papers [21,22,23,24,25,26,27], where numerical methods were developed to study eigenvalues and eigenfunctions, include only a left-sided fractional derivative/derivatives. Consequently, they do not lead to orthogonal systems of eigenfunctions, and in many cases, the spectrum includes a complex part. A similar remark applies to other models, where one-sided fractional derivatives or a mixture of classical and one-sided fractional derivatives were applied to develop rigorous analytical results for eigenvalue problems with Dirichlet or extended Dirichlet boundary conditions, dependent on the fractional order of derivatives [28,29,30,31]. Some results concerning discrete spectrum and orthogonal eigenfunction systems have been studied in papers [32,33], where an integer order part of the differential operator was supplemented with a composition or difference of the left and right fractional derivatives, respectively. Therefore, the motivation of our variational approach, first presented in papers [5,34], is to develop the formulation of a fractional Sturm–Liouville problem with an orthogonal system of eigenfunctions and real eigenvalues.

In this approach, the Sturm–Liouville operator contains both the left and right fractional derivatives [35,36], and equations containing this type of differential operators are known as the fractional Euler–Lagrange equations [37,38,39]. The application of this type of operator causes many difficulties in solving fractional differential equations. The same problem appears in calculating eigenvalues. It should be highlighted that for eigenvalue problems with differential operators, being a composition of the left and right fractional derivatives, only a few exact solutions are available [6,9,11,12,34]. For this reason, many scientists have been working on numerical methods dedicated to the fractional Euler-Lagrange equations [40,41,42,43,44] and/or the fractional Sturm–Liouville problems [45,46,47,48].

In our previous papers, we studied fractional eigenvalue problems with homogeneous Dirichlet [20,49,50], Neumann [51], Robin [52] and mixed [47,48] boundary conditions. In these papers, it was shown that the FSLP with the adequate boundary conditions has a purely discrete spectrum, and the corresponding system of eigenfunctions is orthogonal and complete in a suitable Hilbert space.

In the paper [48], we proposed the numerical method, based on the transformation of the FSLP, subjected to the homogeneous mixed boundary conditions into the equivalent integral eigenvalue problem. The integral operator with a kernel depending on both the form of the fractional differential operator and on boundary conditions was discretized to calculate approximate eigenvalues and eigenvectors. In addition, the presented approach gave us a possibility to approximate eigenfunctions, keeping the orthogonality of the eigenvectors and the approximated eigenfunctions at each step of the algorithm. We used the experimental rate of convergence to control the convergence of the developed numerical scheme, and we obtained a rate of convergence close to $2\alpha$.

In the present paper, the analogous technique is developed for FSLP with Neumann boundary conditions. However, when the homogeneous, fractional Neumann boundary conditions are applied on both boundaries of the domain, an additional integral constraint is necessary. The additional integral condition causes the development of a numerical scheme to become slightly more complex, as compared to the previous case (FSLP with mixed boundary conditions). In particular, two numerical schemes are studied—in the first, we use the same discretization procedure for both integrals, i.e., for an integral determining the operator and for an integral appearing in the definition of an integral constraint. In turn, in the second scheme, which we refer to as a hybrid numerical scheme, two different discretization procedures are allowed. The main results of the paper include two numerical schemes dedicated to FSLP with fractional Neumann boundary conditions, the results of the analysis of the eigenvectors’ orthogonality and examples of numerical solutions.

The paper is structured as follows: Section 2 deals with the fundamentals of FSLP theory. Section 3 describes two types of constructions of the discrete versions of integral FSLP. Section 4 presents results for the numerical solution of FSLP with homogeneous Neumann boundary conditions. Finally, Section 5 concludes the paper.

2. Preliminaries

Let us recall definitions of fractional operators, based on the following books [36,53]. First, the left and right Riemann–Liouville fractional derivatives are defined as:

$$\begin{array}{ccc}\hfill D_{a^{+}}^{\alpha }y\left(x\right)& :=& \frac{d}{dx}{I}_{a^{+}}^{1-\alpha }y\left(x\right),\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill D_{b^{-}}^{\alpha }y\left(x\right)& :=& -\frac{d}{dx}{I}_{b^{-}}^{1-\alpha }y\left(x\right),\hfill \end{array}$$

(2)

where integral operators $I_{a^{+}}^{\alpha }$ and $I_{b^{-}}^{\alpha }$ are the left and right fractional Riemann–Liouville integrals:

$$\begin{array}{ccc}\hfill I_{a^{+}}^{\alpha }y\left(x\right)& :=& \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}\underset{a}{\overset{x}{\int }}\frac{y\left(\tau \right)}{{\left(x-\tau \right)}^{1-\alpha }}d\tau ,x>a,\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill I_{b^{-}}^{\alpha }y\left(x\right)& :=& \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}\underset{x}{\overset{b}{\int }}\frac{y\left(\tau \right)}{{\left(\tau -x\right)}^{1-\alpha }}d\tau ,x<b.\hfill \end{array}$$

(4)

Then, the left and right Caputo fractional derivatives are defined as follows:

$$\begin{array}{ccc}\hfill {}^C{D}_{a^{+}}^{\alpha }y\left(x\right)& :=& D_{a^{+}}^{\alpha }(y\left(x\right)-y\left(a\right)),\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill {}^C{D}_{b^{-}}^{\alpha }y\left(x\right)& :=& D_{b^{-}}^{\alpha }(y\left(x\right)-y\left(b\right)).\hfill \end{array}$$

(6)

In this preliminary part of the paper, we shall report previous results enclosed in papers [34,51], relevant to developing a numerical method of solution of a fractional eigenvalue problem in the case when the solutions’ space is restricted by Neumann-type boundary conditions. Let us quote the general formulation of a regular fractional Sturm–Liouville problem (FSLP) introduced in [5,34].

Definition 1.

[Compare with Definition 5 in [34]. Let $\alpha \in (0,1]$. With the notation:

$${\mathcal{L}}_q:=D_{b-}^{\alpha }p\left(x\right){}^C{D}_{a+}^{\alpha }+q\left(x\right),$$

(7)

consider the fractional Sturm–Liouville Equation (FSLE):

$${\mathcal{L}}_q{y}_{\lambda }\left(x\right)=\lambda w\left(x\right)y_{\lambda }\left(x\right),$$

(8)

where $p\left(x\right)\ne 0,w\left(x\right)>0\forall x\in [a,b]$, functions $p,q,w$ are real-valued continuous functions in $[a,b]$ and boundary conditions are:

$$\begin{array}{ccc}& c_1{y}_{\lambda }\left(a\right)+c_2{I}_{b-}^{1-\alpha }p\left(x\right){}^C{D}_{a+}^{\alpha }{y}_{\lambda }\left(x\right){\mid }_{x=a}=0,& \hfill \end{array}$$

(9)

$$\begin{array}{ccc}& d_1{y}_{\lambda }\left(b\right)+d_2{I}_{b-}^{1-\alpha }p\left(x\right){}^C{D}_{a+}^{\alpha }{y}_{\lambda }\left(x\right){\mid }_{x=b}=0& \hfill \end{array}$$

(10)

with $c_1^2+c_2^2\ne 0$ and $d_1^2+d_2^2\ne 0.$ The problem of finding number λ (eigenvalue) such that the BVP has a non-trivial solution, $y_{\lambda }$ (eigenfunction) will be called the regular fractional Sturm–Liouville eigenvalue problem (FSLP).

In [34], we proved that eigenvalues generated while solving the above Sturm–Liouville problem are real, and eigenfunctions associated to distinct eigenvalues are orthogonal with respect to the following scalar product:

$${\langle f,g\rangle }_w:={\int }_a^{b}w\left(x\right)\overline{f\left(x\right)}g\left(x\right)dx.$$

Let us observe that for $\alpha =1$, we recover the classical Sturm–Liouville problem (CSLP), where operator (7) becomes the second order differential operator, and Equation (8) is the classical Sturm–Liouville equation:

$$-\frac{d}{dx}p\left(x\right)\frac{d{y}_{\lambda }\left(x\right)}{dx}+q\left(x\right)y_{\lambda }\left(x\right)=\lambda w\left(x\right)y_{\lambda }\left(x\right)$$

with boundary conditions appearing as follows:

$$c_1{y}_{\lambda }\left(a\right)+c_2\frac{d{y}_{\lambda }\left(a\right)}{dx}=0,d_1{y}_{\lambda }\left(b\right)+d_2\frac{d{y}_{\lambda }\left(b\right)}{dx}=0.$$

(11)

When we choose $c_1=d_1=0$ in Equation (11), determining the boundary conditions, we arrive at CSLP with homogeneous Neumann boundary conditions:

$$\frac{d{y}_{\lambda }\left(a\right)}{dx}=\frac{d{y}_{\lambda }\left(b\right)}{dx}=0.$$

The same choice of constants in the general fractional boundary conditions (9), (10) yields the fractional version of homogeneous Neumann boundary conditions. In this way, we restrict the general regular FSLP to the case with fractional homogeneous Neumann boundary conditions (14), (15), described in the definition below:

Definition 2.

Let $\alpha \in (0,1]$. With the following notation:

$${\mathcal{L}}_q:=D_{b-}^{\alpha }p\left(x\right){}^C{D}_{a+}^{\alpha }+q\left(x\right),$$

(12)

consider the fractional Sturm–Liouville equation:

$${\mathcal{L}}_q{y}_{\lambda }\left(x\right)=\lambda w\left(x\right)y_{\lambda }\left(x\right),$$

(13)

where $p\left(x\right)\ne 0,w\left(x\right)>0\forall x\in [a,b]$, functions $p,q,w$ are real-valued continuous functions in $[a,b]$ and boundary conditions are:

$$\begin{array}{ccc}& I_{b-}^{1-\alpha }p\left(x\right){}^C{D}_{a+}^{\alpha }{y}_{\lambda }\left(x\right){\mid }_{x=a}=0,& \hfill \end{array}$$

(14)

$$\begin{array}{ccc}& I_{b-}^{1-\alpha }p\left(x\right){}^C{D}_{a+}^{\alpha }{y}_{\lambda }\left(x\right){\mid }_{x=b}=0.& \hfill \end{array}$$

(15)

In the problem of finding number λ (eigenvalue) such that the BVP has a non-trivial solution, $y_{\lambda }$ (eigenfunction) will be called the regular fractional Sturm–Liouville eigenvalue problem with homogeneous fractional Neumann boundary conditions (FSLPN).

We focus the further investigations of numerical solutions to FSLPN, defined above, in the case when $q=0$, i.e., the fractional Sturm–Liouville operator (FSLO) is ${\mathcal{L}}_0$:

$${\mathcal{L}}_0:=D_{b-}^{\alpha }p\left(x\right){}^C{D}_{a+}^{\alpha }$$

(16)

and we restrict fractional order to $\alpha \in (1/2,1]$. These assumptions are motivated by the fact that, in such case, we can apply analytical results obtained in the paper [51]. These results are based on the transformation of the differential FSLPN into the equivalent integral problem. Moreover, they include theorems on purely discrete spectrum of differential and integral fractional Sturm–Liouville operators. The schedule of deriving the analytical spectral result for FSLPN contains the following preliminary steps: construction of the suitable solutions’ space, transformation to the integral FSLPN and proof of equivalence of both types of eigenvalue problems.

First, when we restrict our considerations to continuous solutions, we arrive at the following relation resulting from Equation (8) (for $q=0$) and composition properties of fractional derivatives and integrals:

$${\mathcal{L}}_0\left[y_{\lambda }\left(x\right)-I_{a+}^{\alpha }\frac{1}{p\left(x\right)}{I}_{b-}^{\alpha }\lambda w\left(x\right)y_{\lambda }\left(x\right)\right]=0.$$

(17)

Similar to the case of FSLP subjected to the homogeneous mixed, Dirichlet, and Robin boundary conditions [20,47,48,52], we note that while assuming the continuity of $y_{\lambda }$ and $1\ge \alpha >1/2$, this eigenfunction fulfills the following integral Equation of interval $[a,b]$:

$$y_{\lambda }\left(x\right)-I_{a+}^{\alpha }\frac{1}{p\left(x\right)}{I}_{b-}^{\alpha }\lambda w\left(x\right)y_{\lambda }\left(x\right)=A_1+A_2{I}_{a+}^{\alpha }\frac{{(b-x)}^{\alpha -1}}{p\left(x\right)},$$

(18)

where the function on the right-hand side of Equation (18) is an arbitrary continuous function from the kernel of operator ${\mathcal{L}}_0$. Therefore, in fact, we work on the subspace of continuous functions, defined previously in [51] and described in the definition below.

Definition 3.

Let $C_p[a,b]\subset C[a,b]$ be the subspace of continuous functions given as follows:

$$f\in C_p[a,b]\iff f\left(x\right)=I_{a+}^{\alpha }\frac{1}{p\left(x\right)}{I}_{b-}^{\alpha }\psi \left(x\right)+A_1+A_2{I}_{a+}^{\alpha }\frac{{(b-x)}^{\alpha -1}}{p\left(x\right)},$$

(19)

where $x\in [a,b]$, $\psi \in C[a,b]$, and $A_1,A_2$ are constants.

On the considered function space, $f\left(a\right)=A_1$, therefore, the fractional Caputo derivative ${}^C{D}_{a+}^{\alpha }f$ exists for every point in $[a,b)$. Now, we restrict the defined subspace to functions fulfilling the homogeneous fractional Neumann boundary conditions. In addition, an integral condition results from initial FSLE (8) and assumed boundary conditions. This property of the solutions’ space is also necessary in the classical case $\alpha =1$.

Definition 4.

Let $C_{p,N}[a,b]\subset C_p[a,b]$ be the subspace of continuous functions given in Definition 3, fulfilling homogeneous, fractional Neumann boundary conditions (14), (15) and the integral condition that:

$${\int }_a^{b}w\left(x\right)f\left(x\right)dx=0.$$

(20)

It is easy to verify that constants $A_1$ and $A_2$ in Equation (18) are determined by the homogeneous fractional Neumann boundary conditions (14), (15) and condition (20). Namely, we have:

$$A_1=-\lambda \frac{{\int }_a^{b}w\left(x\right)\left(I_{a+}^{\alpha }\frac{1}{p\left(x\right)}{I}_{b-}^{\alpha }w\left(x\right)y_{\lambda }\left(x\right)\right)dx}{{\int }_a^{b}w\left(\xi \right)y_{\lambda }\left(\xi \right)d\xi },$$

(21)

$$A_2=0.$$

(22)

Thence, we obtain the integral Equation (18) in the exact form:

$$y_{\lambda }\left(x\right)=\lambda T_w{y}_{\lambda }\left(x\right):=\lambda I_{a+}^{\alpha }\frac{1}{p\left(x\right)}{I}_{b-}^{\alpha }w\left(x\right)y_{\lambda }\left(x\right)$$

(23)

$$-\lambda \frac{{\int }_a^{b}w\left(x\right)\left(I_{a+}^{\alpha }\frac{1}{p\left(x\right)}{I}_{b-}^{\alpha }w\left(x\right)y_{\lambda }\left(x\right)\right)dx}{{\int }_a^{b}w\left(\xi \right)y_{\lambda }\left(\xi \right)d\xi }.$$

The integral operator $T_w$, given above as the composition of the left and right Riemann–Liouville integrals, can be expressed as an integral operator with kernel $K_w$ [51]:

$$y_{\lambda }\left(x\right)=\lambda T_w{y}_{\lambda }\left(x\right)=\lambda \underset{a}{\overset{b}{\int }}{K}_w\left(x,s\right)y_{\lambda }\left(s\right)ds,$$

(24)

where kernel $K_w$ is given below:

$$K_w\left(x,s\right)=w\left(s\right)\left(K_1\left(x,s\right)-\frac{{\displaystyle \underset{a}{\overset{b}{\int }}w\left(\xi \right)K_1\left(\xi ,s\right)d\xi }}{{\displaystyle \underset{a}{\overset{b}{\int }}w\left(\xi \right)d\xi }}\right)$$

(25)

and symmetric kernel part $K_1$ is of the form:

$$K_1\left(x,s\right)=\frac{1}{{\mathsf{\Gamma }}^2\left(\alpha \right)}\underset{a}{\overset{min\left(x,s\right)}{\int }}\frac{{\left(x-t\right)}^{\alpha -1}{\left(s-t\right)}^{\alpha -1}}{p\left(t\right)}dt.$$

(26)

The above explicit expressions for kernels $K_1$ and $K_w$ result from the application of the integration rule-change of the order of integration. We note that kernel $K_w$, given in Equation (25), fulfills the condition:

$${\int }_a^b{K}_w(x,s)w\left(x\right)dx=0$$

(27)

which yields the properties:

$$f\in C_{p,N}[a,b]⟹T_{w}f\in C_{p,N}[a,b]$$

(28)

and

$$f\in \widetilde{L_w^2}(a,b)⟹T_{w}f\in C_{p,N}[a,b],$$

(29)

where Hilbert space $\widetilde{L_w^2}(a,b)$ is defined by Equation (33).

We start proving the above properties by observation that in the case when the weight function $w\in C[a,b]$ and the fractional order fulfills $\alpha \in (1/2,1]$ we have for any function $f\in L_w^2(a,b)$ that the results of fractional integration are continuous functions, i.e., $I_{a+}^{\alpha }f\in C[a,b]$ and $I_{b-}^{\alpha }f\in C[a,b]$. Then, we apply composition rules for fractional operators [36,51] and Formula (23) to check Neumann boundary conditions (14), (15) for the respective images of continuous or Lebegue functions:

$$I_{b-}^{1-\alpha }p\left(x\right){}^C{D}_{a+}^{\alpha }{T}_{w}f\left(x\right){\mid }_{x=a}=$$

$$=I_{b-}^{1-\alpha }p\left(x\right){}^C{D}_{a+}^{\alpha }\left(I_{a+}^{\alpha }\frac{1}{p\left(x\right)}{I}_{b-}^{\alpha }w\left(x\right)f\left(x\right)-\frac{{\int }_a^{b}w\left(x\right)\left(I_{a+}^{\alpha }\frac{1}{p\left(x\right)}{I}_{b-}^{\alpha }w\left(x\right)f\left(x\right)\right)dx}{{\int }_a^{b}w\left(\xi \right)y_{\lambda }\left(\xi \right)d\xi }\right){\mid }_{x=a}$$

$$=I_{b-}^{1}w\left(x\right)f\left(x\right)dx{\mid }_{x=a}={\int }_a^{b}w\left(x\right)f\left(x\right)dx=0$$

because, by assumption, function $f\in C_{p,N}[a,b]$ or $f\in \widetilde{L_w^2}(a,b)$. Similarly, we obtain at the endpoint of the interval for continuous or Lebegue function f the following result:

$$I_{b-}^{1-\alpha }p\left(x\right){}^C{D}_{a+}^{\alpha }{T}_{w}f\left(x\right){\mid }_{x=b}=$$

$$=I_{b-}^{1-\alpha }p\left(x\right){}^C{D}_{a+}^{\alpha }\left(I_{a+}^{\alpha }\frac{1}{p\left(x\right)}{I}_{b-}^{\alpha }w\left(x\right)f\left(x\right)-\frac{{\int }_a^{b}w\left(x\right)\left(I_{a+}^{\alpha }\frac{1}{p\left(x\right)}{I}_{b-}^{\alpha }w\left(x\right)f\left(x\right)\right)dx}{{\int }_a^{b}w\left(\xi \right)y_{\lambda }\left(\xi \right)d\xi }\right){\mid }_{x=b}$$

$$=I_{b-}^{1-\alpha }\left(I_{b-}^{\alpha }w\left(x\right)f\left(x\right)\right){\mid }_{x=b}=0$$

as for order $\alpha \in (1/2,1]$ and $w\in C[a,b]$ we have $I_{b-}^{\alpha }wf\in C[a,b]$.

In summary, we observe that images $T_{w}f$ fulfill Neumann boundary conditions (14), (15) for any function from $C_{p,N}$ or $\widetilde{L_w^2}(a,b)$ spaces. Now, let us check the integral condition (20) for image $T_{w}f$. We apply kernel property (27) and obtain for arbitrary $f\in C_{p,N}$ or $f\in \widetilde{L_w^2}(a,b)$ the following integral formula after the change in the order of integration:

$${\int }_a^{b}w\left(x\right)T_{w}f\left(x\right)dx={\int }_a^{b}w\left(x\right)\left({\int }_a^b{K}_w(x,s)f\left(s\right)ds\right)dx$$

$$={\int }_a^{b}f\left(s\right)\left({\int }_a^b{K}_w(x,s)w\left(x\right)dx\right)ds=0.$$

The above calculations show that properties (28), (29) are valid.

The following two lemmas lead to the main result on spectral properties of the considered FSLPN. The first one says that on the constructed solutions’ space operator, $T_w$ is an inverse operator to the fractional Sturm–Liouville operator ${\frac{1}{w}\mathcal{L}}_0$, and it results from work performed in the previous paper [51] (compare relations in Equations (23) and (24), and the corresponding calculations).

Lemma 1.

Let $1\ge \alpha >1/2,\frac{1}{p},w\in C[a,b]$ and $f\in C_{p,N}[a,b]$. Then, the following relations are valid for any $x\in [a,b]$:

$$T_w\frac{1}{w\left(x\right)}{\mathcal{L}}_{0}f\left(x\right)=f\left(x\right),$$

(30)

$$\frac{1}{w\left(x\right)}{\mathcal{L}}_0{T}_{w}f\left(x\right)=f\left(x\right).$$

(31)

The next lemma on the equivalence of the integral and differential FSLPN is a straightforward corollary of Equations (30) and (31).

Lemma 2.

[Compare Lemma 1 in [51]]. Let $1\ge \alpha >1/2,\frac{1}{p},w\in C[a,b]$ and $f\in C_{p,N}[a,b]$.

Then, the following equivalence is valid:

$${\mathcal{L}}_{0}f\left(x\right)=\lambda w\left(x\right)f\left(x\right)⟺T_{w}f\left(x\right)=\frac{1}{\lambda }f\left(x\right).$$

(32)

In the paper [51], we studied properties of the $T_w$ operator and derived the following result on its spectrum. We denote the Hilbert space $\widetilde{L_w^2}(a,b)$ as follows:

$$\widetilde{L_w^2}(a,b):=\{f\in L_w^2(a,b);{\int }_a^{b}w\left(x\right)f\left(x\right)dx=0\}$$

(33)

which means that it is a subspace of $L_w^2(a,b)$ containing functions fulfilling condition (20). Now, let us quote the theorem on the spectrum of integral and differential FSLP with Neumann boundary conditions (14), (15).

Theorem 1.

Let $1\ge \alpha >\frac{1}{2}$ and $\frac{1}{p},w\in C[a,b]$. Then, operator $T_w$, given in (23), (24), has a purely discrete spectrum enclosed in interval $(-1,1)$, with 0 being its only limit point. Eigenfunctions $y_n$ corresponding to the respective eigenvalues belong to the $C_{p,N}[a,b]$-space and form a basis in the $\widetilde{L_w^2}(a,b)$-space (33).

Finally, we recall that the above theorem, together with Lemma 2, yields the result on the spectrum and eigenfunctions for the differential FSLP with homogeneous Neumann type boundary conditions.

Theorem 2.

Let $1\ge \alpha >\frac{1}{2}$ and $\frac{1}{p},w\in C[a,b]$. Then, operator ${\mathcal{L}}_0$, given in (16), considered on the $C_{p,N}[a,b]$-space, has a purely discrete real spectrum with $|{\lambda }_n|\to \infty$. Eigenfunctions $y_n$ form a basis in the $\widetilde{L_w^2}(a,b)$-space (33). For positive function $p>0$, the spectrum is positive, whereas for negative $p<0$, it is negative. Moreover, the following number series is convergent:

$$\sum _{n=1}^{\infty }\frac{1}{{\left({\lambda }_n\right)}^2}<\infty$$

(34)

and the inequality below is fulfilled for certain $M_{+}>0$:

$$\frac{|y_n\left(x\right)|}{|{\lambda }_n|}\le M_{+},x\in [a,b],n\in \mathbf{N}.$$

(35)

Having summarized the approach, developed in [51], to study properties and the spectrum of FSLO (16) on function space subject to the homogeneous fractional Neumann boundary conditions, we are now ready to present the numerical methods of solution to FSLPN.

3. Main Results

In this section, two approaches to the construction of the discrete version of integral FSLPN are presented. We discuss an influence of each of the proposed methods on the ortogonality property of eigenvectors of discrete integral operator.

3.1. Case I—One Numerical Integration Scheme in the Construction of Discrete Integral FSLPN

Let us observe that passing to the discrete form of the integral eigenvalue problem (24), we construct discrete versions of integrals describing the $T_w$-operator itself (24), its kernel (25) and its integral condition (20). First, we shall consider the case when we use the same numerical procedure for integrals in (20) and (24). We divide the interval $[a,b]$ into N subintervals and choose arbitrary points $x_i,i=1,$ …$,N^{\prime }$ according to the applied rule of numerical integration. In particular, for an equidistant partition into N subintervals, we have $N^{\prime }=N$ for rectangular approximation, $N^{\prime }=N+1$ for the trapezoid method while Simpson’s method requires $N^{\prime }=N+1$ (N must be even). In general, the numerical integral is given as:

$${\int }_a^{b}y\left(x\right)dx\approx \sum _{i=1}^{N^{\prime }}{\beta }_{i}y\left(x_i\right)$$

(36)

and the $T_w$-integral operator acts as follows:

$$T_{w}y\left(x_i\right)=\sum _{j=1}^{N^{\prime }}{\beta }_j{K}_w(x_i,x_j)y\left(x_j\right),$$

(37)

where approximation weights ${\beta }_i,i=1,$ …$,N^{\prime }$ depend on the specific choice of points $x_i,i=1,$ …$,N^{\prime }$ corresponding to the applied method of numerical integration. Denoting values of function $y\left(x_i\right)=Y_i$ and vector $Y=(Y_1,Y_2,$ …$,Y_{N^{\prime }}{)}^T$, we rewrite the above Equation to the discrete form, where operator ${\mathbf{T}}_w$ acts on vectors from $N^{\prime }$-dimensional vector space:

$${\left({\mathbf{T}}_{w}Y\right)}_i=\sum _{j=1}^{N^{\prime }}{\beta }_j{K}_w(x_i,x_j)Y_j$$

(38)

subjected to the discrete version of condition (20):

$$\sum _{i=1}^{N^{\prime }}{\beta }_i{w}_i{Y}_i=0.$$

(39)

At this point, we write the discrete version of the kernel, assuming that the integrals on the right-hand side of equality (27) are discretized according to Equation (36):

$$K_w(x_i,x_j)=w_j{K}_1(x_i,x_j)-\frac{w_j{\sum }_{k=1}^{N^{\prime }}{\beta }_k{w}_k{K}_1(x_k,x_j)}{{\sum }_{l=1}^{N^{\prime }}{\beta }_l{w}_l}.$$

(40)

We note that the discrete analog of the kernel property (27) is fulfilled:

$$\sum _{i=1}^{N^{\prime }}{\beta }_i{w}_i{K}_w(x_i,x_j)=0$$

(41)

and that the image of the vector fulfilling the condition (39) obeys the same condition:

$$\sum _{i=1}^{N^{\prime }}{\beta }_i{w}_i{\left({\mathbf{T}}_{w}Y\right)}_i=0.$$

(42)

This property of the discrete integral operator ${\mathbf{T}}_w$ can be easily checked:

$$\sum _{i=1}^{N^{\prime }}{\beta }_i{w}_i{\left({\mathbf{T}}_{w}Y\right)}_i=\sum _{i=1}^{N^{\prime }}{\beta }_i{w}_i\left(\sum _{j=1}^{N^{\prime }}{\beta }_j{K}_w(x_i,x_j)Y_j\right)$$

$$=\sum _{j=1}^{N^{\prime }}{\beta }_j{Y}_j\left(\sum _{i=1}^{N^{\prime }}{\beta }_i{w}_i{K}_w(x_i,x_j)\right)=0.$$

In conclusion, we consider the discrete eigenvalue problem:

$$Y_{\lambda }=\lambda {\mathbf{T}}_w{Y}_{\lambda },{\left(Y_{\lambda }\right)}_i=\lambda \sum _{j=1}^{N^{\prime }}{\beta }_j{K}_w(x_i,x_j){\left(Y_{\lambda }\right)}_j$$

(43)

on the $N^{\prime }$-dimensional space of vectors fulfilling condition (39). In the next step, we shall study the presented discrete eigenvalue problem, where all the discrete analogs of integrals from the fractional integral eigenvalue problem are constructed by using the same numerical integration method. In addition, the scalar product on the vector space is also defined according to the discretization procedure (36) and weight w:

$${\langle F,G\rangle }_w:=\sum _{i=1}^{N^{\prime }}{\beta }_i{w}_i{F}_i{G}_i.$$

(44)

The following lemma describes the orthogonality property of eigenvectors of discrete operator ${\mathbf{T}}_w$.

Lemma 3.

Eigenvectors corresponding to distinct eigenvalues of the discrete FSLP (43) are orthogonal with respect to scalar product (44), i.e.,

$${\langle Y_{\lambda },Y_{\rho }\rangle }_w=0\lambda \ne \rho .$$

(45)

Proof. 

In the proof, we apply condition (39) and the symmetricity of kernel part $K_1$:

$$\frac{1}{\rho }{\langle Y_{\lambda },Y_{\rho }\rangle }_w={\langle Y_{\lambda },{\mathbf{T}}_w{Y}_{\rho }\rangle }_w=\sum _{i=1}^{N^{\prime }}{\beta }_i{w}_i{\left(Y_{\lambda }\right)}_i{\left({\mathbf{T}}_w{Y}_{\rho }\right)}_i=$$

$$=\sum _{i=1}^{N^{\prime }}{\beta }_i{w}_i{\left(Y_{\lambda }\right)}_i\left(\sum _{j=1}^{N^{\prime }}{\beta }_j{K}_w(x_i,x_j){\left(Y_{\rho }\right)}_j\right)$$

$$=\sum _{i=1}^{N^{\prime }}{\beta }_i{w}_i{\left(Y_{\lambda }\right)}_i\left[\sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j{\left(Y_{\rho }\right)}_j\left(K_1(x_i,x_j)-\frac{{\sum }_{k=1}^{N^{\prime }}{\beta }_k{w}_k{K}_1(x_k,x_j)}{{\sum }_{l=1}^{N^{\prime }}{\beta }_l{w}_l}\right)\right]$$

$$=\sum _{j=1}^{N^{\prime }}\sum _{i=1}^{N^{\prime }}{\beta }_i{w}_i{\left(Y_{\lambda }\right)}_i{\beta }_j{w}_j{\left(Y_{\rho }\right)}_j{K}_1(x_i,x_j)$$

$$=\sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j{\left(Y_{\rho }\right)}_j\left(\sum _{i=1}^{N^{\prime }}{\beta }_i{w}_i{\left(Y_{\lambda }\right)}_i{K}_1(x_j,x_i)\right)$$

$$=\sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j{\left(Y_{\rho }\right)}_j\left[\sum _{i=1}^{N^{\prime }}{\beta }_i{\left(Y_{\lambda }\right)}_i\left(w_i{K}_1(x_j,x_i)-\frac{w_i{\sum }_{k=1}^{N^{\prime }}{\beta }_k{w}_k{K}_1(x_k,x_i)}{{\sum }_{l=1}^{N^{\prime }}{\beta }_l{w}_l}\right)\right]$$

$$=\sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j{\left(Y_{\rho }\right)}_j\left(\sum _{i=1}^{N^{\prime }}{\beta }_i{\left(Y_{\lambda }\right)}_i{K}_w(x_j,x_i)\right)$$

$$=\sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j{\left(Y_{\rho }\right)}_j{\left({\mathbf{T}}_w{Y}_{\lambda }\right)}_j={\langle Y_{\rho },{\mathbf{T}}_w{Y}_{\lambda }\rangle }_w=\frac{1}{\lambda }{\langle Y_{\lambda },Y_{\rho }\rangle }_w.$$

Finally, we obtain:

$$\left(\frac{1}{\rho }-\frac{1}{\lambda }\right){\langle Y_{\lambda },Y_{\rho }\rangle }_w=0$$

and as the eigenvalues are distinct, we end the proof:

$${\langle Y_{\lambda },Y_{\rho }\rangle }_w=0.$$

□

3.2. Case II—Hybrid Numerical Integration Scheme in the Construction of Discrete Integral FSLPN

Now, we shall investigate an influence of introducing two numerical integration schemes in the construction of the discrete version of integral FSLPN. Again, we divide the interval into N subintervals and choose points $S_1=\{{\tilde{x}}_1,$ …$,{\tilde{x}}_{N_1}\}$ according to the rule of numerical integration which will be applied in the construction of the ${\mathbf{T}}_w$ operator and the set of points (some may coincide with the previous ones) $S_2=\{{\overline{x}}_1,$ …$,{\overline{x}}_{N_2}\}$ corresponding to numerical integration for kernel (25). In addition, we assume $S_1\subseteq S_2$ and denote the respective weights as $\tilde{\beta }=({\tilde{\beta }}_1,$ …$,{\tilde{\beta }}_{N_1})$ and $\overline{\gamma }=({\overline{\gamma }}_1,$ …$,{\overline{\gamma }}_{N_2})$. Then, the numerical versions of integral over interval $[a,b]$ are:

$${\int }_a^{b}y\left(x\right)dx\approx \sum _{i=1}^{N_1}{\tilde{\beta }}_{i}y\left({\tilde{x}}_i\right),$$

(46)

$${\int }_a^{b}y\left(x\right)dx\approx \sum _{i=1}^{N_2}{\overline{\gamma }}_{i}y\left({\overline{x}}_i\right).$$

(47)

The discrete form of integral FSLP looks as follows:

$$y\left({\overline{x}}_i\right)=\lambda {\mathbf{T}}_{w}y\left({\overline{x}}_i\right){\mathbf{T}}_{w}y\left({\overline{x}}_i\right)=\sum _{j=1}^{N_1}{\tilde{\beta }}_j{K}_w({\overline{x}}_i,{\tilde{x}}_j)y\left({\tilde{x}}_j\right),i=1,\cdots ,N_2$$

(48)

where kernel is ($i=1,$ …$,N_2$ and $j=1,$ …$,N_1$)

$$K_w({\overline{x}}_i,{\tilde{x}}_j)={\tilde{w}}_j{K}_1({\overline{x}}_i,{\tilde{x}}_j)-\frac{{\tilde{w}}_j{\sum }_{k=1}^{N_2}{\overline{\gamma }}_k{\overline{w}}_k{K}_1({\overline{x}}_k,{\tilde{x}}_j)}{{\sum }_{l=1}^{N_2}{\overline{\gamma }}_l{\overline{w}}_l}.$$

(49)

In the case $S_1=S_2$, we have $N_1=N_2$. Therefore, the matrix of operator ${\mathbf{T}}_w$ is a quadratic one, and set (48) simultaneously yields the eigenvalues and eigenvectors of the discrete version of the integral FSLPN. In turn, when we work under assumption $S_1\subset S_2$, we note that the matrix of operator ${\mathbf{T}}_w$ is a rectangular one, and the set of Equation (48) should be split into two parts. The first is:

$$y\left({\overline{x}}_i\right)=\lambda {\mathbf{T}}_{w}y\left({\overline{x}}_i\right){\mathbf{T}}_{w}y\left({\overline{x}}_i\right)=\sum _{j=1}^{N_1}{\tilde{\beta }}_j{K}_w({\overline{x}}_i,{\tilde{x}}_j)y\left({\tilde{x}}_j\right),i=1,\cdots ,N_1$$

(50)

and we call it a reduced discrete integral FSLPN. It provides $N_1$ non-zero eigenvalues of discrete integral FSLPN and $N_1$ eigenvectors up to the first $N_1$ coordinates. The remaining equations of set (48) allow us to calculate eigenvector coordinates for $i=N_1+1,$ …$,N_2$. We reformulate numerical integration schemes by joining all the chosen points in set $S=S_1\cup S_2=S_2$, which means:

$$x_i={\tilde{x}}_i={\overline{x}}_i{x}_i\in S_1\wedge x_i={\tilde{x}}_i{x}_i\notin S_1.$$

Then, we construct the respective extended weights as follows:

$${\beta }_i={\overline{\beta }}_i{x}_i\in S_1\wedge {\beta }_i=0x_i\notin S_1,$$

$${\gamma }_i={\tilde{\gamma }}_i{x}_i\in S_2\wedge {\gamma }_i=0x_i\notin S_2.$$

In general, the numerical integration rules (46) and (47) now become:

$${\int }_a^{b}y\left(x\right)dx\approx \sum _{i=1}^{N_1}{\tilde{\beta }}_{i}y\left({\tilde{x}}_i\right)=\sum _{i=1}^{N^{\prime }}{\beta }_{i}y\left(x_i\right),$$

(51)

$${\int }_a^{b}y\left(x\right)dx\approx \sum _{i=1}^{N_2}{\overline{\gamma }}_{i}y\left({\overline{x}}_i\right)=\sum _{i=1}^{N^{\prime }}{\gamma }_{i}y\left(x_i\right).$$

(52)

The $T_w$-integral operator, rewritten by using numerical integration rule (51), is:

$${\mathbf{T}}_{w}y\left(x_i\right)=\sum _{j=1}^{N^{\prime }}{\beta }_j{K}_w(x_i,x_j)y\left(x_j\right),$$

(53)

where weights ${\beta }_i,i=1,$ …$,N^{\prime }$ are described above. The discrete extended version of kernel appears as follows:

$$K_w(x_i,x_j)=w_j{K}_1(x_i,x_j)-\frac{w_j{\sum }_{k=1}^{N^{\prime }}{\gamma }_k{w}_k{K}_1(x_k,x_j)}{{\sum }_{l=1}^{N^{\prime }}{\gamma }_l{w}_l},$$

(54)

where we assume that the integrals on the right-hand side of equality (25) are discretized according to rule (52). Similar to the previous scheme, the discrete analog of kernel property (27) is fulfilled:

$$\sum _{i=1}^{N^{\prime }}{\gamma }_i{w}_i{K}_w(x_i,x_j)=0j=1,\cdots ,N^{\prime }.$$

(55)

This property of kernel leads to the $N^{\prime }$-dimensional vector space of solutions obeying the numerical version of condition (20) in the form of:

$$\sum _{i=1}^{N^{\prime }}{\gamma }_i{w}_i{Y}_i=0.$$

(56)

Similar to the scheme discussed previously, the image of the vector obeying condition (56) obeys the same condition:

$$\sum _{i=1}^{N^{\prime }}{\gamma }_i{w}_i{\left({\mathbf{T}}_{w}Y\right)}_i=0$$

(57)

and this property of the discrete integral operator ${\mathbf{T}}_w$ is a straightforward corollary of property (55):

$$\sum _{i=1}^{N^{\prime }}{\gamma }_i{w}_i{\left({\mathbf{T}}_{w}Y\right)}_i=\sum _{i=1}^{N^{\prime }}{\gamma }_i{w}_i\left(\sum _{j=1}^{N^{\prime }}{\beta }_j{K}_w(x_i,x_j)Y_j\right)$$

$$=\sum _{j=1}^{N^{\prime }}{\beta }_j{Y}_j\left(\sum _{i=1}^{N^{\prime }}{\gamma }_i{w}_i{K}_w(x_i,x_j)\right)=0.$$

In conclusion, in this section, we consider some properties of solutions of discrete FSLPN:

$$Y_{\lambda }=\lambda {\mathbf{T}}_w{Y}_{\lambda },{\left(Y_{\lambda }\right)}_i=\lambda \sum _{j=1}^{N^{\prime }}{\beta }_j{K}_w(x_i,x_j){\left(Y_{\lambda }\right)}_j$$

(58)

on the $N^{\prime }$-dimensional space of vectors fulfilling condition (56), equipped with the scalar product defined according to the numerical integration rule (52):

$${\langle F,G\rangle }_w:=\sum _{i=1}^{N^{\prime }}{\gamma }_i{w}_i{F}_i{G}_i.$$

(59)

The following lemma describes the control of orthogonality breaking for the eigenvectors of discrete operator ${\mathbf{T}}_w$.

Lemma 4.

The scalar product (59) of eigenvectors corresponding to distinct non-zero eigenvalues obeys the following inequality:

$$\left|\left(\frac{1}{\rho }-\frac{1}{\lambda }\right){\langle Y_{\lambda },Y_{\rho }\rangle }_w\right|\le 2\left|\right|w\left|\right|\sqrt{\sum _{j=1}^{N^{\prime }}{\beta }_j^2}\sqrt{\sum _{j=1}^{N^{\prime }}{\left(M_{err,j}\right)}^2},$$

(60)

where we denote errors of numerical integration generated by the respective sets of weights $\beta ,\gamma$:

$$\begin{array}{ccc}& Er{r}_{\beta }\left(y_{\lambda }^{ap}(·)K_w(x_j,·)\right)=\left|{\sum }_{i=1}^{N^{\prime }}{\beta }_i{\left(Y_{\lambda }\right)}_i{K}_w(x_j,x_i)-{\int }_a^b{y}_{\lambda }^{ap}\left(x\right)K_w(x_j,x)dx\right|,& \hfill \\ & Er{r}_{\gamma }\left(y_{\lambda }^{ap}(·)K_w(x_j,·)\right)=\left|{\sum }_{i=1}^{N^{\prime }}{\gamma }_i{\left(Y_{\lambda }\right)}_i{K}_w(x_j,x_i)-{\int }_a^b{y}_{\lambda }^{ap}\left(x\right)K_w(x_j,x)dx\right|,& \hfill \end{array}$$

with $y_{\lambda }^{ap}$ being an approximate eigenfunction constructed by using the coordinates of eigenvector $Y_{\lambda }$ and

$$M_{err,j}:=\underset{\beta ,\gamma }{max}\{Er{r}_{\beta }\left(y_{\lambda }^{ap}(·)K_w(x_j,·)\right),Er{r}_{\gamma }\left(y_{\lambda }^{ap}(·)K_w(x_j,·)\right)\}.$$

(61)

Proof. 

In the first part of the proof, we apply condition (56) and the symmetricity of kernel $K_1$:

$$\frac{1}{\rho }{\langle Y_{\lambda },Y_{\rho }\rangle }_w={\langle Y_{\lambda },{\mathbf{T}}_w{Y}_{\rho }\rangle }_w=\sum _{i=1}^{N^{\prime }}{\gamma }_i{w}_i{\left(Y_{\lambda }\right)}_i{\left({\mathbf{T}}_w{Y}_{\rho }\right)}_i=$$

$$=\sum _{i=1}^{N^{\prime }}{\gamma }_i{w}_i{\left(Y_{\lambda }\right)}_i\left(\sum _{j=1}^{N^{\prime }}{\beta }_j{K}_w(x_i,x_j){\left(Y_{\rho }\right)}_j\right)$$

$$=\sum _{i=1}^{N^{\prime }}{\gamma }_i{w}_i{\left(Y_{\lambda }\right)}_i\left[\sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j{\left(Y_{\rho }\right)}_j\left(K_1(x_i,x_j)-\frac{{\sum }_{k=1}^{N^{\prime }}{\gamma }_k{w}_k{K}_1(x_k,x_j)}{{\sum }_{l=1}^{N^{\prime }}{\gamma }_l{w}_l}\right)\right]$$

$$=\sum _{j=1}^{N^{\prime }}\sum _{i=1}^{N^{\prime }}{\gamma }_i{w}_i{\left(Y_{\lambda }\right)}_i{\beta }_j{w}_j{\left(Y_{\rho }\right)}_j{K}_1(x_i,x_j)$$

$$=\sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j{\left(Y_{\rho }\right)}_j\left(\sum _{i=1}^{N^{\prime }}{\gamma }_i{w}_i{\left(Y_{\lambda }\right)}_i{K}_1(x_j,x_i)\right)$$

$$=\sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j{\left(Y_{\rho }\right)}_j\left[\sum _{i=1}^{N^{\prime }}{\gamma }_i{\left(Y_{\lambda }\right)}_i\left(w_i{K}_1(x_j,x_i)-\frac{w_i{\sum }_{k=1}^{N^{\prime }}{\gamma }_k{w}_k{K}_1(x_k,x_i)}{{\sum }_{l=1}^{N^{\prime }}{\gamma }_l{w}_l}\right)\right]$$

$$=\sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j{\left(Y_{\rho }\right)}_j\left(\sum _{i=1}^{N^{\prime }}{\gamma }_i{\left(Y_{\lambda }\right)}_i{K}_w(x_j,x_i)\right)$$

$$=\sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j{\left(Y_{\rho }\right)}_j{\left({\mathbf{T}}_w{Y}_{\lambda }\right)}_j+\sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j{\left(Y_{\rho }\right)}_j\left(\sum _{i=1}^{N^{\prime }}({\gamma }_i-{\beta }_i){\left(Y_{\lambda }\right)}_i{K}_w(x_j,x_i)\right)=$$

$$={\langle Y_{\rho },{\mathbf{T}}_w{Y}_{\lambda }\rangle }_w+\sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j{\left(Y_{\rho }\right)}_j\left(\sum _{i=1}^{N^{\prime }}({\gamma }_i-{\beta }_i){\left(Y_{\lambda }\right)}_i{K}_w(x_j,x_i)\right)=$$

$$=\frac{1}{\lambda }{\langle Y_{\lambda },Y_{\rho }\rangle }_w+\sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j{\left(Y_{\rho }\right)}_j\left(\sum _{i=1}^{N^{\prime }}({\gamma }_i-{\beta }_i){\left(Y_{\lambda }\right)}_i{K}_w(x_j,x_i)\right).$$

Finally, we obtain:

$$\left(\frac{1}{\rho }-\frac{1}{\lambda }\right){\langle Y_{\lambda },Y_{\rho }\rangle }_w=\sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j{\left(Y_{\rho }\right)}_j\left(\sum _{i=1}^{N^{\prime }}({\gamma }_i-{\beta }_i){\left(Y_{\lambda }\right)}_i{K}_w(x_j,x_i)\right).$$

(62)

Now, we shall estimate the absolute value of the expression on the left-hand side of equality (62). In the calculations below, we use the notation: $\left|\right|w\left|\right|$ for the supremum norm of weight function w, $Er{r}_{\beta }(·)$ for error in the numerical integration generated by rule (51) and $Er{r}_{\gamma }(·)$ for error generated by rule (52). In addition, we normalize coordinates of eigenvectors by condition: $|{\left(Y_{\lambda }\right)}_i|\le 1,i=1,$ …$,N^{\prime }$ and apply the Cauchy–Schwarz inequality:

$$\left|\left(\frac{1}{\rho }-\frac{1}{\lambda }\right){\langle Y_{\lambda },Y_{\rho }\rangle }_w\right|\le$$

$$\le \sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j\left|{\left(Y_{\rho }\right)}_j\right|\left(\left|\sum _{i=1}^{N^{\prime }}({\gamma }_i-{\beta }_i){\left(Y_{\lambda }\right)}_i{K}_w(x_j,x_i)\right|\right)$$

$$\le \sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j\left|{\left(Y_{\rho }\right)}_j\right|\times$$

$$\times \left(\left|\sum _{i=1}^{N^{\prime }}{\gamma }_i{\left(Y_{\lambda }\right)}_i{K}_w(x_j,x_i)\mp {\int }_a^b{y}_{\lambda }^{ap}\left(x\right)K_w(x_j,x)dx-\sum _{i=1}^{N^{\prime }}{\beta }_i{\left(Y_{\lambda }\right)}_i{K}_w(x_j,x_i)\right|\right)$$

$$\le \sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j\left|{\left(Y_{\rho }\right)}_j\right|\left(Er{r}_{\gamma }(y_{\lambda }^{ap}(·)K_w(x_j,·)+Er{r}_{\beta }(y_{\lambda }^{ap}(·)K_w(x_j,·)\right)$$

$$\le 2\sum _{j=1}^{N^{\prime }}{\beta }_j{w}_j{M}_{err,j}\le 2\left|\right|w\left|\right|\sqrt{{\sum }_{j=1}^{N^{\prime }}{\beta }_j^2}\sqrt{{\sum }_{j=1}^{N^{\prime }}{\left(M_{err,j}\right)}^2}.$$

□

4. Examples of Numerical Solution

Now, we shall present some numerical results of solution FSLP with homogeneous Neumann boundary conditions. We choose a very simple quadrature rule—the midpoint rectangular rule because the numerical evaluations of kernel values $K_w$ are very time-consuming operations. Thus, we obtain quadrature nodes ${\tilde{x}}_i=a+(i-0.5)\Delta x$, for $i=1,$ …$,N$ and $\Delta x=\frac{b-a}{N}$, and we obtain the following system of N linear algebraic equations corresponding to the reduced discrete integral FSLPN (50):

$$y\left({\tilde{x}}_i\right)=\lambda {\mathbf{T}}_{w}y\left({\tilde{x}}_i\right){\mathbf{T}}_{w}y\left({\tilde{x}}_i\right)=\frac{b-a}{N}\sum _{j=1}^N{K}_w({\tilde{x}}_i,{\tilde{x}}_j)y\left({\tilde{x}}_j\right),i=1,\cdots ,N.$$

(63)

We write the above system of equations in the following matrix form:

$$Y=\lambda {\mathbf{T}}_{w}Y\mathrm{or}{\left({\mathbf{T}}_w\right)}^{-1}Y=\lambda Y,$$

(64)

where

$$Y=\left[\begin{array}{c}{y}_1\hfill \\ y_2\hfill \\ ⋮\hfill \\ y_N\hfill \end{array}\right],{\mathbf{T}}_w=\frac{b-a}{N}\left[\begin{array}{cccc}{\left(K_w\right)}_{1,1}\hfill & {\left(K_w\right)}_{1,2}& \cdots & {\left(K_w\right)}_{1,N}\\ {\left(K_w\right)}_{2,1}\hfill & {\left(K_w\right)}_{2,2}& & {\left(K_w\right)}_{2,N}\\ ⋮\hfill & & \ddots & ⋮\\ {\left(K_w\right)}_{N,1}\hfill & {\left(K_w\right)}_{N,2}& \cdots & {\left(K_w\right)}_{N,N}\end{array}\right]$$

(65)

while $y_i=y\left({\tilde{x}}_i\right)$ and ${\left(K_w\right)}_{i,j}=K_w\left({\tilde{x}}_i,{\tilde{x}}_j\right),$ where

$$K_w\left({\tilde{x}}_i,{\tilde{x}}_j\right)=w\left({\tilde{x}}_j\right)\left(K_1\left({\tilde{x}}_i,{\tilde{x}}_j\right)-\frac{{\int }_a^{b}w\left(\xi \right)K_1\left(\xi ,{\tilde{x}}_j\right)d\xi }{{\int }_a^{b}w\left(\xi \right)d\xi }\right)$$

(66)

with the symmetric kernel part:

$$K_1\left({\tilde{x}}_i,{\tilde{x}}_j\right)=\frac{1}{{\mathsf{\Gamma }}^2\left(\alpha \right)}\underset{a}{\overset{min\left({\tilde{x}}_i,{\tilde{x}}_j\right)}{\int }}\frac{{\left({\tilde{x}}_i-t\right)}^{\alpha -1}{\left({\tilde{x}}_j-t\right)}^{\alpha -1}}{p\left(t\right)}dt.$$

(67)

For function $w\left(x\right)=1$, Equation (66) can be simplified to the form:

$$K_w\left({\tilde{x}}_i,{\tilde{x}}_j\right)=K_1\left({\tilde{x}}_i,{\tilde{x}}_j\right)-\frac{1}{b-a}\underset{a}{\overset{b}{\int }}{K}_1\left(\xi ,{\tilde{x}}_j\right)d\xi .$$

(68)

The discrete/matrix eigenvalue problem (64) can be solved by using mathematical software. The resulting solution is the set of eigenvalues ${\lambda }_{\left(k\right)}$, for $k=1,$ …$,N$ and the set of eigenvectors $Y_{{\lambda }_{\left(k\right)}}$, for $k=1,$ …$,N$, that correspond to the eigenvalues ${\lambda }_{\left(k\right)}$, satisfying:

$$\left({\left({\mathbf{T}}_w\right)}^{-1}-{\lambda }_{\left(k\right)}\mathbf{I}\right)Y_{{\lambda }_{\left(k\right)}}=0.$$

(69)

In order to obtain the numerical solution at all points in the considered interval $\left[a,b\right]$, we construct the approximate eigenfunctions, for example, using the step function $\chi$:

$$y_{{\lambda }_{\left(k\right)}}^{ap}\left(x\right)=\sum _{j=1}^N{\left({\mathbf{Y}}_{{\lambda }_{\left(k\right)}}\right)}_j{\chi }_{\left[{\tilde{x}}_j-e,{\tilde{x}}_j+e\right)}\left(x\right),e:=\frac{b-a}{2N}.$$

(70)

Now, we report on two examples of numerical calculations of eigenvalues and eigenfunctions to verify the proposed numerical method. In both examples, we consider the interval of calculations $[0,1]$ and the order of derivatives $\alpha \in \{1,0.8,0.6\}$. Representative results for the test problem are collected and presented in the form of graphs and tables. The values of $er{r}_{\lambda }$, presented in tables, have been calculated for the fixed parameters $\alpha$ and variable values of N utilizing the following formula [48]:

$$er{c}_{\lambda }\left(N,\alpha ,k\right)={log}_2\frac{{\lambda }_{\left(k\right)}^{\left(N,\alpha \right)}-{\lambda }_{\left(k\right)}^{\left(N/2,\alpha \right)}}{{\lambda }_{\left(k\right)}^{\left(2N,\alpha \right)}-{\lambda }_{\left(k\right)}^{\left(N,\alpha \right)}}.$$

(71)

Moreover, the approximate eigenfunctions were normalized by:

$${\int }_a^{b}w\left(x\right)y_{{\lambda }_{\left(k\right)}}^{ap}\left(x\right)y_{{\lambda }_{\left(k\right)}}^{ap}\left(x\right)dx=1,k=1,\cdots ,N.$$

(72)

4.1. Example I

The first example is devoted to the fractional equivalent of the classical harmonic oscillator problem with $p\left(x\right)=1$, $q\left(x\right)=0$ and $w\left(x\right)=1$. Figure 1 shows graphs of the approximate eigenfunctions corresponding to the first four eigenvalues for the considered problem. The calculations that are presented on the plot have been performed for $N=4000$. In

Table 1. Numerical values of the first 8 eigenvalues and the experimental rates of convergence $er{c}_{\lambda }$ for $\alpha \in \{1,0.8,0.6\}$, $p\left(x\right)=1$, $q\left(x\right)=0$ and $w\left(x\right)=1$.

		$\mathit{\alpha }=1$		$\mathit{\alpha }=0.8$		$\mathit{\alpha }=0.6$
$\mathit{k}$	$\mathit{N}$	${\mathit{\lambda }}_{\left(\mathit{k}\right)}$	${\mathit{erc}}_{\mathit{\lambda }}$	${\mathit{\lambda }}_{\left(\mathit{k}\right)}$	${\mathit{erc}}_{\mathit{\lambda }}$	${\mathit{\lambda }}_{\left(\mathit{k}\right)}$	${\mathit{erc}}_{\mathit{\lambda }}$
1	250	9.8694745	-	7.8556424	-	5.4750947	-
	500	9.8695719	1.997	7.8571537	1.610	5.5029938	1.191
	1000	9.8695963	2.000	7.8576488	1.607	5.5152125	1.197
	2000	9.8696024	2.024	7.8578113	1.608	5.5205439	1.199
	4000	9.8696039	-	7.8578646	-	5.5228666	-
2	250	39.476340	-	20.327286	-	10.008656	-
	500	39.477898	1.998	20.337158	1.602	10.102136	1.182
	1000	39.478288	2.007	20.340410	1.602	10.143348	1.192
	2000	39.478385	2.015	20.341481	1.601	10.161385	1.197
	4000	39.478409	-	20.341834	-	10.169253	-
3	250	88.815920	-	38.881625	-	16.167499	-
	500	88.823810	2.000	38.917945	1.603	16.412925	1.170
	1000	88.825782	2.000	38.929901	1.603	16.522008	1.187
	2000	88.826275	1.991	38.933836	1.603	16.569919	1.195
	4000	88.826399	-	38.935131	-	16.590853	-
4	250	157.88042	-	59.665008	-	21.512600	-
	500	157.90536	2.001	59.749754	1.599	21.948949	1.159
	1000	157.91159	1.998	59.777723	1.601	22.144364	1.182
	2000	157.91315	2.000	59.786945	1.601	22.230487	1.192
	4000	157.91354	-	59.789984	-	22.268175	-
5	250	246.65895	-	85.292837	-	27.925213	-
	500	246.71982	2.000	85.466870	1.600	28.665274	1.147
	1000	246.73504	2.002	85.524297	1.602	28.999570	1.177
	2000	246.73884	2.000	85.543221	1.602	29.147466	1.190
	4000	246.73979	-	85.549456	-	29.212295	-
6	250	355.13747	-	112.54371	-	33.623633	-
	500	355.26368	2.000	112.84505	1.596	34.701683	1.135
	1000	355.29524	2.000	112.94470	1.600	35.192608	1.171
	2000	355.30313	2.002	112.97757	1.599	35.410598	1.187
	4000	355.30510	-	112.98842	-	35.506312	-
7	250	483.29886	-	144.00145	-	40.125471	-
	500	483.53266	2.000	144.49694	1.596	41.671012	1.122
	1000	483.59113	2.001	144.66082	1.600	42.381118	1.165
	2000	483.60574	1.997	144.71487	1.601	42.697703	1.185
	4000	483.60940	-	144.73269	-	42.836953	-
8	250	631.12288	-	176.75038	-	45.983903	-
	500	631.52170	2.000	177.49404	1.594	48.024089	1.110
	1000	631.62143	2.000	177.74045	1.598	48.969362	1.160
	2000	631.64637	2.001	177.82183	1.600	49.392428	1.182
	4000	631.65260	-	177.84867	-	49.578838	-

, we present the numerical values of the first eight eigenvalues and the experimental rate of convergence $er{c}_{\lambda }$ of numerical calculations of the k-th eigenvalue. The table contains numerical results obtained for orders $\alpha \in \{1,0.8,0.6\}$ and different values of $N\in \{250,500,1000,2000,4000\}$.

4.2. Example II

The second example contains numerical results obtained for FSLPN with functions: $p\left(x\right)=2x^2+1$, $q\left(x\right)=0$ and $w\left(x\right)=cos\left(4\pi x\right)+2$. Figure 2 shows graphs of the approximate eigenfunctions corresponding to the first four eigenvalues for the considered problem. In this case, the calculations presented on the plot have also been performed for $N=4000$. In

Table 2. Numerical values of the first 8 eigenvalues and the experimental rates of convergence $er{c}_{\lambda }$ for $\alpha \in \{1,0.8,0.6\}$, $p\left(x\right)=2x^2+1$, $q\left(x\right)=0$ and $w\left(x\right)=cos\left(4\pi x\right)+2$.

		$\mathit{\alpha }=1$		$\mathit{\alpha }=0.8$		$\mathit{\alpha }=0.6$
$\mathit{k}$	$\mathit{N}$	${\mathit{\lambda }}_{\left(\mathit{k}\right)}$	${\mathit{erc}}_{\mathit{\lambda }}$	${\mathit{\lambda }}_{\left(\mathit{k}\right)}$	${\mathit{erc}}_{\mathit{\lambda }}$	${\mathit{\lambda }}_{\left(\mathit{k}\right)}$	${\mathit{erc}}_{\mathit{\lambda }}$
1	250	7.0687342	-	6.2956345	-	4.6051127	-
	500	7.0687782	2.000	6.2967500	1.622	4.6361847	1.195
	1000	7.0687892	1.974	6.2971124	1.620	4.6497610	1.199
	2000	7.0687920	2.000	6.2972303	1.615	4.6556764	1.200
	4000	7.0687927	-	6.2972688	-	4.6582511	-
2	250	24.644959	-	12.771392	-	6.4333679	-
	500	24.646018	1.999	12.777978	1.602	6.4923883	1.181
	1000	24.646283	2.005	12.780147	1.603	6.5184265	1.192
	2000	24.646349	1.957	12.780861	1.597	6.5298259	1.197
	4000	24.646366	-	12.781097	-	6.5347998	-
3	250	70.274276	-	30.075510	-	12.136357	-
	500	70.279914	1.999	30.101699	1.604	12.331176	1.174
	1000	70.281324	2.002	30.110313	1.604	12.417505	1.189
	2000	70.281676	2.000	30.113146	1.604	12.455371	1.195
	4000	70.281764	-	30.114078	-	12.471906	-
4	250	127.27968	-	46.851483	-	16.228206	-
	500	127.29583	1.999	46.912533	1.602	16.596768	1.171
	1000	127.29987	2.000	46.932645	1.603	16.760495	1.187
	2000	127.30088	2.014	46.939266	1.603	16.832387	1.195
	4000	127.30113	-	46.941446	-	16.863793	-
5	250	197.81960	-	67.232260	-	20.906261	-
	500	197.88521	1.999	67.418158	1.600	21.637418	1.143
	1000	197.90162	2.001	67.479487	1.602	21.968596	1.175
	2000	197.90572	2.007	67.499697	1.602	22.115285	1.189
	4000	197.90674	-	67.506354	-	22.179616	-
6	250	285.02794	-	89.803850	-	25.836658	-
	500	285.14300	2.000	90.081339	1.598	26.887072	1.144
	1000	285.17177	2.001	90.173017	1.600	27.362312	1.175
	2000	285.17896	1.998	90.203255	1.601	27.572716	1.189
	4000	285.18076	-	90.213224	-	27.664975	-
7	250	381.52715	-	112.94393	-	30.374162	-
	500	381.75312	2.000	113.43053	1.598	31.911832	1.117
	1000	381.80962	2.001	113.59131	1.601	32.620555	1.164
	2000	381.82374	1.996	113.64431	1.601	32.936757	1.185
	4000	381.82728	-	113.66178	-	33.075861	-
8	250	495.02316	-	136.38509	-	34.838979	-
	500	495.40389	2.000	137.06049	1.593	36.579470	1.111
	1000	495.49909	2.000	137.28436	1.598	37.385104	1.159
	2000	495.52289	2.000	137.35829	1.600	37.745787	1.182
	4000	495.52884	-	137.38267	-	37.904751	-

, we present the numerical values of the first eight eigenvalues and the experimental rate of convergence $er{c}_{\lambda }$ of numerical calculations of the k-th eigenvalue. The table contains the numerical results obtained for orders $\alpha \in \{1,0.8,0.6\}$ and different values of $N\in \{250,500,1000,2000,4000\}$.

5. Conclusions

In the paper, the FSLP with the fractional Neumann boundary conditions was studied. The considered homogeneous Neumann boundary conditions require an additional integral constraint on the solutions’ space. The additional integral condition causes the considered eigenvalue problem to become more subtle and complex than the FSLP with homogeneous Dirichlet, mixed or Robin conditions. This fact was a premise for examining the particular problem with assumption $q=0$. Such a problem (the differential FSLPN) was transformed to the equivalent integral one on a suitable function space. This transformation was based on results presented in the paper [51]. In the main part of the paper, the construction of the discrete version of the integral FSLPN was developed and studied. Furthermore, the orthogonality property of teh eigenvectors of the discrete integral operator was analyzed. These considerations include two particular cases of the construction of the discrete integral FSLPN. The first case is devoted to the discrete integral FSLPN received by utilizing a single numerical integration scheme. Based on such an assumption, the orthogonality (with respect to the adequate scalar product) of eigenvectors corresponding to distinct eigenvalues of the considered discrete FSLPN was proved. The second case covers the situation when the discrete integral FSLPN was derived by applying a hybrid numerical scheme (two different numerical integration schemes: one in the construction of the discrete integral FSLO and another for the integral condition restricting the solutions’ space). In this case, the orthogonality of the eigenvectors of the constructed discrete FSLPN has not been proved. However, the control of orthogonality breaking for eigenvectors of the discrete integral operator was established in Lemma 4. In the final part of the paper, two examples of numerical calculations of the eigenvalues and eigenvectors were reported. The performed calculations show that the experimental rate of convergence depends on the fractional order $\alpha$ and is close to $2\alpha$.