Accurate Eigenvalues for the Sturm-Liouville Problems, Involving Generalized and Periodic Ones

In the paper, the eigenvalues of Sturm-Liouville problems (SLPs), generalized SLPs and periodic SLPs are solved. First, we propose a new method to transform the SLP with mixed boundary conditions to a generalized SLP for a transformed variable, for which the Dirichlet boundary conditions occur on two-side, but the coe ﬃ cients are nonlinear functions of eigenvalue. To computing the eigenvalue and eigenfunction, we further recast the transformed system to an initial value problem for a new variable. In terms of the relative residual of two consecutive terminal values of the new variable a nonlinear equation is solved for seeking the eigenvalue by the ﬁctitious time integration method (FTIM), which monotonically converges to the exact eigenvalue. We solve a numerically characteristic equation by the half-interval method (HIM) and a derivative-free iterative scheme LHL (Liu, Hong & Li, 2021) to achieve high precision eigenvalues. Next, the generalized SLP is transformed to a new one, so that the Dirichlet boundary condition happens on the right-end. By using the boundary shape function method and the uniqueness condition of the transformed variable, a deﬁnite initial value problem is derived for the new variable. To match the right-end Dirichlet boundary condition a numerically characteristic equation is obtained and solved by the HIM and LHL. Finally, new techniques for solving the periodic SLPs with three types periodic boundary conditions are proposed, which preserve the periodic boundary conditions with the aids of boundary shape functions. Three iterative algorithms are developed, which converge quickly. All the proposed iterative algorithms are identiﬁed by testing some examples.


Introduction
During 1836-37, Sturm and Liouville created a new mathematical analysis of second-order ordinary differential equation (ODE) equipped with the Sturm-Liouville boundary conditions (Lutzen, 1984). It is nowadays known as the Sturm-Liouville theory to deal with the Sturm-Liouville problem (SLP): hu(a) − p(a)u (a) = 0, where p(x) > 0 and w(x) > 0 are strictly positive functions of x ∈ [a, b], and λ is a constant parameter to be determined. Eqs. (1) and (2) is a regular SLP, which allows non-trivial solutions only for certain values of λ, known as the eigenvalues.
The Sturm-Liouville problem is very important in partial differential equations, vibrations in continuum mechanics, heat conduction problems, Schrödinger equation, the transport of micro wave, and the Lax pair of KdV equation, etc. For example, as an application of the Sturm-Liouville problem, describes the longitudinal wave motion of a rod when the Young's modulus E(x) is varying, with the cross section area A(x), a constant density ρ of mass, and the axial displacement u(x, t). Let u(x, t) = e iωt y(x) with ω a natural vibration frequency; Eq. (3) reduces to Eq. (1) with p(x) = E(x)A(x), q(x) = 0, w(x) = ρA(x), λ = ω 2 .
In the design of engineering structure, knowing a first few frequencies ω is utmost important.
Besides some limited cases, there exists no general solution of the Sturm-Liouville problem. The Rayleigh minimal principle was often employed to solve the Sturm-Liouville problem, which is carried out in a large space and leaves a finite number of free parameters to match the minimality condition, but it is effective only for a few lower-order eigenvalues (Gould, 1995;Hinton & Schaefer, 1997;Ghelardoni1997). In general, for higher-order eigenvalues the accuracy is lost very fast.
As to be proved in Theorem 1 below, the SLP with mixed boundary conditions can be transformed to a generalized Sturm-Liouville problem with the coefficients depending on λ nonlinearly. In the existent literature, the three problems of SLP, generalized SLP and periodic SLP were tackled separately. We will propose a unified treatment of these three SLPs by designing different boundary shape functions and transforming them to the corresponding initial value problems, and then some solvers of nonlinear equation are used to solve different target equations for computing the eigenvalues.
We arrange the paper as follows. In Section 2, we describe the SLP which with the Liouville transformation is recast to a canonical form equipped with a single potential function in the unit interval, and then we further transform it to a generalized SLP with the Dirichlet boundary conditions happened on two-side. With the aid of a simple transformation made in Section 3, we derive an initial value problem which together with the fictitious time integration method (FTIM) is used to determine the eigenvalue iteratively, and two numerical examples are given to test the FTIM. In Section 4, we propose a normalization condition for the transformed variable to ascertain the terminal value of the new variable explicitly, such that a definite initial value problem is derived for a new variable. To precisely match the right-end Dirichlet boundary condition, a numerically characteristic equation provides a medium to obtain the eigenvalues, which is solved by the half-interval method (HIM) and a derivative-free iterative scheme (LHL), and very accurate eigenvalues are revealed. In Section 5, we solve the generalized SLP which is transformed to a new generalized SLP with a simplified Dirichlet boundary condition on the right-end. By developing the boundary shape function method and a definite initial value problem method, we then apply the HIM and LHL to solve the numerically characteristic equation. In Section 6, we discuss the periodic SLP by developing the boundary shape function methods. The iterative algorithms are developed in Section 7 to determine the eigenvalues. The conclusions are drawn in Section 8.

A generalized Sturm-Liouville Problem with Two-Side Dirichlet Boundary Conditions
The Liouville transformation (Fu, Wang & Wei, 2015) is given by which are used to transform Eqs. (1) and (2) to the following canonical form of the SLP, endowing with the mixed boundary conditions: where The original SLP becomes a new SLP in the unit interval with the same eigenvalues multiplied by a constant factor c 2 , and only the potential function Q(t) is required.
Therefore, we consider a specific SLP by where q(x, λ) is a linear function of λ. It is different from the Liouville transformation in Eq. (5) that a new one in terms of u(x) and u (x) is derived below, which can transform the above SLP with mixed boundary conditions to a new generalized SLP with the Dirichlet boundary conditions happened on two-side.
Expressing u(x) and u (x) in terms of v(x) and v (x) by using Eqs. (27) and (28), we can obtain the inverse transformations by where D(x, λ) is defined by Eq. (24). Further differentiating Eq.

Transforming to an Initial Value Problem
With the aid of Theorem 1, we derive new methods in Sections 3 and 4 to solve Eqs. (11) and (12), and give examples to test the proposed methods.
Theorem 2 For a new variable y(x) ∈ C 2 [0, 1], the transformed variable v(x) is given by which satisfies the boundary conditions in Eq. (14).

Fictitious Time Integration Method
For solving a nonlinear equation: we can transform it to a first-order differential equation (Liu & Atluri, 2008): where ξ is a fictitious time variable and ν is a constant parameter. Integrating Eq. (34) by the forward Euler scheme, yields where ξ k = k∆ξ, ∆ξ > 0 is an increment of fictitious time, and λ(k) denotes the eigenvalue at the kth iteration. Eq. (35) is the fictitious time integration method (FTIM) developed by Liu and Atluri (2008) to solve nonlinear equation (33). Now, we explore the FTIM (Liu & Atluri, 2008) to determine λ in Eqs. (11) and (12) with the aids of Theorems 1 and 2. It follows from Eqs. (13) and (32) that y(x) replacing v(x), where y(1) is an unknwon constant to be determined by the iteration method. With the given initial values y(0) and y (0), if y can be integarted from Eq. (36), Theorem 2 guarantees that v(x) satisfies the Dirichlet boundary conditions in Eq. (14). Then, by Eqs. (25) and (26), u(x) satisfies the mixed boundary conditions in Eq. (12).
Let λ(k) be the eigenvalue λ and y k (1) be the terminal value of y(1) in Eq. (36) at the kth iteration. We set which is using the relative residual of the iterative sequence y k (1) as a driving force in Eq. (35), so that we can enforce the iterative sequence of λ(k) tending to its limit point as an approximation to the exact eigenvalue λ with F(λ) = 0.
For the SLP endowed with the mixed boundary conditions, we can transform it to the Dirichlet type generalized SLP as shown in Theorem 1. Two examples are given below to demonstrate the usefulness of Theorem 1 and the FTIM.

Example 1
To explore the efficiency and accuracy of the proposed iterative algorithm, we consider where λ k , k = 0, 1, 2, . . . are solved from the characteristic equation: while the kth order eigenfunction is given by With µ 1 = 1, µ 2 = −2, in the iterative algorithm we take y(0) = 1, y (0) = 0, y 0 (1) = 10, 1 = 10 −10 , ∆ξ = 0.001 and N = 500 and iteratively determine λ by Eqs. (35) and (37), whose results are listed in Table 1 to show λ(0), the computed λ, the exact λ, the absolute error of eigenvalue (EE), ν, the number of iterations (NI), the absolute error of right boundary condition (ERBC) and the maximum error (ME) of eigenfunction. Notice that the absolute error of left boundary condition is zero, so that it is not listed in the
We fix y(0) = 1, y (0) = 0, y 0 (1) = 10, 1 = 10 −10 , ∆ξ = 0.001 and N = 500 and iteratively determine λ by Eqs. (35) and (37), whose results are shown in Table 2 for λ(0), the computed λ, the exact λ, EE, ν, NI, and ERBC. The above two examples as shown in Tables 1 and 2 show that the FTIM is convergent very fast with small numbers of iterations. The parameter ∆ξ must be small enough to avoid over-shooting, while the value of the parameter ν can be chosen as large as possible to accelerate the convergence speed. However, the accuracy of eigenvalues can be further improved by other methods given below.

Terminal Value y(1) of the New Variable y(x)
In view of Eqs. (13) and (14), if v(x) is a transformed eigenfunction, then γv(x), γ 0 is also a transformed eigenfunction. For the uniqueness of v(x) and then u(x) by Eq. (29), we can consider a normalization condition as follows: where A 0 is a given nonzero constant. We fix A 0 = 1.
Theorem 3 If Eq. (45) is imposed for the uniqueness of v(x), then the terminal value y(1) in Eq. (32) is given by Proof. Differentiating Eq. (32) with respect to x and inserting x = 0, yields Therefore, we come to a definite second-order initial value problem from Eqs. (36) and (46) for the new variable y(x): which is a nonhomogeneous ODE, such that the given initial values can be taken as y(0) = 0 and y (0) = 0 without leading to y(x) = 0. Hence, Eq. (49) reduces to where we have taken A 0 = 1. Notice the difference between Eqs. (36) and (50). Eq. (36) exists an unknown constant y(1), which is disappeared in Eq. (50).
Up to here we have presented a new formulation of the SLP. When the original SLP in Eqs. (11) and (12) exist two unknowns: an unknown initial value u(0) or u (0) and the eigenvalue λ, the present formulation in Eq. (50) only has one unknown λ to be determined below.
The curve of v(1; λ) vs. λ is an eigenvalue curve. We can adjust the value of λ until v(1; λ) satisfies |v(1; λ)| < 2 , where 2 is a given tolerance of the error to mismatch the right-end Dirichlet boundary condition in Eq. (14): which is a numerically characteristic equation to determine λ. We can apply the half-interval method to solve Eq. (52); emanating from y(0) = y (0) = 0, we repeatedly integrate Eq. (50) with N steps from x = 0 to x = 1 to obtain y(1) and then v(1; λ) = y(1) + 1 is available for each λ. In the original SLP of Eqs. (11) and (12) the target equation a 2 u(1) + b 2 u (1) = 0 upon comparing to Eq. (52) is more complicated to determine the eigenvalue. In the present formulation of Eq. (50), the target equation y(1) + 1 = 0 is simple to determine the eigenvalue.
Like Eq. (42) for Example 1, Eq. (52) can be used to determine the eigenvalue λ. However, when Eq. (42) with its explicit form can be effectively solved by the Newton iteration method, Eq. (52) has to be solved by other methods, since it is an implicit form nonlinear equation of λ. The half-interval method (HIM) is described as follows. If λ d is a desired eigenvalue, we can detect the eigenvalue curve and decide two initial values λ 0 and λ with λ 0 < λ d < λ, such that v(1; λ 0 )v(1; λ) < 0. Then, we repeat the following processes until convergence: Endif N u is the upper bound of the number of iterations, which is set to be N u = 100. Besides the first two integrations, twice integrations of Eq. (50) in each iteration are required in the HIM.

A Derivative-Free Iterative Scheme
As to be shown below by numerical examples, the HIM is convergent slowly. Liu, Hong & Lee (2021) have proposed an iterative scheme: where in which x * is a simple root of the nonlinear equation f (x) = 0. The iterative scheme (53) has third-order convergence, which is shortened as LHL.
Like the half-interval method the first step used in the LHL is choosing two initial guesses x 0 and x 2 such that f (x 0 ) f (x 2 ) < 0 to render x * ∈ (x 0 , x 2 ). Then, we take x 1 = (x 0 + x 2 )/2. As the approximations of α and β in Eq. (54), we can use the technique of finite differences: In summary, the procedures of the LHL are given as follows: (i) Give the initial guesses of x 0 and x 2 to render f (x 0 ) f (x 2 ) < 0. (ii) Compute α and β by Eq. (55). (iii) For n=0,1,. . . , doing , until | f (x n )| < 3 .
We will apply the LHL to solve Eq. (52), which is a derivative-free iterative scheme. We fix the convergence criterion 1 for the FTIM, 2 for the HIM and 3 for the LHL.

Examples 1 and 2
For example 1, we fix 2 = 3 = 10 −15 and N = 5000 and iteratively determine λ by applying the half-interval method (HIM) and LHL to Eq. (52), and the results are listed in Table 3 to show the computed λ, the exact λ, EE, NI a , NI b , ERBC and ME. Since the convergence criterion and N used in the HIM and LHL are the same, they obtained the same results, but the numbers of iterations for HIM denoted as NI a are much larger than that for LHL denoted as NI b . When the CPU time spent in the LHL is about 4.28 s, the HIM needs 43.38 s. The LHL is much saving CPU time than the HIM. For example 2, we fix 2 = 3 = 10 −15 and N = 5000, and the results are listed in Table 4 to show the computed λ, the exact λ, EE, NI a , NI b , and ERBC.  Tables 3 and 4 to Tables 1 and 2, the improvement of the accuracy in every aspect is apparent with about three orders. For the purpose to achieve the high precision in Tables 3 and 4, N is increased to N = 5000. The NIs in the half-interval method (HIM) is increased and the CPU time for integrating Eq. (50) with 2(NI+1) times is increased, but in the FTIM it needs NI times integrations of the ODE with a less CPU time. Although the CPU time is increased for the HIM, but this is deserved because of the increase of accuracy with three orders than that for the FTIM.

Example 3
We consider a singular SLP (El-Gamel & El-Hady, 2013;Bailey, Everitt & Zettl, 1991): which is already in the form of Theorem 1, and thus v(x) is replaced by u(x) in Sections 4.1 and 4.2.
We fix 2 = 3 = 10 −15 , ∆x = 1/5000 and N = 5000, and the results are listed in Table 5 to show the computed λ, the exact λ, EE, NI a , and NI b . Notice that ERBC is zero. The accuracy is competitive to that in (El-Gamel & El-Hady, 2013), but is more accurate than that in (Bailey, Everitt & Zettl, 1991) about four orders, where the code SLEIGN2 was used.  Vol. 14, No. 4;2022 To further test the accuracy, we consider (El-Gamel & El-Hady, 2013) We fix 2 = 3 = 10 −15 and ∆x = 1/N = 1/8000, and the results are tabulated in Table 6 to show the computed λ, the exact λ, EE, NI a , and NI b . The accuracy is slightly better than (El-Gamel & El-Hady, 2013) in the first three eigenvalues.

Generalized Sturm-Liouville Problems
Let us consider Eq. (11) with the λ-dependent mixed boundary conditions: which is a generalized Sturm-Liouville problem.

A Main Result
Theorem 4 The generalized Sturm-Liouville problem (11) and (60) with a 2 (λ)b 2 (λ) 0 can be transformed to a new one with the Dirichlet boundary condition on the right-end: where in which the primes in q (x, λ) denotes the differential with respect to x.
Differentiating Eq. (65) with respect to x and inserting Eq. (11) for u (x), we have From Eqs. (65) and (66), we can express u(x) and u (x) in terms of w(x) and w (x) as follows: where D(x, λ) is defined in Eq. (64).

Journal of Mathematics Research
It follows from Eqs. (61), (64) and (72) that If we can solve y with the initial values y(0) and y (0) being given, Theorem 5 guarantees that w(x) satisfies the boundary conditions in Eq. (62). However, y(1) in G(x, λ) of Eq. (71) is unknown, whose exact value can be determined below.

A Definite Initial Value Problem
In view of Eqs. (61) and (62), if w(x) is a transformed eigenfunction, then γw(x), γ 0 is also a transformed eigenfunction. For the uniqueness of w(x) and then u(x) by Eq. (67), we can consider a normalization condition as follows: where A 0 is a given nonzero constant.
Theorem 6 If Eq. (74) is imposed for the uniqueness of w(x), then the right-end value y(1) in Eq. (72) is given by Proof. Differentiating Eq. (72) with respect to x and inserting x = 0, yields which immediately leads to Eq. (75) by using w (0) = A 0 in Eq. (74). 2 Therefore, we come to a definite second-order ODE from Eqs. (73), (71) and (75) for the new variable y(x): where and G (x, λ) = 0 was taken into account. Upon giving the initial values y(0) and y (0) and fixing A 0 = 1, Eq. (77) is a definite initial value problem with only unknown eigenvalue λ involved.
We can integrate Eq. (77) to obtain y(1); hence, inserting x = 1 into Eq. (72), yields where y(1) is obtained from the integration of Eq. (77) with the given initial conditions y(0) and y (0), which is no longer the one in Eq. (75).
To match the right-end Dirichlet boundary condition in Eq. (62), we can derive which is a numerically characteristic equation to determine λ. We can adjust the value of λ until w(1; λ) satisfies |w(1; λ)| < 2 , by applying the HIM or LHL to solve Eq. (80) to determine the eigenvalue λ.

Example 7
Consider (Chanane, 2005) u where the exact eigenvalues can be obtained from (Chanane, 2005) by In order to raise the accuracy, we fix 2 = 3 = 10 −15 and ∆x = 1/5000, and iteratively determine λ by applying the HIM and LHL to Eq. (80), and the results are listed in Table 9 to show the computed λ, the exact λ, EE, NI a , NI b , and ERBC.

Periodic Sturm-Liouville Problem
Here, we start from the periodic Sturm-Liouville problem (Andrew, 1989): and consider three types periodic boundary conditions: case (c) : For case (a), A 0 and B 0 are not given, such that three unknown values λ, y(0) and y (0) are to be determined. For case (b), A 0 is not given and B 0 is given, two unknown values λ and y(0) to be determined. For case (c), A 0 and B 0 are given, one unknown value λ to be determined. Obviously, the case (a) is more difficult than the cases (b) and (c).
Theorem 8 For any functions z 1 (x), z 2 (x) ∈ C 1 [0, ], we set the boundary shape functions to be where s 0 0. If then y 1 (x) and y 2 (x) satisfy Eq. (99) for case (b) of the periodic Sturm-Liouville problem.

Iterative Algorithms
Based on the boundary shape functions in Theorems 7-9, we can develop iterative algorithms to determine the eigenvalue and eigenfunction of Eqs. (89)-(92), namely the boundary shape function methods.

Conclusions
A new result is transformed the Sturm-Liouville problem with mixed boundary conditions into a generalized Sturm-Liouville eigenvalue problem with the Dirichlet boundary conditions happened on two sides, which greatly enhanced the convergence speed of the iterative algorithm to determine the eigenvalues. Combining to the fictitious time integration method (FTIM), an iterative algorithm was presented by transforming the problem to an initial value problem for the new variable. The eigenvalue and a terminal value of the new variable are determined iteratively by using the FTIM, from which the iterative sequence of eigenvalues are monotonically converged to its limit point to approximate the true eigenvalue when the relative residual reduced to zero. We proposed a normalization condition for the transformed variable and derived the terminal value of the new variable explicitly, such that we can obtain a definite initial value problem for the new variable. In order to precisely meet the right-end Dirichlet boundary condition, a numerically characteristic equation was provided which was solved by the half-interval method (HIM) and a derivative-free iterative method LHL to obtain the eigenvalues with high precision. For the generalized Sturm-Liouville problem it needs to solve two unknowns: a missing initial value and the eigenvalue. The new formulation presented in the paper overcame this difficulty by transforming the generalized Sturm-Liouville problem to a definite initial value problem with the initial values given, and the Dirichlet boundary condition on the right-end acting as a target equation. Then, the uniqueness condition of the eigenfunction is an extra requirement to ascertain y(1) in Eq. (75), such that only the eigenvalues are unknown to be determined. The sequence of iterative eigenvalues generated from the numerically characteristic equation converged very fast to the true eigenvalues by using the HIM and LHL. For the periodic Sturm-Liouville problems, we considered three types periodic boundary conditions and utilized the boundary shape functions methods to transform them into initial value problems involved unknown terminal values of the new variables. Three iterative algorithms were developed to compute the eigenvalues, whose high efficiency and accuracy can be appreciated.
In the whole paper the boundary shape function methods and the initial value problem methods are developed to simplify the three types Sturm-Liouville problems, involving only unknown eigenvalues to be determined. The Dirichlet type target equations are solved by the HIM and LHL. Both methods are effectively to determine the eigenvalues; however, we found that the LHL is convergent faster than the HIM. Many numerical examples assured the high-performance of the proposed iterative algorithms through a few iterations and with high accuracy of the obtained eigenvalues and eigenfunctions to preserve the mixed boundary conditions or periodic boundary conditions.
Overall, the relative errors of the eigenvalues obtained by the HIM and LHL are located between 10 −13 and 10 −15 as shown in Tables 3-9, which are more accurate than that presented in the literature. Owing to the high-accuracy of the proposed methods, they might be applied to solve more difficult nonlinear Sturm-Liouville problems and inverse Sturm-Liouville problems in the future.