COMPUTING EIGENPAIRS OF TWO-PARAMETER STURM-LIOUVILLE SYSTEMS USING THE BIVARIATE SINC-GAUSS FORMULA

. The use of sampling methods in computing eigenpairs of two-parameter boundary value problems is extremely rare. As far as we know, there are only two studies up to now using the bivariate version of the classical and regularized sampling series. These series have a slow convergence rate. In this paper, we use the bivariate sinc-Gauss sampling formula that was pro-posed in [6] to construct a new sampling method to compute eigenpairs of a two-parameter Sturm-Liouville system. The convergence rate of this method will be of exponential order, i.e. O (e − δN / √ N ) where δ is a positive number and N is the number of terms in the bivariate sinc-Gaussian formula. We estimate the amplitude error associated to this formula, which gives us the possibility to establish the rigorous error analysis of this method. Numerical illustrative examples are presented to demonstrate our method in comparison with the results of the bivariate classical sampling method.


1.
Introduction. The issue of computing eigenvalues of one-parameter eigenvalue problems using sinc methods has attracted many researchers. During the period 1996-2018, six sampling methods have been developed to compute the eigenvalues of boundary value problems of various types. These are the classical sinc (1996) [13], the regularized sinc (2005) [14], the sinc-Gaussian (2008) [2], the Hermite (2012) [3], the Hermite-Gauss (2016) [4], and the generalized sinc-Gaussian method (2018) [7]. The use of sampling methods in computing eigenpairs of two-parameter boundary value problems is extremely rare. As far as we know, there are only two studies, cf. [1,15]. In [1], the authors used the bivariate classical sampling series of Whittaker, Kotelnikov and Shannon (WKS) to find a representation for the eigencurves of a two-parameter Sturm-Liouville eigenvalue system. The convergence rate of this method is of order O(ln(N )/ √ N . The authors of [15] used the regularized sampling method to compute the eigenpairs of a two-parameter Sturm-Liouville eigenvalue problem with three-point boundary conditions, without studying the error analysis.
The convergence rate of the regularized sampling series is not much better than that of the classical sampling method and it still is of polynomial order.
By a theorem on analytic parameter dependence, the function is an entire function in λ and µ, cf. e.g. [11]. The couples (λ 2 n , µ 2 n ) are the eigenpairs of system (1.1)-(1.3) if and only if (λ n , µ n ) are the common zeros of the equations (1.6), [8, p. 106]. Those zeros cannot be computed exactly except for extremely rare cases. If {y r (·, λ n , µ n )} r=1,2 is a corresponding set of simultaneous solutions of (1.1)-(1.3), then the product 2 r=1 y r (·, λ n , µ n ) is an eigenfunction of this system corresponding to the eigenpair (λ 2 n , µ 2 n ), and it is unique up to a multiplicative factor. Under the condition (1.5), the set of eigenfunctions is complete in L 2 [0, b] 2 with respect to the weight function δ(x 1 , x 2 ), cf. [8,Theorem 10.6.1].
The rest of the paper has been organized as follows: The next section is devoted to briefly describe the bivariate sinc-Gauss sampling formula. Since alternative samples will be used in our sampling formula, the amplitude error appears in our method. For this reason, we will derive estimates for the amplitude error associated with the bivariate sinc-Gauss sampling formula, which gives us the possibility to establish the rigorous error analysis of this method. The method is constructed in Section 3 and 4, where we also establish a rigorous error analysis associated with this method. Section 5 deals with illustrative examples and comparisons. Lastly, Section 6 concludes the paper.
2. Bivariate sinc-Gauss sampling formula. This section is devoted to describe briefly the bivariate sinc-Gauss sampling formula and investigate the amplitude error associated with it. Let γ = (γ 1 , γ 2 ), γ r > 0, r = 1, 2. The Bernstein space B γ,∞ (R 2 ) is the class of entire functions of two variables satisfying the following growth condition which belong to L ∞ (R 2 ) when restricted to R 2 . For h r ∈ (0, π/γ r ], set δ r := (π − γ r h r )/2 and let E 2 be the class of all entire functions on C 2 . For the Bernstein space B 2 γ,∞ (R 2 ), in [6] we defined a localization operator G h,N : The sinc function is defined by This operator is generalized in [5] including samples from the function and its partial derivatives. Let us mention here that formula (2.1) is defined for wider classes than the Bernstein space B γ,∞ (R 2 ), cf. [6], but here this space is sufficient to treat our problem. There, we estimated the absolute error |f (z) − G h,N [f ](z)| where f ∈ B γ,∞ (R 2 ) and bounds of exponential order were found. Since the eigenpairs of where δ := min{δ 1 , δ 2 } and A δ,N is the real-valued function defined by it does not affect the convergence rate of the error bound in (2.2). Formula (2.1) has a convergence rate of exponential order. In Table 1, we compare between three methods with the same numbers of samples (2N + 1) 2 Thus, the bivariate sinc-Gaussian sampling will give us higher accuracy results when we use it in computing eigenpairs of a two-parameter Sturm-Liouville system with the additional cost that the function is approximated on a smaller domain or using additional samples. It is approximating the function from the Bernstein space B γ,∞ (R 2 ) using only a finite number of samples of the function. However, sometimes these samples cannot be computed explicitly. This is why the amplitude error associated with formula (2.1) appears. This amplitude error arises when the exact values f (k 1 h 1 , k 2 h 2 ) of (2.1) are replaced by close approximate ones. We assume that the approximated samplesf (k 1 h 1 , k 2 h 2 ) are close to the original The amplitude error associated with In the following theorem, we will show a bound for |G h, on the real domain, which will be used in the investigation of the error analysis of our method in Section 4.
Proof. Using condition (2.4), we obtain for Here we have used the fact |sinc πh −1 r x r − k r π | ≤ 1 for all x r ∈ R, and the index Z N (x r ) is defined as follows (2.7) Combining Mills' Ratio inequality, cf. [17], with (2.7) and (2.6) implies (2.5).
3. The method. This section is devoted to the construction of our method. The main concept of this approximation is to split the two simultaneous equations in (1.6) into two parts. The first is known and the second is unknown, but it belongs to the Bernstein space B γ,∞ (R 2 ). Therefore, we approximate the unknown parts of both equations (1.6) using our bivariate sinc-Gauss sampling formula (2.1). We distinguish between the two cases α r = 0 and α r = 0 because the representation of 3.1. The case α r = 0. Let us start with the following result which we will use in the sequel.
Lemma 3.1. If a, b > 0 and λ, µ ∈ C, then we have from which the assertion follows immediately.
To use the Liouville transformation, denote by
In the following theorem, we will show that the unknown part U r belongs to the Bernstein space B γr,∞ (R 2 ). That gives us the possibility to approximate U r via the bivariate sinc-Gauss sampling formula (2.1).
Proof. Applying inequality (3.1) implies It is easy to see that for all τ r ∈ [0, b] where ρ r is defined in (3.2). Combining (3.10) and (3.9), using the fact where γ r1 and γ r2 are defined above. It follows from (3.11) that U r belongs to the Bernstein space B γr,∞ (R 2 ).
Since U r ∈ B γr,∞ (R 2 ), we can approximate it using the bivariate sinc-Gauss sampling formula (2.1), i.e. U r (λ, µ) ≈ G h,N [U r ](λ, µ) where the samples are given by N (λ, µ) and β r ∈ [0, π). Unfortunately, the samples U r (k 1 h 1 , k 2 h 2 ) cannot be determined explicitly in the general case. That is why the amplitude error usually appears. LetŨ r (k 1 h 1 , k 2 h 2 ) be the approximation of the samples U r (k 1 h 1 , k 2 h 2 ) when the solution y r (b, k 1 h 1 , k 2 h 2 ) and its derivative y (b, k 1 h 1 , k 2 h 2 ) are computed numerically at the nodes (k 1 h 1 , k 2 h 2 ), k ∈ Z 2 N (λ, µ). Now, let us define the following interesting functioñ where K r is defined in (3.8) and G h,N [Ũ r ] is the bivariate sinc-Gauss sampling formula (2.1) which is constructed by the approximated samplesŨ r (k 1 h 1 , k 2 h 2 ). The simultaneous functions∆ r,h,N (λ, µ) are determined explicitly and will be very close to the functions ∆ r (λ, µ) which are defined in (1.6), as we will see in the next result. Therefore, the zeros of the simultaneous equations∆ r,h,N (λ, µ) will be very close to the desired zeros, which are precisely the eigenpairs of the system (1.1)-(1.3), of ∆ r (λ, µ).
Since U r,αr ∈ B γr,∞ (R 2 ), we can approximate it using the bivariate sinc-Gauss sampling formula (2.1). Here, we complete the method in the same way as we have done in the last case. Remark 1. If α 1 = 0 and α 2 = 0, the solutions y 1 and y 2 will be as given in (3.4) and (3.14), respectively. The cases α 1 = 0 and α 2 = 0 are similar. In these cases, the method will be completed as indicated above.
Denote by J ∆ the Jacobian matrix If (λ * , µ * ) is a zero of both equations in (1.6), then the determinant |J ∆ (λ * , µ * )| is nonzero, cf. [16,Lemma 4.1]. This fact will be used in the proof of the following theorem.

6.
Conclusions. This work is devoted to constructing a new sampling method to compute eigenpairs of a two-parameter Sturm-Liouville eigenvalue system with separate boundary conditions. This method is built by using the bivariate sinc-Gauss sampling formula which was established by the authors in 2016. Due to the exponential order of convergence, the method described here, i.e. the bivariate sinc-Gauss method, is superior to comparable methods for computing the eigenpairs of the boundary value system. The accuracy of the method increases without additional cost when the parameter N is fixed and h is decreasing, except that the function is approximated on a smaller domain. Examples are provided to illustrate the effectiveness of the approximation. This method can be applied to compute eigenpairs for various types of two-parameter boundary eigenvalue systems. This will be studied in future works.