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

Authors: M. M. Tharwat; S. M. Al-Harbi
Journal: Abstract and Applied Analysis
 (2013)

---

## 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.


---

## Body

## 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σ2 be the Paley-Wiener space of all L2(ℝ), entire functions of exponential type σ. Assume that f(t)∈PWσ2⊂PW2σ2. Then f(t) can be reconstructed via the Hermite-type sampling series as
(1)f(t)=∑n=-∞∞‍[sin(σt-nπ)σf(nπσ)Sn2(t)+f′(nπσ)sin(σt-nπ)σSn(t)],
where Sn(t) is the sequences of sinc functions as follows:
(2)Sn(t):={sin(σt-nπ)(σt-nπ),t≠nπσ,1,t=nπσ.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)RN(f)(t):=f(t)-fN(t),N∈ℤ+,t∈ℝ,
where fN(t) is the truncated series as follows:
(4)fN(t)=∑|n|≤N‍[sin(σt-nπ)σf(nπσ)Sn2(t)+f′(nπσ)sin(σt-nπ)σSn(t)].
It is proved in [22] that if f(t)∈PWσ2 and f(t) is sufficiently smooth in the sense that there exists k∈ℤ+ such that tkf(t)∈L2(ℝ), then for t∈ℝ, |t|<Nπ/σ, we have
(5)|RN(f)(t)|≤TN,k,σ(t)≔ξk,σEk|sinσt|23(N+1)k(1(Nπ-σt)3/2+1(Nπ+σt)3/2)+ξk,σ(σEk+kEk-1)|sinσt|2σ(N+1)k×(1Nπ-σt+1Nπ+σt),
where the constants Ek and ξk,σ are given by
(6)Ek:=∫-∞∞‍|tkf(t)|2dt,ξk,σ:=σk+1/2πk+11-4-k.
The amplitude error occurs when approximate samples are used instead of the exact ones, which we cannot compute. It is defined to be
(7)𝒜(ɛ,f)(t)=∑n=-∞∞‍[{f(nπσ)-f~(nπσ)}Sn2(t)+{f′(nπσ)-f′~(nπσ)}sin(σt-nπ)σSn(t)],iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiit∈ℝ,
where f~(nπ/σ) and f′~(nπ/σ) are approximate samples of f(nπ/σ) and f′(nπ/σ), respectively. Let us assume that the differences ɛn:=f(nπ/σ)-f~(nπ/σ), ɛn′:=f′(nπ/σ)-f′~(nπ/σ), and n∈ℤ are bounded by a positive number ɛ, that is, |ɛn|,|ɛn′|≤ɛ. If f(t)∈PWσ2 satisfies the natural decay conditions
(8)|ɛn|≤|f(nπσ)|,|ɛn′|≤|f′(nπσ)|,(9)|f(t)|≤Mf|t|ν+1,t∈ℝ-{0},0<ν≤1, then for 0<ɛ≤min{π/σ,σ/π,1/e}, we have, [22],
(10)∥𝒜(ɛ,f)∥∞≤4e1/4σ(ν+1){3e(1+σ)+((πσ)A+Mf)ρ(ɛ)+(σ+2+log(2))Mf((πσ)A+Mf)}ɛlog(1ɛ),
where(11)A:=3σπ(|f(0)|+Mf(σπ)ν),ρ(ɛ):=γ+10log(1ɛ),
and γ:=limn→∞[∑k=1n‍1/k-logn]≅0.577216 is the Euler-Mascheroni constant.The classical [23] sampling theorem of Whittaker, Kotelnikov, and Shannon (WKS) for f∈PWσ2 is the series representation as follows:
(12)f(t)=∑n=-∞∞‍f(nπσ)Sn(t),t∈ℝ,
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] f∈PWσ2 if and only if there is g(·)∈L2(-σ,σ) such that
(13)f(t)=12π∫-σσ‍g(x)eixtdx.
Therefore, f′(t)∈PWσ2, that is, f′(t) also has an expansion of the form (12). However, f′(t) can also be obtained by term-by-term differentiation formula of (12) as follows:
(14)f′(t)=∑n=-∞∞‍f(nπσ)Sn′(t),
see [23, page 52] for convergence. Thus, the use of Hermite interpolations will not cost any additional computational efforts since the samples f(nπ/σ) will be used to compute both f(t) and f′(t) according to (12) and (14), respectively.Now, we consider the following differential equations:(15)ℓ(y)≔-y′′(x,μ)+q(x)y(x,μ)=μ2y(x,μ),iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiix∈[0,1],
with the following boundary conditions:
(16)a1y(0,μ)+a2y′(0,μ)=μ2(a1′y(0,μ)+a2′y′(0,μ)),(17)b1y(1,μ)+b2y′(1,μ)=μ2(b1′y(1,μ)+b2′y′(1,μ)),
where μ is a complex spectral parameter, q(·) is assumed to be real valued and continuous on [0,1], and ai,bi,ai′,bi′∈ℝ, i=0,1 satisfying
(18)((a1′,a2′)=(0,0)ora1a2′-a1′a2>0),((b1′,b2′)=(0,0)orb1′b2-b1b2′>0).
The eigenvalue problem (15)–(17) will be denoted by Π(q,a,b,a′,b′) when (a1′,a2′)≠(0,0)≠(b1′,b2′). It is a Sturm-Liouville problem when the eigenparameter μ appears linearly in both boundary conditions. The classical problem when a1′=a2′=b1′=b2′=0, which we denote by Π(q,a,b,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 Π(q,a,b,0,0) numerically. In [36], see also [37], Annaby and Tharwat proved that Π(q,a,b,a′,b′) 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 a1′=a2′=0, (b1′,b2′)≠(0,0) or equivalently when (a1′,a2′)≠(0,0) and b1′=b2′=0. These problems will be denoted by Π(q,a,b,0,b′), Π(q,a,b,a′,0), respectively. The aim of the present work is to compute the eigenvalues of Π(q,a,b,a′,b′), Π(q,a,b,0,b′), and Π(q,a,b,a′,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 Π(q,a,b,a′,b′) and Π(q,a,b,0,b′) with error estimates, respectively. The last section involves some illustrative examples.


## 2. Treatment ofΠ(q,a,b,a′,b′)

In this section, we derive approximate values of the eigenvalues ofΠ(q,a,b,a′,b′). Let y(·,μ) denote the solution of (15) satisfying the following initial conditions:
(19)y(0,μ)=a2-a2′μ2,y′(0,μ)=a1′μ2-a1.
Thus, y(·,μ) satisfies the boundary condition (16). The eigenvalues of the problem Π(q,a,b,a′,b′) are the zeros of the function as follows:
(20)Δ(μ):=(b1′μ2-b1)y(1,μ)+(b2′μ2-b2)y′(1,μ).
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 y(x,μ) satisfies Volterra integral equation as follows:
(21)y(x,μ)=(a2-a2′μ2)cosμx-(a1-a1′μ2)sinμxμ+T[y](x,μ),
where T is the Volterra operator defined by
(22)T[y](x,μ)=∫0x‍sinμ(x-t)μq(t)y(t,μ)dt.
Differentiating (21), we get
(23)y′(x,μ)=(a2′μ2-a2)μsinμx+(a1′μ2-a1)cosμx+T~[y](x,μ),
where T~ is the Volterra operator
(24)T~[y](x,μ)=∫0x‍cosμ(x-t)q(t)y(t,μ)dt.
Define f(·,μ) and g(·,μ) to be
(25)f(x,μ):=T[y](x,μ),g(x,μ):=T~[y](x,μ).
In the following, we will make use of the estimates [44] as follows:
(26)|cosz|≤e|ℑz|,|sinzz|≤c01+|z|e|ℑz|,
where c0 is some constant (we may take c0≃1.72). For convenience, we define the constants by
(27)τ≔∫01‍|q(t)|dt,c1≔|a2|+c0|a1|,c2≔|a2′|+c0|a1′|,c3≔c0τ,c4:=expc3,c5:=max{c1,c2,|b1|+|b2|τ,|b1′|+|b2′|τ}.
From (21) and (25), we get
(28)f(x,μ)=∫0x‍sinμ(x-t)μq(t)[sinμtμ(a2-a2′μ2)cosμt-(a1-a1′μ2)sinμtμ]dt+∫0x‍sinμ(x-t)μq(t)f(t,μ)dt.Lemma 1.
For0≤x≤1, μ∈ℂ, the following estimates hold:
(29)f(x,μ)≤c3c4(c1+c2|μ|2)1+|μ|e|ℑμ|x,(30)g(x,μ)≤τc3c4(c1+c2|μ|2)1+|μ|e|ℑμ|x.Proof.
We dividef(·,μ) into two parts f1(·,μ) and f2(·,μ) and estimate each of them. Indeed, for x∈[0,1] and μ∈ℂ we have
(31)|f1(x,μ)|=|∫0x‍sinμ(x-t)μq(t)[sinμtμ(a2-a2′μ2)cosμt-(a1-a1′μ2)sinμtμ]dt∫0x‍sinμ(x-t)μ|≤e|ℑμ|x∫0x‍|q(t)|c0(x-t)1+|μ|(x-t)×[c0t1+|μ|t|a2|+|a2′||μ|2+(|a1|+|a1′||μ|2)×c0t1+|μ|t]dt≤e|ℑμ|xc0x1+|μ|x∫0x‍|q(t)|[|a2|+|a2′||μ|2+(|a1|+|a1′||μ|2)c0t]dt≤e|ℑμ|xc01+|μ|∫01‍|q(t)|[|a2|+|a2′||μ|2+(|a1|+|a1′||μ|2)c0t]dt.
Moreover, 0≤x≤1, μ∈ℂ,
(32)|f2(x,μ)|=|∫0x‍sinμ(x-t)μq(t)f(t,μ)dt|≤∫0x‍c0(x-t)1+|μ|(x-t)e|ℑμ|(x-t)×|q(t)||f(t,μ)|dt≤c0e|ℑμ|x∫0x‍e-|ℑμ|t|q(t)||f(t,μ)|dt.
Combining (31) and (32), we obtain 0≤x≤1, μ∈ℂ,
(33)|f(x,μ)|≤e|ℑμ|xc01+|μ|∫01‍|q(t)|[|a2|+|a2′||μ|2+(|a1|+|a1′||μ|2)c0t]dt+c0e|ℑμ|x∫0x‍e-|ℑμ|t|q(t)||f(t,μ)|dt.
Applying Gronwall's inequality, cf. for example, [34, page 51], yields μ∈ℂ,
(34)e-|ℑμ|x|f(x,μ)|≤[c01+|μ|∫01|q(t)|[|a2|+|a2′||μ|2+(|a1|+|a1′||μ|2)c0t]dtc01+|μ|∫01|q(t)|[|a2|+|a2′||μ|2]×exp(c0∫0x‍|q(t)|dt)≤[c01+|μ|∫01|q(t)|[|a2|+|a2′||μ|2+(|a1|+|a1′||μ|2)c0t]dtc01+|μ|∫01|q(t)|[|a2|+|a2′||μ|2]×exp(c0∫01‍|q(t)|dt),
from which we get
(35)|f(x,μ)|≤e|ℑμ|x[∑c0[|a2|+|a2′||μ|2+(|a1|+|a1′||μ|2)c0]1+|μ|×∫01|q(t)|dt∑c0[|a2|+|a2′||μ|2+(|a1|+|a1′||μ|2)c0]1+|μ|]exp(c0∫01‍|q(t)|dt)=c3c4(c1+c2|μ|2)1+|μ|e|ℑμ|x.
Then from (25) and (29), we obtain the estimate (30).Now we splitΔ(μ) into two parts via
(36)Δ(μ)=𝒢(μ)+𝒮(μ),
where 𝒢(μ) is known part
(37)𝒢(μ)=(b1′μ2-b1)[(a2-a2′μ2)cosμ-(a1-a1′μ2)sinμμ]+(b2′μ2-b2)[(a2′μ2-a2)μsinμ+(a1′μ2-a1)cosμ],
and 𝒮(μ) is unknown part
(38)𝒮(μ)=(b1′μ2-b1)f(1,μ)+(b2′μ2-b2)g(1,μ).
Then, from Lemma 1, we have the following lemma.Lemma 2.
The function𝒮(μ) is entire in μ and the following estimate holds:
(39)|𝒮(μ)|≤c3c4c5(1+|μ|2)21+|μ|e|ℑμ|.Proof.
Since(40)𝒮(μ)≤(|b1′||μ|2+|b1|)|f(1,μ)|+(|b2′||μ|2+|b2|)|g(1,μ)|,
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 m∈ℤ+, m>4 be fixed. Let ℱθ,m(μ) be the function
(41)ℱθ,m(μ):=(sinθμθμ)m𝒮(μ),μ∈ℂ.
More specifications on m,θ will be given later on. Then we have the next lemma.Lemma 3.
ℱ
θ
,
m
(
μ
) is an entire function of μ which satisfies the estimates
(42)|ℱθ,m(μ)|≤c3c4c5c0m(1+|μ|2)2(1+θ|μ|)m+1e|ℑμ|(1+mθ).
Moreover, μm-4ℱθ,m(μ)∈L2(ℝ) and
(43)Em-4(ℱθ,m)=∫-∞∞‍|μm-4ℱθ,m(μ)|2dμ≤2c3c4c5c0mν0,
where
(44)ν0≔((m(4m2-1)Γ[2m+2]θ2m+1((m(2m-1)+4θ2)Γ[2m+2])-1((m(2m-1)+4θ2)Γ[2m+2]+144m(4m2-1)θ4×(280θ4Γ[2m-7]+20θ2Γ[2m-5]+Γ[2m-3]))×(m(4m2-1)Γ[2m+2]θ2m+1((m(2m-1)+4θ2)Γ[2m+2])-1)1/2.Proof.
Since𝒮(μ) is entire, then also ℱθ,m(μ) is entire in μ. Combining the estimates |sinz/z|≤(c0/(1+|z|))e|ℑz| and (39), we obtain
(45)|ℱθ,m(μ)|≤(c01+θ|μ|)me|ℑμ|mθ·c3c4c5(1+|μ|2)21+|μ|e|ℑμ|,μ∈ℂ,
leading to (42). Therefore, we get
(46)|μm-4ℱθ,m(μ)|≤c3c4c5c0m|μ|m-4(1+|μ|2)2(1+θ|μ|)m+1,μ∈ℝ.
That is, μm-4ℱθ,m(μ)∈L2(ℝ). Moreover, we get
(47)∫-∞∞‍|μm-4ℱθ,m(μ)|2dμ≤c32c42c52c02m×∫-∞∞‍|μ|2m-8(1+|μ|2)4(1+θ|μ|)2m+2dμ=2c32c42c52c02mν02.What we have just proved is thatℱθ,m(μ) belongs to the Paley-Wiener space PWσ2 with σ=1+mθ. Since ℱθ,m(μ)∈PWσ2⊂PW2σ2, then we can reconstruct the functions ℱθ,m(μ) via the following sampling formula:
(48)ℱθ,m(μ)=∑n=-∞∞‍[sin(σμ-nπ)σℱθ,m(nπσ)Sn2(μ)+ℱθ,m′(nπσ)sin(σμ-nπ)σSn(μ)].
Let N∈ℤ+, N>m and approximate ℱθ,m(μ) by its truncated series ℱθ,m,N(μ), where
(49)ℱθ,m,N(μ):=∑n=-NN‍[sin(σμ-nπ)σℱθ,m(nπσ)Sn2(μ)+ℱθ,m′(nπσ)sin(σμ-nπ)σSn(μ)].
Since all eigenvalues are real, then from now on we restrict ourselves to μ∈ℝ. Since μm-4ℱθ,m(μ)∈L2(ℝ), the truncation error, cf. (5), is given for |μ|<Nπ/σ by
(50)|ℱθ,m(μ)-ℱθ,m,N(μ)|≤TN,m-4,σ(μ),
where
(51)TN,m-4,σ(μ)≔ξm-4,σEm-4|sinσμ|23(N+1)m-4(1(Nπ-σμ)3/2+1(Nπ+σμ)3/2)+ξm-4,σ(σEm-4+(m-4)Em-5)|sinσμ|2σ(N+1)m-4×(1Nπ-σμ+1Nπ+σμ).
The samples {ℱθ,m(nπ/σ)}n=-NN and {ℱθ,m′(nπ/σ)}n=-NN, in general, are not known explicitly. So we approximate them by solving numerically 8N+4 initial value problems at the nodes {nπ/σ}n=-NN. Let {ℱ~θ,m(nπ/σ)}n=-NN and let {ℱ~θ,m′(nπ/σ)}n=-NN be the approximations of the samples of {ℱθ,m(nπ/σ)}n=-NN and {ℱθ,m′(nπ/σ)}n=-NN, respectively. Now we define ℱ~θ,m,N(μ), which approximates ℱθ,m,N(μ) as
(52)ℱ~θ,m,N(μ)≔∑n=-NN‍[(nπσ)sin(σμ-nπ)σℱ~θ,m(nπσ)Sn2(μ)+ℱ~θ,m′(nπσ)sin(σμ-nπ)σSn(μ)],iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiN>m.
Using standard methods for solving initial problems, we may assume that for |n|<N,
(53)|ℱθ,m(nπσ)-ℱ~θ,m(nπσ)|<ɛ,|ℱθ,m′(nπσ)-ℱ~θ,m′(nπσ)|<ɛ,
for a sufficiently small ɛ. From (42) we can see that ℱθ,m(μ) satisfies the condition (9) when m>4 and therefore whenever 0<ɛ≤min{π/σ,σ/π,1/e} we have
(54)|ℱθ,m,N(μ)-ℱ~θ,m,N(μ)|≤𝒜(ɛ),μ∈ℝ,
where there is a positive constant Mℱθ,m for which, cf. (10),
(55)𝒜(ɛ)≔2e1/4σ{∫3e(1+σ)+(πσA+Mℱθ,m)ρ(ɛ)+(σ+2+log(2))Mℱθ,m∫}ɛlog(1ɛ).
Here
(56)A:=3σπ(|ℱθ,m(0)|+σπMℱθ,m),ρ(ɛ):=γ+10log(1ɛ).
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,
(57)Δ(μ*)=𝒢(μ*)+(sinθμ*θμ*)-mℱθ,m(μ*)=0.
Then it follows that
(58)𝒢(μ*)+(sinθμ*θμ*)-mℱ~θ,m,N(μ*)=(sinθμ*θμ*)-mℱ~θ,m,N(μ*)-(sinθμ*θμ*)-mℱθ,m(μ*)=[(sinθμ*θμ*)-mℱ~θ,m,N(μ*)-(sinθμ*θμ*)-mℱθ,m,N(μ*)]+[(sinθμ*θμ*)-mℱθ,m,N(μ*)-(sinθμ*θμ*)-mℱθ,m(μ*)]
and so
(59)|𝒢(μ*)+(sinθμ*θμ*)-mℱ~θ,m,N(μ*)|≤|sinθμ*θμ*|-m(TN,m-4,σ(μ*)+𝒜(ɛ)).
Since 𝒢(μ*)+(sinθμ*/θμ*)-mℱ~θ,m,N(μ*) is given and |sinθμ*/θμ*|-m(TN,m-4,σ(μ*)+𝒜(ɛ)) has computable upper bound, we can define an enclosure for μ* by solving the following system of inequalities:
(60)-|sinθμ*θμ*|-m(TN,m-4,σ(μ*)+𝒜(ɛ))≤𝒢(μ*)+(sinθμ*θμ*)-mℱ~θ,m,N(μ*)≤|sinθμ*θμ*|-m(TN,m-4,σ(μ*)+𝒜(ɛ)).
Its solution is an interval containing μ*, and over which the graph 𝒢(μ*)+(sinθμ*/θμ*)-mℱ~θ,m,N(μ*) is squeezed between the graphs as follows:
(61)-|sinθμ*θμ*|-m(TN,m-4,σ(μ*)+𝒜(ɛ)),|sinθμ*θμ*|-m(TN,m-4,σ(μ*)+𝒜(ɛ)).
Using the fact that
(62)ℱ~θ,m,N(μ)⟶ℱθ,m(μ)
uniformly over any compact set and since μ* is a simple root, we obtain the following for large N and sufficiently small ɛ:
(63)∂∂μ(𝒢(μ)+(sinθμθμ)-mℱ~θ,m,N(μ))≠0
in a neighborhood of μ*. Hence, the graph of 𝒢(μ)+(sinθμ/θμ)-mℱ~θ,m,N(μ) intersects the graphs -|sinθμ/θμ|-m(TN,m-4,σ(μ)+𝒜(ɛ)) and |sinθμ/θμ|-m(TN,m-4,σ(μ)+𝒜(ɛ)) at two points with abscissae a-(μ*,N,ɛ)≤a+(μ*,N,ɛ) and the solution of the system of inequalities (60) is the interval
(64)Iɛ,N:=[a-(μ*,N,ɛ),a+(μ*,N,ɛ)]
and in particular μ*∈Iɛ,N. 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 N0∈ℤ+ and sufficiently small ɛ such that μ*∈Iɛ,N  for N>N0. Moreover, we get
(65)[a-(μ*,N,ɛ),a+(μ*,N,ɛ)]⟶{μ*}iiiiiiiiiiiiiiiiiiiiasN⟶∞,ɛ⟶0.Proof.
Since all eigenvalues ofΠ(q,a,b,a′,b′) are simple, then for large N and sufficiently small ɛ we have (∂/∂μ)(𝒢(μ)+(sinθμ/θμ)-mℱ~θ,m,N(μ))>0, in a neighborhood of μ*. Choose N0 such that
(66)𝒢(μ)+(sinθμθμ)-mℱ~θ,m,N0(μ)=±|sinθμθμ|-m(TN0,m-4,σ(μ)+𝒜(ɛ))
has two distinct solutions which we denote by a-(μ*,N0,ɛ)≤a+(μ*,N0,ɛ). The decay of TN,m-4,σ(μ)→0 as N→∞ and 𝒜(ɛ)→0 as ɛ→0 will ensure the existence of the solutions a-(μ*,N,ɛ) and a+(μ*,N,ɛ) as N→∞ and ɛ→0. For the second point we recall that ℱ~θ,m,N(μ)→ℱθ,m(μ) as N→∞ and as ɛ→0. Hence, by taking the limit we obtain
(67)𝒢(a+(μ*,∞,0))+(sinθμ*θμ*)-mℱθ,m(a+(μ*,∞,0))=0,𝒢(a-(μ*,∞,0))+(sinθμ*θμ*)-mℱθ,m(a-(μ*,∞,0))=0.
That is, Δ(a+)=Δ(a-)=0. This leads us to conclude that a+=a-=μ*, since μ* is a simple root.LetΔ~N(μ):=𝒢(μ)+(sinθμ/θμ)-mℱ~θ,m,N(μ). Then (50) and (54) imply
(68)|Δ(μ)-Δ~N(μ)|≤|sinθμθμ|-m(TN,m-4,σ(μ)+𝒜(ɛ)),iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii|μ|<Nπσ
and θ is chosen sufficiently small for which |θμ|<π. Therefore, θ,m must be chosen so that for |μ|<Nπ/σ(69)m>4,θ∈(0,1),|θμ|<π.
Let μ* be an eigenvalue and let μN be its approximation. Thus, Δ(μ*)=0 and Δ~N(μN)=0. From (68) we have |Δ~N(μ*)|≤|sinθμ*/θμ*|-m(TN,m-4,σ(μ*)+𝒜(ɛ)). Now we estimate the error |μ*-μN| for an eigenvalue μ*.Theorem 5.
Letμ*2 be an eigenvalue of Π(q,a,b,a′,b′). For sufficient large N we have the following estimate:
(70)|μ*-μN|<|sinθμNθμN|-mTN,m-4,σ(μN)+𝒜(ɛ)infζ∈Iɛ,N|Δ′(ζ)|.Proof.
SinceΔ(μN)-Δ~N(μN)=Δ(μN)-Δ(μ*), then from (68) and after replacing μ by μN, we obtain
(71)|Δ(μN)-Δ(μ*)|≤|sinθμNθμN|-m(TN,m-4,σ(μN)+𝒜(ɛ)).
Using the mean value theorem yields that for some ζ∈Jɛ,N:=[min(μ*,μN),max(μ*,μN)],
(72)|(μ*-μN)Δ′(ζ)|≤|sinθμNθμN|-m(TN,m-4,σ(μN)+𝒜(ɛ)),ζ∈Jɛ,N⊂Iɛ,N.
Since the eigenvalues are simple, then for sufficiently large Ninfζ∈Iɛ,N|Δ′(ζ)|>0 and we get (70).


## 3. The Case ofΠ(q,a,b,0,b′)

This section includes briefly a treatment similarly to that of the previous section for the eigenvalue problemΠ(q,a,b,0,b′) introduced in Section 1. Notice that condition (18) implies that the analysis of problem Π(q,a,b,0,b′) is not included in that of Π(q,a,b,a′,b′). Let ψ(·,μ) denote the solution of (15) satisfying the following initial conditions:
(73)ψ(0,μ)=a2,ψ′(0,μ)=-a1.
Thus, ψ(·,μ) satisfies the boundary condition (16). The eigenvalues of the problem Π(q,a,b,0,b′) are the zeros of the function as follows:
(74)Ω(μ):=(b1′μ2-b1)ψ(1,μ)+(b2′μ2-b2)ψ′(1,μ).
Recall that Π(q,a,b,a′,b′) has denumerable set of real and simple eigenvalues, cf. [38]. Using the method of variation of constants, the solution ψ(x,μ) satisfies Volterra integral equation as follows:
(75)ψ(x,μ)=a2cosμx-a1sinμxμ+T[ψ](x,μ),
where T is the Volterra operator defined in (22). Differentiating (75), we get
(76)ψ′(x,μ)=-a2μsinμx-a1cosμx+T~[ψ](x,μ),
where T~ is the Volterra operator defined in (24). Define h1(·,μ) and h2(·,μ) to be
(77)h1(x,μ):=T[ψ](x,μ),h2(x,μ):=T~[ψ](x,μ).
As in the preceding section we split Ω(μ) into
(78)Ω(μ):=𝒦(μ)+𝒰(μ),
where 𝒦(μ) is the known part
(79)𝒦(μ)=(b1′μ2-b1)[a2cosμ-a1sinμμ]+(b2′μ2-b2)[-a2μsinμ-a1cosμ],
and 𝒰(μ) is the unknown one
(80)𝒰(μ):=(b1′μ2-b1)h1(1,μ)+(b2′μ2-b2)h2(1,μ).
Then, as in the previous section, 𝒰(μ) is entire in μ for each x∈[0,1] for which
(81)|𝒰(μ)|≤c1c3c4c6(1+|μ|2)1+|μ|e|ℑμ|,μ∈ℂ,
where c6:=max{|b1|+|b2|τ,|b1′|+|b2′|τ}.Letθ∈(0,1) and let m be as in the previous section, but m>2. Define ℛm,θ(μ) to be
(82)ℛm,θ(μ)=(sinθμθμ)m𝒰(μ),μ∈ℂ.
Hence,
(83)|ℛm,θ(μ)|≤c0mc1c3c4c6(1+|μ|2)(1+θ|μ|)m+1e|ℑμ|(1+mθ),μ∈ℂ,
and μm-2ℛm,θ(μ)∈L2(ℝ) with
(84)Em-2(ℛm,θ)=∫-∞∞‍|μm-2ℛm,θ(μ)|2dμ≤2c0mc1c3c4c6ω0,
where
(85)ω0:=12θ2Γ[2m-3]+Γ[2m-1]θ2m-1Γ[2(m+1)].
Thus, ℛm,θ(μ) belongs to the Paley-Wiener space PWσ2 with σ=1+mθ. Since ℛθ,m(μ)∈PWσ2⊂PW2σ2, then we can reconstruct the functions ℛθ,m(μ) via the following sampling formula:
(86)ℛθ,m(μ)=∑n=-∞∞‍[sin(σμ-nπ)σℛθ,m(nπσ)Sn2(μ)+ℛθ,m′(nπσ)sin(σμ-nπ)σSn(μ)].
Let N∈ℤ+, N>m, and approximate ℛθ,m(μ) by its truncated series ℛθ,m,N(μ), where
(87)ℛθ,m,N(μ)≔∑n=-NN‍[sin(σμ-nπ)σsin(σμ-nπ)σℛθ,m(nπσ)Sn2(μ)+ℛθ,m′(nπσ)sin(σμ-nπ)σSn(μ)].
Since all eigenvalues are real, then from now on we restrict ourselves to μ∈ℝ. Since μm-2ℛθ,m(μ)∈L2(ℝ), the truncation error, cf. (5), is given for |μ|<Nπ/σ by
(88)|ℛθ,m(μ)-ℛθ,m,N(μ)|≤TN,m-2,σ(μ),
where
(89)TN,m-2,σ(μ)≔ξm-2,σEm-2|sinσμ|23(N+1)m-2(1(Nπ-σμ)3/2+1(Nπ+σμ)3/2)+ξm-2,σ(σEm-2+(m-2)Em-3)|sinσμ|2σ(N+1)m-2×(1Nπ-σμ+1Nπ+σμ).
The samples {ℛθ,m(nπ/σ)}n=-NN and {ℛθ,m′(nπ/σ)}n=-NN, in general, are not known explicitly. So we approximate them by solving numerically 4N+2 initial value problems at the nodes {nπ/σ}n=-NN. Let {ℛ~θ,m(nπ/σ)}n=-NN and let {ℛ~θ,m′(nπ/σ)}n=-NN be the approximations of the samples of {ℛθ,m(nπ/σ)}n=-NN and {ℛθ,m′(nπ/σ)}n=-NN, respectively. Now we define ℛ~θ,m,N(μ), which approximates ℛθ,m,N(μ)(90)ℛ~θ,m,N(μ)≔∑n=-NN‍[sin(σμ-nπ)σℛ~θ,m(nπσ)Sn2(μ)+ℛ~θ,m′(nπσ)sin(σμ-nπ)σSn(μ)],iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiN>m.Using standard methods for solving initial problems, we may assume that for|n|<N(91)|ℛθ,m(nπσ)-ℛ~θ,m(nπσ)|<ɛ,|ℛθ,m′(nπσ)-ℛ~θ,m′(nπσ)|<ɛ,
for a sufficiently small ɛ. From (83) we can see that ℛθ,m(μ) satisfies the condition (9) when m>2 and therefore whenever 0<ɛ≤min{π/σ,σ/π,1/e} we have
(92)|ℛθ,m,N(μ)-ℛ~θ,m,N(μ)|≤𝒜(ɛ),μ∈ℝ,
where there is a positive constant Mℛθ,m for which, cf. (10), and
(93)𝒜(ɛ)≔2e1/4σ{3e(1+σ)+(πσA+Mℛθ,m)ρ(ɛ)+(σ+2+log(2))Mℛθ,m3e(1+σ)+(πσA+Mℛθ,m)}ɛlog(1ɛ).
Here
(94)A:=3σπ(|ℛθ,m(0)|+σπMℛθ,m),ρ(ɛ):=γ+10log(1ɛ).
As in the above section, we have the following lemma.Lemma 6.
For any eigenvalueμ*2 of the problem Π(q,a,b,0,b′), we can find N0∈ℤ+ and sufficiently small ɛ such that μ*∈ℐɛ,N for N>N0, where
(95)ℐɛ,N:=[b-(μ*,N,ɛ),ab+(μ*,N,ɛ)],b-,b+ are the solutions of the inequalities
(96)-|sinθμθμ|-m(TN,m-2,σ(μ)+𝒜(ɛ))≤Ω~N(μ)≤|sinθμθμ|-m(TN,m-2,σ(μ)+𝒜(ɛ)).
Moreover, we get
(97)[b-(μ*,N,ɛ),b+(μ*,N,ɛ)]⟶{μ*}asN⟶∞,ɛ⟶0.LetΩ~N(μ):=𝒦(μ)+(sinθμ/θμ)-mℛ~θ,m,N(μ). Then (88) and (92) imply
(98)|Ω(μ)-Ω~N(μ)|≤|sinθμθμ|-m(TN,m-2,σ(μ)+𝒜(ɛ)),|μ|<Nπσ
and θ is chosen sufficiently small for which |θμ|<π. Therefore, θ,m must be chosen so that for |μ|<Nπ/σ(99)m>2,θ∈(0,1),|θμ|<π.
Let μ* be an eigenvalue and μN be its approximation. Thus Ω(μ*)=0 and Ω~N(μN)=0. From (98) we have |Ω~N(μ*)|≤|sinθμ*/θμ*|-m(TN,m-2,σ(μ*)+𝒜(ɛ)). Now we estimate the error |μ*-μN| for an eigenvalue μ*. Finally we have the following estimate.Theorem 7.
Letμ*2 be an eigenvalue of the problem Π(q,a,b,0,b′). For sufficient large N we have the following estimate
(100)|μ*-μN|<|sinθμNθμN|-mTN,m-2,σ(μN)+𝒜(ɛ)infζ∈ℐɛ,N|Ω′(ζ)|.


## 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 m 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 a±(μ) and b±(μ) are defined by
(101)a±(μ)=Δ~N(μ)±|sinθμθμ|-m(TN,m-3,σ(μ)+𝒜(ɛ)),iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii|μ|<Nπσ,(102)b±(μ)=Ω~N(μ)±|sinθμθμ|-m(TN,m-2,σ(μ)+𝒜(ɛ)),iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii|μ|<Nπσ,
respectively. Recall also that the enclosure intervals Iɛ,N:=[a-,a+] and ℐɛ,N:=[b-,b+] are determined by solving
(103)a±(μ)=0,|μ|<Nπσ,(104)b±(μ)=0,|μ|<Nπσ.
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]
(105)-y′′(x,μ)-y(x,μ)=μ2y(x,μ),0≤x≤1,y′(0,μ)=0,y(1,μ)+y′(1,μ)=μ2(2y(1,μ)+y′(1,μ)),
is a special case of the problem treated in the previous section with a1=0, a2=b1=b2=b2′=1, b1′=2 and q(x)=-1. The characteristic function is
(106)Ω(μ)=(-1+2μ2)cos[1+μ2]-(-1+μ2)1+μ2sin[1+μ2].
The function𝒦(μ) will be
(107)𝒦(μ)=(-1+2μ2)cos[μ]-μ(-1+μ2)sin[μ].
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 m, θ and ɛ. By exact we mean the zeros of Ω(μ) computed by Mathematica.
Figures1 and 2 illustrate the enclosure intervals dominating μ1 for N=15, m=10, θ=1/5, and ɛ=10-5 and ɛ=10-10 respectively. The middle curve represents Ω(μ), while the upper and lower curves represent the curves of b+(μ),b-(μ), respectively. We notice that when ɛ=10-10, the two curves are almost identical. Similarly, Figures 3 and 4 illustrate the enclosure intervals dominating μ2 for N=15, m=10, θ=1/5, and ɛ=10-5andɛ=10-10 respectively.Table 1
WithN=15, the approximation μk,N and the exact solution μk for different choices of m and θ.
μ
k
μ
1
μ
2
μ
3
μ
4
Exactμk
1.0493258679653497
3.5207214555369464
6.505146961583527
9.578576417519093
μ
k
,
N
m=6
θ = 1/9
1.0493258487568435
3.5207214557864277
6.5051469611825405
9.578576417536349
θ = 1/12
1.0493291624256957
3.5207212860735546
6.5051469673141336
9.578576393056661
m=10
θ = 1/5
1.0493258679653554
3.5207214555369624
6.505146961583516
9.57857641751921
θ = 1/8
1.049325865990155
3.520721455899672
6.505146961464727
9.578576417559768Table 2
Absolute error|μk-μk,N|.
μ
k
μ
1
μ
2
μ
3
μ
4
m
=
6
θ = 1/9
1.92085 × 10−8
2.49481 × 10−10
4.00987 × 10−10
1.72555 × 10−11
θ = 1/12
3.29446 × 10−6
1.69463 × 10−7
5.73061 × 10−9
2.44624 × 10−8
m
=
10
θ = 1/5
5.77316 × 10−15
1.59872 × 10−14
1.15463 × 10−14
1.1724 × 10−13
θ = 1/8
1.97519 × 10−9
3.62725 × 10−10
1.188 × 10−10
4.0675 × 10−11Table 3
ForN=15, m=10, and θ=1/5, the exact solution μk are all inside the interval [b-, b+] for different values of ε.
μ
k
μ
1
μ
2
μ
3
μ
4
Exactμk
1.0493258679653497
3.5207214555369464
6.505146961583527
9.578576417519093
ℐ
ε
,
N, ε=10-5
[1.04294069, 1.05557896]
[3.51981844, 3.52162396]
[6.50375768, 6.50653831]
[9.55222305, 9.60712093]
ℐ
ε
,
N, ε=10-10
[1.04932561, 1.049326118]
[3.52072141, 3.52072149]
[6.50514690, 6.50514702]
[9.57857529, 9.57857754]
E
8
(
ℛ
θ
,
m
)
=
4
.
51845
×
10
8, E7(ℛθ,m)=2.29709×105, ν=1, Mℛθ,m=4.55609×104.Figure 1
The enclosure interval dominatingμ1 for m=10, N=15, θ=1/5, and ɛ=10-5.Figure 2
The enclosure interval dominatingμ1 for m=10, N=15, θ=1/5, and ɛ=10-10.Figure 3
The enclosure interval dominatingμ2 for m=10, N=15, θ=1/5, and ɛ=10-5.Figure 4
The enclosure interval dominatingμ2 for m=10, N=15, θ=1/5, and ɛ=10-10.Example 2.
The boundary value problem [45]
(108)-y′′(x,μ)-y(x,μ)=μ2y(x,μ)0≤x≤1,y(0,μ)=μ2y′(0,μ),y′(1,μ)=μ2y(1,μ),
is a special case of problem Π(q,a,b,a′,b′), when q(x)=-1, a2=a1′=b1=b2′=0, and a1=a2′=b1′=b2=1. Here the characteristic function is
(109)Δ(μ)=(1-μ4)cosμ2+1-(2μ2+μ4)sinμ2+1μ2+1.
After computing 𝒢(μ), we obtain
(110)𝒢(μ)=(1+μ2)((1-μ2)cosμ-μsinμ).
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
WithN=40, the approximation μk,N and the exact solution μk for different choices of θ.
μ
k
μ
1
μ
2
μ
3
μ
4
Exactμk
0.4828692021748484
1.966318052350425
4.827089429919572
7.919684444168381
μ
k
,
N
m=8
θ = 1/32
0.48286920221045176
1.96631805234574
4.827089429919605
7.919684444168366
θ = 1/35
0.4828692337692527
1.966318047624416
4.8270894299720776
7.91968444416245Table 5
Absolute error|μk-μk,N|.
μ
k
μ
1
μ
2
μ
3
μ
4
m
=
8
θ = 1/32
3.56034 × 10−11
4.68492 × 10−12
3.28626 × 10−14
1.5099 × 10−14
θ = 1/35
3.15944 × 10−8
4.72601 × 10−9
5.25047 × 10−11
5.93126 × 10−12Table 6
The approximationμk,N and the exact solution μk for N=40, m=14 and θ=1/26.
μ
k
exactμk
μ
k
,
N
|
μ
k
-
μ
k
,
N
|
μ
1
0.4828692021748484698568637
0.4828692021748484678442680
2.012596 × 10−18
μ
2
1.966318052350424642326091
1.9663180523504246423320204
5.93 × 10−21
μ
3
4.8270894299195722717631337
4.8270894299195722717463715
1.6762 × 10−20
μ
4
7.9196844441683813942255769
7.9196844441683813942260057
4.29 × 10−22Table 7
ForN=40, m=14 and θ=1/26, the exact solution μk are all inside the interval [a-, a+] for different values of ε.
μ
k
μ
1
μ
2
μ
3
μ
4
Exactμk
0.4828692021748484
1.966318052350425
4.827089429919572
7.919684444168381
ℐ
ε
,
N,, ε = 10−5
[0.47918888, 0.48651557]
[1.96592879, 1.96670680]
[4.82707252, 4.82710633]
[7.91968171, 7.919687175]
ℐ
ε
,
N, ε = 10−10
[0.48284084, 0.48289756]
[1.96631794, 1.96631815]
[4.82708919, 4.82708966]
[7.919684437, 7.919684450]
E
10
(
ℱ
θ
,
m
)
=
2
.
83057
×
10
18, E9(ℱθ,m)=1.12829×1014, ν=1, Mℱθ,m=1.57716×107.Figure 5
The enclosure interval dominatingμ1 for m=14, N=40, θ=1/26, and ɛ=10-5.Figure 6
The enclosure interval dominatingμ1 for m=14, N=40, θ=1/26, and ɛ=10-10.Figure 7
The enclosure interval dominatingμ4 for m=14, N=40, θ=1/26, and ɛ=10-5.Figure 8
The enclosure interval dominatingμ4 for m=14, N=40, θ=1/26, and ɛ=10-10.



---

*Source: 412028-2013-11-07.xml*