A corrected spectral method for Sturm-Liouville problems with unbounded potential at one endpoint

TL;DR

This paper introduces a spectral matrix method using Legendre polynomials for accurately approximating eigenvalues of Sturm-Liouville problems with unbounded potentials at one endpoint, including error analysis and numerical validation.

Contribution

It develops a corrected spectral method tailored for Sturm-Liouville problems with unbounded potentials, enhancing eigenvalue approximation accuracy.

Findings

The method effectively approximates eigenvalues near singularities.

Error analysis guides the correction procedures.

Numerical experiments confirm improved accuracy and efficiency.

Abstract

In this paper, we shall derive a spectral matrix method for the approximation of the eigenvalues of (weakly) regular and singular Sturm-Liouville problems in normal form with an unbounded potential at the left endpoint. The method is obtained by using a Galerkin approach with an approximation of the eigenfunctions given by suitable combinations of Legendre polynomials. We will study the errors in the eigenvalue estimates for problems with unsmooth eigenfunctions in proximity of the left endpoint. The results of this analysis will be then used conveniently to determine low-cost and effective procedures for the computation of corrected numerical eigenvalues. Finally, we shall present and discuss the results of several numerical experiments which confirm the effectiveness 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	Name	${\xi }_n$	${\eta }_n$	${\theta }_n$
$y\left(\pm 1\right)=0$	Dirichlet–Dirichlet	1	0	$-1$
$y^{\prime }\left(\pm 1\right)=0$	Neumann–Neumann	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}$	Dirichlet-Neumann	$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}$	Neumann-Dirichlet	$1$	$-\frac{\left(2n+3\right)}{{\left(n+2\right)}^2}$	$-{\left(\frac{n+1}{n+2}\right)}^2$
$\begin{array}{c}y\left(-1\right)=-y^{\prime }\left(-1\right)\\ y\left(1\right)=y^{\prime }\left(1\right)\end{array}$	$\begin{array}{c}\text{Robin-Robin}\\ \text{symmetric}\end{array}$	$1$	$0$	$-\frac{n\left(n+1\right)-2}{\left(n+2\right)\left(n+3\right)-2}$
$\begin{array}{c}y\left(-1\right)=0\\ y\left(1\right)=y^{\prime }\left(1\right)\end{array}$	$\begin{array}{c}\text{Dirichlet-Robin}\\ n\ge 1\end{array}$	$1$	$\frac{2n+3}{{\left(n+2\right)}^2-2}$	$-\frac{{\left(n+1\right)}^2-2}{{\left(n+2\right)}^2-2}$

Table 2. Table 2: Order of convergence for the weakly regular problems ( 1 )-( 98 ).

$\gamma =0.25,p=5$
$N$	$\delta {\lambda }_{5,N}$	order	$\delta {\lambda }_{10,N}$	order	$\delta {\lambda }_{20,N}$	order
$49$	$9.9201E-08$	$5.006$	$1.1937E-07$	$5.002$	$1.2280E-07$	$4.981$
$99$	$3.0866E-09$	$5.001$	$3.7263E-09$	$4.997$	$3.8890E-09$	$4.988$
$199$	$9.6399E-11$	$5.000$	$1.1670E-10$	$5.121$	$1.2255E-10$	$4.865$
$399$	$3.0127E-12$	–	$3.3538E-12$	–	$4.2064E-12$	–
$\gamma =0.5,p=4$
$N$	$\delta {\lambda }_{5,N}$	order	$\delta {\lambda }_{10,N}$	order	$\delta {\lambda }_{20,N}$	order
$49$	$2.1098E-05$	$4.003$	$3.0250E-05$	$3.999$	$3.2895E-05$	$3.981$
$99$	$1.3159E-06$	$4.001$	$1.8917E-06$	$4.000$	$2.0828E-06$	$3.999$
$199$	$8.2192E-08$	$4.000$	$1.1819E-07$	$4.000$	$1.3031E-07$	$4.000$
$399$	$5.1362E-09$	–	$7.3861E-09$	–	$8.1446E-09$	–
$\gamma =0.75,p=3$
$N$	$\delta {\lambda }_{5,N}$	order	$\delta {\lambda }_{10,N}$	order	$\delta {\lambda }_{20,N}$	order
$49$	$1.9714E-03$	$2.999$	$5.1330E-03$	$2.996$	$7.5944E-03$	$2.981$
$99$	$2.4665E-04$	$3.000$	$6.4360E-04$	$3.000$	$9.6156E-04$	$2.998$
$199$	$3.0833E-05$	$3.000$	$8.0475E-05$	$3.000$	$1.2036E-04$	$3.000$
$399$	$3.8541E-06$	–	$1.0060E-05$	–	$1.5048E-05$	–

Table 3. Table 3: Numerical approximations of λ 15 subscript 𝜆 15 \lambda_{15} for problems ( 1 )-( 100 ).

$\gamma$	MATSLISE2 [16]	SLEDGE [19]	SLEIGN2 [3]	${\mu }_{15}^{\left(3000\right)}$
$0.40$	$523.9182398826$	$523.9182763992$	$523.9182711601$	$523.9182763990$
$0.65$	$528.1784764701$	$528.1830115554$	$528.1791421864$	$528.1830147149$
$0.90$	$551.0460761225$	$551.9133020569$	$550.5488188848$	$552.2447514722$

Table 4. Table 4: Order of convergence for problems ( 1 )-( 104 ).

$\alpha =1/8,p=\sqrt{1+4\alpha }=2\sqrt{1.5}\approx 2.450$
$N$	$\delta {\lambda }_{5,N}$	order	$\delta {\lambda }_{10,N}$	order	$\delta {\lambda }_{20,N}$	order
$49$	$1.4443E-04$	$2.448$	$6.2160E-04$	$2.445$	$2.8090E-03$	$2.431$
$99$	$2.6461E-05$	$2.449$	$1.1412E-04$	$2.449$	$5.2076E-04$	$2.448$
$199$	$4.8448E-06$	$2.449$	$2.0900E-05$	$2.449$	$9.5467E-05$	$2.449$
$399$	$8.8697E-07$	–	$3.8264E-06$	–	$1.7481E-05$	–
$\alpha =1/2,p=\sqrt{1+4\alpha }=2\sqrt{3}\approx 3.464$
$N$	$\delta {\lambda }_{5,N}$	order	$\delta {\lambda }_{10,N}$	order	$\delta {\lambda }_{20,N}$	order
$49$	$8.4050E-05$	$3.463$	$4.0854E-04$	$3.459$	$2.4019E-03$	$3.438$
$99$	$7.6244E-06$	$3.464$	$3.7163E-05$	$3.463$	$2.2161E-04$	$3.462$
$199$	$6.9096E-07$	$3.464$	$3.3691E-06$	$3.464$	$2.0117E-05$	$3.464$
$399$	$6.2613E-08$	–	$3.0530E-07$	–	$1.8233E-06$	–
$\alpha =1,p=2\sqrt{1+4\alpha }=2\sqrt{5}\approx 4.472$
$N$	$\delta {\lambda }_{5,N}$	order	$\delta {\lambda }_{10,N}$	order	$\delta {\lambda }_{20,N}$	order
$49$	$8.6382E-06$	$4.470$	$4.5493E-05$	$4.465$	$3.2872E-04$	$4.425$
$99$	$3.8980E-07$	$4.472$	$2.0601E-06$	$4.471$	$1.5299E-05$	$4.469$
$199$	$1.7566E-08$	$4.472$	$9.2871E-08$	$4.472$	$6.9082E-07$	$4.472$
$399$	$7.9135E-10$	–	$4.1846E-09$	–	$3.1136E-08$	–

Equations

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

$${\alpha }_{a}y(a)+{\beta }_a{y}^{\prime }(a)$$

$${\alpha }_{b}y(b)+{\beta }_b{y}^{\prime }(b)$$

$$(a,b)\equiv (-1,1),$$

$$q(x)=f(x)+\frac{g(x)}{(1+x{)}^{\gamma }},$$

$$\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}$$

$$q_{mn}$$

$$⟨v,u^{\prime \prime }⟩-⟨u,v^{\prime \prime }⟩={\left[u^{\prime }(x)v(x)-u(x)v^{\prime }(x)\right]}_{-1}^1,$$

$$⟨{\mathcal{R}}_m,{\mathcal{R}}_n^{\prime \prime }⟩=⟨{\mathcal{R}}_n,{\mathcal{R}}_m^{\prime \prime }⟩+{\left[{\mathcal{R}}_n^{\prime }(x){\mathcal{R}}_m(x)-{\mathcal{R}}_n(x){\mathcal{R}}_m^{\prime }(x)\right]}_{-1}^1=⟨{\mathcal{R}}_n,{\mathcal{R}}_m^{\prime \prime }⟩.$$

$${\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.$$

$$V_{n}=\left(\begin{array}[]{r|r|r}\alpha_{a}-\frac{n(n+1)}{2}\beta_{a}&-\alpha_{a}+\frac{(n+1)(n+2)}{2}\beta_{a}&\alpha_{a}-\frac{(n+2)(n+3)}{2}\beta_{a}\\ &&\vspace{-2mm}\hfil\\ \hline\cr&&\vspace{-2mm}\hfil\\ \alpha_{b}+\frac{n(n+1)}{2}\beta_{b}&\alpha_{b}+\frac{(n+1)(n+2)}{2}\beta_{b}&\alpha_{b}+\frac{(n+2)(n+3)}{2}\beta_{b}\end{array}\right).$$

$$\displaystyle\left\{\begin{array}[]{r}\left(\alpha_{a}-\frac{n(n+1)}{2}\beta_{a}\right)\xi_{n}+\left(\alpha_{a}-\frac{(n+2)(n+3)}{2}\beta_{a}\right)\theta_{n}=0,\\ \eta_{n}=0,\end{array}\right.$$

$${\xi }_n$$

$${\eta }_n$$

$${\theta }_n$$

$$\text{mbox}det\left(V_n(:,[13])\right)=(2n+3)({\alpha }_a{\beta }_b+{\alpha }_b{\beta }_a)\ne 0,\text{mbox}foreachn.$$

$${\xi }_n$$

$${\eta }_n$$

$${\theta }_n$$

$$\mbox{for each $n$ let $\nu_{n}$ be such that }\|\left(\xi_{n},\eta_{n},\theta_{n}\right)^{T}\|_{\infty}=1\mbox{ and }\xi_{n}\geq 0.$$

$${\xi }_n$$

$${\theta }_n$$

$${\eta }_n$$

$$a_{nn}$$

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

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

Full text

1 Introduction

The direct Sturm-Liouville problem (SLP) in normal form with separated boundary conditions is given by

[TABLE]

where the potential $q,$ the domain $(a,b)$ and the coefficients $\alpha_{a},\beta_{a},\alpha_{b},\beta_{b}$ represent the data of the problem while the unknowns are the eigenvalues $\lambda$ and the corresponding eigenfunctions $y.$ It is surely a classical problem that has been extensively studied both from the theoretical and from the numerical point of views. Many numerical schemes are nowadays available for its solution which can be subdivided into two main families: matrix methods and shooting techniques [20]. A number of well-established numerical codes that are able to solve regular problems as well as many singular ones had been developed and are freely available for the scientific community. Among them we surely mention the MATSLISE2 [16], the SLEDGE [19] and the SLEIGN2 [3] codes, but there are many others. In spite of this, we think that the approach of deriving matrix schemes by using spectral methods, instead of finite difference or element ones, deserves further insights and this is the topic of the present paper. In particular, we will consider the case of a bounded domain which, without loss of generality, we assume to be

[TABLE]

and study problems with

[TABLE]

where $f$ and $g$ are analytic functions inside and on a Bernstein ellipse containing $[-1,1],$ and ${\gamma}\geq 0.$ Problems of this type have many applications in physics (cf., for example, [5, 6, 24, 25]). We recall that if $q\in L_{1}(-1,1)$ then the problem is regular (sometimes called weakly regular if ${\gamma}\in(0,1)$ and $g(-1)\neq 0$) otherwise it is singular. More precisely, the singular left endpoint is of type [2, 9, 15, 23]

Limit-Circle (LC) and nonoscillatory if ${\gamma}\in[1,2)$ and $g(-1)\neq 0$ or ${\gamma}=2$ and $g(-1)\in[-1/4,3/4)\setminus\{0\};$ 2. 2.

Limit-Circle and oscillatory if ${\gamma}=2$ and $g(-1)<-1/4$ or ${\gamma}>2$ and $g(-1)<0;$ 3. 3.

Limit-Point (LP) and nonoscillatory if ${\gamma}=2$ and $g(-1)\geq 3/4$ or ${\gamma}>2$ and $g(-1)>0.$

In the sequel, we will always exclude the case of an oscillatory endpoint which is definitely much more difficult to be treated numerically and we will assume ${\gamma}\in[0,2]$ leaving the generalization to larger values of ${\gamma}$ to future investigation. Concerning the boundary condition at $x=-1,$ which is required in the LC case, we will consider the Friedrichs one, namely in this context the Dirichlet condition. This means that for any $\lambda$ we select the principal solution of the equation which is frequently the most significant in the applications, [18, 20]. We recall that the Dirichlet condition is the only possible one in the LP case.

Under these assumptions, it is known that the eigenvalues of (1)–(5) are real, simple and that they 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)–(5).

Talking about numerical methods based on shooting techniques, the approach used most frequently for solving this type of problems is based on the selection of a layer, namely (1) is usually solved over $(-1+\epsilon,1).$ In particular, suitable algorithms have been studied for an adaptive selection of $\epsilon$ and for the computation of the condition to be imposed at $x=-1+\epsilon$ (see, for example, [3, 17] for further details).

Concerning classical matrix methods, it is clear that they can be employed if based on a discretization of the differential equation and of the boundary conditions which do not require the evaluation of $q,$ or even of its derivative, at the left endpoint. This is the case, for example, of the classical three-point formula. The main difficulty with this approach may be that of an order reduction caused by the fact that the derivative, of suitable order, of the eigenfunctions may be unbounded as $x$ approaches $-1.$

In this paper, we shall derive a matrix method by using a spectral Galerkin approach based on Legendre polynomials. Before proceeding, it must be said that if the problem is subject to Dirichlet boundary conditions at both endpoints then schemes based on spectral (collocation) methods which use orthogonal polynomials or sinc functions are already available in the literature (see, for example, [7, 10, 12, 14]). These are surely effective methods to be employed whenever the eigenfunctions are sufficiently smooth. Our main purpose is therefore that of treating general boundary conditions and to get accurate approximations of the eigenvalues even in the case where the eigenfunctions are not so regular.

The remaining part of this paper is organized as follows. In Section 2, we describe the approach, derive the generalized eigenvalue problem which discretize the continuous one and discuss how the entries of the matrices involved can be computed efficiently. Section 3 is devoted to the analysis of the errors in the resulting eigenvalue approximations for problems with unbounded potential at the left endpoint. Moreover, in the same section we shall derive low cost and effective procedures for the computation of corrected numerical eigenvalues. Finally, the results of several numerical experiments are reported and discussed in Section 4.

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

[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 (9) 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. The same property holds for $A_{N}$ thanks to the well-known Green’s identity

[TABLE]

by which one gets, for each ${\cal{R}}_{n}$ and ${\cal{R}}_{m}\in{\cal{S}}_{N},$

[TABLE]

Clearly the effectiveness of the procedure is strictly connected to the choice of the basis functions. The main criterion we have considered is the computational cost of the method which is essentially determined by the calculation of the coefficient matrices and by the solution of (10). This suggests to use suitable combinations of the classical Legendre polynomials as described in the next subsection.

2.1 Basis functions

As done in [21], we look for a basis function ${\cal{R}}_{n}$ of the following form

[TABLE]

where ${\cal{P}}_{j}$ is the Legendre polynomial of degree $j$ for which it is known that [1]

[TABLE]

Therefore, with some computations one gets that ${\cal{R}}_{n}\in{\cal S}_{N},$ see (6)-(7), if and only if $\left(\xi_{n},\eta_{n},\theta_{n}\right)^{T}$ belongs to the kernel of

[TABLE]

We must now distinguish the following two possibilities:

$\alpha_{a}\beta_{b}+\alpha_{b}\beta_{a}=0,$ i.e. problems subject to symmetric BCs. In this case it is natural to set $\eta_{n}=0$ so that ${\cal{R}}_{n}$ is an even or an odd function, depending on $n,$ which implies that if (2) holds true then (3) is verified automatically. In this way, one obtains the following system of equations

[TABLE]

whose general solution can be written as

[TABLE]

where $\nu_{n}\neq 0$ is a free parameter; 2. 2.

$\alpha_{a}\beta_{b}+\alpha_{b}\beta_{a}\neq 0.$ From the previous considerations, one deduces that $\eta_{n}$ must be different from zero. Using the Matlab notation, this is confirmed by the fact that

[TABLE]

Hence, if we let as before $\nu_{n}\neq 0$ be a free parameter then we got

[TABLE]

Concerning the parameters $\left\{\nu_{n}\right\}_{n\in{\mathbb{N}}_{0}}$ we decided to establish a criterion for their selection in order to obtain basis functions independent of the scaling of $(\alpha_{a},\beta_{a})$ and/or $(\alpha_{b},\beta_{b}).$ Now, a natural choice would have been that of choosing $\nu_{n}$ so that $\|{\cal{R}}_{n}\|_{2}=1$ for each $n$ since we are working in $L_{2}.$ Nevertheless, we preferred not to proceed in this way to avoid the computation of square roots at least at this level. As an alternative to $\|{\cal{R}}_{n}\|_{2}$ we considered its uniform norm which is not known in closed form but it satisfies $\|{\cal{R}}_{n}\|_{\infty}\leq 3\|\left(\xi_{n},\eta_{n},\theta_{n}\right)^{T}\|_{\infty}$. We thus applied the following criterion

[TABLE]

The so-obtained coefficients for the four problems subject to natural BCs and for two general ones are listed in Table 1. The unspecified values for the last BCs, of Dirichlet-Robin type, are $\xi_{0}=2/3,\,\eta_{0}=1$ and $\theta_{0}=1/3.$

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

[TABLE]

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

In this section we are going to show that the entries of $A_{N}$ and $B_{N}$ in (11)-(12) can be determined analytically thanks predominantly to the orthogonality of the Legendre polynomials with respect to the standard inner product.

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 (15). Consequently, $A_{N}=A_{N}^{T}$ is diagonal with diagonal entries

[TABLE]

where for the last two equalities we used (14) and (16). For later convenience, 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.3 The matrix $Q_{N}.$

Let us consider first of all the case of a regular problem with $\gamma<1.$ From (13), 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 (32) and

[TABLE]

We recall that the Legendre polynomials obey the recurrence relation

[TABLE]

This allows to prove the following result.

Proposition 2.1

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

[TABLE]

with $\hat{{\mathbf{q}}}_{-1}$ the zero sequence. If we define the linear tridiagonal operator ${\mathbf{z}}\in\ell_{\infty}\mapsto\mathcal{H}\,{\mathbf{z}}\in\ell_{\infty}$ where

[TABLE]

then we get

[TABLE]

Proof: see [13, Propositions 1,2] with $\alpha=0.$ $\Box$

The immediate consequence of this proposition is that the recurrence in (46) permits to determine the entire matrix $\hat{Q}_{N}$ once $\hat{q}_{m0}=\hat{f}_{m0}+\hat{g}_{m0}$ have been computed for each $m=0,1,\ldots,2N+2.$ In fact, these values are sufficient to determine $\hat{q}_{m1}$ for each $m=0,1,\ldots,2N+1$ by using (46) with $n=0$ and the fact that $\mathcal{H}$ is tridiagonal. At this point, the application of (46) with $n=1$ allows to compute $\hat{q}_{m2}$ for each $m=0,1,\ldots,2N,$ and so forth. Clearly, in the actual implementation, the symmetry of $\hat{Q}_{N}$ is taken into account.

In particular, the coefficients $\hat{f}_{m0}$ in (36) decay exponentially as $m$ increases due to the assumption that $f$ is analytical inside and over a Bernstein ellipses containing $[-1,1].$ In addition, we can use the routine legcoeffs of the well-established open-source software package Chebfun [8] to determine the numerically significant values. Thus if $L+1$ is the length of the vector of Legendre coefficients of $f$ provided by such routine, we approximate upto machine precision the symmetric matrix $\hat{F}_{N}$ with its banded portion with bandwidth $2L+1.$

Remark 2.1

It is important to stress that if $q$ is analytical, i.e. if $g\equiv 0,$ then the generalized eigenvalue problem (10) which discretize the SLP involves only sparse matrices. In particular, $B_{N}$ is symmetric positive definite and pentadiagonal while $A_{N}+Q_{N}=A_{N}+R_{N}^{T}\hat{F}_{N}R_{N}$ is symmetric with bandwith $2L+5.$

Concerning the computation of the required entries of $\hat{{\mathbf{g}}}_{0}=\left(\hat{g}_{10},\hat{g}_{20},\ldots\right)^{T},$ see (37), we applied arguments similar to the ones used in [13, Propositions 2,3]. In detail, recalling that by assumption $g$ is analytical inside and over a Bernstein ellipse containing $[-1,1]$ too, the operator $g(\mathcal{H})=\sum_{\ell=0}^{+\infty}\frac{\langle g,{\cal{P}}_{\ell}\rangle}{\langle{\cal{P}}_{\ell},{\cal{P}}_{\ell}\rangle}\,{\cal{P}}_{\ell}\left(\mathcal{H}\right)$ is well defined. Consequently, [13, Proposition 2] and [11, 16.4 formula (2)] allow to get that if $\gamma<1$ then

[TABLE]

where $(t)_{\ell}$ is the Pochhammer symbol. We observe that $\hat{g}_{m}$ coincides with $\hat{g}_{m0}$ if $g(x)\equiv 1$. In the actual implementation, we proceed as follows: we get a polynomial approximation of $g$ by transforming it in a Chebfun function, which is accurate up to machine precision, then we apply the previous formula to compute the first $2N+2$ entries of $\hat{{\mathbf{g}}}_{0}.$

Let us now discuss the case of singular problems, namely how we determine $Q_{N}$ if $\gamma\in[1,2].$ We recall that the corresponding Friedrichs boundary condition at the singular endpoint is $y(-1)=0.$ The basis functions have therefore a root at $x=-1$ and we need to highlight this fact. This is done in the following proposition.

Proposition 2.2

If $\beta_{a}=0$ and ${\cal{P}}_{\ell}^{(0,1)}$ is the Jacobi polynomial of degree $\ell$ with weighting function $\omega(x)=(1+x),$ then ${\cal{R}}_{n}$ in (15) can be written as

[TABLE]

Proof: The first equality is evident with ${\cal{U}}_{n}\in\Pi_{n+1}$ since ${\cal{R}}_{n}\in\Pi_{n+2}.$ In addition

[TABLE]

This implies that there exist suitable $\tilde{\xi}_{n}$ and $\tilde{\theta}_{n}$ such that

[TABLE]

It remains to prove that $\tilde{\theta}_{n}=\theta_{n}$ and $\tilde{\xi}_{n}=\xi_{n}.$ The former equality is an immediate consequence of the fact that ${\cal{P}}_{n+2}$ and ${\cal{P}}_{n+1}^{(0,1)}$ have the same leading coefficient, [1]. Concerning the latter equality, it is then trivial if the problem is subject to Dirichlet-Dirichlet BCs (see Table 1 and recall that ${\cal{P}}_{\ell}^{(0,1)}(1)=1$ for each $\ell$). On the other hand, if $\beta_{b}\neq 0$ then from (21) and (23) with $(\alpha_{a},\beta_{a})=(1,0)$ we get

[TABLE]

In addition, it is known that $\frac{d}{dx}{\cal{P}}_{n}^{(0,1)}(1)=n(n+2)/2.$ With this information, one verifies that the polynomial at the right hand-side of (48) satisfies the BC at $x=1$ and this completes the proof. $\Box$

The matrix $Q_{N}$ can be therefore written as

[TABLE]

see (32) and (35), where

[TABLE]

Now with an approach similar to the one considered in Proposition 2.1 and in the subsequent paragraph, which is essentially based on the recurrence relation for the Jacobi polynomials ${\cal{P}}_{n}^{(0,1)}$ [1, 11], we obtain that if we know the values of $\tilde{g}_{m0}$ for each $m=0,1,\ldots,2N+1$ then we can determine the remaining required values recursively. Moreover, if we let

[TABLE]

then we get, see [11, 16.4 formula (2)],

[TABLE]

where

[TABLE]

Remark 2.2

If ${\gamma}=1$ then $\tilde{g}_{m}=0$ for each $m\geq 1.$ Therefore, in this case, $Q_{N}$ can be approximated up to machine precision with a banded matrix where the bandwidth depends on the number of numerically significant coefficients of the Legendre-Fourier series expansion of $g$ and $f.$ In particular, $Q_{N}$ is tridiagonal for the Boyd equation for which $g(x)\equiv g(-1)$ and $f\equiv 0.$

3 Error analysis and computation of corrected numerical eigenvalues.

In this section, we shall study the behavior of the error in the resulting numerical eigenvalues as $N$ increases and for a fixed index.

As usual, this is related to the regularity of the solution, namely, in this context, to the regularity of the eigenfunctions. In particular, if $\gamma=0$ or $g(x)\equiv 0,$ then $q,$ and consequently $y,$ belongs to $C^{\infty}[-1,1].$ In this case, it is well-known that the errors in the approximations provided by a spectral method decay exponentially.

Problems which require a deeper analysis are therefore those for which $q$ is unbounded at the left endpoint. We must observe that from (34)–(37) and (49)–(56) one deduces that the spectral method we have derived is well defined for each ${\gamma}\in(0,3),$ with $y(-1)=0$ if the left endpoint is singular. Nevertheless, we shall consider only the case ${\gamma}\in(0,2]$ with $g(-1)\neq 0,$ namely only problems for which $x=-1$ is a regular singular endpoint. The generalization to essential singularities will be the topic of future research.

In this context, the results we are going to present are not only interesting from the theoretical point of view but they will also provide very simple, economical and effective techniques for the computation of corrected 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]

By construction of $z_{N}$ and of $\lambda^{(N)},$ see (8)-(9), we can write

[TABLE]

On the other hand, it is evident that

[TABLE]

In addition

[TABLE]

since $z_{N}^{\prime\prime}\in\Pi_{N-1}$ and $y-y_{N}\in\Pi_{N-1}^{\bot}.$ From these formulas we get

[TABLE]

being $\varepsilon_{N}=\langle z_{N},y-y_{N}\rangle/\langle z_{N},y\rangle.$ Therefore, an analysis of the behavior of the coefficients $c_{n}$ in (63) as $n$ increases is required and the following result constitutes a first step.

Proposition 3.1

If $n$ is sufficiently larger than the index of the eigenvalue, ${\gamma}\in(0,2]$ and $g(-1)\neq 0$ then

[TABLE]

Proof: It is evident that

[TABLE]

We recall that we assumed $f$ in (5) to be analytical inside and over a Bernstein ellipse containing $[-1,1]$ and this implies that its Legendre coefficients decay exponentially. Thus, $\langle{\cal{R}}_{n},(\lambda-f)y\rangle$ becomes negligible with respect to $\langle{\cal{R}}_{n},(1+x)^{-{\gamma}}gy\rangle$ as $n$ increases and this complete the proof. $\Box$

From (64), by using similar arguments, one deduces that the main contribution to the error in the eigenvalue approximation is given by

[TABLE]

We recall the following asymptotic estimate [22].

Proposition 3.2

Let $\psi\in C^{\infty}(-1,1)\bigcap C[-1,1]$ have the expansion

[TABLE]

with $\psi(-1)\neq 0,\psi_{0}=1$ and $\sigma_{0}=0<\sigma_{1}<\ldots<\sigma_{L}<1.$ If $s\in(-1,+\infty)\setminus{\mathbb{N}}_{0}$ and if $n$ is sufficiently large then

[TABLE]

We can now prove the following result which concerns (weakly) regular problems not subject to the Dirichlet boundary condition at the left endpoint.

Theorem 3.1

If $y(-1)g(-1)\neq 0,$ ${\gamma}\in(0,1)$ and if $N$ is sufficiently larger than the index of the eigenvalue then

[TABLE]

where

[TABLE]

Proof: From (15), (67) and the previous Proposition with $s=-\gamma,$ $\psi=g\,z_{N}$ or $\psi=gy,$ and $L=0$ we obtain

[TABLE]

Now, we can rewrite the last equation in (47) as follows

[TABLE]

Therefore, by using the following expansion of the ratio of two gamma functions

[TABLE]

with $z=n+1/2,$ we get

[TABLE]

and, consequently,

[TABLE]

We recall that if $n$ is sufficiently large then $\xi_{n}=1,$ see (24). In addition, by using (18)–(29) it is possible to verify with some computations that if $\beta_{a}\neq 0$ then

[TABLE]

Hence, see (31) and (71),

[TABLE]

Therefore

[TABLE]

which is the statement of the theorem. $\Box$

This result immediately suggests a very simple formula for the correction of the numerical eigenvalue. First of all, we assume the numerical and the exact eigenfunctions have been normalized so that

[TABLE]

By using the orthogonality of the Legendre polynomials, this permits the estimates $y(-1)\approx z_{N}(-1),$ $\langle z_{N},y\rangle\approx 1,$ and consequently, see (70),

[TABLE]

In addition, we observe that the term $\lambda^{(N)}\varepsilon_{N}$ in (64) can be of some relevance if $N$ is not so much larger than the index of the eigenvalue. By virtue of (65), (69)-(69) and of (73), we therefore decided to consider the following approximation of $\varepsilon_{N}$ which can be computed with a very low cost

[TABLE]

where

[TABLE]

It is worth recalling that $\langle z_{N},{\cal{R}}_{N}\rangle=b_{N,N-2}\zeta_{N-2,N}+b_{N,N-1}\zeta_{N-1,N},$ and $\langle z_{N},{\cal{R}}_{N+1}\rangle=b_{N+1,N-1}\zeta_{N-1,N}$ (see (10),(12) and (33)).

All these arguments lead to the following formula for the correction of the numerical eigenvalue to be used when the BC at the left endpoint is not of Dirichlet type

[TABLE]

The main steps of the procedure for the computation of such $\lambda^{(N)}$ and $\mu^{(N)}$ are summarized in Algorithm 1. With respect to the notation we have used so far, we add the further index $k$ which represents the index of the eigenvalue. Its value belongs to $\left\{1,\ldots,M\right\}$ being $M$ the number of smallest eigenvalues requested.

Let us now consider problems subject to the Dirichlet boundary condition at the left endpoint with ${\gamma}\in(0,2).$ In this case, Proposition 3.2 with $s=-{\gamma}$ is not directly applicable to the inner products in (67). In fact, we need $s>-1$ and $\psi(-1)\neq 0.$ Nevertheless, it is evident that since $z_{N}(-1)=y(-1)=0$ we can write

[TABLE]

with $\hat{z}_{N}(-1)=z_{N}^{\prime}(-1)\neq 0$ and $\hat{y}(-1)=y^{\prime}(-1)\neq 0.$ Consequently

[TABLE]

Now, an analysis of the behavior of $\hat{y}$ in proximity of $x=-1$ is required for the application of (69)-(69) and a Frobenius-type method provides the following result.

Lemma 3.1

If $y(-1)=0,$ $g(-1)\neq 0$ and ${\gamma}\in(0,2)$ then an eigenfunction $y$ admits the following expansion as $x\rightarrow-1^{+}:$

[TABLE]

Here

[TABLE]

i.e.

[TABLE]

Remark 3.1

It must be observed that the term with the summation in (78) represents the truncation of a fractional power series expansion at $x=-1$ of the solution of (1) with

[TABLE]

In fact, if one sets $t=(1+x)^{1-{\gamma}/2}$ and $y(x)=u(t)$ then one gets

[TABLE]

A solution of this equation subject to $u(0)=0$ is proportional to [4]

[TABLE]

Here $J_{\nu}$ is the Bessel function of the first kind and $\,{}_{0}F_{1}$ a confluent hypergeometric limit function. Therefore, the solution of the initial value problem (1)-(81) is

[TABLE]

Remark 3.2

For the special value ${\gamma}=1,$ like the Boyd equation, it is possible to verify that a solution of (1) subject to $y(-1)=0$ admits a (classical) power series expansion at $x=-1$. Moreover, the coefficients of the expansion of $y$ in (63) decays exponentially and, see Remark 2.2, the matrix $Q_{N}$ can be approximated up to machine precision with a banded one with bandwidth independent of $N.$ From these observations, we deduce that if ${\gamma}=1$ then the errors in the approximation of the eigenvalues decay exponentially with respect to $N.$

Problems to be studied are therefore those with ${\gamma}\in(0,2)\setminus\left\{1\right\}$ and this is done in the next theorem.

Theorem 3.2

If $y(-1)=0,$ $g(-1)\neq 0$ and ${\gamma}\in(0,2)\setminus\left\{1\right\}$ then

[TABLE]

where $L$ is defined in (79) and, see (80),

[TABLE]

Proof: From the premises of this theorem, Propositions 2.2 and 3.2, (48), (57), (72) with $z=n+1$, (82), and (24)-(29) we obtain that if $n\geq N$ with $N$ sufficiently large then

[TABLE]

Similarly, by using also (83) and the previous lemma, we get

[TABLE]

with

[TABLE]

Hence, by considering (31) and (65) we obtain

[TABLE]

so that, see (64),

[TABLE]

The statement follows by using an integral estimate. $\Box$

Let us now discuss how one can use this result for the correction of the numerical eigenvalues. Following the idea used for problems with $y(-1)\neq 0,$ we consider the normalization specified in (74), with $z^{\prime}_{N}(-1)>0$ and $y^{\prime}(-1)>0,$ by which we get $\langle z_{N},y\rangle\approx 1.$ On the other hand, the estimate $y^{\prime}(-1)\approx z_{N}^{\prime}(-1)$ turns out to be rather poor. In fact, we have just established that if $n$ is sufficiently large then $c_{n}=O\left((n+3/2)^{-p/2-1}\right)$ and it is possible to verify that

[TABLE]

Therefore, $c_{n}{\cal{R}}_{n}^{\prime}(-1)=O\left((n+3/2)^{-p/2}\right)$ approaches zero rather slowly when $p\rightarrow 2^{+},$ i.e. ${\gamma}\rightarrow 2^{-}$. For example, if ${\gamma}=1.9$ then $p/2=6/5.$ By considering (8) and (63), the following approximation

[TABLE]

turns out to be more appropriate. Now, from (84) we obtain

[TABLE]

and, consequently, $y^{\prime}(-1)\approx z_{N}^{\prime}(-1)-d_{N}y^{\prime}(-1),$ i.e. $y^{\prime}(-1)\approx z_{N}^{\prime}(-1)/(1+d_{N})$.

Summarizing, if the problem is subject to the Dirichlet BC at the singular endpoint and if ${\gamma}\in(0,2)\setminus\left\{1\right\}$ then we correct the numerical eigenvalues as follows

[TABLE]

where $\bar{\varepsilon}_{N}$ is defined in (75) with $\bar{c}_{N}$ and $\bar{c}_{N+1}$ given by (84) with $n=N,N+1,$ respectively. The main steps of the procedure for their computation are listed in Algorithm 2.

The final error analysis we are going to present concerns problems with ${\gamma}=2$ and $g(-1)>0.$ In this case, the application of the Frobenius method allows to state that the exact eigenfunction satisfies

[TABLE]

where $\chi$ is a free parameter while $\varrho$ is the positive root of the indicial equation $\varrho^{2}-\varrho-g(-1),$ i.e.

[TABLE]

For instance, if $g(x)\equiv g(-1)$ and $f(x)\equiv 0$ then a solution of (1) subject to $y(-1)=0$ is proportional to [4]

[TABLE]

We shall proceed by assuming the exact eigenfunction has been normalized so that

[TABLE]

Concerning the corresponding numerical eigenfunction, we assume

[TABLE]

being

[TABLE]

In particular, the second formula in (90) state that $\hat{z}_{N}(-1)$ and $\hat{y}_{N}(-1)$ are infinitesimal of the same order as $N$ increases.

With these notations, we can prove the following result.

Theorem 3.3

If ${\gamma}=2,$ $y(-1)=0,$ $\varrho\notin{\mathbb{N}},$ see (89), and if $N$ is sufficiently larger then the index of the eigenvalue then

[TABLE]

where $\kappa$ is the limit value in (90) and

[TABLE]

Proof: We begin by studying the asymptotic behavior of the coefficients $c_{n}$ with an approach similar to the one used in the proof of the previous theorems. From (65), (69)-(69) and (88) we get

[TABLE]

In particular, for the last estimate we used the fact that $g(-1)=\varrho(\varrho-1)$ and (31).

Let us now consider $\langle{\cal{R}}_{n},(1+x)^{-2}gz_{N}\rangle=\langle{\cal{U}}_{n},g\hat{z}_{N}\rangle.$ If we let $g(x)=g(-1)+(1+x)\tilde{g}(x),$ with $\tilde{g}\in C^{\infty}([-1,1]),$ then we obtain

[TABLE]

with a remainder that decreases exponentially. In addition, if one uses [11, 16.4, formula(20)] then one gets $\langle{\cal{P}}_{m}^{(0,1)},{\cal{P}}_{n}^{(0,1)}\rangle=P_{m}^{(0,1)}(-1)\int_{-1}^{1}{\cal{P}}_{n}^{(0,1)}(x)dx=(-1)^{n}\,2/(n+1)$ for each $m\leq n.$ Clearly, this implies $\langle{\cal{U}}_{m},{\cal{U}}_{n}\rangle={\cal{U}}_{m}(-1)\int_{-1}^{1}{\cal{U}}_{n}(x)dx,$ for each $m<n$ and consequently

[TABLE]

Therefore, the first estimate (92) follows from (66), with $\langle z_{N},y\rangle\approx 1,$ and from the application of an integral estimate. It must be said that this result is only a starting point and it is not particularly useful both from the theoretical point of view and for the derivation of a correction formula. This is because $\hat{z}_{N}(-1)$ approaches zero as $N$ increases and if we don’t know its infinitesimal order then we don’t know the order of convergence of $\lambda-\lambda^{(N)}.$ Concerning the computation of a corrected numerical eigenvalue we need an estimate of $\chi$. Let us therefore discuss the approximations in (93). By using the formula obtained for $c_{n}$ and by observing that ${\cal{U}}_{n}(-1)=(-1)^{n}(2n+3)+O(1),$ one gets

[TABLE]

where $\kappa$ is defined in (90). This implies $\hat{z}_{N}(-1)=O((N+1)^{2-2\varrho})$ and consequently

[TABLE]

Moreover

[TABLE]

which, with a simple substitution, completes the proof of (93). $\Box$

We must now spend few words concerning the limit value $\kappa.$ It is evident that intuitively one would expect $\kappa=1,$ namely $\hat{z}_{N}(-1)\approx\hat{y}_{N}(-1).$ Nevertheless, the results of several numerical experiments we have conducted by considering different $g,$ $f$ and different boundary conditions at $x=1$ indicate that this is not the case. In particular, such tests lead us to the following assumption

[TABLE]

Currently, this is an experimental result but we didn’t a find a problem for which it doesn’t work. By using it, we get that we can correct the numerical eigenvalues with the following very simple formula

[TABLE]

where, once again, $\bar{\varepsilon}_{N}$ is defined in (75) with

[TABLE]

The procedure for their computation is sketched in Algorithm 3.

4 Numerical tests

The method described and the algorithms for the a posteriori correction were implemented in Matlab (ver.R2017a). In particular, we solved the arising generalized eigenvalue problem (10) by using the eigs function with option “SM” for getting the ones of smallest magnitude. In addition, as we already said in Section 2.3, routines included in the open-source Chebfun package [8] were conveniently used to determine the Fourier-Legendre coefficients of the functions $f$ and $g$ required for the computation of the coefficient matrix $Q_{N}.$

For each of the three types of problems we have studied in the previous section, which give rise to the three algorithms we have sketched, we now present the results obtained for two different $g,$ $f$ and boundary conditions. In several tests, we needed an accurate estimate of the exact eigenvalues for the evaluation of the errors in the uncorrected and/or corrected numerical ones. In this regard, we decided to consider as “exact” those provided by a well-established routine or, alternatively, the corrected ones obtained with $N$ very large. Further details will be given in the sequel. Before proceeding, we must say that when talking about relative errors we usually refer to

[TABLE]

with $\bar{\lambda}_{k}$ the $k$th reference eigenvalue used.

Let us begin with (weakly) regular problems not subject to the Dirichlet boundary condition at $x=-1$ (see Algorithm 1). The first results we present confirm that the error in the uncorrected numerical eigenvalues behaves like $O((N+1)^{-p})$ where $p=6-4{\gamma},$ see Theorem 3.1. In particular, we considered the problems with the following potentials and BCs

[TABLE]

and we used the classical formula

[TABLE]

for the numerical estimate of the order of convergence. The results we got for the eigenvalues of index $k=5,10,20,$ listed in Table 2, surely confirm the statement of the theorem previously mentioned.

Concerning the effectiveness and utility of the application of the a posteriori correction, we applied Algorithm 1 for solving the problems with

[TABLE]

For the computation of corresponding reference eigenvalues, we first of all tried to use the well-established and general-purpose codes MATSLISE2 [16], SLEDGE [19], and SLEIGN2 [3] with a tolerance for the absolute and/or relative error equal to $10^{-13}.$ In particular, for the SLEIGN2 routine we set the input parameter ${\tt NCA}$ equal to two which indicate that the left endpoint is weakly regular. The approximations we obtained for $\lambda_{15}$ are listed in Table 3. As one can see, the number of significant digits for which such estimates agrees decreases as ${\gamma}$ approaches one. Indeed, this fact is underlined in the documentation of these three softwares and all our tests indicate that it is more relevant if $y(-1)\neq 0.$ We therefore decided to use as reference eigenvalues $\lambda_{k}\approx\bar{\lambda}_{k}\equiv\mu_{k}^{(N_{T})}$ with $N_{T}\gg N$ and, in particular, for the results shown in Figure 1, we set $N_{T}=3000$ (the values of $\mu_{15}^{(3000)}$ are listed in the last column of Table 3). In more details, in the three subplots at the top of Figure 1, the resulting relative errors in the approximation of the fifteenth eigenvalue are plotted versus $N$ with $N$ ranging from $40$ to $320.$ For the subplots at the bottom, instead, we fixed $N=80$ and we depict the errors for the index $k$ ranging from $1$ to $30.$ 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. As one can see, the correction improves noticeably the accuracy of the numerical eigenvalues. As a matter of fact, see the subplots at the top of Figure 1, it results always

[TABLE]

with $k=15.$ On the other hand, from the subplots at the bottom one deduces that for $N=80$ and $1\leq k\leq 30$ the gain resulting from the correction is at least of two significant digits for each eigenvalue but it is frequently much larger.

The next examples regards problems with ${\gamma}\in(1,2)$ subject to $y(-1)=0$ at the singular endpoint which is of type limit-circle. In this case, we applied Algorithm 2 for getting approximations of the eigenvalues and we used as reference “exact” ones, i.e. $\bar{\lambda}_{k},$ those provided by MATSLISE2 [16] with a tolerance for the absolute and/or relative errors equal to $10^{-13}.$ It must be said that this choice was motivated only by the fact that MATSLISE2 is a Matlab code and that the results we are going to present would have been essentially the same if we had decided to use one of the other two previously mentioned codes. With the same notation used for the second example, in Figure 2 we represent the errors for the problems with

[TABLE]

while in Figure 3 those corresponding to

[TABLE]

For both these examples, we observe that the spectral matrix method already provides accurate approximation for the smallest $\gamma$’s and the application of the correction further improves such estimates. In particular, see the subplots at the bottom left, the relative errors in the first fifty uncorrected numerical eigenvalues, determined with $N=128,$ are smaller than $10^{-10}$ while those in the corresponding corrected eigenvalues are smaller than $10^{-13},$ i.e. smaller than the tolerance used for the computation of the reference ones. This is the reason for which they are not depicted. Concerning the results obtained for ${\gamma}\in[1.5,2),$ the advantage arising from the application of the correction is undeniable since (101), with $k=25,$ holds almost always.

Finally, we consider problems with ${\gamma}=2$ and $g(-1)>0.$ In Table 4, we list the experimental orders of convergence, see (99), for the problems with

[TABLE]

The results obtained are in perfect agreement with the statement of Theorem 3.3. In fact, $g(x)=\alpha\cos(4\pi x)$ and $p=2\sqrt{1+4g(-1)}=2\sqrt{1+4\alpha},$ see (94).

Regarding Algorithm 3, we applied it for solving the problems with

[TABLE]

The corresponding errors with respect to the eigenvalue estimates provided by MATSLISE2 are represented in Figure 4. As one can see, the comments we have done for the previous examples concerning the effectiveness of the a posteriori correction surely apply even to this last one. In addition, in support of (91)-(95), in Figure 5 some graphs of $\hat{z}_{N}(x)$ and of $\hat{y}_{N}(x)$ (obtained as partial sum of $\hat{z}_{8000}(x)$) in proximity of $x=-1$ are reported. In the same figure, some ratios $\hat{z}_{N}(-1)/\hat{y}_{N}(-1)$ are also depicted which show that such values approach $(2\rho-1)/\rho^{2}$ as $N$ increases as stated in (95) (see also (89)).

Acknowledgments

The author is very indebted to Prof. Paolo Ghelardoni for helpful discussions and suggestions.