Centrosymmetric Matrices in the Sinc Collocation Method for Sturm-Liouville Problems

Recently, we used the Sinc collocation method with the double exponential transformation to compute eigenvalues for singular Sturm-Liouville problems. In this work, we show that the computation complexity of the eigenvalues of such a differential eigenvalue problem can be considerably reduced when its operator commutes with the parity operator. In this case, the matrices resulting from the Sinc collocation method are centrosymmetric. Utilizing well known properties of centrosymmetric matrices, we transform the problem of solving one large eigensystem into solving two smaller eigensystems. We show that only 1/(N+1) of all components need to be computed and stored in order to obtain all eigenvalues, where (2N+1) corresponds to the dimension of the eigensystem. We applied our result to the Schr\"odinger equation with the anharmonic potential and the numerical results section clearly illustrates the substantial gain in efficiency and accuracy when using the proposed algorithm.


Introduction
In science and engineering, differential eigenvalue problems occur abundantly. Differential eigenvalue problems can arise when partial differential equations are solved using the method of separation of variables. Consequently, they also play an important role in Sturm-Liouville (SL) differential eigenvalue problems [1]. For example, the solution of the wave equation can be expressed as the sum of standing waves. The frequencies of these standing waves are precisely the eigenvalues of its corresponding Sturm-Liouville problem. Similarly, in quantum mechanics, the energy eigenvalues associated with a Hamiltonian operator are modelled using the time-independent Schrödinger equation which is in fact a special case of a Sturm-Liouville differential eigenvalue problem.
Recently, collocation and spectral methods have shown great promise for solving singular Sturm-Liouville differential eigenvalue problems [2,3]. More specifically, the Sinc collocation method (SCM) [4][5][6] has been shown to yield exponential convergence. During the last three decades the SCM has been used extensively to solve many problems in numerical analysis. The applications include numerical integration, linear and non-linear ordinary differential equations, partial differential equations, interpolation and approximations to functions [7,8]. The SCM applied to Sturm-Liouville problems consists of expanding the solution of a SL problem using a basis of Sinc functions. By evaluating the resulting approximation at the Sinc collocation points separated by a fixed mesh size h, one obtains a matrix eigenvalue problem or generalized matrix eigenvalue problem for which the eigenvalues are approximations to the eigenvalues of the SL operator. In [9], we used the double exponential Sinc collocation method (DESCM) to compute the eigenvalues of singular Sturm-Liouville boundary value problems. The DESCM leads to a generalized eigenvalue problem where the matrices are symmetric and positive-definite. In addition, we demonstrate that the convergence of the DESCM is of the rate O N 5/2 log(N ) 2 e −κN/ log(N ) for some κ > 0 as N → ∞, where 2N + 1 is the dimension of the resulting generalized eigenvalue system. The DESCM was also applied successfully to the Schrödinger equation with the anharmonic oscillators [10].
Using the eigenspectrum properties of symmetric centrosymmetric matrices presented in [12], we apply the DESCM algorithm to Sturm-Liouville eigenvalue problems and demonstrate that solving the resulting generalized eigensystem of dimension (2N +1)×(2N +1) is equivalent to solving the two smaller eigensystems of dimension N × N and (N + 1) × (N + 1). Moreover, we also demonstrate that only 1 N +1 of all components need to be stored at every iteration in order to obtain all generalized eigenvalues. To illustrate the gain in efficiency obtained by this method, we apply the DESCM method to the time independent Schrödinger equation with an anharmonic potential. Furthermore, it is worth mentioning that research concerning inverse eigenvalue problems where the matrices are assumed centrosymmetric has been the subject of much research recently [23,28]. Consequently, the combination of these results and our findings could lead to a general approach for solving inverse Sturm-Liouville problems.
All calculations are performed using the programming language Julia and all the codes are available upon request.

Definitions and basic properties
The sinc function valid for all z ∈ C is defined by the following expression: For j ∈ Z and h a positive number, we define the Sinc function S(j, h)(x) by: The Sinc function defined in (2) form an interpolatory set of functions with the discrete orthogonality property: where δ j,k is the Kronecker delta function.
Definition 2.1. [7] Given any function v defined everywhere on the real line and any h > 0, the symmetric truncated Sinc expansion of v is defined by the following series: where v j,h = v(jh).
The Sturm-Liouville (SL) equation in Liouville form is defined as follows: where −∞ ≤ a < b ≤ ∞. Moreover, we assume that the function q(x) is non-negative and the weight function ρ(x) is positive. The values λ are known as the eigenvalues of the SL equation.
In [34], we apply the DESCM to obtain an approximation to the eigenvalues λ of equation (5). We initially applied Eggert et al.'s transformation to equation (5) since it was shown that the proposed change of variable results in a symmetric discretized system when using the Sinc collocation method [35]. The proposed change of variable is of the form [35, Defintion 2.1]: where φ −1 (x) is a conformal map of a simply connected domain in the complex plane with boundary points Applying the change of variable (6) into equation (5), one obtains [35]: where:q To implement the double exponential transformation, we use a conformal mapping φ(x) such that the solution to equation (7) decays double exponentially. In other words, we need to find a function φ(x) such that: for some positive constants A, B, γ. Examples of such mappings are given in [34,36].
Applying the SCM method, we obtain the following generalized eigenvalue problem: where the vectors v and C M (v, h) are given by: and µ are approximations of the eigenvalues λ of equation (7). For more details on the application of the SCM, we refer the readers to [9].
As in [37], we let δ (l) j,k denote the l th Sinc differentiation matrix with unit mesh size: The entries A j,k of the (2N + 1) × (2N + 1) matrix A are then given by: and the entries D 2 j,k of the (2N + 1) × (2N + 1) diagonal matrix D 2 are given by: As previously mentioned, Eggert et al.'s transformation leads to the matrices A and D 2 to be symmetric and positive definite. However, as will be illustrated in the next section, given certain parity assumptions, these matrices yield even more symmetry.
3 Centrosymmetric properties of the matrices A and D 2 In this section, we present some properties of the matrix A and D 2 that will be beneficial in the computation of their eigenvalues. The matrices A and D 2 are symmetric positive definite matrices when equation (7) is discretized using the Sinc collocation method. Additionally, given certain parity assumptions on the functions q(x), φ(x) and ρ(x) in equation (7), the matrices A and D 2 will also be centrosymmetric.
where f (x) is a well defined function being acted upon by J .
An operator B is said to commute with parity operator J if it satisfies the following relation: Equivalently, we can say that the the commutator between B and J is zero, that is: where J is an exchange matrix of dimension (2N + 1) × (2N + 1). Writing equation (19) in a component form, we have the following relation: We now present the following Theorem establishing the connection between symmetries of the Sturm-Liouville operator and its resulting matrix approximation.
Theorem 3.5. Let L denote the operator of the transformed Sturm-Liouville problem in equation (7): If the commutator [L, J ] = 0, where J is the parity operator, then the matrices A and D 2 defined by equations (13) and (14)  If φ(x) is an odd function, then φ ′ (x) is even, φ ′′ (x) is odd and φ ′′′ (x) is even. From this and equation (8), it follows thatq(x) is even.
In order to show that the resulting matrices A and D 2 are centrosymmetric, we demonstrate that both these matrices satisfy equation (20). Before doing so, it is important to notice that the l th Sinc differentiation matrices defined in equation (12) have the following symmetric properties: Hence, the l th Sinc differentiation matrices are centrosymmetric if l is even. It is worth noting that when l is odd, the Sinc differentiation matrices are skew-centrosymmetric [27]. Consequently, investigating the form for the components of the matrix A in equation (13), we obtain: Similarly, investigating the form for the components of the matrix D 2 in equation (14), we obtain: Both matrices A and D 2 satisfy equation (20). From this it follows that A and D 2 are centrosymmetric.
Theorem 3.5 illustrates that Sinc basis functions preserve the parity property of the Sturm-Liouville operator when discretized. Hence, when the matrices A and D 2 are symmetric centrosymmetric positive definite matrices, we can utilize these symmetries when solving for their generalized eigenvalues. In [12], Cantoni et al. proved several properties of symmetric centrosymmetric matrices. In the following, we will utilize some of these properties to facilitate our task of obtaining approximations to the generalized eigenvalues of the matrices A and D 2 . The following lemma will demonstrate the internal block structure of symmetric centrosymmetric matrices.
where S, C are matrices of size N × N , J is the exchange matrix of size N × N , x is a column vector of length N and h is a scalar. In addition, S T = S and C T = JCJ.
The next lemma simplifies the calculation needed to solve for these eigenvalues.
Since our problem consists of solving a generalized eigenvalue problem where one matrix is a full symmetric centrosymmetric and the other is a diagonal centrosymmetric matrix, we propose the following Theorem.
Theorem 3.9. Let H and W be square symmetric centrosymmetric matrices of the same size, such that: then solving the generalized eigenvalue problem det(H − λW) = 0 is equivalent to solving the two smaller generalized eigenvalue problems: Proof This proof relies on the unitary transformation matrix presented in [12, Lemma 3]: where I is the identity matrix and J is the exchange matrix.
It is easy to verify that: where V is the matrix in lemma 3.7.
This result is analogous for the matrix W with a change in notation. Hence: from which the result follows.
Theorem 3.9 is very useful when N is large since it is less costly to solve two symmetric generalized eigensystems of dimensions N × N and (N + 1) × (N + 1) rather than one symmetric eigensystem of dimension (2N + 1) × (2N + 1). Additionally, Lemma 3.6 also has large ramifications when it comes to saving storage space. As is discussed in [37], the l th Sinc differentiation matrices are symmetric toeplitz matrices. Therefore, for a symmetric toeplitz matrix of dimension (2N + 1) × (2N + 1), only 2N + 1 elements need to be stored. Investigating the definition of the matrix A in equation (13), we can see that A is defined as the sum of a symmetric toeplitz matrix and a diagonal matrix. Moreover, from Lemma 3.6 and Theorem 3.9, using only the antidiagonal and anti-upper triangular half of matrix C, the vector x, the scalar h, the diagonal and lower triangular half of the matrix S, the vector w and the scalar w, we can create all the elements needed to solve for the generalized eigenvalues of the matrices A and D 2 . Hence, the ratio of elements needed to be computed and stored at each iteration N in order to solve for these eigenvalues is given by: Thus, only 1 N +1 of the entries need to be generated and stored at every iteration to obtain all of the generalized eigenvalues.
In the following section, we will illustrate the gain in efficiency of these results by applying the DESCM to the Schrödinger equation with an anharmonic oscillator.

The anharmonic oscillator
The time independent Schrödinger equation given by: In equation (34), the Hamiltonian is given by the following linear operator: where V (x) is the potential energy function and E is the energy eigenvalue of the hamiltonian operator H. In our case, we are treating the anharmonic oscillator potential V (x) defined by: In [10], we successfully applied the DESCM to time independent Schrödinger equation with an anharmonic potential. As we can see, the time independent Schrödinger equation (34) where E are approximations of the energy eigenvalues E.
The matrices A and D 2 defined by equation (13) and (14) are given by: where:Ṽ and: Since the anharmonic potential V (x) defined in (35) is an even function, the function ρ(x) = 1 is an even function and the conformal map φ(x) = sinh(x) is an odd function, we know that Theorem 3.5 applies. Hence, the matrices A and D 2 are symmetric centrosymmetric.

Numerical Discussion
In this section, we test the computational efficiency of the results obtained in Theorem 3.9. All calculations are performed using the programming language Julia in double precision. The eigenvalue solvers in Julia use the linear algebra package LAPACK.
In [40], Chaudhuri et al. presented several potentials with known analytic solutions for energy levels calculated using supersymmetric quantum mechanics, namely: Figure 1 presents the absolute error between our approximation and the exact values given in (40). The absolute error is defined by: The optimal mesh size obtained in [10]: where W (x) is the Lambert W function, is used in the calculation.
As can be seen from Figure 1, using the centrosymmetric property improves the convergence rate of the DESCM significantly.

Conclusion
Sturm-Liouville eigenvalue problems are abundant in scientific and engineering problems. In certain applications, these problems possess a symmetry structure which results in the Sturm-Liouville operator to be commutative with the parity operator. As was proven in Theorem 3.5, applying the DESCM will preserve this symmetry and results in a generalized eigenvalue problem where the matrices are symmetric centrosymmetric. The centrosymmetric property leads to a substantial reduction in the computational cost when computing the eigenvalues by splitting the original eigenvalue problem of dimension (2N + 1) × (2N + 1) into two smaller generalized eigensystems of dimension N × N and (N + 1) × (N + 1). Moreover, due to the internal block structure of the matrices obtained using the DESCM, we have shown that only 1 N + 1 of all entries need to be computed and stored at every iteration in order to find all of their eigenvalues. Numerical results are presented for the time independent Schrödinger equation (34) with an anharmonic oscillator potential (35). Four exact potentials with known eigenvalues are tested and the results clearly demonstrated the reduction in complexity and increase in convergence.

Tables and Figures
Absolute Error for energy level 0