Approximation of Eigenvalues of Sturm-Liouville Problems by Using Hermite Interpolation

Abstract

Eigenvalue problems with eigenparameter appearing in the boundary conditions usually have complicated characteristic determinant where zeros cannot be explicitly computed. In this paper, we use the derivative sampling theorem “Hermite interpolations” to compute approximate values of the eigenvalues of Sturm-Liouville problems with eigenvalue parameter in one or two boundary conditions. We use recently derived estimates for the truncation and amplitude errors to compute error bounds. Also, using computable error bounds, we obtain eigenvalue enclosures. Also numerical examples, which are given at the end of the paper, give comparisons with the classical sinc method and explain that the Hermite interpolations method gives remarkably better results.

Introduction

The mathematical modeling of many practical problems in mechanics and other areas of mathematical physics requires solutions of boundary value problems (see, [1–7]) and fractional differential equations (see, [8–13]). It is well known that many topics in mathematical physics require the investigation of the eigenvalues and eigenfunctions of Sturm-Liouville type boundary value problems. The literature on computing eigenvalues of various types of Sturm-Liouville problems is little and we refer to [14–17].

Let 𝜎>0 and let PW²_𝜎 be the Paley-Wiener space of all 𝐿²(ℝ), entire functions of exponential type 𝜎. Assume that 𝑓(𝑡)∈PW²_𝜎⊂PW²_2𝜎. Then 𝑓(𝑡) can be reconstructed via the Hermite-type sampling series as 𝑓(𝑡)=^∞_𝑛=−∞𝑓𝑛𝜋𝜎𝑆²_𝑛(𝑡)+𝑓^′𝑛𝜋𝜎sin(𝜎𝑡−𝑛𝜋)𝜎𝑆_𝑛(𝑡), (1) where 𝑆_𝑛(𝑡) is the sequences of sinc functions as follows: 𝑆_𝑛(𝑡)()∶=sin𝜎𝑡−𝑛𝜋(𝜎𝑡−𝑛𝜋),𝑡≠𝑛𝜋𝜎,1,𝑡=𝑛𝜋𝜎. (2)

Series (1) converges absolutely and uniformly on ℝ, cf. [18–21]. Sometimes, series (1) is called the derivative sampling theorem. Our task is to use (1) to compute eigenvalues of Sturm-Liouville problems with eigenvalue parameter in boundary conditions numerically. This approach is a fully new technique that uses the recently obtained estimates for the truncation and amplitude errors associated with (1), cf. [22]. Both types of errors normally appear in numerical techniques that use interpolation procedures. In the following we summarize these estimates. The truncation error associated with (1) is defined to be 𝑅_𝑁𝑓(𝑡)∶=𝑓(𝑡)−𝑓_𝑁(𝑡),𝑁∈ℤ⁺,𝑡∈ℝ, (3) where 𝑓_𝑁(𝑡) is the truncated series as follows: 𝑓_𝑁(𝑡)=_{|𝑛|≤𝑁}𝑓𝑛𝜋𝜎𝑆²_𝑛(𝑡)+𝑓^′𝑛𝜋𝜎sin(𝜎𝑡−𝑛𝜋)𝜎𝑆_𝑛(𝑡). (4) It is proved in [22] that if 𝑓(𝑡)∈PW²_𝜎 and 𝑓(𝑡) is sufficiently smooth in the sense that there exists 𝑘∈ℤ⁺ such that 𝑡^𝑘𝑓(𝑡)∈𝐿²(ℝ), then for 𝑡∈ℝ, |𝑡|<𝑁𝜋/𝜎, we have 𝑅_𝑁𝑓(𝑡)≤𝑇_{𝑁,𝑘,𝜎}(𝑡)𝜉∶=^𝑘,𝜎𝐸_𝑘|sin𝜎𝑡|²√3(𝑁+1)^𝑘1(𝑁𝜋−𝜎𝑡)^3/2+1(𝑁𝜋+𝜎𝑡)^3/2+𝜉_𝑘,𝜎𝜎𝐸_𝑘+𝑘𝐸_𝑘−1|sin𝜎𝑡|²𝜎(𝑁+1)^𝑘×1√+1𝑁𝜋−𝜎𝑡√,𝑁𝜋+𝜎𝑡 (5) where the constants 𝐸_𝑘 and 𝜉_𝑘,𝜎 are given by 𝐸_𝑘∶=^∞_−∞𝑡^𝑘𝑓(𝑡)²𝑑𝑡,𝜉_𝑘,𝜎𝜎∶=^𝑘+1/2𝜋^𝑘+11−4^−𝑘. (6) The amplitude error occurs when approximate samples are used instead of the exact ones, which we cannot compute. It is defined to be 𝒜(𝑡)=𝜀,𝑓_{∞𝑛=−∞}𝑓𝑛𝜋𝜎−𝑓𝑛𝜋𝜎𝑆²_𝑛(𝑡)+𝑓^′𝑛𝜋𝜎−𝑓^′𝑛𝜋𝜎sin(𝜎𝑡−𝑛𝜋)𝜎𝑆_𝑛(𝑡),𝑡∈ℝ, (7) where 𝑓(𝑛𝜋/𝜎) and 𝑓^′(𝑛𝜋/𝜎) are approximate samples of 𝑓(𝑛𝜋/𝜎) and 𝑓^′(𝑛𝜋/𝜎), respectively. Let us assume that the differences 𝜀_𝑛∶=𝑓(𝑛𝜋/𝜎)−𝑓(𝑛𝜋/𝜎), 𝜀^′_𝑛∶=𝑓^′(𝑛𝜋/𝜎)−𝑓^′(𝑛𝜋/𝜎), and 𝑛∈ℤ are bounded by a positive number 𝜀, that is, |𝜀_𝑛|,|𝜀^′_𝑛|≤𝜀. If 𝑓(𝑡)∈PW²_𝜎 satisfies the natural decay conditions 𝜀_𝑛≤𝑓𝑛𝜋𝜎,𝜀^′_𝑛≤𝑓^′𝑛𝜋𝜎,𝑓(𝑡≤𝑀)^𝑓|𝑡|^𝜈+1,𝑡∈ℝ−{0}, (8)(9)0<𝜈≤1, then for 0<𝜀≤min{𝜋/𝜎,𝜎/𝜋,1/√𝑒}, we have, [22], 𝒜𝜀,𝑓_∞≤4𝑒^1/4𝜎()√𝜈+1()+𝜋3𝑒1+𝜎𝜎𝐴+𝑀_𝑓𝜌(𝜀)+𝜎+2+log(2)𝑀_𝑓1𝜀log𝜀, (10) where𝐴∶=3𝜎𝜋𝑓(0)+𝑀_𝑓𝜎𝜋^𝜈,𝜌(𝜀)1∶=𝛾+10log𝜀, (11) and 𝛾∶=lim_𝑛→∞[∑^𝑛_𝑘=11/𝑘−log𝑛]≅0.577216 is the Euler-Mascheroni constant.

The classical [23] sampling theorem of Whittaker, Kotelnikov, and Shannon (WKS) for 𝑓∈PW²_𝜎 is the series representation as follows: 𝑓(𝑡)=^∞_𝑛=−∞𝑓𝑛𝜋𝜎𝑆_𝑛(𝑡),𝑡∈ℝ, (12) where the convergence is absolute and uniform on ℝ and it is uniform on compact sets of ℂ cf. [23–25]. Series (12), which is of Lagrange interpolation type, has been used to compute eigenvalues of second-order eigenvalue problems, see for example, [17, 26–29]. The use of (12) in numerical analysis is known as the sinc method established by Stenger et al., cf. [30–32]. The aim of this paper is to investigate the possibilities of using Hermite interpolations rather than Lagrange interpolations, to compute the eigenvalues numerically. Notice that, due to Paley-Wiener's theorem [33] 𝑓∈PW²_𝜎 if and only if there is 𝑔(⋅)∈𝐿²(−𝜎,𝜎) such that 𝑓(𝑡)=1√2𝜋^𝜎−𝜎𝑔(𝑥)𝑒^𝑖𝑥𝑡𝑑𝑥. (13) Therefore, 𝑓^′(𝑡)∈PW²_𝜎, that is, 𝑓^′(𝑡) also has an expansion of the form (12). However, 𝑓^′(𝑡) can also be obtained by term-by-term differentiation formula of (12) as follows: 𝑓^′(𝑡)=^∞_𝑛=−∞𝑓𝑛𝜋𝜎𝑆^′_𝑛(𝑡), (14) see [23, page 52] for convergence. Thus, the use of Hermite interpolations will not cost any additional computational efforts since the samples 𝑓(𝑛𝜋/𝜎) will be used to compute both 𝑓(𝑡) and 𝑓^′(𝑡) according to (12) and (14), respectively.

Now, we consider the following differential equations: ℓ𝑦∶=−𝑦^′′𝑥,𝜇+𝑞(𝑥)𝑦𝑥,𝜇=𝜇²𝑦,[],𝑥,𝜇𝑥∈0,1 (15) with the following boundary conditions: 𝑎₁𝑦0,𝜇+𝑎₂𝑦^′0,𝜇=𝜇²𝑎^′₁𝑦0,𝜇+𝑎^′₂𝑦^′𝑏0,𝜇,₁𝑦1,𝜇+𝑏₂𝑦^′1,𝜇=𝜇²𝑏^′₁𝑦1,𝜇+𝑏^′₂𝑦^′1,𝜇, (16)(17) where 𝜇 is a complex spectral parameter, 𝑞(⋅) is assumed to be real valued and continuous on [0,1], and 𝑎_𝑖,𝑏_𝑖,𝑎^′_𝑖,𝑏^′_𝑖∈ℝ, 𝑖=0,1 satisfying 𝑎^′₁,𝑎^′₂=(0,0)or𝑎₁𝑎^′₂−𝑎^′₁𝑎₂,𝑏>0_′1,𝑏^′₂=(0,0)or𝑏^′₁𝑏₂−𝑏₁𝑏^′₂>0. (18) The eigenvalue problem (15)–(17) will be denoted by Π(𝑞,𝑎,𝑏,𝑎^′,𝑏^′) when (𝑎^′₁,𝑎^′₂)≠(0,0)≠(𝑏^′₁,𝑏^′₂). It is a Sturm-Liouville problem when the eigenparameter 𝜇 appears linearly in both boundary conditions. The classical problem when 𝑎^′₁=𝑎^′₂=𝑏^′₁=𝑏^′₂=0, which we denote by Π(𝑞,𝑎,𝑏,0,0) has a countable set of real and simple eigenvalues with ∞ as the only possible limit point, [34, 35]. In [14], the authors used Hermite-type sampling series (1) to compute the eigenvalues of problem Π(𝑞,𝑎,𝑏,0,0) numerically. In [36], see also [37], Annaby and Tharwat proved that Π(𝑞,𝑎,𝑏,𝑎^′,𝑏^′) has a denumerable set of real and simple eigenvalues with ∞ as the limit point using techniques similar of those established in [38–40], where also sampling theorems have been established. Similar results are established in [38] for the problem when the eigenparameter appears in one condition, that is, when 𝑎^′₁=𝑎^′₂=0, (𝑏^′₁,𝑏^′₂)≠(0,0) or equivalently when (𝑎^′₁,𝑎^′₂)≠(0,0) and 𝑏^′₁=𝑏^′₂=0. These problems will be denoted by Π(𝑞,𝑎,𝑏,0,𝑏^′), Π(𝑞,𝑎,𝑏,𝑎^′,0), respectively. The aim of the present work is to compute the eigenvalues of Π(𝑞,𝑎,𝑏,𝑎^′,𝑏^′), Π(𝑞,𝑎,𝑏,0,𝑏^′), and Π(𝑞,𝑎,𝑏,𝑎^′,0) numerically by the Hermite interpolations with an error analysis. This method is based on sampling theorem, Hermite interpolations, but applied to regularized functions. Hence, avoiding any (multiple) integration and keeping the number of terms in the Cardinal series manageable. It has been demonstrated that the method is capable of delivering higher order estimates of the eigenvalues at a very low cost, see [41–43]. In Sections 2 and 3 we derive the Hermite interpolation technique to compute the eigenvalues of Π(𝑞,𝑎,𝑏,𝑎^′,𝑏^′) and Π(𝑞,𝑎,𝑏,0,𝑏^′) with error estimates, respectively. The last section involves some illustrative examples.

2. Treatment of Π(𝑞,𝑎,𝑏,𝑎^′,𝑏^′)

In this section, we derive approximate values of the eigenvalues of Π(𝑞,𝑎,𝑏,𝑎^′,𝑏^′). Let 𝑦(⋅,𝜇) denote the solution of (15) satisfying the following initial conditions: 𝑦0,𝜇=𝑎₂−𝑎^′₂𝜇²,𝑦^′0,𝜇=𝑎^′₁𝜇²−𝑎₁. (19) Thus, 𝑦(⋅,𝜇) satisfies the boundary condition (16). The eigenvalues of the problem Π(𝑞,𝑎,𝑏,𝑎^′,𝑏^′) are the zeros of the function as follows: Δ𝜇∶=𝑏^′₁𝜇²−𝑏₁𝑦1,𝜇+𝑏^′₂𝜇²−𝑏₂𝑦^′1,𝜇. (20) These zeros are real and simple. The function Δ(𝜇) is an entire function of 𝜇. We aim to approximate Δ(𝜇) and hence its zeros, that is, the eigenvalues by the use of the Hermite Interpolation. The idea is to split Δ(𝜇) into two parts, one is known and the other is unknown, but lies in a Paley-Wiener space. Then we approximate the unknown part to get the approximate Δ(𝜇) and then compute the approximate zeros. Using the method of variation of constants, the solution 𝑦(𝑥,𝜇) satisfies Volterra integral equation as follows: 𝑦𝑥,𝜇=𝑎₂−𝑎^′₂𝜇²−𝑎cos𝜇𝑥₁−𝑎^′₁𝜇²sin𝜇𝑥𝜇+𝑇𝑦𝑥,𝜇, (21) where 𝑇 is the Volterra operator defined by 𝑇𝑦𝑥,𝜇=^𝑥₀sin𝜇(𝑥−𝑡)𝜇𝑞(𝑡)𝑦𝑡,𝜇𝑑𝑡. (22) Differentiating (21), we get 𝑦^′𝑥,𝜇=𝑎^′₂𝜇²−𝑎₂+𝑎𝜇sin𝜇𝑥_′1𝜇²−𝑎₁cos𝜇𝑥+𝑇𝑦𝑥,𝜇, (23) where 𝑇 is the Volterra operator 𝑇𝑦𝑥,𝜇=^𝑥₀cos𝜇(𝑥−𝑡)𝑞(𝑡)𝑦𝑡,𝜇𝑑𝑡. (24) Define 𝑓(⋅,𝜇) and 𝑔(⋅,𝜇) to be 𝑓𝑥,𝜇∶=𝑇𝑦𝑥,𝜇,𝑔𝑥,𝜇∶=𝑇𝑦𝑥,𝜇. (25) In the following, we will make use of the estimates [44] as follows: |cos𝑧|≤𝑒^|ℑ𝑧|,sin𝑧𝑧≤𝑐₀|𝑧|𝑒1+^|ℑ𝑧|, (26) where 𝑐₀ is some constant (we may take 𝑐₀≃1.72). For convenience, we define the constants by 𝜏∶=¹₀𝑞(𝑡)𝑑𝑡,𝑐₁∶=𝑎₂+𝑐₀𝑎₁,𝑐₂∶=𝑎^′₂+𝑐₀𝑎^′₁,𝑐₃∶=𝑐₀𝑐𝜏,₄∶=exp𝑐₃,𝑐₅∶=max𝑐₁,𝑐₂,𝑏₁+𝑏₂𝜏,𝑏^′₁+𝑏^′₂𝜏. (27) From (21) and (25), we get 𝑓=𝑥,𝜇_𝑥0sin𝜇(𝑥−𝑡)𝜇𝑞(𝑡)𝑎₂−𝑎^′₂𝜇²−𝑎cos𝜇𝑡₁−𝑎^′₁𝜇²sin𝜇𝑡𝜇+𝑑𝑡_𝑥0sin𝜇(𝑥−𝑡)𝜇𝑞(𝑡)𝑓𝑡,𝜇𝑑𝑡. (28)

Lemma 1. For 0≤𝑥≤1, 𝜇∈ℂ, the following estimates hold: 𝑓≤𝑐𝑥,𝜇³𝑐₄𝑐₁+𝑐₂𝜇²𝜇𝑒1+^{|ℑ𝜇|𝑥},𝑔𝑥,𝜇≤𝜏𝑐₃𝑐₄𝑐₁+𝑐₂𝜇²𝜇𝑒1+^{|ℑ𝜇|𝑥}. (29)(30)

Proof. We divide 𝑓(⋅,𝜇) into two parts 𝑓₁(⋅,𝜇) and 𝑓₂(⋅,𝜇) and estimate each of them. Indeed, for 𝑥∈[0,1] and 𝜇∈ℂ we have 𝑓₁=𝑥,𝜇_𝑥0sin𝜇(𝑥−𝑡)𝜇𝑞(𝑡)𝑎₂−𝑎^′₂𝜇²−𝑎cos𝜇𝑡₁−𝑎^′₁𝜇²sin𝜇𝑡𝜇𝑑𝑡≤𝑒^{|ℑ𝜇|𝑥𝑥}₀𝑞(𝑡𝑐)⁰(𝑥−𝑡)𝜇)×𝑎1+(𝑥−𝑡₂+𝑎^′₂𝜇²+𝑎₁+𝑎^′₁𝜇²×𝑐₀𝑡1+𝜇𝑡𝑑𝑡≤𝑒^{|ℑ𝜇|𝑥}𝑐₀𝑥1+𝜇𝑥^𝑥₀𝑞(𝑡)𝑎₂+𝑎^′₂𝜇²+𝑎₁+𝑎^′₁𝜇²𝑐₀𝑡𝑑𝑡≤𝑒^{|ℑ𝜇|𝑥}𝑐₀1+𝜇¹₀𝑞(𝑡)𝑎₂+𝑎^′₂𝜇²+𝑎₁+𝑎^′₁𝜇²𝑐₀𝑡𝑑𝑡. (31) Moreover, 0≤𝑥≤1, 𝜇∈ℂ, 𝑓₂𝑥,𝜇=^𝑥₀sin𝜇(𝑥−𝑡)𝜇𝑞(𝑡)𝑓𝑡,𝜇𝑑𝑡≤^𝑥₀𝑐₀(𝑥−𝑡)𝜇)𝑒1+(𝑥−𝑡^{|ℑ𝜇|(𝑥−𝑡)}×𝑞(𝑡)𝑓𝑡,𝜇𝑑𝑡≤𝑐₀𝑒^{|ℑ𝜇|𝑥𝑥}₀𝑒^{−|ℑ𝜇|𝑡}𝑞(𝑡)𝑓𝑡,𝜇𝑑𝑡. (32) Combining (31) and (32), we obtain 0≤𝑥≤1, 𝜇∈ℂ, 𝑓𝑥,𝜇≤𝑒^{|ℑ𝜇|𝑥}𝑐₀1+𝜇¹₀𝑞(𝑡)𝑎₂+𝑎^′₂𝜇²+𝑎₁+𝑎^′₁𝜇²𝑐₀𝑡𝑑𝑡+𝑐₀𝑒^{|ℑ𝜇|𝑥𝑥}₀𝑒^{−|ℑ𝜇|𝑡}𝑞(𝑡)𝑓𝑡,𝜇𝑑𝑡. (33) Applying Gronwall's inequality, cf. for example, [34, page 51], yields 𝜇∈ℂ, 𝑒^{−|ℑ𝜇|𝑥}𝑓≤𝑐𝑥,𝜇₀1+𝜇¹₀𝑞(𝑡)𝑎₂+𝑎^′₂𝜇²+𝑎₁+𝑎^′₁𝜇²𝑐₀𝑡𝑐𝑑𝑡×exp₀^𝑥₀𝑞(𝑡≤𝑐)𝑑𝑡₀1+𝜇¹₀𝑞(𝑡)𝑎₂+𝑎^′₂𝜇²+𝑎₁+𝑎^′₁𝜇²𝑐₀𝑡𝑐𝑑𝑡×exp₀¹₀𝑞(𝑡)𝑑𝑡, (34) from which we get 𝑓𝑥,𝜇≤𝑒^{|ℑ𝜇|𝑥}𝑐₀𝑎₂+𝑎^′₂𝜇²+𝑎₁+𝑎^′₁𝜇²𝑐₀1+𝜇×¹₀𝑞(𝑡)𝑑𝑡exp𝑐₀¹₀𝑞(𝑡=𝑐)𝑑𝑡₃𝑐₄𝑐₁+𝑐₂𝜇²𝜇𝑒1+^{|ℑ𝜇|𝑥}. (35) Then from (25) and (29), we obtain the estimate (30).

Now we split Δ(𝜇) into two parts via Δ𝜇=𝒢𝜇+𝒮𝜇, (36) where 𝒢(𝜇) is known part 𝒢𝜇=𝑏^′₁𝜇²−𝑏₁𝑎₂−𝑎^′₂𝜇²cos𝜇−𝑎₁−𝑎^′₁𝜇²sin𝜇𝜇+𝑏^′₂𝜇²−𝑏₂𝑎^′₂𝜇²−𝑎₂𝜇sin𝜇+𝑎^′₁𝜇²−𝑎₁cos𝜇, (37) and 𝒮(𝜇) is unknown part 𝒮𝜇=𝑏^′₁𝜇²−𝑏₁𝑓1,𝜇+𝑏^′₂𝜇²−𝑏₂𝑔1,𝜇. (38) Then, from Lemma 1, we have the following lemma.

Lemma 2. The function 𝒮(𝜇) is entire in 𝜇 and the following estimate holds: 𝒮𝜇≤𝑐₃𝑐₄𝑐₅1+𝜇²²𝜇𝑒1+^|ℑ𝜇|. (39)

Proof. Since 𝒮𝜇≤𝑏^′₁𝜇²+𝑏₁𝑓1,𝜇+𝑏^′₂𝜇²+𝑏₂𝑔1,𝜇, (40) then from (29) and (30) we get (39).

The analyticity of 𝒮(𝜇) and estimate (39) are not adequate to prove that 𝒮(𝜇) lies in a Paley-Wiener space. To solve this problem, we will multiply 𝒮(𝜇) by a regularization factor. Let 𝜃∈(0,1) and let 𝑚∈ℤ⁺, 𝑚>4 be fixed. Let ℱ_𝜃,𝑚(𝜇) be the function ℱ_𝜃,𝑚𝜇∶=sin𝜃𝜇𝜃𝜇^𝑚𝒮𝜇,𝜇∈ℂ. (41) More specifications on 𝑚,𝜃 will be given later on. Then we have the next lemma.

Lemma 3. ℱ_𝜃,𝑚(𝜇) is an entire function of 𝜇 which satisfies the estimates ℱ_𝜃,𝑚𝜇≤𝑐₃𝑐₄𝑐₅𝑐^𝑚₀1+𝜇²²1+𝜃𝜇^𝑚+1𝑒^{|ℑ𝜇|(1+𝑚𝜃)}. (42) Moreover, 𝜇^𝑚−4ℱ_𝜃,𝑚(𝜇)∈𝐿²(ℝ) and 𝐸_𝑚−4ℱ_𝜃,𝑚=^∞_−∞𝜇^𝑚−4ℱ_𝜃,𝑚𝜇²𝑑𝜇≤√2𝑐₃𝑐₄𝑐₅𝑐^𝑚₀𝜈₀, (43) where 𝜈₀∶=𝑚(2𝑚−1)+4𝜃²Γ[2𝑚+2]+144𝑚4𝑚²−1𝜃⁴×280𝜃⁴Γ[2𝑚−7]+20𝜃²Γ[][]×𝑚2𝑚−5+Γ2𝑚−34𝑚²−1Γ[2𝑚+2]𝜃^{2𝑚+1−11/2}. (44)

Proof. Since 𝒮(𝜇) is entire, then also ℱ_𝜃,𝑚(𝜇) is entire in 𝜇. Combining the estimates |sin𝑧/𝑧|≤(𝑐₀/(1+|𝑧|))𝑒^|ℑ𝑧| and (39), we obtain ℱ_𝜃,𝑚𝜇≤𝑐₀1+𝜃𝜇^𝑚𝑒^{|ℑ𝜇|𝑚𝜃}⋅𝑐₃𝑐₄𝑐₅1+𝜇²²𝜇𝑒1+^|ℑ𝜇|,𝜇∈ℂ, (45) leading to (42). Therefore, we get 𝜇^𝑚−4ℱ_𝜃,𝑚𝜇≤𝑐₃𝑐₄𝑐₅𝑐^𝑚₀𝜇^𝑚−41+𝜇²²1+𝜃𝜇^𝑚+1,𝜇∈ℝ. (46) That is, 𝜇^𝑚−4ℱ_𝜃,𝑚(𝜇)∈𝐿²(ℝ). Moreover, we get ^∞_−∞𝜇^𝑚−4ℱ_𝜃,𝑚𝜇²𝑑𝜇≤𝑐²₃𝑐²₄𝑐²₅𝑐₀^2𝑚×^∞_−∞𝜇^2𝑚−81+𝜇²⁴1+𝜃𝜇^2𝑚+2𝑑𝜇=2𝑐²₃𝑐²₄𝑐²₅𝑐₀^2𝑚𝜈²₀. (47)

What we have just proved is that ℱ_𝜃,𝑚(𝜇) belongs to the Paley-Wiener space PW²_𝜎 with 𝜎=1+𝑚𝜃. Since ℱ_𝜃,𝑚(𝜇)∈PW²_𝜎⊂PW²_2𝜎, then we can reconstruct the functions ℱ_𝜃,𝑚(𝜇) via the following sampling formula: ℱ_𝜃,𝑚𝜇=^∞_𝑛=−∞ℱ_𝜃,𝑚𝑛𝜋𝜎𝑆²_𝑛𝜇+ℱ^′_𝜃,𝑚𝑛𝜋𝜎sin𝜎𝜇−𝑛𝜋𝜎𝑆_𝑛𝜇. (48) Let 𝑁∈ℤ⁺, 𝑁>𝑚 and approximate ℱ_𝜃,𝑚(𝜇) by its truncated series ℱ_{𝜃,𝑚,𝑁}(𝜇), where ℱ_{𝜃,𝑚,𝑁}𝜇∶=^𝑁_𝑛=−𝑁ℱ_𝜃,𝑚𝑛𝜋𝜎𝑆²_𝑛𝜇+ℱ^′_𝜃,𝑚𝑛𝜋𝜎sin𝜎𝜇−𝑛𝜋𝜎𝑆_𝑛𝜇. (49) Since all eigenvalues are real, then from now on we restrict ourselves to 𝜇∈ℝ. Since 𝜇^𝑚−4ℱ_𝜃,𝑚(𝜇)∈𝐿²(ℝ), the truncation error, cf. (5), is given for |𝜇|<𝑁𝜋/𝜎 by ℱ_𝜃,𝑚𝜇−ℱ_{𝜃,𝑚,𝑁}𝜇≤𝑇_{𝑁,𝑚−4,𝜎}𝜇, (50) where 𝑇_{𝑁,𝑚−4,𝜎}𝜇𝜉∶=^{𝑚−4,𝜎}𝐸_𝑚−4sin𝜎𝜇²√3(𝑁+1)^𝑚−41𝑁𝜋−𝜎𝜇^3/2+1𝑁𝜋+𝜎𝜇^3/2+𝜉_{𝑚−4,𝜎}𝜎𝐸_𝑚−4+(𝑚−4)𝐸_𝑚−5sin𝜎𝜇²𝜎(𝑁+1)^𝑚−4×1+1𝑁𝜋−𝜎𝜇.𝑁𝜋+𝜎𝜇 (51) The samples {ℱ_𝜃,𝑚(𝑛𝜋/𝜎)}^𝑁_𝑛=−𝑁 and {ℱ^′_𝜃,𝑚(𝑛𝜋/𝜎)}^𝑁_𝑛=−𝑁, in general, are not known explicitly. So we approximate them by solving numerically 8𝑁+4 initial value problems at the nodes {𝑛𝜋/𝜎}^𝑁_𝑛=−𝑁. Let {ℱ_𝜃,𝑚(𝑛𝜋/𝜎)}^𝑁_𝑛=−𝑁 and let {ℱ^′_𝜃,𝑚(𝑛𝜋/𝜎)}^𝑁_𝑛=−𝑁 be the approximations of the samples of {ℱ_𝜃,𝑚(𝑛𝜋/𝜎)}^𝑁_𝑛=−𝑁 and {ℱ^′_𝜃,𝑚(𝑛𝜋/𝜎)}^𝑁_𝑛=−𝑁, respectively. Now we define ℱ_{𝜃,𝑚,𝑁}(𝜇), which approximates ℱ_{𝜃,𝑚,𝑁}(𝜇) as ℱ_{𝜃,𝑚,𝑁}𝜇∶=^𝑁_𝑛=−𝑁ℱ_𝜃,𝑚𝑛𝜋𝜎𝑆²_𝑛𝜇+ℱ^′_𝜃,𝑚𝑛𝜋𝜎sin𝜎𝜇−𝑛𝜋𝜎𝑆_𝑛𝜇,𝑁>𝑚. (52) Using standard methods for solving initial problems, we may assume that for |𝑛|<𝑁, ℱ_𝜃,𝑚𝑛𝜋𝜎−ℱ_𝜃,𝑚𝑛𝜋𝜎ℱ<𝜀,_′𝜃,𝑚𝑛𝜋𝜎−ℱ^′_𝜃,𝑚𝑛𝜋𝜎<𝜀, (53) for a sufficiently small 𝜀. From (42) we can see that ℱ_𝜃,𝑚(𝜇) satisfies the condition (9) when 𝑚>4 and therefore whenever 0<𝜀≤min{𝜋/𝜎,𝜎/𝜋,1/√𝑒} we have ℱ_{𝜃,𝑚,𝑁}𝜇−ℱ_{𝜃,𝑚,𝑁}𝜇≤𝒜(𝜀),𝜇∈ℝ, (54) where there is a positive constant 𝑀_ℱ𝜃,𝑚 for which, cf. (10), 𝒜(𝜀)∶=2𝑒^1/4𝜎√3𝑒(1+𝜎)+𝜋𝜎𝐴+𝑀_ℱ𝜃,𝑚𝜌(𝜀)+𝜎+2+log(2)𝑀_ℱ𝜃,𝑚1𝜀log𝜀. (55) Here 𝐴∶=3𝜎𝜋ℱ_𝜃,𝑚(0+𝜎)𝜋𝑀_ℱ𝜃,𝑚,𝜌(𝜀)1∶=𝛾+10log𝜀. (56) In the following we use the technique of [26], where only truncation error analysis is considered to determine enclosure intervals for the eigenvalues, see also [41]. Let 𝜇^∗2 be an eigenvalue; that is, Δ𝜇^∗=𝒢𝜇^∗+sin𝜃𝜇^∗𝜃𝜇^∗−𝑚ℱ_𝜃,𝑚𝜇^∗=0. (57) Then it follows that 𝒢𝜇^∗+sin𝜃𝜇^∗𝜃𝜇^∗−𝑚ℱ_{𝜃,𝑚,𝑁}𝜇^∗=sin𝜃𝜇^∗𝜃𝜇^∗−𝑚ℱ_{𝜃,𝑚,𝑁}𝜇^∗−sin𝜃𝜇^∗𝜃𝜇^∗−𝑚ℱ_𝜃,𝑚𝜇^∗=sin𝜃𝜇^∗𝜃𝜇^∗−𝑚ℱ_{𝜃,𝑚,𝑁}𝜇^∗−sin𝜃𝜇^∗𝜃𝜇^∗−𝑚ℱ_{𝜃,𝑚,𝑁}𝜇^∗+sin𝜃𝜇^∗𝜃𝜇^∗−𝑚ℱ_{𝜃,𝑚,𝑁}𝜇^∗−sin𝜃𝜇^∗𝜃𝜇^∗−𝑚ℱ_𝜃,𝑚𝜇^∗ (58) and so 𝒢𝜇^∗+sin𝜃𝜇^∗𝜃𝜇^∗−𝑚ℱ_{𝜃,𝑚,𝑁}𝜇^∗≤sin𝜃𝜇^∗𝜃𝜇^∗−𝑚𝑇_{𝑁,𝑚−4,𝜎}𝜇^∗+𝒜(𝜀). (59) Since 𝒢(𝜇^∗)+(sin𝜃𝜇^∗/𝜃𝜇^∗)^−𝑚ℱ_{𝜃,𝑚,𝑁}(𝜇^∗) is given and |sin𝜃𝜇^∗/𝜃𝜇^∗|^−𝑚(𝑇_{𝑁,𝑚−4,𝜎}(𝜇^∗)+𝒜(𝜀)) has computable upper bound, we can define an enclosure for 𝜇^∗ by solving the following system of inequalities: −sin𝜃𝜇^∗𝜃𝜇^∗−𝑚𝑇_{𝑁,𝑚−4,𝜎}𝜇^∗(𝜀)𝜇+𝒜≤𝒢^∗+sin𝜃𝜇^∗𝜃𝜇^∗−𝑚ℱ_{𝜃,𝑚,𝑁}𝜇^∗≤sin𝜃𝜇^∗𝜃𝜇^∗−𝑚𝑇_{𝑁,𝑚−4,𝜎}𝜇^∗+𝒜(𝜀). (60) Its solution is an interval containing 𝜇^∗, and over which the graph 𝒢(𝜇^∗)+(sin𝜃𝜇^∗/𝜃𝜇^∗)^−𝑚ℱ_{𝜃,𝑚,𝑁}(𝜇^∗) is squeezed between the graphs as follows: −sin𝜃𝜇^∗𝜃𝜇^∗−𝑚𝑇_{𝑁,𝑚−4,𝜎}𝜇^∗+𝒜(𝜀),sin𝜃𝜇^∗𝜃𝜇^∗−𝑚𝑇_{𝑁,𝑚−4,𝜎}𝜇^∗+𝒜(𝜀). (61) Using the fact that ℱ_{𝜃,𝑚,𝑁}𝜇⟶ℱ_𝜃,𝑚𝜇 (62) uniformly over any compact set and since 𝜇^∗ is a simple root, we obtain the following for large 𝑁 and sufficiently small 𝜀: 𝜕𝒢𝜇+𝜕𝜇sin𝜃𝜇𝜃𝜇^−𝑚ℱ_{𝜃,𝑚,𝑁}𝜇≠0 (63) in a neighborhood of 𝜇^∗. Hence, the graph of 𝒢(𝜇)+(sin𝜃𝜇/𝜃𝜇)^−𝑚ℱ_{𝜃,𝑚,𝑁}(𝜇) intersects the graphs −|sin𝜃𝜇/𝜃𝜇|^−𝑚(𝑇_{𝑁,𝑚−4,𝜎}(𝜇)+𝒜(𝜀)) and |sin𝜃𝜇/𝜃𝜇|^−𝑚(𝑇_{𝑁,𝑚−4,𝜎}(𝜇)+𝒜(𝜀)) at two points with abscissae 𝑎₋(𝜇^∗,𝑁,𝜀)≤𝑎₊(𝜇^∗,𝑁,𝜀) and the solution of the system of inequalities (60) is the interval 𝐼_𝜀,𝑁∶=𝑎₋𝜇^∗,𝑁,𝜀,𝑎₊𝜇^∗,𝑁,𝜀 (64) and in particular 𝜇^∗∈𝐼_𝜀,𝑁. Summarizing the above discussion, we arrive at the following lemma which is similar to that of [26].

Lemma 4. For any eigenvalue 𝜇^∗2, we can find 𝑁₀∈ℤ⁺ and sufficiently small 𝜀 such that 𝜇^∗∈𝐼_𝜀,𝑁  for 𝑁>𝑁₀. Moreover, we get 𝑎₋𝜇^∗,𝑁,𝜀,𝑎₊𝜇^∗,𝑁,𝜀⟶𝜇^∗𝑎𝑠𝑁⟶∞,𝜀⟶0. (65)

Proof. Since all eigenvalues of Π(𝑞,𝑎,𝑏,𝑎^′,𝑏^′) are simple, then for large 𝑁 and sufficiently small 𝜀 we have (𝜕/𝜕𝜇)(𝒢(𝜇)+(sin𝜃𝜇/𝜃𝜇)^−𝑚ℱ_{𝜃,𝑚,𝑁}(𝜇))>0, in a neighborhood of 𝜇^∗. Choose 𝑁₀ such that 𝒢𝜇+sin𝜃𝜇𝜃𝜇^−𝑚ℱ_{𝜃,𝑚,𝑁0}𝜇=±sin𝜃𝜇𝜃𝜇^−𝑚𝑇_{𝑁0,𝑚−4,𝜎}𝜇+𝒜(𝜀) (66) has two distinct solutions which we denote by 𝑎₋(𝜇^∗,𝑁₀,𝜀)≤𝑎₊(𝜇^∗,𝑁₀,𝜀). The decay of 𝑇_{𝑁,𝑚−4,𝜎}(𝜇)→0 as 𝑁→∞ and 𝒜(𝜀)→0 as 𝜀→0 will ensure the existence of the solutions 𝑎₋(𝜇^∗,𝑁,𝜀) and 𝑎₊(𝜇^∗,𝑁,𝜀) as 𝑁→∞ and 𝜀→0. For the second point we recall that ℱ_{𝜃,𝑚,𝑁}(𝜇)→ℱ_𝜃,𝑚(𝜇) as 𝑁→∞ and as 𝜀→0. Hence, by taking the limit we obtain 𝒢𝑎₊𝜇^∗,∞,0+sin𝜃𝜇^∗𝜃𝜇^∗−𝑚ℱ_𝜃,𝑚𝑎₊𝜇^∗𝒢𝑎,∞,0=0,₋𝜇^∗,∞,0+sin𝜃𝜇^∗𝜃𝜇^∗−𝑚ℱ_𝜃,𝑚𝑎₋𝜇^∗,∞,0=0. (67) That is, Δ(𝑎₊)=Δ(𝑎₋)=0. This leads us to conclude that 𝑎₊=𝑎₋=𝜇^∗, since 𝜇^∗ is a simple root.

Let Δ_𝑁(𝜇)∶=𝒢(𝜇)+(sin𝜃𝜇/𝜃𝜇)^−𝑚ℱ_{𝜃,𝑚,𝑁}(𝜇). Then (50) and (54) imply Δ𝜇−Δ_𝑁𝜇≤sin𝜃𝜇𝜃𝜇^−𝑚𝑇_{𝑁,𝑚−4,𝜎}𝜇(𝜀),𝜇<+𝒜𝑁𝜋𝜎 (68) and 𝜃 is chosen sufficiently small for which |𝜃𝜇|<𝜋. Therefore, 𝜃,𝑚 must be chosen so that for |𝜇|<𝑁𝜋/𝜎𝑚>4,𝜃∈(0,1),𝜃𝜇<𝜋. (69) Let 𝜇^∗ be an eigenvalue and let 𝜇_𝑁 be its approximation. Thus, Δ(𝜇^∗)=0 and Δ_𝑁(𝜇_𝑁)=0. From (68) we have |Δ_𝑁(𝜇^∗)|≤|sin𝜃𝜇^∗/𝜃𝜇^∗|^−𝑚(𝑇_{𝑁,𝑚−4,𝜎}(𝜇^∗)+𝒜(𝜀)). Now we estimate the error |𝜇^∗−𝜇_𝑁| for an eigenvalue 𝜇^∗.

Theorem 5. Let 𝜇^∗2 be an eigenvalue of Π(𝑞,𝑎,𝑏,𝑎^′,𝑏^′). For sufficient large 𝑁 we have the following estimate: 𝜇^∗−𝜇_𝑁<sin𝜃𝜇_𝑁𝜃𝜇_𝑁^−𝑚𝑇_{𝑁,𝑚−4,𝜎}𝜇_𝑁+𝒜(𝜀)inf_{𝜁∈𝐼𝜀,𝑁}Δ^′(𝜁.) (70)

Proof. Since Δ(𝜇_𝑁)−Δ_𝑁(𝜇_𝑁)=Δ(𝜇_𝑁)−Δ(𝜇^∗), then from (68) and after replacing 𝜇 by 𝜇_𝑁, we obtain Δ𝜇_𝑁−Δ𝜇^∗≤sin𝜃𝜇_𝑁𝜃𝜇_𝑁^−𝑚𝑇_{𝑁,𝑚−4,𝜎}𝜇_𝑁+𝒜(𝜀). (71) Using the mean value theorem yields that for some 𝜁∈𝐽_𝜀,𝑁∶=[min(𝜇^∗,𝜇_𝑁),max(𝜇^∗,𝜇_𝑁)], 𝜇^∗−𝜇_𝑁Δ^′(𝜁≤)sin𝜃𝜇_𝑁𝜃𝜇_𝑁^−𝑚𝑇_{𝑁,𝑚−4,𝜎}𝜇_𝑁+𝒜(𝜀),𝜁∈𝐽_𝜀,𝑁⊂𝐼_𝜀,𝑁. (72) Since the eigenvalues are simple, then for sufficiently large 𝑁inf_{𝜁∈𝐼𝜀,𝑁}|Δ^′(𝜁)|>0 and we get (70).

3. The Case of Π(𝑞,𝑎,𝑏,0,𝑏^′)

This section includes briefly a treatment similarly to that of the previous section for the eigenvalue problem Π(𝑞,𝑎,𝑏,0,𝑏^′) introduced in Section 1. Notice that condition (18) implies that the analysis of problem Π(𝑞,𝑎,𝑏,0,𝑏^′) is not included in that of Π(𝑞,𝑎,𝑏,𝑎^′,𝑏^′). Let 𝜓(⋅,𝜇) denote the solution of (15) satisfying the following initial conditions: 𝜓0,𝜇=𝑎₂,𝜓^′0,𝜇=−𝑎₁. (73) Thus, 𝜓(⋅,𝜇) satisfies the boundary condition (16). The eigenvalues of the problem Π(𝑞,𝑎,𝑏,0,𝑏^′) are the zeros of the function as follows: Ω𝜇∶=𝑏^′₁𝜇²−𝑏₁𝜓1,𝜇+𝑏^′₂𝜇²−𝑏₂𝜓^′1,𝜇. (74) Recall that Π(𝑞,𝑎,𝑏,𝑎^′,𝑏^′) has denumerable set of real and simple eigenvalues, cf. [38]. Using the method of variation of constants, the solution 𝜓(𝑥,𝜇) satisfies Volterra integral equation as follows: 𝜓𝑥,𝜇=𝑎₂cos𝜇𝑥−𝑎₁sin𝜇𝑥𝜇+𝑇𝜓𝑥,𝜇, (75) where 𝑇 is the Volterra operator defined in (22). Differentiating (75), we get 𝜓^′𝑥,𝜇=−𝑎₂𝜇sin𝜇𝑥−𝑎₁cos𝜇𝑥+𝑇𝜓𝑥,𝜇, (76) where 𝑇 is the Volterra operator defined in (24). Define ℎ₁(⋅,𝜇) and ℎ₂(⋅,𝜇) to be ℎ₁𝑥,𝜇∶=𝑇𝜓𝑥,𝜇,ℎ₂𝑥,𝜇∶=𝑇𝜓𝑥,𝜇. (77) As in the preceding section we split Ω(𝜇) into Ω𝜇∶=𝒦𝜇+𝒰𝜇, (78) where 𝒦(𝜇) is the known part 𝒦𝜇=𝑏^′₁𝜇²−𝑏₁𝑎₂cos𝜇−𝑎₁sin𝜇𝜇+𝑏^′₂𝜇²−𝑏₂−𝑎₂𝜇sin𝜇−𝑎₁cos𝜇, (79) and 𝒰(𝜇) is the unknown one 𝒰𝜇∶=𝑏^′₁𝜇²−𝑏₁ℎ₁1,𝜇+𝑏^′₂𝜇²−𝑏₂ℎ₂1,𝜇. (80) Then, as in the previous section, 𝒰(𝜇) is entire in 𝜇 for each 𝑥∈[0,1] for which 𝒰𝜇≤𝑐₁𝑐₃𝑐₄𝑐₆1+𝜇²𝜇𝑒1+^|ℑ𝜇|,𝜇∈ℂ, (81) where 𝑐₆∶=max{|𝑏₁|+|𝑏₂|𝜏,|𝑏^′₁|+|𝑏^′₂|𝜏}.

Let 𝜃∈(0,1) and let 𝑚 be as in the previous section, but 𝑚>2. Define ℛ_𝑚,𝜃(𝜇) to be ℛ_𝑚,𝜃𝜇=sin𝜃𝜇𝜃𝜇^𝑚𝒰𝜇,𝜇∈ℂ. (82) Hence, ℛ_𝑚,𝜃𝜇≤𝑐^𝑚₀𝑐₁𝑐₃𝑐₄𝑐₆1+𝜇²1+𝜃𝜇^𝑚+1𝑒^{|ℑ𝜇|(1+𝑚𝜃)},𝜇∈ℂ, (83) and 𝜇^𝑚−2ℛ_𝑚,𝜃(𝜇)∈𝐿²(ℝ) with 𝐸_𝑚−2ℛ_𝑚,𝜃=^∞_−∞𝜇^𝑚−2ℛ_𝑚,𝜃𝜇²≤√𝑑𝜇2𝑐^𝑚₀𝑐₁𝑐₃𝑐₄𝑐₆𝜔₀, (84) where 𝜔₀∶=12𝜃²Γ[2𝑚−3]+Γ[2𝑚−1]𝜃^2𝑚−1Γ[2(.𝑚+1)] (85) Thus, ℛ_𝑚,𝜃(𝜇) belongs to the Paley-Wiener space PW²_𝜎 with 𝜎=1+𝑚𝜃. Since ℛ_𝜃,𝑚(𝜇)∈PW²_𝜎⊂PW²_2𝜎, then we can reconstruct the functions ℛ_𝜃,𝑚(𝜇) via the following sampling formula: ℛ_𝜃,𝑚𝜇=^∞_𝑛=−∞ℛ_𝜃,𝑚𝑛𝜋𝜎𝑆²_𝑛𝜇+ℛ^′_𝜃,𝑚𝑛𝜋𝜎sin𝜎𝜇−𝑛𝜋𝜎𝑆_𝑛𝜇. (86) Let 𝑁∈ℤ⁺, 𝑁>𝑚, and approximate ℛ_𝜃,𝑚(𝜇) by its truncated series ℛ_{𝜃,𝑚,𝑁}(𝜇), where ℛ_{𝜃,𝑚,𝑁}𝜇∶=^𝑁_𝑛=−𝑁ℛ_𝜃,𝑚𝑛𝜋𝜎𝑆²_𝑛𝜇+ℛ^′_𝜃,𝑚𝑛𝜋𝜎sin𝜎𝜇−𝑛𝜋𝜎𝑆_𝑛𝜇. (87) Since all eigenvalues are real, then from now on we restrict ourselves to 𝜇∈ℝ. Since 𝜇^𝑚−2ℛ_𝜃,𝑚(𝜇)∈𝐿²(ℝ), the truncation error, cf. (5), is given for |𝜇|<𝑁𝜋/𝜎 by ℛ_𝜃,𝑚𝜇−ℛ_{𝜃,𝑚,𝑁}𝜇≤𝑇_{𝑁,𝑚−2,𝜎}𝜇, (88) where 𝑇_{𝑁,𝑚−2,𝜎}𝜇𝜉∶=^{𝑚−2,𝜎}𝐸_𝑚−2sin𝜎𝜇²√3(𝑁+1)^𝑚−21𝑁𝜋−𝜎𝜇^3/2+1𝑁𝜋+𝜎𝜇^3/2+𝜉_{𝑚−2,𝜎}𝜎𝐸_𝑚−2+(𝑚−2)𝐸_𝑚−3sin𝜎𝜇²𝜎(𝑁+1)^𝑚−2×1+1𝑁𝜋−𝜎𝜇.𝑁𝜋+𝜎𝜇 (89) The samples {ℛ_𝜃,𝑚(𝑛𝜋/𝜎)}^𝑁_𝑛=−𝑁 and {ℛ^′_𝜃,𝑚(𝑛𝜋/𝜎)}^𝑁_𝑛=−𝑁, in general, are not known explicitly. So we approximate them by solving numerically 4𝑁+2 initial value problems at the nodes {𝑛𝜋/𝜎}^𝑁_𝑛=−𝑁. Let {ℛ_𝜃,𝑚(𝑛𝜋/𝜎)}^𝑁_𝑛=−𝑁 and let {ℛ^′_𝜃,𝑚(𝑛𝜋/𝜎)}^𝑁_𝑛=−𝑁 be the approximations of the samples of {ℛ_𝜃,𝑚(𝑛𝜋/𝜎)}^𝑁_𝑛=−𝑁 and {ℛ^′_𝜃,𝑚(𝑛𝜋/𝜎)}^𝑁_𝑛=−𝑁, respectively. Now we define ℛ_{𝜃,𝑚,𝑁}(𝜇), which approximates ℛ_{𝜃,𝑚,𝑁}(𝜇)ℛ_{𝜃,𝑚,𝑁}𝜇∶=^𝑁_𝑛=−𝑁ℛ_𝜃,𝑚𝑛𝜋𝜎𝑆²_𝑛𝜇+ℛ^′_𝜃,𝑚𝑛𝜋𝜎sin𝜎𝜇−𝑛𝜋𝜎𝑆_𝑛𝜇,𝑁>𝑚. (90)

Using standard methods for solving initial problems, we may assume that for |𝑛|<𝑁ℛ_𝜃,𝑚𝑛𝜋𝜎−ℛ_𝜃,𝑚𝑛𝜋𝜎ℛ<𝜀,_′𝜃,𝑚𝑛𝜋𝜎−ℛ^′_𝜃,𝑚𝑛𝜋𝜎<𝜀, (91) for a sufficiently small 𝜀. From (83) we can see that ℛ_𝜃,𝑚(𝜇) satisfies the condition (9) when 𝑚>2 and therefore whenever 0<𝜀≤min{𝜋/𝜎,𝜎/𝜋,1/√𝑒} we have ℛ_{𝜃,𝑚,𝑁}𝜇−ℛ_{𝜃,𝑚,𝑁}𝜇≤𝒜(𝜀),𝜇∈ℝ, (92) where there is a positive constant 𝑀_ℛ𝜃,𝑚 for which, cf. (10), and 𝒜(𝜀)∶=2𝑒^1/4𝜎√3𝑒(1+𝜎)+𝜋𝜎𝐴+𝑀_ℛ𝜃,𝑚𝜌(𝜀)+𝜎+2+log(2)𝑀_ℛ𝜃,𝑚1𝜀log𝜀. (93) Here 𝐴∶=3𝜎𝜋ℛ_𝜃,𝑚(0+𝜎)𝜋𝑀_ℛ𝜃,𝑚,𝜌(𝜀)1∶=𝛾+10log𝜀. (94) As in the above section, we have the following lemma.

Lemma 6. For any eigenvalue 𝜇^∗2 of the problem Π(𝑞,𝑎,𝑏,0,𝑏^′), we can find 𝑁₀∈ℤ⁺ and sufficiently small 𝜀 such that 𝜇^∗∈ℐ_𝜀,𝑁 for 𝑁>𝑁₀, where ℐ_𝜀,𝑁∶=𝑏₋𝜇^∗,𝑁,𝜀,𝑎𝑏₊𝜇^∗,𝑁,𝜀, (95)𝑏₋,𝑏₊ are the solutions of the inequalities −sin𝜃𝜇𝜃𝜇^−𝑚𝑇_{𝑁,𝑚−2,𝜎}𝜇(𝜀)≤Ω+𝒜_𝑁𝜇≤sin𝜃𝜇𝜃𝜇^−𝑚𝑇_{𝑁,𝑚−2,𝜎}𝜇+𝒜(𝜀). (96) Moreover, we get 𝑏₋𝜇^∗,𝑁,𝜀,𝑏₊𝜇^∗,𝑁,𝜀⟶𝜇^∗as𝑁⟶∞,𝜀⟶0. (97)

Let Ω_𝑁(𝜇)∶=𝒦(𝜇)+(sin𝜃𝜇/𝜃𝜇)^−𝑚ℛ_{𝜃,𝑚,𝑁}(𝜇). Then (88) and (92) imply Ω𝜇−Ω_𝑁𝜇≤sin𝜃𝜇𝜃𝜇^−𝑚𝑇_{𝑁,𝑚−2,𝜎}𝜇+𝒜(𝜀),𝜇<𝑁𝜋𝜎 (98) and 𝜃 is chosen sufficiently small for which |𝜃𝜇|<𝜋. Therefore, 𝜃,𝑚 must be chosen so that for |𝜇|<𝑁𝜋/𝜎𝑚>2,𝜃∈(0,1),𝜃𝜇<𝜋. (99) Let 𝜇^∗ be an eigenvalue and 𝜇_𝑁 be its approximation. Thus Ω(𝜇^∗)=0 and Ω_𝑁(𝜇_𝑁)=0. From (98) we have |Ω_𝑁(𝜇^∗)|≤|sin𝜃𝜇^∗/𝜃𝜇^∗|^−𝑚(𝑇_{𝑁,𝑚−2,𝜎}(𝜇^∗)+𝒜(𝜀)). Now we estimate the error |𝜇^∗−𝜇_𝑁| for an eigenvalue 𝜇^∗. Finally we have the following estimate.

Theorem 7. Let 𝜇^∗2 be an eigenvalue of the problem Π(𝑞,𝑎,𝑏,0,𝑏^′). For sufficient large 𝑁 we have the following estimate 𝜇^∗−𝜇_𝑁<sin𝜃𝜇_𝑁𝜃𝜇_𝑁^−𝑚𝑇_{𝑁,𝑚−2,𝜎}𝜇_𝑁+𝒜(𝜀)inf_{𝜁∈ℐ𝜀,𝑁}Ω^′(𝜁.) (100)

4. Numerical Examples

This section includes two detailed worked examples illustrating the above technique. Examples 1 and 2 computed in [27, 45] with the classical sinc method, where only truncation error analysis is considered, respectively. It is clearly seen that our new method (Hermite interpolations) gives remarkably better results than in [27, 45], see also [41–43]. We indicate in these examples the effect of the amplitude error in the method by determining enclosure intervals for different values of 𝜀. We also indicate the effect of the parameters 𝑚 and 𝜃 by several choices. Each example is exhibited via figures that accurately illustrate the procedure near to some of the approximated eigenvalues. More explanations are given below. Recall that 𝑎_±(𝜇) and 𝑏_±(𝜇) are defined by 𝑎_±𝜇=Δ_𝑁𝜇±sin𝜃𝜇𝜃𝜇^−𝑚𝑇_{𝑁,𝑚−3,𝜎}𝜇(𝜀),𝜇<+𝒜𝑁𝜋𝜎,𝑏_±𝜇=Ω_𝑁𝜇±sin𝜃𝜇𝜃𝜇^−𝑚𝑇_{𝑁,𝑚−2,𝜎}𝜇(𝜀),𝜇<+𝒜𝑁𝜋𝜎, (101)(102) respectively. Recall also that the enclosure intervals 𝐼_𝜀,𝑁∶=[𝑎₋,𝑎₊] and ℐ_𝜀,𝑁∶=[𝑏₋,𝑏₊] are determined by solving 𝑎_±𝜇=0,𝜇<𝑁𝜋𝜎,𝑏_±𝜇=0,𝜇<𝑁𝜋𝜎. (103)(104) respectively. We would like to mention that Mathematica has been used to obtain the exact values for the three examples where eigenvalues cannot be computed concretely. Mathematica is also used in rounding the exact eigenvalues, which are square roots.

Example 1. The boundary value problem [27] −𝑦^′′𝑥,𝜇−𝑦𝑥,𝜇=𝜇²𝑦𝑦𝑥,𝜇,0≤𝑥≤1,_′𝑦0,𝜇=0,1,𝜇+𝑦^′1,𝜇=𝜇²2𝑦1,𝜇+𝑦^′1,𝜇, (105) is a special case of the problem treated in the previous section with 𝑎₁=0, 𝑎₂=𝑏₁=𝑏₂=𝑏^′₂=1, 𝑏^′₁=2 and 𝑞(𝑥)=−1. The characteristic function is Ω𝜇=−1+2𝜇²cos1+𝜇²−−1+𝜇²1+𝜇²sin1+𝜇². (106)
The function 𝒦(𝜇) will be 𝒦𝜇=−1+2𝜇²cos𝜇−𝜇−1+𝜇²sin𝜇. (107) As is clearly seen, the eigenvalues cannot be computed explicitly. Tables 1, 2, and 3 indicate the application of our technique to this problem and the effect of 𝑚, 𝜃 and 𝜀. By exact we mean the zeros of Ω(𝜇) computed by Mathematica.
Figures 1 and 2 illustrate the enclosure intervals dominating 𝜇₁ for 𝑁=15, 𝑚=10, 𝜃=1/5, and 𝜀=10⁻⁵ and 𝜀=10⁻¹⁰ respectively. The middle curve represents Ω(𝜇), while the upper and lower curves represent the curves of 𝑏₊(𝜇),𝑏₋(𝜇), respectively. We notice that when 𝜀=10⁻¹⁰, the two curves are almost identical. Similarly, Figures 3 and 4 illustrate the enclosure intervals dominating 𝜇₂ for 𝑁=15, 𝑚=10, 𝜃=1/5, and 𝜀=10⁻⁵and𝜀=10⁻¹⁰ respectively.

Table 1 

With 𝑁=15, the approximation 𝜇_𝑘,𝑁 and the exact solution 𝜇_𝑘 for different choices of 𝑚 and 𝜃.

Table 2 

Absolute error |𝜇_𝑘−𝜇_𝑘,𝑁|.

Table 3 

For 𝑁=15, 𝑚=10, and 𝜃=1/5, the exact solution 𝜇_𝑘 are all inside the interval [𝑏₋, 𝑏₊] for different values of ε.

Figure 1 

The enclosure interval dominating 𝜇1 for 𝑚=10, 𝑁=15, 𝜃=1/5, and 𝜀=10−5.

Figure 2 

The enclosure interval dominating 𝜇1 for 𝑚=10, 𝑁=15, 𝜃=1/5, and 𝜀=10−10.

Figure 3 

The enclosure interval dominating 𝜇2 for 𝑚=10, 𝑁=15, 𝜃=1/5, and 𝜀=10−5.

Figure 4 

The enclosure interval dominating 𝜇2 for 𝑚=10, 𝑁=15, 𝜃=1/5, and 𝜀=10−10.

Example 2. The boundary value problem [45] −𝑦^′′𝑥,𝜇−𝑦𝑥,𝜇=𝜇²𝑦𝑦𝑥,𝜇0≤𝑥≤1,0,𝜇=𝜇²𝑦^′0,𝜇,𝑦^′1,𝜇=𝜇²𝑦1,𝜇, (108) is a special case of problem Π(𝑞,𝑎,𝑏,𝑎^′,𝑏^′), when 𝑞(𝑥)=−1, 𝑎₂=𝑎^′₁=𝑏₁=𝑏^′₂=0, and 𝑎₁=𝑎^′₂=𝑏^′₁=𝑏₂=1. Here the characteristic function is Δ𝜇=1−𝜇⁴cos𝜇²+1−2𝜇²+𝜇⁴sin𝜇²+1𝜇²+1. (109) After computing 𝒢(𝜇), we obtain 𝒢𝜇=1+𝜇²1−𝜇²cos𝜇−𝜇sin𝜇. (110) As is clearly seen, the eigenvalues cannot be computed explicitly. As in the previous example, Figures 5, 6, 7, and 8 illustrate the results of Tables 4, 5, 6, and 7.

Table 4 

With 𝑁=40, the approximation 𝜇_𝑘,𝑁 and the exact solution 𝜇_𝑘 for different choices of θ.

Table 5 

Absolute error |𝜇_𝑘−𝜇_𝑘,𝑁|.

Table 6 

The approximation 𝜇_𝑘,𝑁 and the exact solution 𝜇_𝑘 for 𝑁=40, 𝑚=14 and 𝜃=1/26.

Table 7 

For 𝑁=40, 𝑚=14 and 𝜃=1/26, the exact solution 𝜇_𝑘 are all inside the interval [𝑎₋, 𝑎₊] for different values of ε.

Figure 5 

The enclosure interval dominating 𝜇1 for 𝑚=14, 𝑁=40, 𝜃=1/26, and 𝜀=10−5.

Figure 6 

The enclosure interval dominating 𝜇1 for 𝑚=14, 𝑁=40, 𝜃=1/26, and 𝜀=10−10.

Figure 7 

The enclosure interval dominating 𝜇4 for 𝑚=14, 𝑁=40, 𝜃=1/26, and 𝜀=10−5.

Figure 8 

The enclosure interval dominating 𝜇4 for 𝑚=14, 𝑁=40, 𝜃=1/26, and 𝜀=10−10.

Acknowledgment

This research has been supported by a grant from the Institute of Scientific Research at Umm AL-Qura University, Saudi Arabia.