Super-Exponentially Convergent Parallel Algorithm for a Fractional Eigenvalue Problem of Jacobi-Type

TL;DR

This paper introduces a super-exponentially convergent parallel algorithm for solving fractional Jacobi-type eigenvalue problems, transforming a nonlinear problem into a sequence of linear problems with high convergence speed.

Contribution

It proposes a novel recursive algorithm using truncated functional discrete schemes that achieves super-exponential convergence for fractional eigenvalue problems.

Findings

Algorithm converges super-exponentially as m increases.

Eigenpairs can be computed in parallel for all indices.

Numerical examples validate the theoretical convergence rates.

Abstract

A new algorithm for eigenvalue problems for the fractional Jacobi type ODE is proposed. The algorithm is based on piecewise approximation of the coefficients of the differential equation with subsequent recursive procedure adapted from some homotopy considerations. As a result, the eigenvalue problem (which is in fact nonlinear) is replaced by a sequence of linear boundary value problems (besides the first one) with a singular linear operator called the exact functional discrete scheme (EFDS). A finite subsequence of $m$ terms, called truncated functional discrete scheme (TFDS), is the basis for our algorithm. The approach provides an super-exponential convergence rate as $m \to \infty$. The eigenpairs can be computed in parallel for all given indexes. The algorithm is based on some recurrence procedures including the basic arithmetical operations with the coefficients of some expansions only. This is an exact symbolic algorithm (ESA) for $m=\infty$ and a truncated symbolic algorithm (TSA) for a finite $m$. Numerical examples are presented to support the theory.

Figures

Tables

Table 1. Table 1 : Approximations λ 20 n subscript 20 𝜆 𝑛 \overset{20}{\lambda}_{n} to the eigenvalues λ n subscript 𝜆 𝑛 \lambda_{n} and the values of the ratio r n subscript 𝑟 𝑛 r_{n} of the geometric progression for n = 0 , 1 , 2 , 3 , 4 , 10 𝑛 0 1 2 3 4 10 n=0,1,2,3,4,10 .

$n$	$r_n$	${\stackrel{20}{\lambda }}_n$
0	1.505	0.630053891717269391713596782178
1	0.602	2.30514376782605437183729542346
2	0.446	4.56663095405447807274309207123
3	0.370	7.26852041545097231248033107978
4	0.324	10.3587519350445765031940934458
10	0.208	35.2508805155975526382623041970

Table 2. Table 2 : Corrections λ n ( m ) , n = 0 , 10 formulae-sequence subscript superscript 𝜆 𝑚 𝑛 𝑛 0 10 \lambda^{(m)}_{n},n=0,10 , and the norms of corrections to the corresponding eigenfunctions u n ( m ) , n = 0 , 10 formulae-sequence subscript superscript 𝑢 𝑚 𝑛 𝑛 0 10 u^{(m)}_{n},n=0,10 , for the FD-method of the ranks m = 0 , 1 , … ​ 10 , 20 𝑚 0 1 … 10 20 m=0,1,...10,20 .

$m$	${\lambda }_0^{\left(m\right)}$	$\Vert u_0^{\left(m\right)}\Vert$	${\lambda }_{10}^{\left(m\right)}$	$\Vert u_{10}^{\left(m\right)}\Vert$
0	0.6647	1.373	35.25	0.3627
1	$-0.03095$	0.06253	$-0.0002053$	0.009806
2	$-0.003648$	0.001742	0.00009424	0.00007673
3	$-0.00001679$	0.00004165	$-1.207\cdot {10}^{-8}$	0.000002101
4	0.000001138	0.000002067	$6.993\cdot {10}^{-10}$	$2.059\cdot {10}^{-8}$
5	$9.723\cdot {10}^{-8}$	$4.556\cdot {10}^{-8}$	$-7.034\cdot {10}^{-13}$	$5.656\cdot {10}^{-10}$
6	$-1.319\cdot {10}^{-10}$	$5.386\cdot {10}^{-9}$	$1.597\cdot {10}^{-14}$	$5.657\cdot {10}^{-12}$
7	$-2.950\cdot {10}^{-10}$	$1.856\cdot {10}^{-10}$	$-4.932\cdot {10}^{-17}$	$1.553\cdot {10}^{-13}$
8	$-4.443\cdot {10}^{-12}$	$1.438\cdot {10}^{-11}$	$7.277\cdot {10}^{-19}$	$1.558\cdot {10}^{-15}$
9	$7.229\cdot {10}^{-13}$	$7.175\cdot {10}^{-13}$	$-6.316\cdot {10}^{-21}$	$4.277\cdot {10}^{-17}$
10	$2.848\cdot {10}^{-14}$	$3.451\cdot {10}^{-14}$	$2.910\cdot {10}^{-23}$	$4.293\cdot {10}^{-19}$
20	$-3.979\cdot {10}^{-27}$	$2.621\cdot {10}^{-26}$	$-3.694\cdot {10}^{-45}$	$6.793\cdot {10}^{-35}$

Table 3. Table 3 : Approximations λ 𝑚 n subscript 𝑚 𝜆 𝑛 \overset{m}{\lambda}_{n} of the ranks m = 20 , 30 𝑚 20 30 m=20,30 to the eigenvalues λ n subscript 𝜆 𝑛 \lambda_{n} and the values of the ratio r n subscript 𝑟 𝑛 r_{n} of the geometric progression for n = 0 , 1 , 2 , 3 , 4 , 10 𝑛 0 1 2 3 4 10 n=0,1,2,3,4,10 .

$n$	$r_n$	${\stackrel{20}{\lambda }}_n$	${\stackrel{30}{\lambda }}_n$
0	1.202	0.15067840864298071809615545712	0.15067840864298071808633595280
1	1.093	1.43411319565086825800583685557	1.43411319565086825801565633056
2	0.687	3.36697390018171063527745195043	3.36697390018171063527745197976
3	0.543	5.81865643939502492320620997309	5.81865643939502492320620997309
4	0.463	8.69578418818622939490618135946	8.69578418818622939490618135946
10	0.286	32.6418532333094312351786092007	32.6418532333094312351786092007

Table 4. Table 4 : Corrections λ n ( m ) , n = 0 , 1 , 2 , 3 , 4 formulae-sequence subscript superscript 𝜆 𝑚 𝑛 𝑛 0 1 2 3 4 \lambda^{(m)}_{n},n=0,1,2,3,4 for the FD-method of the ranks m = 0 , 1 , … , 10 , 20 , 30 𝑚 0 1 … 10 20 30 m=0,1,...,10,20,30 .

$m$	${\lambda }_0^{\left(m\right)}$	${\lambda }_1^{\left(m\right)}$	${\lambda }_2^{\left(m\right)}$	${\lambda }_3^{\left(m\right)}$	${\lambda }_4^{\left(m\right)}$
0	$0.07396$	$1.294$	$3.236$	$5.691$	$8.569$
1	$0.08152$	$0.1389$	$0.1297$	$0.1269$	$0.1260$
2	$-0.005109$	$0.001165$	$0.001380$	$0.0004626$	$0.0002623$
3	$0.0003151$	$-0.0003367$	$0.7655\cdot {10}^{-5}$	$0.9807\cdot {10}^{-5}$	$0.2132\cdot {10}^{-5}$
4	$-0.5669\cdot {10}^{-5}$	$0.8132\cdot {10}^{-5}$	$-0.2832\cdot {10}^{-5}$	$2.742\cdot {10}^{-7}$	$2.891\cdot {10}^{-9}$
5	$-0.1551\cdot {10}^{-5}$	$0.1676\cdot {10}^{-5}$	$-1.164\cdot {10}^{-7}$	$-9.865\cdot {10}^{-9}$	$1.484\cdot {10}^{-9}$
6	$1.943\cdot {10}^{-7}$	$-1.979\cdot {10}^{-7}$	$4.191\cdot {10}^{-9}$	$-5.421\cdot {10}^{-10}$	$-1.394\cdot {10}^{-11}$
7	$-4.321\cdot {10}^{-9}$	$3.925\cdot {10}^{-9}$	$4.053\cdot {10}^{-10}$	$-7.920\cdot {10}^{-12}$	$-1.403\cdot {10}^{-12}$
8	$-1.618\cdot {10}^{-9}$	$1.619\cdot {10}^{-9}$	$-5.889\cdot {10}^{-13}$	$1.278\cdot {10}^{-13}$	$-2.140\cdot {10}^{-14}$
9	$2.199\cdot {10}^{-10}$	$-2.186\cdot {10}^{-10}$	$-1.253\cdot {10}^{-12}$	$1.466\cdot {10}^{-14}$	$-7.319\cdot {10}^{-17}$
10	$-4.351\cdot {10}^{-12}$	$4.382\cdot {10}^{-12}$	$-3.167\cdot {10}^{-14}$	$5.091\cdot {10}^{-16}$	$1.464\cdot {10}^{-18}$
20	$-1.116\cdot {10}^{-20}$	$1.116\cdot {10}^{-20}$	$2.576\cdot {10}^{-26}$	$-2.037\cdot {10}^{-30}$	$7.184\cdot {10}^{-34}$
30	$8.072\cdot {10}^{-30}$	$-8.072\cdot {10}^{-30}$	$7.176\cdot {10}^{-38}$	$1.316\cdot {10}^{-44}$	$4.848\cdot {10}^{-50}$

Equations

$${}_{-1}{D}_x^{s}v(x)=\frac{1}{\Gamma (1-s)}\frac{d}{dx}{\int }_{-1}^x\frac{v(y)}{(x-y{)}^s}dy,x\in \left(-1,1\right),$$

$${}_x{D}_1^{s}v(x)=-\frac{1}{\Gamma (1-s)}\frac{d}{dx}{\int }_x^1\frac{v(y)}{(y-x{)}^s}dy,x\in \left(-1,1\right),$$

$$\begin{array}{rl}& {}^{+}{L}_{\alpha ,\beta }^{2s}u(x)+\left(q(x)-\lambda \right)u(x)=0,x\in \left(-1,1\right),s\in \left(0,1\right),\alpha \ge 0,\beta >-1,\\ & \left|u(-1)\right|<\infty ,\left|u(1)\right|<c{\delta }_{\alpha ,0},c=\text{const}<\infty \end{array}$$

$${}^{+}{L}_{\alpha ,\beta }^{2s}u(x)={\omega }^{(\alpha ,-\beta )}{}_{-1}{D}_x^s\left\{{\omega }^{(-\alpha +s,\beta +s)}{}_x{D}_1^{s}u(x)\right\},x\in \left(-1,1\right)$$

$${\omega }^{(\alpha ,\beta )}={\omega }^{(\alpha ,\beta )}(x)=(1-x{)}^{\alpha }(1+x{)}^{\beta }.$$

$${\left|{\omega }^{(-\alpha +s,\beta +s)}{}_x{D}_1^s\left\{(1-x{)}^{\alpha }{P}_n^{(\alpha ,\beta )}(x)\right\}\right|}_{x=1}<\infty .$$

$$u_n(x)=\underset{j=0}{\sum ^{\infty }}{u}_n^{(j)}(x),{\lambda }_n=\underset{j=0}{\sum ^{\infty }}{\lambda }_n^{(j)}$$

$${\stackrel{m}{u}}_n(x)=\underset{j=0}{\sum ^m}{u}_n^{(j)}(x),{\stackrel{m}{\lambda }}_n=\underset{j=0}{\sum ^m}{\lambda }_n^{(j)}$$

$$\begin{array}{rl}& {}^{+}{L}_{\alpha ,\beta }^{2s}{u}_n^{(j+1)}(x)-{\lambda }_n^{(0)}{u}_n^{(j+1)}(x)=F_n^{(j+1)}(x),x\in \left(-1,1\right),\\ & s\in \left(0,1\right),\alpha \ge 0,\beta >-1,\\ & \left|u_n^{(j+1)}(-1)\right|<\infty ,\left|u_n^{(j+1)}(1)\right|\le c{\delta }_{\alpha ,0},c=const<\infty ,\end{array}$$

$$F_n^{(j+1)}(x)=\sum _{p=0}^j{\lambda }_n^{(j+1-p)}{u}_n^{(p)}(x)-q(x)u_n^{(j)}(x),j=0,1,2,...,n=0,1,2,...$$

$${\left|{\omega }^{(-\alpha +s,\beta +s)}{}_x{D}_1^s\left\{u_n^{(j+1)}(x)\right\}\right|}_{x=1}<\infty .$$

$$\begin{array}{rl}& {}^{+}{L}_{\alpha ,\beta }^{2s}{u}_n^{(0)}(x)-{\lambda }_n^{(0)}{u}_n^{(0)}(x)=0,x\in \left(-1,1\right),s\in \left(0,1\right),\alpha \ge 0,\beta >-1\\ & \left|u_n^{(0)}(-1)\right|<\infty ,\left|u_n^{(0)}(1)\right|\le c{\delta }_{\alpha ,0},c=const<\infty \end{array}$$

$${\left|{\omega }^{(-\alpha +s,\beta +s)}{}_x{D}_1^s\left\{u_n^{(0)}(x)\right\}\right|}_{x=1}<\infty .$$

$${}^{+}{J}_n^{(-\alpha ,\beta )}(x)=(1-x{)}^{\alpha }{P}_n^{(\alpha ,\beta )}(x),$$

$${\int }_{-1}^1{F}_n^{(j+1)}(x)u_n^{(0)}(x){\omega }^{(-\alpha ,\beta )}(x)dx=0.$$

$${\int }_{-1}^1{u}_n^{(0)}(x)u_n^{(j+1)}(x){\omega }^{(-\alpha ,\beta )}(x)dx=0.$$

$${\int }_{-1}^1{P}_n^{(\alpha ,\beta )}(x)P_k^{(\alpha ,\beta )}(x){\omega }^{(\alpha ,\beta )}(x)dx={\gamma }_n^{(\alpha ,\beta )}{\delta }_{n,k},$$

$${\gamma }_n^{(\alpha ,\beta )}=\frac{2^{\alpha +\beta +1}\Gamma (n+\alpha +1)\Gamma (n+\beta +1)}{(2n+\alpha +\beta +1)n!\Gamma (n+\alpha +\beta +1)}.$$

$$\begin{array}{rl}& u_n^{(0)}(x)={\left({\gamma }_n^{(\alpha ,\beta )}\right)}^{-1/2}\cdot {}^{+}{J}_n^{(-\alpha ,\beta )}(x),\\ & {\lambda }_n^{(0)}=\frac{\Gamma (n+\alpha +1)\Gamma (n+\beta +s+1)}{\Gamma (n+\alpha -s+1)\Gamma (n+\beta +1)},n=0,1,2,...\end{array}$$

$${\lambda }_n^{(j+1)}={\int }_{-1}^{1}q(x)u_n^{(0)}(x)u_n^{(j)}(x){\omega }^{(-\alpha ,\beta )}(x)dx.$$

$$u_n^{(j+1)}(x)=\underset{p=0,p\ne n}{\sum ^{\infty }}{\int }_{-1}^1{F}_n^{\left(j+1\right)}\left(\xi \right)u_p^{(0)}\left(\xi \right){\omega }^{(-\alpha ,\beta )}\left(\xi \right)d\xi \frac{u_p^{(0)}(x)}{{\lambda }_n^{(0)}-{\lambda }_p^{(0)}}.$$

$$\begin{array}{rl}& \Vert u_n^{(j+1)}\Vert ={\left(\underset{p=0,p\ne n}{\sum ^{\infty }}{\left({\int }_{-1}^1{F}_n^{\left(j+1\right)}\left(\xi \right)u_p^{(0)}\left(\xi \right){\omega }^{(-\alpha ,\beta )}\left(\xi \right)d\xi {\left({\lambda }_n^{(0)}-{\lambda }_p^{(0)}\right)}^{-1}\right)}^2\right)}^{1/2}\\ & \le \max \left\{{\left({\lambda }_n^{(0)}-{\lambda }_{n-1}^{(0)}\right)}^{-1},{\left({\lambda }_{n+1}^{(0)}-{\lambda }_n^{(0)}\right)}^{-1}\right\}\Vert F_n^{\left(j+1\right)}\Vert =M_n\Vert F_n^{\left(j+1\right)}\Vert ,\end{array}$$

$$\Vert u_n^{(j+1)}\Vert ={\Vert u_n^{(j+1)}\Vert }_{L_{\omega }^2}={\left({\int }_{-1}^1{\left(u_n^{(j+1)}\left(x\right)\right)}^2{\omega }^{(-\alpha ,\beta )}\left(x\right)dx\right)}^{1/2},\Vert u_n^{(0)}\Vert =1,$$

$$\begin{array}{rl}& M_n=\max \left\{\frac{1}{{\lambda }_n^{(0)}-{\lambda }_{n-1}^{(0)}},\frac{1}{{\lambda }_{n+1}^{(0)}-{\lambda }_n^{(0)}}\right\}\\ & =\frac{\Gamma \left(n+\beta +1\right)\Gamma \left(n+\alpha +1-s\right)}{s\left(2n+\beta +\alpha \right)\Gamma \left(n+\alpha \right)\Gamma \left(n+\beta +s\right)}\\ & \times \max \left\{1,\frac{\left(n+\beta +1\right)\left(2n+\beta +\alpha \right)\left(n+\alpha +1-s\right)}{\left(2n+\beta +\alpha +2\right)\left(n+\alpha \right)\left(n+\beta +s\right)}\right\}.\end{array}$$

$$\frac{\Gamma \left(z+\gamma \right)}{\Gamma \left(z+\eta \right)}=z^{\gamma -\eta }\left[1+\frac{\left(\gamma -\eta \right)\left(\gamma +\eta -1\right)}{2z}+o\left({\left|z\right|}^{-2}\right)\right]\text{as}z\to \infty$$

$$\begin{split}&\underset{n\to\infty}{\lim}M_{n}=\frac{1}{2s}\cdot\left\{\begin{array}[]{ll}\infty,&\;s\in(0,1/2),\\ 1,&\;s=1/2,\\ 0,&\;s\in(1/2,1).\end{array}\right.\end{split}$$

$$\Vert u_n^{(j+1)}\Vert \le M_n{\Vert q\Vert }_{\infty }\left\{\sum _{l=1}^j\Vert u_n^{(j-l)}\Vert \Vert u_n^{(l)}\Vert +\Vert u_n^{(j)}\Vert \right\},\left|{\lambda }_n^{(j+1)}\right|\le {\Vert q\Vert }_{\infty }\Vert u_n^{(j)}\Vert ,$$

$$\Vert u_n^{(j+1)}\Vert \le M_n{\Vert q\Vert }_{\infty }\sum _{l=0}^j\Vert u_n^{(j-l)}\Vert \Vert u_n^{(l)}\Vert .$$

$${\overline{U}}_{j+1}=\sum _{l=0}^j{\overline{U}}_{j-l}{\overline{U}}_l.$$

$${\overline{U}}_{j+1}=\frac{\left(2j+2\right)!}{\left(j+1\right)!\left(j+2\right)!}=4^{j+1}2\frac{\left(2j+1\right)!!}{\left(2j+4\right)!!}.$$

Full text

{topmatter}

1 Introduction

Using various definitions of fractional derivatives one can consider, e.g., differential operators of fractional order, boundary and eigenvalue problems for these operators, as well as various approximation methods for them (see e.g. [14, 8]).

The eigenvalue problem (EVP) is the problem of finding eigenvalues (frequencies) and eigenfunctions (vibration shapes) or so called eigenpairs and plays an important role in various applications concerned with vibrations and wave processes [3, 19]. Often it is needed to compute a great number (hundreds of thousands) of eigenvalues and eigenfunctions including eigenpairs with great indexes (see e.g. [19, p. 273]). A very efficient approach to solve such problems represent methods based on perturbation and homotopy ideas [2, 4], such as the FD-method [15, 16, 17, 7] (compare with the Adomian decomposition method [1, 20]). These approaches allow in a natural way the use of computer algebra tools (see e.g. [13]).

In [6] a spectral approximation of Jacobi-type fractional differential equations (FDEs) is considered. Chen, Shen and Wang defined a new class of generalized Jacobi functions (GJFs), which are the eigenfunctions of some fractional differential operator and can serve as natural basis functions for properly designed spectral methods for FDEs. The spectral approximation results for these GJFs in weighted Sobolev spaces involving fractional derivatives are established and 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. One of the important drawbacks of these methods is the following: their accuracy decreaces with the growth of the eigenvalue index.

In order to find numerically the higher eigenvalues we propose a new approach described below which we will refer to as the FD-method (following [15, 16, 17, 7, 11]). This approach has been applied also to EVPs with multiple eigenvalues [11]. We would like to emphasize the following features of our approach, which differ from many other methods: 1) all eigenpairs can be computed in parallel, 2)The convergence rate increases as the index of the eigenpair increases, 3) the approximation called the truncated discrete scheme (TFDS) is derived from an exact discrete representation of the solution which we call the exact discrete scheme (EFDS), 4) the computational algorithm operates with the coefficients of some expansions of the initial data and uses the basic arithmetical operations only so that the computer algebra tools are the natural medium of choice, 5) under some assumptions it is possible to pass to the limit at the truncated sum of the series and to obtain the exact solution in a closed form [18]. This is the reason, why we call our algorithm symbolic.

The article is organized as follows. First, we describe the FD-method for the fractional Jacobi-type differential operator. In Section 3 we prove the super-exponential convergence rate of our method for the case when the potential function is approximated by the constant zero. Section 4 deals with the reformulation of the algorithm in the case of a polynomial potential as a recursive procedure for the coefficients of the representation of corrections of the FD-method through the generalized Jacobi functions. Although these corrections are the solutions of some BVPs with a singular differential operator, this procedure does not use any BVP-solver but some recursions for the expansion coefficients with basic arithmetical operations only. We call such procedures symbolic algorithms. Note that similar algorithms for some special potential functions were proposed in [7]. In Section 5 numerical examples are given to support the theoretical results.

2 Algorithm of the FD-method for a fractional Sturm-Liouville-type operator with an Jacobi-type weight function

Let us define the left-side and right-side Riemann-Liouville fractional derivatives of order $s\in\left(0,1\right)$

[TABLE]

[TABLE]

where $\Gamma\left(\alpha\right)$ is the gamma-function. Let us consider the following eigenvalue problem of Sturm-Liouville-type

[TABLE]

for the fractional Jacobi-type differential operator

[TABLE]

with a piecewise smooth potential coefficient $q(x)$ and the Jacobi-type weight function

[TABLE]

In (3) $\delta_{n,k}$ denotes the Kronecker delta. A solution of the eigenvalue problem (3) consists of an eigenvalue $\lambda$ and of a corresponding eigenfunction $u(x)$, in other words, of the eigenpair $\lambda,u(x)$. If the set of eigenpairs is countable, we identify each pair by a natural number $n=1,2,...$. According to [6] (see Corollary 3.1) the Sturm-Liouville singular differential operator (4) is selfadjoint provided that $s\in\left(0,1\right)$. For $s-1<\alpha<0$ the second boundary condition in (3) is given as follows:

[TABLE]

We apply to problem (3) the FD-method (see e.g. [7] for details) which is a combination of perturbation of the original differential operator by a parameter dependent operator (embedding) and then a “trip” along this parameter from a “simple” problem to the original one (homotopy). The homotopy idea was exploited in various ways e.g. in [2, 4]. The FD-method consists of the following steps: 1) the function $q(x)$ (the potential) is approximated by a piece-wise constant function $\overline{q}(x)$; 2) problem (3) is embedded into the parametric family of problems with the potential $\overline{q}(x)+t(q(x)-\overline{q}(x))$ (the perturbation of the operator); 3) the solution of the perturbed family is represented as a power series in $t$ with the coefficients which are solutions of a recursive sequence of problems with the potential $\overline{q}(x)$; 4) by setting $t=1$ one obtains the series representations of the solution of the original problem. Below we consider the simplest case of the approximation of $q(x)$ by the constant $\overline{q}(x)\equiv 0$. Note that in this case the idea of the FD-method is close to the idea of the Adomian method [1, 20].

The exact solution of the eigenvalue problem is then represented by the series

[TABLE]

provided that these series converge. The sufficient conditions for the convergence of the series (6) will be presented later in Section 3 (see (31)). The approximate solution to problem (3) is represented by a pair of corresponding truncated series

[TABLE]

which is called an approximation of rank $m$. The summands of series (6), (7) are the solutions of the following recursive sequence of problems:

[TABLE]

where

[TABLE]

and $\delta_{\alpha,\beta}$ is the Kronecker delta.

For $s-1<\alpha<0$ the second boundary condition in (8) is given as follows:

[TABLE]

The initial approximation $u_{n}^{(0)}(x),$ $\lambda_{n}^{(0)}$ is the solution of the so called base problem, that is,

[TABLE]

and for $s-1<\alpha<0$ the second boundary condition in (11) is given by

[TABLE]

We call the representation of the exact solution by the series (6), i.e., by the sequences $\{u_{n}^{(j)},j=0,1,...\}$, and we call $\{\lambda_{n}^{(j)},j=0,1,...\}$ we call the exact discrete scheme (EDisS) and the approximate representation by the corresponding truncated series the truncated discrete scheme (TDisS). The idea in broad sense is related to the idea of exact and truncated difference schemes (EDS and TDS) from [10] but the algorithms are rather different.

The solution of (11) is the generalized Jacobi function

[TABLE]

where $P_{n}^{(\alpha,\beta)}(x)$ is the classical Jacobi polynomial. This is a generalization of the eigenfunctions of the classical Jacobi differential operator (see (17) and Theorem 3.2 in [6]). Since the operator in (8) is singular these problems are solvable under the solvability condition

[TABLE]

The solution of (11) can be found up to a constant whose value is determined by the orthogonality condition

[TABLE]

The polynomials $P_{n}^{(\alpha,\beta)}(x)$ are orthogonal with the Jacobi weight function $\omega^{(\alpha,\beta)}(x)$ (see [12]), i.e.

[TABLE]

where

[TABLE]

The normalized solution of problem (11) is

[TABLE]

Using the solvability condition and the orthogonality property (14) we obtain from (8), (9) the following formula for the eigenvalue corrections

[TABLE]

We consider the solution $u_{n}^{(j+1)}(x)$ of problem (8), (9) satisfying the orthogonality condition (14) as an element of the Hilbert space $L_{w}^{2}\left[-1,1\right]$ with the scalar product $(u,v)=\int_{-1}^{1}\omega^{(-\alpha,\beta)}(x)u(x)v(x)dx$. Now, using (14) we obtain the following representation for the eigenfunction corrections

[TABLE]

\theoremstyle

definition

Remark 2.1

Note that (6) can be considered as an exact discrete scheme for the boundary value problem (3) where we have switched from a nonlinear problem to the sequence $u_{n}^{(j)}(x)$ of solutions of linear problems. The truncated series (7) can be viewed as an algorithmically realizable truncated discrete scheme.

3 Convergence of the FD-method with $\overline{q}(x)\equiv 0$

In this section we find sufficient convergence conditions of the FD-method as well as the accuracy estimates.

From (19) we obtain

[TABLE]

where

[TABLE]

and

[TABLE]

Here we have $M_{n}=1$ if coincidentally $s=1/2$ and $\alpha=\beta$ or $s=1/2$ and $\alpha=\beta+1$.

The asymptotic formula from [22], namely,

[TABLE]

together with (22) yields

[TABLE]

From (20) and (18) we obtain

[TABLE]

where $\left\|q\right\|_{\infty}=\mathop{\max}\limits_{x\in\left[-1,1\right]}\left|q(x)\right|$. The estimates (25) lead to

[TABLE]

Substituting here $U_{j}=\left(\left\|q\right\|_{\infty}M_{n}\right)^{-j}\left\|u_{n}^{\left(j\right)}\right\|$ and $U_{0}=\left\|u_{n}^{\left(0\right)}\right\|=1$ and replacing the new variables by the majorant variables subject to $U_{j}\leq\overline{U}_{j}$ and $\overline{U}_{0}=U_{0}=1$, we come to the majorant equation

[TABLE]

The solution of this convolution type equation is (see e.g. [23, p. 159-161,210] and [21])

[TABLE]

Returning to the old variables we obtain the following estimate for the solution of (26)

[TABLE]

and then from (25) the next estimate for the eigenvalue corrections

[TABLE]

The last part of inequalities (29) and (30) was obtained using the reflections like those from the proof of the Wallis formula (see e.g. [9, p. 344]).

From (29) and (30) follows the next assertion.

Theorem 3.1

Under the condition

[TABLE]

the FD-method converges super-exponentially and his accuracy is characterized by the estimates

[TABLE]

and

[TABLE]

Remark 3.2

The relation (24) shows that for a fixed $s\in(1/2,1)$ there exists some $n_{0}$ such that for all $n\geq n_{0}$ condition (31) of the theorem is fulfilled and the FD-method converges super-exponentially. The error estimates (32) and (33) show that the accuracy of our method increases as we increase the eigenvalue number. This remarkable property can be lost for problem (3) with a fixed $0<s<1/2$ or if coincidentally $s=1/2$ and $\alpha=\beta$ or $s=1/2$ and $\alpha=\beta+1$. For these cases condition (31) for a given $n$ can be fulfilled for $\left\|q\right\|_{\infty}$ small enough. If this is not the case then the trivial variant of the FD-method with $\overline{q}(x)\equiv 0$ is divergent and one should apply the variant with a piecewise constant $\overline{q}(x)$ (see [15, 16]).

4 Symbolic algorithm of the FD-method with $\overline{q}(x)\equiv 0$

A new algorithmic realization of the FD-method for problem (3) with $\bar{q}(x)\equiv 0,\alpha=\beta=0$ and with the potential $q(x)=x^{2}$ was proposed in [7]. It was shown that the corrections to eigenfunctions of the FD-method are linear combinations of Legendre polynomials with a number of summands depending on the degree of the potential polynomial and on the correction number $j$. The coefficients of these linear combinations can be represented recursively through the corresponding coefficients computed at previous steps using the basic arithmetical operations only. The approximations of the eigenvalues are represented through these coefficients too. Thus, the method uses these arithmetical operations and does not involve solutions of any supplementary BVPs and computation of any integrals unlike the traditional implementations. In this sense our algorithm is a symbolic one since it operates with the expansion coefficients only.

Further we extend these results and describe such algorithm for problem (3) with the potential

[TABLE]

where $c_{l}$, $l=0,1,...,r$ are real constant coefficients. In this case the FD-method is exactly realizable in the sense that the corrections to the eigenfunctions can be explicitly represented as linear combinations of the generalized Jacobi functions ${}^{+}J_{n}^{(-\alpha,\beta)}(x)$ with coefficients which can be explicitly represented as polynomials of $c_{l}$, $l=0,1,...,r$ with rational coefficients.

The solution $u_{n}^{(0)}(x)$ of the base problem (17) (up to an arbitrary multiplicative constant) can be represented as $u_{n}^{(0)}(x)=a_{n}^{(0)}{}^{+}J_{n}^{(-\alpha,\beta)}(x)$, $a_{n}^{(0)}\neq 0$ (see Lemma 4.1 below). Note that we can renounce here the normalizing of the eigenfunctions of the base problem which we have used above to prove the convergence of the FD-method. For the orthonormal eigenfunctions $\left\{u_{n}^{(0)}(x)\right\}_{n=1}^{\infty}$ in accordance with (17) we have $a_{n}^{(0)}=\left(\gamma_{n}^{(\alpha,\beta)}\right)^{-1/2}$.

Further we use the following recurrence relation for the Jacobi polynomials (see, e.g., [5, Section 10.8])

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with $k=\max(n-r,0),\max(n-r,0)+1,...,n+r$ and $r=2,3,...$. The $(j+1)$-th step of the FD-method consists of solving the BVP (8), (9) with

[TABLE]

Using the solution of the base problem (17) and considering the properties of the classical Jacobi polynomials $P_{n}^{(\alpha,\beta)}(x)$, of the generalized Jacobi functions ${}^{+}J_{k}^{(-\alpha,\beta)}(x)$, the properties of the fractional Sturm-Liouville type operator ${}^{+}L_{\alpha,\beta}^{2s}u(x)$ as well as the relations (35) and (36) we come to the following assertion.

Lemma 4.1

The solution of problem (8) with the right-hand sides (36) can be represented by

[TABLE]

where $a_{n}^{(j+1)}=0$, $j=0,1,2,...$, and $a_{n}^{(0)}\neq 0$.

Let us substitute (37) into (36) and use (35) with the aim to represent (36) through the generalized Jacobi functions ${}^{+}J_{k}^{(-\alpha,\beta)}(x)$ only. Then, by changing the summation order we arrive at

[TABLE]

where $\left[\sigma\right]$ is the entire part of $\sigma$. Substituting (37) into (14) and comparing the coefficients in the front of the generalized Jacobi functions ${}^{+}J_{k}^{(-\alpha,\beta)}(x),$ where $k$ is an integer in $\left[\max(0,\;n-r(j+1)),n+r(j+1)\right],$ we obtain the following formulas for the eigenvalue corrections

[TABLE]

[TABLE]

[TABLE]

as well as the following recursive representations for the coefficients in (37):

[TABLE]

with $m=n+rj+1,n+rj+2,...,n+r(j+1),\;\;m\neq n,\;\;j=0,1,2,...$,

[TABLE]

with $m=\max(0,n-rj),\max(0,n-rj)+1,...,n+rj,\;\;m\neq n,\;\;j=0,1,2,...$,

[TABLE]

with $m=\max(0,n-r(j+1)),\max(0,n-r(j+1))+1,...,\max(0,n-rj-1),\;\;m\neq n,\;\;j=0,1,2,...$, and

[TABLE]

Formulas (37)-(44) represent the symbolic algorithm of the FD-method for problem (3) with the polynomial potential (34).

5 Numerical examples

Example 5.1

We consider the eigenvalue problem (3) with $\alpha=1/2,\beta=0,$ $s=3/4$, and with potential (34), where $r=3$, $c_{3}=1/4$ and $c_{l}=0$, $l=0,1,2$. In this case the sufficient convergence condition $r_{n}<1$ (see (31)) is fulfilled for $n\geq 1$ (see Table 1). We have $\left\|q\right\|_{\infty}=1/4$ and in (22)

[TABLE]

The computer algebra system Maple 17 was used, where the corrections (37) and (39) to the eigenpairs were computed exactly as mathematical expressions, i.e., we had no rounding errors. Below we give the eigenvalue corrections (39) and the coefficients of (37) for some first steps of the FD-method:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The approximations $\overset{m}{\lambda}_{n}$ of rank $m$ to the eigenvalues $\lambda_{n}$ for $n=0,1,2,3,4,10$ using the FD-method of rank $m=20$ are given in Table 1. The behavior of the corrections $\lambda_{n}^{(m)},\;u_{n}^{(m)}(x)$, $m=0,1,...,10,20$ for the eigenpairs $\lambda_{n},$ $u_{n}(x),$ $n=0,10,$ is illustrated in Table 2, which contains the corrections to these eigenvalues as well as the norms (21) of the eigenfunction corrections $||u_{n}^{(m)}||$. One can observe that the convergence of the method improves when the eigenpair number $n$ increases (see Remark 3.1). Table 2 shows that the method converges for $n=0$ too, i.e., the conditions of Theorem 3.1 are rough and can be improved.

Example 5.2

Let us consider the eigenvalue problem (3) in the case $s-1<\alpha<0$ with the second boundary condition (5), namely with $\alpha=-1/8$, $\beta=-1/2$, $s=3/4$ and with potential (34), where $r=3$, $c_{l}=1/12$, $l=0,1,2,3$. In (22) we have

[TABLE]

and $\left\|q\right\|_{\infty}=\mathop{\max}\limits_{x\in\left[-1,1\right]}\left|q(x)\right|=1/3$. Table 3 shows that the sufficient convergence condition (31) is fulfilled, i.e., $r_{n}<1$ for $n\geq 2$.

As in Example 5.1 the computer algebra system Maple 17 was used, where the corrections to the eigenpairs (37) and (39) were computed exactly as mathematical expressions. The approximations $\overset{m}{\lambda}_{n}$ of rank $m$ to the eigenvalues $\lambda_{n}$ for $n=0,1,2,3,4,10$ using the FD-method of ranks $m=20$ and $m=30$ are given in Table 3. The normalized solution of problem (11) with the second boundary condition (12) is given by (17). The behavior of the corrections $\lambda_{n}^{(m)}$, $m=0,1,...,10,20,30,$ for the eigenvalues $\lambda_{n},$ $n=0,1,2,3,4,$ is illustrated in Table 4. One can observe that the method converges for $n=0,1$ too, i.e. the conditions of Theorem 3.1 are rough and can be improved. Besides we can see that the convergence of the method improves when the eigenpair number $n$ increases (see Remark 3.1).

Example 5.3

Let us consider problem (3) with $\alpha=\beta=0,\,\,s=3/4$ and with the potential $q(x)=(\text{sgn}(x)+1)/2$. The potential is not polynomial but we can use the general idea of the symbolic algorithms above representing the potential as a series with respect to the base problem solution. Then we represent the solution $u_{n}^{(j+1)}(x)$ of the problems (8) by a series of the same type and obtain some recurrence relations for the series coefficients like (40)-(43) but with $\infty$ as the upper bound for $m$. By restricting $m$ by some finite $N$ we have an approximate symbolic algorithm which we describe below in detail.

The solution of the base problem in this case is

[TABLE]

where $P_{n}(x)$ is the Legendre polynomial. We look for the solution of problem (8), (9) in the form

[TABLE]

which satisfies the condition

[TABLE]

Let us find for the right-hand side of (8) the following expansion like (similar to 46)

[TABLE]

To find the coefficients of this expansion we use the formula

[TABLE]

with

[TABLE]

(see, e.g., [5], Volume 1, Section 3.12, formula (15)). Then we obtain

[TABLE]

Now, from (8), (9) and (48), we obtain the following exact symbolic algorithm

[TABLE]

The truncated sums with $N$ summands represent a practically realizable symbolic algorithm.

The computations for the case $n=0,N=64$, using Maple with 32 significant digits, resulted in

[TABLE]

[TABLE]

[TABLE]

The correction to the eigenvalue in the 16-th step (i.e. the 16-th summand of the series) is

[TABLE]

and all corrections with odd numbers are equal to zero. It is easy to obtain analytically the corrections $\lambda_{0}^{(0)}=\frac{3\sqrt{2}}{8\pi}\Gamma^{2}\left(\frac{3}{4}\right),\,\,u_{0}^{(0)}(x)=1,\;\lambda_{0}^{(1)}=1.$

Acknowledgment. We want to thank the anonymous referees for their helpful remarks which contributed to improvement of the article.