Weakly regular Sturm-Liouville problems: a corrected spectral matrix method

TL;DR

This paper introduces a spectral matrix method for solving weakly regular Sturm-Liouville eigenproblems with unbounded potentials, providing accurate eigenvalue approximations and validated through numerical experiments.

Contribution

It proposes a novel Galerkin spectral matrix approach with error analysis and a correction formula for eigenvalues in weakly regular Sturm-Liouville problems.

Findings

Convergence analysis confirms the method's accuracy.

The correction formula improves eigenvalue estimates.

Numerical experiments validate the effectiveness of the approach.

Abstract

In this paper, we consider weakly regular Sturm-Liouville eigenproblems with unbounded potential at both endpoints of the domain. We propose a Galerkin spectral matrix method for its solution and we study the error in the eigenvalue approximations it provides. The result of the convergence analysis is then used to derive a low-cost and very effective formula for the computation of corrected numerical eigenvalues. Finally, we present and discuss the results of several numerical experiments which confirm the validity of the approach.

Tables

Table 1. Table 1: Coefficients ξ n , η n subscript 𝜉 𝑛 subscript 𝜂 𝑛 \xi_{n},\eta_{n} and θ n subscript 𝜃 𝑛 \theta_{n} for some BCs.

BCs		${\xi }_n$	${\eta }_n$	${\theta }_n$
$y\left(\pm 1\right)=0$	$n\ge 0$	$1$	$0$	$-1$
$y^{\prime }\left(\pm 1\right)=0$	$n\ge 0$	$1$	$0$	$-\frac{n\left(n+1\right)}{\left(n+2\right)\left(n+3\right)}$
$\begin{array}{c}y\left(-1\right)=0\\ y^{\prime }\left(1\right)=0\end{array}$	$n\ge 0$	$1$	$\frac{\left(2n+3\right)}{{\left(n+2\right)}^2}$	$-{\left(\frac{n+1}{n+2}\right)}^2$
$\begin{array}{c}{y}^{\prime }\left(-1\right)=0\\ y\left(1\right)=0\end{array}$	$n\ge 0$	$1$	$-\frac{\left(2n+3\right)}{{\left(n+2\right)}^2}$	$-{\left(\frac{n+1}{n+2}\right)}^2$
$\begin{array}{c}{y}^{\prime }\left(-1\right)=y\left(-1\right)\\ y^{\prime }\left(1\right)=\mathrm{\hspace{0.17em}0}\end{array}$	$n\ge 0$	$1$	$\frac{2\left(2n+3\right)}{{\left(n+2\right)}^2\left(2+\left(n+1\right)\left(n+3\right)\right)}$	$-\frac{{\left(n+1\right)}^2\left(2+n\left(n+2\right)\right)}{{\left(n+2\right)}^2\left(2+\left(n+1\right)\left(n+3\right)\right)}$
	$n=0$	$\frac{2}{3}$	$1$	$\frac{1}{3}$
$y^{\prime }\left(\pm 1\right)=y\left(\pm 1\right)$
	$n\ge 1$	$1$	$\frac{4\left(2n+3\right)}{\left(n+1\right){\left(n+2\right)}^2\left(n+3\right)-4}$	$-\frac{n{\left(n+1\right)}^2\left(n+2\right)+4}{\left(n+1\right){\left(n+2\right)}^2\left(n+3\right)-4}$

Table 2. Table 2: Order of convergence for problems with potential ( 58 ).

$y^{\prime }\left(-1\right)=y\left(1\right)=0,p=p_L=6-4\times \left(1/4\right)=5$
$N$	$\delta {\lambda }_{5,N}$	order	$\delta {\lambda }_{10,N}$	order	$\delta {\lambda }_{20,N}$	order
$49$	$4.4416E-06$	$5.002$	$5.5319E-06$	$4.993$	$4.8732E-06$	$4.969$
$99$	$1.3859E-07$	$5.001$	$1.7368E-07$	$4.999$	$1.5556E-07$	$4.996$
$199$	$4.3295E-09$	$4.999$	$5.4312E-09$	$4.997$	$4.8743E-09$	$5.000$
$399$	$1.3534E-10$	–	$1.7002E-10$	–	$1.5234E-10$	–
$y\left(-1\right)=y^{\prime }\left(1\right)=0,p=p_R=6-4\times \left(3/4\right)=3$
$N$	$\delta {\lambda }_{5,N}$	order	$\delta {\lambda }_{10,N}$	order	$\delta {\lambda }_{20,N}$	order
$49$	$2.1678E-03$	$3.002$	$7.9981E-03$	$2.999$	$1.2424E-02$	$2.985$
$99$	$2.7065E-04$	$3.000$	$1.0005E-03$	$3.000$	$1.5693E-03$	$2.999$
$199$	$3.3820E-05$	$3.000$	$1.2505E-04$	$3.000$	$1.9636E-04$	$3.000$
$399$	$4.2273E-06$	–	$1.5631E-05$	–	$2.4547E-05$	–
$y^{\prime }\left(\pm 1\right)=0,p=p_R=3$
$N$	$\delta {\lambda }_{5,N}$	order	$\delta {\lambda }_{10,N}$	order	$\delta {\lambda }_{20,N}$	order
$49$	$1.7360E-03$	$3.001$	$7.4699E-03$	$3.003$	$1.2385E-02$	$2.989$
$99$	$2.1688E-04$	$3.000$	$9.3168E-04$	$3.001$	$1.5602E-03$	$2.999$
$199$	$2.7105E-05$	$3.000$	$1.1637E-04$	$3.000$	$1.9510E-04$	$3.000$
$399$	$3.3881E-06$	–	$1.4543E-05$	–	$2.4386E-05$	–
$y^{\prime }\left(\pm 1\right)=y\left(\pm 1\right),p=p_R=3$
$N$	$\delta {\lambda }_{5,N}$	order	$\delta {\lambda }_{10,N}$	order	$\delta {\lambda }_{20,N}$	order
$49$	$2.0295E-03$	$3.002$	$8.1569E-03$	$3.004$	$1.2697E-02$	$2.990$
$99$	$2.5338E-04$	$3.000$	$1.0168E-03$	$3.001$	$1.5983E-03$	$3.000$
$199$	$3.1663E-05$	$3.000$	$1.2698E-04$	$3.000$	$1.9983E-04$	$3.000$
$399$	$3.9576E-06$	–	$1.5869E-05$	–	$2.4977E-05$	–

Table 3. Table 3: Order of convergence for problems with potential ( 59 ).

$y^{\prime }\left(-1\right)=y\left(1\right)=0,p=p_L=6-4\times \left(1/2\right)=4$
$N$	$\delta {\lambda }_{5,N}$	order	$\delta {\lambda }_{10,N}$	order	$\delta {\lambda }_{20,N}$	order
$49$	$6.1520E-05$	$3.998$	$7.5495E-05$	$3.993$	$6.9334E-05$	$3.974$
$99$	$3.8510E-06$	$3.999$	$4.7406E-06$	$3.999$	$4.4108E-06$	$3.997$
$199$	$2.4080E-07$	$4.000$	$2.9654E-07$	$4.000$	$2.7632E-07$	$4.000$
$399$	$1.5052E-08$	–	$1.8538E-08$	–	$1.7275E-08$	–
$y\left(-1\right)=y^{\prime }\left(1\right)=0,p=p_R=6-4\times \left(7/8\right)=2.5$
$N$	$\delta {\lambda }_{5,N}$	order	$\delta {\lambda }_{10,N}$	order	$\delta {\lambda }_{20,N}$	order
$49$	$6.9840E-03$	$2.503$	$2.4576E-02$	$2.502$	$4.1088E-02$	$2.491$
$99$	$1.2317E-03$	$2.501$	$4.3385E-03$	$2.501$	$7.3081E-03$	$2.500$
$199$	$2.1760E-04$	$2.500$	$7.6648E-04$	$2.500$	$1.2921E-03$	$2.500$
$399$	$3.8460E-05$	–	$1.3547E-04$	–	$2.2839E-04$	–
$y^{\prime }\left(\pm 1\right)=0,p=p_R=2.5$
$N$	$\delta {\lambda }_{5,N}$	order	$\delta {\lambda }_{10,N}$	order	$\delta {\lambda }_{20,N}$	order
$49$	$5.6379E-03$	$2.507$	$2.3091E-02$	$2.510$	$4.0904E-02$	$2.497$
$99$	$9.9209E-04$	$2.502$	$4.0545E-03$	$2.503$	$7.2462E-03$	$2.501$
$199$	$1.7508E-04$	$2.501$	$7.1512E-04$	$2.501$	$1.2797E-03$	$2.501$
$399$	$3.0932E-05$	–	$1.2633E-04$	–	$2.2612E-04$	–
$y^{\prime }\left(\pm 1\right)=y\left(\pm 1\right),p=p_R=2.5$
$N$	$\delta {\lambda }_{5,N}$	order	$\delta {\lambda }_{10,N}$	order	$\delta {\lambda }_{20,N}$	order
$49$	$6.4759E-03$	$2.506$	$2.5025E-02$	$2.510$	$4.2056E-02$	$2.498$
$99$	$1.1403E-03$	$2.502$	$4.3929E-03$	$2.503$	$7.4452E-03$	$2.502$
$199$	$2.0131E-04$	$2.501$	$7.7477E-04$	$2.501$	$1.3146E-03$	$2.501$
$399$	$3.5571E-05$	–	$1.3686E-04$	–	$2.3228E-04$	–

Equations

$$-y^{\prime \prime }(x)+q(x)y(x)$$

$${\alpha }_{L}y(-1)+{\beta }_L{y}^{\prime }(-1)$$

$${\alpha }_{R}y(-1)+{\beta }_R{y}^{\prime }(-1)$$

$$q(x)=\sum _{i=1}^S\frac{g_i(x)}{(1-x{)}^{{\beta }_i}(1+x{)}^{{\gamma }_i}},{\beta }_i,{\gamma }_i<1,i=1,\dots ,S,$$

$$\left\{{\lambda }_1<{\lambda }_2<{\lambda }_3<\dots \right\}$$

$${\mathcal{S}}_N$$

$$\equiv$$

$$z_N(x)=\sum _{n=0}^{N-1}{\zeta }_{n,N}{\mathcal{R}}_n(x)\approx y(x)$$

$$\sum _{n=0}^{N-1}⟨{\mathcal{R}}_m,-{\mathcal{R}}_n^{\prime \prime }+(q-{\lambda }^{(N)}){\mathcal{R}}_n⟩{\zeta }_{n,N}=0,\text{mbox}foreachm=0,\dots ,N-1.$$

$$⟨u,v⟩={\int }_{-1}^{1}u(x)v(x)dx,u,v\in L_2([-1,1]),$$

$$\left(A_{N}+Q_{N}\right)\mbox{\boldmath$\zeta$}_{N}=\lambda^{(N)}B_{N}\mbox{\boldmath$\zeta$}_{N}$$

$$A_N=\left(a_{mn}\right),B_N=\left(b_{mn}\right),Q_N=\left(q_{mn}\right),m,n=0,\dots ,N-1,$$

$$a_{mn}=-⟨{\mathcal{R}}_m,{\mathcal{R}}_n^{\prime \prime }⟩,b_{mn}=⟨{\mathcal{R}}_m,{\mathcal{R}}_n⟩,q_{mn}=⟨{\mathcal{R}}_m,q{\mathcal{R}}_n⟩.$$

$${\mathcal{R}}_n(x)={\xi }_n{\mathcal{P}}_n(x)+{\eta }_n{\mathcal{P}}_{n+1}(x)+{\theta }_n{\mathcal{P}}_{n+2}(x)$$

$${\mathcal{P}}_j(1)=(-1{)}^j{\mathcal{P}}_j(-1)=1,{\mathcal{P}}_j^{\prime }(1)=(-1{)}^{j-1}{\mathcal{P}}_j^{\prime }(-1)=j(j+1)/2$$

$${\xi }_n$$

$${\theta }_n$$

$${\eta }_n$$

$$a_{nn}$$

$$a_{nn}=4\left(n+\frac{3}{2}\right)\left(1+O\left(n^{-1}\right)\right),n\gg 1.$$

$${\hat{b}}_n=⟨{\mathcal{P}}_n,{\mathcal{P}}_n⟩=2/(2n+1),$$

$$\hat{B}_{N}=\left(\begin{array}[]{ccc}\hat{b}_{0}\\ &\ddots\\ &&\hat{b}_{N+1}\end{array}\right),\quad R_{N}=\left(\begin{array}[]{ccccc}\xi_{0}\\ \eta_{0}&\ddots\\ \theta_{0}&\ddots&\ddots\\ &\ddots&\ddots&\xi_{N-1}\\ &&\ddots&\eta_{N-1}\\ &&&\theta_{N-1}\end{array}\right),$$

$$B_N=R_N^T{\hat{B}}_N{R}_N.$$

$$Q_N$$

$${\hat{q}}_{mn}=⟨{\mathcal{P}}_m,q{\mathcal{P}}_n⟩=\sum _{i=1}^S⟨{\mathcal{P}}_m,g_i{\mathcal{P}}_n{⟩}_{({\beta }_i,{\gamma }_i)}\equiv \sum _{i=1}^S{\hat{q}}_{mn}^{(i)},$$

$$⟨u,v{⟩}_{(\beta ,\gamma )}={\int }_{-1}^1\frac{u(x)v(x)}{(1-x{)}^{\beta }(1+x{)}^{\gamma }}dx,\beta ,\gamma <1.$$

$${\hat{\mathbf{q}}}_n$$

$$\mathcal{H}$$

$$\hat{{\mathbf{q}}}_{n+1}=\frac{2n+1}{n+1}\,\mathcal{H}\hat{{\mathbf{q}}}_{n}-\frac{n}{n+1}\,\hat{{\mathbf{q}}}_{n-1},\qquad n\geq 0.\qquad\mbox{\hfill$\Box$}\\$$

$$\hat{{\mathbf{q}}}_{0}=g_{1}(\mathcal{H})\left(\begin{array}[]{c}\hat{q}_{0}^{(1)}\\ \hat{q}_{1}^{(1)}\\ \vdots\end{array}\right),\quad\hat{q}_{m}^{(1)}=\langle{\cal{P}}_{m},{\cal{P}}_{0}\rangle_{({\beta}_{1},{\gamma}_{1})}=\langle{\cal{P}}_{m},1\rangle_{({\beta}_{1},{\gamma}_{1})}.$$

$${\hat{q}}_m^{(1)}=\frac{(-1{)}^m{2}^{1-{\gamma }_1}({\gamma }_1{)}_m}{(1-{\gamma }_1{)}_{m+1}}$$

$${\hat{q}}_m^{(1)}=\frac{2^{1-{\beta }_1}({\beta }_1{)}_m}{(1-{\beta }_1{)}_{m+1}};$$

$$\hat{q}_{m}^{(1)}=\alpha\,_{3}F_{2}\left(\begin{array}[]{c}-m,\,\,1+m,\,\,1-{\beta}_{1}\\ 1,\,\,2-{\beta}_{1}-{\gamma}_{1}\end{array};1\right)$$

Full text

Keywords: Sturm-Liouville eigenproblems, spectral matrix methods, Legendre polynomials, acceleration of convergence.

MSC: 65L15, 65L60, 65L70, 65B99

1 Introduction

Recently, the author studied a corrected spectral matrix method for solving weakly regular and singular Sturm-Liouville problems defined over the bounded domain $(-1,1)$ with an unbounded potential at the left endpoint, [7]. The numerical results provided by such technique are definitely satisfactory for weakly regular problems. This suggested to study a generalization of the method for the approximation of the eigenvalues and of the eigenfunctions of problems of the following type

[TABLE]

where the potential $q$ is given by

[TABLE]

with functions $g_{i}$ at the numerators that are analytical inside and on a Bernstein ellipse containing $[-1,1].$ In the literature, problems of this type with $q$ unbounded at least at one endpoint are sometimes called weakly regular and it is well known that their spectrum is composed by real and simple eigenvalues which can be ordered as an increasing sequence tending to infinity. We will number them starting from index $k=1,$ i.e. we will call

[TABLE]

the exact spectrum of (1)–(4).

In [7], we considered the case $S=2$ with $\beta_{1}=\gamma_{1}=\beta_{2}=0,$ namely problems with a potential of the form $q(x)=g_{1}(x)+g_{2}(x)/(1+x)^{{\gamma}_{2}},$ and a special algorithm for ${\gamma}_{2}\in(0,1)$ and $y(-1)\neq 0$ was derived. As remarked in the same paper, the obtained results appear to be more reliable than those given by well-known and well-established general-purpose codes based on shooting techniques like, for example, the MATSLISE2 [5, 6], the SLEDGE [8] and the SLEIGN2 [1] ones. A possible explanation is that the common basic idea in them implemented is essentially the selection of suitable layers. In particular, if $q$ is unbounded at both endpoints, then the approach is that of solving a suitable problem over $(-1+\varepsilon_{1},1-\varepsilon_{2})$ with $\varepsilon_{1}$ and $\varepsilon_{2}$ small positive values automatically selected, [9]. As indicated in their documentation, this may cause a loss of accuracy and from all our tests we deduced that this may be more relevant if the problem is not subject to the Dirichlet condition at the endpoint where $q$ is unbounded. As an example, in the following table we list some numerical eigenvalues for the problem with $q(x)=10/(1-x^{2})^{3/4}$ subject to $y^{\prime}(\pm 1)=0$ that we computed by using such codes with a tolerance equal to $10^{-13}.$

[TABLE]

These considerations justify the interest in generalizing the method proposed in [7] and the outline of this paper is the following. In Section 2, we recall the basic facts concerning the spectral Legendre-Galerkin matrix method introduced in [7] and we discuss the computation of the coefficient matrix that corresponds to a potential $q$ of the form in (4). An analysis of the error in the numerical eigenvalues with respect to the generalized eigenvalue problem size is carried out in Section 3. In addition, in the same section, we derive a low cost and effective procedure for an a posteriori correction of the numerical eigenvalues. Finally, in Section 4 we report and discuss the results of some numerical experiments.

2 Spectral Legendre-Galerkin method

Let $\Pi_{N+1}$ be the space of polynomials of maximum degree $N+1,$ for a fixed $N\in{\mathbb{N}},$ and let

[TABLE]

We look for an approximation of an eigenfunction $y$ of the following type

[TABLE]

where the coefficients $\zeta_{n,N}$ and the numerical eigenvalue $\lambda^{(N)}$ are determined by imposing, see (1),

[TABLE]

Here $\langle\cdot,\cdot\rangle$ is the standard inner product in $L_{2}([-1,1]),$ i.e.

[TABLE]

which is naturally suggested by the Liouville normal form of the SLP we are studying. We can write (8) as the following generalized eigenvalue problem

[TABLE]

where \mbox{\boldmath\zeta}_{N}=\left(\zeta_{0N},\ldots,\zeta_{N-1,N}\right)^{T},

[TABLE]

with

[TABLE]

The matrices $B_{N}$ and $Q_{N}$ are clearly symmetric and the same property holds for $A_{N}$ thanks to the well-known Green’s identity, [7].

The basis function ${\cal{R}}_{n}$ is chosen as follows [10]

[TABLE]

where ${\cal{P}}_{j}$ is the Legendre polynomial of degree $j$ and the three coefficients $\xi_{n},\eta_{n}$ and $\theta_{n}$ are such that ${\cal{R}}_{n}$ verifies the boundary conditions (BCs), see (5)-6). The complete discussion of the computation of such coefficients can be found in [7] where we used the fact that

[TABLE]

and we decided to use basis functions that verify $\|{\cal{R}}_{n}\|_{\infty}\leq 3$ for each $n\in{\mathbb{N}}_{0}.$ This is obtained by imposing $\|\left(\xi_{n},\eta_{n},\theta_{n}\right)^{T}\|_{\infty}=1$ and $\xi_{n}\geq 0.$ Here we simply list in Table 1 the three coefficients of the linear combination in (12) for the four problems subject to natural BCs and for two general ones of Robin type.

For later reference, it is important to underline the fact that, as soon as $n$ is sufficiently large, we always got

[TABLE]

More precisely, if the BCs are symmetric, i.e. if $\alpha_{L}\beta_{R}+\alpha_{R}\beta_{L}=0,$ then we always set $\eta_{n}=0$ so that ${\cal{R}}_{n}$ is an even or an odd function if $n$ is even or odd, respectively.

2.1 The matrices $A_{N}$ and $B_{N}$

In this section, we recall the results obtained in [7] about the entries of $A_{N}$ and $B_{N}$ in (10)-(11).

Concerning the first matrix, one immediately gets that $a_{mn}=0$ for each $m>n$ since ${\cal{R}}_{m}$ is orthogonal to any polynomial in $\Pi_{m-1},$ see (12). Consequently, $A_{N}=A_{N}^{T}$ is diagonal with diagonal entries

[TABLE]

see (13). We remark that, independently of the BCs, $a_{nn}$ satisfies

[TABLE]

Regarding $B_{N},$ it is not too difficult to verify that it is pentadiagonal. In more detail, if we let

[TABLE]

[TABLE]

then we get

[TABLE]

2.2 The matrix $Q_{N}.$

From (10)-(11), one obtains that $Q_{N}$ admits a factorization similar to the one just given for $B_{N}.$ Specifically

[TABLE]

where $R_{N}$ is defined in (22),

[TABLE]

being

[TABLE]

As done in [4, 7], it is possible to prove the following result by using the well-known recurrence relation of the Legendre polynomials.

Proposition 2.1

Let $q\in L_{1}([-1,1])$ and

[TABLE]

If we define the linear tridiagonal operator ${\mathbf{z}}\in\ell_{\infty}\mapsto\mathcal{H}\,{\mathbf{z}}\in\ell_{\infty}$ where

[TABLE]

and we let $\hat{{\mathbf{q}}}_{-1}$ be the zero sequence, then we get

[TABLE]

Now, the structure of $\mathcal{H}$ and (35) permit to determine the entire matrix $\hat{Q}_{N}$ once $\hat{q}_{m0}$ have been computed for each $m=0,1,\ldots,2N+2$ (see [4, 7] for the details). We shall proceed by discussing how we determine these values for problems with $S=1$ since the generalization is simple, see (25). In this regard, we observe that (26), [3, 16.4 formula (2)] and arguments similar to the ones used in the proof of [4, Proposition 2] allow to get that

[TABLE]

Let us assume for the moment that we have computed the values of $\hat{q}_{m}^{(1)}$ for each $m=0,\ldots,L$ with $L$ sufficiently large. Then, recalling that by assumption $g_{1}$ is analytical inside and over a Bernstein ellipse containing $[-1,1]$, we proceed in this way. We get a polynomial approximation of $g_{1}$ by transforming it in a Chebfun function [2], which is accurate up to machine precision, and then we apply the previous formula to compute the first $2N+2$ entries of $\hat{{\mathbf{q}}}_{0}.$

Concerning the computation of $\hat{q}_{m}^{(1)}$ we have to distinguish the following cases:

${\beta}_{1}={\gamma}_{1}=0:$ it is evident that $\hat{q}_{0}^{(1)}=2$ and $\hat{q}_{m}^{(1)}=0$ for each $m>0;$ 2. 2.

${\beta}_{1}=0,{\gamma}_{1}\neq 0:$ as discussed in [7] it results

[TABLE]

where $(t)_{\ell}$ is the Pochhammer symbol; 3. 3.

${\beta}_{1}\neq 0,{\gamma}_{1}=0:$ with similar computations one gets

[TABLE] 4. 4.

${\beta}_{1}{\gamma}_{1}\neq 0:$ by using [3, 16.2, formula(6)] we get

[TABLE]

where

[TABLE]

It is clear that the formulas in (37) or (38) allow to compute $\hat{q}_{m}^{(1)}$ easily. For example, if ${\beta}_{1}\neq 0$ and ${\gamma}_{1}=0$ then from (38) one gets

[TABLE]

so that it is possible to proceed recursively. On the other hand, if ${\beta}_{1}{\gamma}_{1}\neq 0$ then the computation of the Gauss hypergeometric function at right hand-side of (39) can be costly and ill-conditioned. We thus preferred to find alternative expressions. In particular, if ${\beta}_{1}={\gamma}_{1}\neq 0$ then the application of the following Whipple sum

[TABLE]

with $a=-m,\,c=1-{\gamma}_{1}$ and $e=1$ gives

[TABLE]

Therefore, if $m$ is odd then $\hat{q}_{m}^{(1)}=0$ (this was indeed already evident from its definition in (36) with ${\beta}_{1}={\gamma}_{1}$). On the other hand, if $m$ is even then

[TABLE]

Remark 2.1

If $q(x)=q(-x)$ and $N$ is even then the matrix $\hat{Q}_{N}$ in (24) is permutation similar to a $2\times 2$ block diagonal matrix with diagonal blocks of size $N/2+1.$ In addition, if the coefficients of the BCs in (2)-(3) verify $\alpha_{L}\beta_{R}+\alpha_{R}\beta_{R}=0$ then the SLP is called symmetric and, as we recalled after (19), $\eta_{n}=0$ for each $n.$ Therefore, the matrices $A_{N},$ $Q_{N}$ and $B_{N}$ are permutation similar to $2\times 2$ block diagonal matrices too, see the paragraph before (20) and (22)–(24). This implies that we can split the generalized eigenvalue problem (9) into two ones of halved size.

It remains to discuss how it is possible to avoid the evaluation of the Gauss hypergeometric function for the computation of $\hat{q}_{m}^{(1)}$ if ${\beta}_{1}{\gamma}_{1}\neq 0$ and ${\beta}_{1}\neq{\gamma}_{1}.$ In this case, even though alternative strategies are possible, we decided to write ${\cal{P}}_{m}\equiv{\cal{P}}^{(0,0)}_{m}$ as a linear combination of $\left\{{\cal{P}}_{\ell}^{(0,-{\gamma}_{1})}\right\}_{\ell=0}^{m}$ where ${\cal{P}}_{\ell}^{(0,-{\gamma}_{1})}$ is the Jacobi polynomial of degree $\ell$ with weighting function $\omega(x)=(1+x)^{-{\gamma}_{1}}.$ In other words, first of all we write

[TABLE]

Then we determine $\hat{q}_{m}^{(1)}$ as follows

[TABLE]

Now, by using formulas in [3, 16.4] we obtain

[TABLE]

[TABLE]

Hence

[TABLE]

where $\alpha$ is defined in (40),

[TABLE]

It is not too difficult to verify that

[TABLE]

being “$\circ$” the Hadamard product and

[TABLE]

i.e. $T$ and $U$ are a lower triangular Toeplitz and an Hankel matrix, respectively. Clearly, a suitable truncation of the vectors and matrices in (42) is operated depending on the number of values of $\hat{q}_{m}^{(1)}$ that we actually need. In addition, we compute $\nu_{\ell},$ $t_{r}$ and $u_{r}$ recursively. Finally, it is worth mentioning that an algorithm similar to the one described in [12] can be used for the matrix-vector product in (42).

3 Error analysis and computation of corrected numerical eigenvalues.

We now study the behavior of the error in the resulting numerical eigenvalues as $N$ increases and for a fixed index. In particular, we consider weakly regular problems with a potential of the type specified in (4) which is unbounded at least at one endpoint. The analysis that we are going to present will be also used to derive a very effective and efficient procedure for an a posteriori correction of the numerical eigenvalues.

Let $\lambda^{(N)}$ be the approximation of the exact eigenvalue $\lambda$ as $N$ increases and let $y$ be the corresponding exact eigenfunction having the following expansion

[TABLE]

The following first result can be proved by using arguments similar to the ones considered in [7].

Theorem 3.1

If $N$ is sufficiently larger than the index of the eigenvalue then, see (7) and (20),

[TABLE]

where

[TABLE]

The asymptotic estimate that we are going to prove in the next theorem is fundamental for proceeding.

Theorem 3.2

Let $\psi\in C^{\infty}(-1,1)\bigcap C^{1}[-1,1].$ If $\psi(-1)\,\psi(1)\neq 0,$ ${\beta},{\gamma}<1$ and if $n$ is sufficiently large then

[TABLE]

where

[TABLE]

being

[TABLE]

Proof: Recalling the definition of ${\cal{R}}_{n}$ in (12), let us consider first of all $\langle{\cal{P}}_{n},\psi\rangle_{({\beta},{\gamma})}.$ In this regard, if we use the results proved in [11], (36)–(38) and we assume that $n$ is sufficiently large then we get

[TABLE]

Now, it is known that the ratio of two gamma functions satisfies

[TABLE]

Therefore, if we use it with $z=n+1/2$ then we obtain

[TABLE]

with

[TABLE]

This implies that to determine an estimate for $\langle{\cal{R}}_{n},\psi\rangle_{({\beta},{\gamma})}$, we have to study these terms

[TABLE]

We recall that if $n$ is sufficiently large then $\xi_{n}=1,$ see (14). In addition, by using the formulas in [7], see also (15)-(19), it is possible to verify with some computations that

if ${\cal{R}}_{n}(-1)=(-1)^{n}(\xi_{n}-\eta_{n}+\theta_{n})\neq 0$ then

[TABLE] 2. 2.

if ${\cal{R}}_{n}(1)=\xi_{n}+\eta_{n}+\theta_{n}\neq 0$ then

[TABLE] 3. 3.

$\xi_{n}-\theta_{n}=2\left(1+O\left(n^{-1}\right)\right).$

Therefore, see (48) and (49),

[TABLE]

The statement follows by collecting all these partial results. $\Box$

It must be underlined that $\hat{\omega}(\delta_{0},\delta_{1})=0$ if $-\delta_{0}\in{\mathbb{N}}_{0},$ see (48). This implies that one or both the terms at the right hand-side of (46) can be zero. For our purposes, this does not constitute a problem since in the convergence analysis that we are going to prove such terms are surely negligible with respect to the others.

We need the following notation to proceed: for each $i=1,\ldots,S,$ let

[TABLE]

i.e. let $r_{i}$ and $\ell_{i}$ be the multiplicities of $x=1$ and $x=-1,$ as zeros of $g_{i},$ respectively. We are now ready for proving the following theorem.

Theorem 3.3** (Convergence)**

Let assume the potential in (4) is unbounded at least at one endpoint and let consider the following subsets of $\left\{1,2,\ldots,S\right\}$

[TABLE]

If $N$ is sufficiently larger than the index of the eigenvalue then

[TABLE]

where, see (47)-(48) and (50),

[TABLE]

being $\inf\emptyset=+\infty$ by convention.

Proof: From the definition of $\kappa_{\pm},$ (7) and (43) we get

[TABLE]

so that recalling (44)-(45) we must determine an estimate for

[TABLE]

To this end, we apply (46) with $({\beta},{\gamma})=(\hat{{\beta}}_{i},\hat{{\gamma}}_{i})$ and $\psi=\hat{g}_{i}\hat{y},$ or with $({\beta},{\gamma})=(\hat{{\beta}}_{j},\hat{{\gamma}}_{j})$ and $\psi=\hat{g}_{j}\hat{z}_{N}.$ In this way, recalling also (21), we obtain

[TABLE]

where

[TABLE]

We now use the following integral estimates with suitable $\bar{p}>0:$

[TABLE]

In particular, we apply the first one with $\bar{p}=(p(\hat{{\gamma}}_{i})+p(\hat{{\gamma}}_{j}))/2,$ or $\bar{p}=(p(\hat{{\beta}}_{i})+p(\hat{{\beta}}_{j}))/2$ and the second estimate with $\bar{p}=(p(\hat{{\gamma}}_{i})+p(\hat{{\beta}}_{j}))/2$ or with $\bar{p}=(p(\hat{{\beta}}_{i})+p(\hat{{\gamma}}_{j}))/2.$ This leads to, see (45),

[TABLE]

Therefore

[TABLE]

and the statement follows by observing that the principal term of such summation behaves like $O\left((N+1)^{-p}\right)$ where $p$ is defined in (51). $\Box$

As done in [7], we now discuss how we can use (56) to improve the accuracy of the numerical eigenvalues. The approach is that of considering the following normalization for the numerical and the exact eigenfunctions, see (52)-(53),

[TABLE]

By using the orthogonality of the Legendre polynomials, the estimates $\langle z_{N},y\rangle\approx 1$ and $\hat{y}(\pm 1)\approx\hat{z}_{N}(\pm 1)$ follow and consequently the next formula for the computation of corrected numerical eigenvalues

[TABLE]

where $\hat{\Delta}_{ij}$ is obtained from $\bar{\Delta}_{ij}$ via the substitutions $\hat{y}(\pm 1)\rightarrow\hat{z}_{N}(\pm 1),$ see (52)–(3). These are done by using $z_{N}(\pm 1)$ or $z_{N}^{\prime}(\pm 1).$ For example, if $y(-1)=0$ and $y(1)\neq 0$ then $y(x)=(1+x)\hat{y}(x)$ and $z_{N}(x)=(1+x)\hat{z}_{N}(x).$ Consequently

[TABLE]

Finally, we must say that in the actual implementation we do not consider the correction terms corresponding to values of $r_{i}$ and of $\ell_{i}$ in (50) that are greater than one since their contributions are surely irrelevant with respect to the others. The advantage is that for the computation of $\mu^{(N)}$ we do not need to evaluate derivatives of $g_{i}$ at the endpoints of order greater than one. In this way, altogether, the cost for the application of (57) is very low since it is essentially given by the evaluations of $g_{i}(\pm 1),$ eventually of $g_{i}^{\prime}(\pm 1)$ and of $\hat{z}_{N}(\pm 1)$ which is simple because the values of ${\cal{R}}_{n}(\pm 1)$ or of ${\cal{R}}_{n}^{\prime}(\pm 1)$ are known (see (12) and (13)).

4 Numerical tests

The method described was implemented in Matlab (ver.R2017a). In particular, we used routines included in the open-source Chebfun package [2] for the computation of the matrix $Q_{N}$ and we solved the generalized eigenvalue problem (9) by using the eigs function, with option “SM” for getting the ones of smallest magnitude.

The first results that we present confirm the statement of Theorem 3.3. In particular, we considered the problems with one of the following potentials

[TABLE]

subject to one of the next four BCs

[TABLE]

In addition, we used the classical formula

[TABLE]

for the numerical estimate of the order of convergence (the lower index $k$ denotes the index of the eigenvalue). The results we got for the eigenvalues of index $k=5,10,20,$ are listed in Table 2 for the first potential and in Table 3 for the second one. As one can see, such results are in perfect agreement with the statement of the theorem previously mentioned.

Concerning the problems with $q$ defined in (59) subject to the first or to the last BCs in (60), we applied the a posteriori correction, namely we computed also $\mu_{k}^{(N)}$ defined in (57). In addition, for these problems and for the subsequent ones, we evaluated the relative errors

[TABLE]

by using as reference “exact” eigenvalue the values of $\bar{\lambda}_{k}\equiv\mu_{k}^{(Nt)}$ with $Nt\gg N\geq k.$ As discussed in the introduction, this choice was motivated by the fact that the accuracy of the numerical approximations of $\lambda_{k}$ provided by the MATSLISE2 [5], the SLEDGE [8] and the SLEIGN2 [1] codes may be not sufficient for our purposes. The resulting relative errors (61) have been reported in Figure 1. In more details, in the subplots at the top of such figure, the relative errors in the approximation of the fifteenth eigenvalue are plotted versus $N$ with $N$ ranging from $50$ to $400.$ For the subplots at the bottom, instead, we fixed $N=100$ and we depict the errors for the index $k$ ranging from $1$ to $25.$ The legend of each graphic and of the subsequent ones is dashed line and solid line for the errors in the uncorrected numerical eigenvalues and in the corrected ones, respectively. These results show that the a posteriori correction is very effective from many point of views. In fact:

• •

from the subplots on the bottom one deduces that for $N=100$ and $1\leq k\leq 25$ the gain resulting from the correction is larger than two significant digits;

• •

the two subplots on the top show that $\mu_{15}^{(N)}$ is always more accurate than $\lambda_{15}^{(2N)};$

• •

the error in the corrected numerical eigenvalues decreases much faster with respect to $N$ than the error in the uncorrected ones. Concerning this point, we used a least-square approach to evaluate numerically the order of convergence $p_{\mu}$ such that

[TABLE]

For these examples, we obtained $p_{\mu}\approx 7$ for the problem subject to Neumann-Dirichlet BCs and $p_{\mu}\approx 5$ for the one subject to the unsimmetric Robin-Robin BCs.

In the following final tests, we compare the performances of our correction technique with those of the classical Richardson extrapolation given by

[TABLE]

where $p$ is specified in the convergence theorem. In particular, we considered symmetric problems with potentials

[TABLE]

and problems with the following not symmetric $q$’s

[TABLE]

The results obtained for some BCs have been reported in Figures 2 and 3 (the errors corresponding to $\rho_{k}^{(N)}$ are depicted in dotted lines). As one can see, the Richardson extrapolation requires $N$ much larger than $k$ to improve the eigenvalue approximation, say $N$ not smaller than $4k+1.$ If this is not the case then it may deteriorates drastically the accuracy of $\lambda_{k}^{(N)}.$ Furthermore, the improvement that we get with our low-cost method is undeniably larger than that obtained with Richardson.