A Boundary Shape Function Method for Computing Eigenvalues and Eigenfunctions of Sturm–Liouville Problems

: In the paper, we transform the general Sturm–Liouville problem (SLP) into two canonical forms: one with the homogeneous Dirichlet boundary conditions and another with the homogeneous Neumann boundary conditions. A boundary shape function method (BSFM) was constructed to solve the SLPs of these two canonical forms. Owing to the property of the boundary shape function, we could transform the SLPs into an initial value problem for the new variable with initial values that were given deﬁnitely. Meanwhile, the terminal value at the right boundary could be entirely determined by using a given normalization condition for the uniqueness of the eigenfunction. In such a manner, we could directly determine the eigenvalues as the intersection points of an eigenvalue curve to the zero line, which was a horizontal line in the plane consisting of the zero values of the target function with respect to the eigen-parameter. We employed a more delicate tuning technique or the ﬁctitious time integration method to solve an implicit algebraic equation for the eigenvalue curve. We could integrate the Sturm–Liouville equation using the given initial values to obtain the associated eigenfunction when the eigenvalue was obtained. Eight numerical examples revealed a great advantage of the BSFM, which easily obtained eigenvalues and eigenfunctions with the desired accuracy.


Introduction
During 1836-1837, Sturm and Liouville created a new subject in the mathematical analysis of second-order ordinary differential equations (ODEs) with the homogeneous Sturm-Liouville boundary conditions. It is known as the Sturm-Liouville theory nowadays, and deals with the Sturm-Liouville problem (SLP): hu(a) − p(a)u (a) = 0, where p(x) > 0 and ω(x) > 0 are strictly positive functions of x in a finite interval [a, b] and λ is a constant parameter. The two-point boundary value problem (1) and (2) is a regular SLP that allows nontrivial solutions only for specific values of λ, which are known as the eigenvalues. Only in a few cases can the eigenvalues be obtained analytically by solving algebraic equations; however, we do not even have explicit algebraic equations to solve λ for most cases. SLPs are essential in partial differential equations, vibrations in continuum mechanics, heat conduction problems, the Schrödinger equation, the transport of microwaves, the Lax pair of the Korteweg-De Vries equation, and fractional derivative problems [1][2][3][4], etc. Many numerical method, such as the residuals-collocation method [5], Chebyshev collocation method [6,7], shooting method [8], exponentially fitted Numerov method [9], differential quadrature method [10,11], Sinc-Galerkin method [12], Adomian decomposition method [13], homotopy analysis method [14], and variational iteration method [15] have been developed in the past several decades to solve SLPs. These methods have been widely applied and efficient and accurate codes have been developed. On the other hand, the modified shooting methods were developed by Ghelardoni and Gheri [16] and Liu [17,18]. Compared with the corrected Numerov's method in recent years, competitive results have been obtained by applying symmetric boundary value methods [19,20] and using high-order difference schemes [21]. This paper's main focus was to develop the boundary shape function (BSF) to solve SLPs by automatically satisfying the boundary conditions. First, the idea of the boundary shape function proposed by Liu and Chang [22] was applied to find the periodic solution of nonlinear jerk equations. Then, some problems such as the optimal control problems of nonlinear Duffing oscillators and the boundary value problem (BVP) with different boundary conditions were addressed by the BSF. New methods based on the BSF were developed to address the SLPs that were easy to formulate and implement. There are various numerical methods to approximate the eigenvalue and the corresponding eigenfunction. The present technique used the BSF to transform the SLP to the initial value problem (IVP), which was drastically different from the Lie-group shooting method developed by Liu [17]. This new method did not need to find the missing initial values via the shooting technique; the difference will be pointed out below.
The paper is arranged as follows. Two main theorems are proven in Section 2 for the SLP with the homogeneous Sturm-Liouville boundary conditions. Here, we developed two novel transformation methods to derive two canonical forms that transformed the general SLP to those with homogeneous Dirichlet and Neumann boundary conditions. In Section 3, we transform the SLP with the Dirichlet boundary conditions to an IVP for a new variable, which guaranteed that the Dirichlet boundary conditions for y(x) could be satisfied automatically. In the transformed IVP for z(x), an unknown constant z(b) appeared that is yet to be determined. In Section 4, we derive z(b) in terms of the initial values for z(x) and a normalized value A 0 given in y (a) = A 0 for the uniqueness of the eigenfunction y(x), which could avoid the iteration to determine z(b). The numerical algorithm based on the boundary shape function method (BSFM) is developed in Section 5. Numerical examples with the homogeneous Dirichlet boundary conditions are given in Section 6 that verify the new results and algorithms. In Section 7, the BSFM is derived to solve the SLP with the homogeneous Neumann boundary conditions. In Section 8, several examples are given of homogeneous Sturm-Liouville boundary conditions. Finally, some conclusions are presented in Section 9.

Sturm-Liouville Boundary Conditions
It is well known that the eigenvalues satisfy the relation: and the eigenfunction corresponding to an eigenvalue λ n has n intersections with the x-axis inside the interval x ∈ (a, b). Moreover, since for any given eigenvalue it is possible to compute an infinite number of eigenfunctions, in general, a normalization condition in the form b a u 2 (x)dx = 1 is required for the uniqueness of u(x). The computation of λ and u(x) traditionally has been quite a challenging task. We will show that inexpensive computing that evaluates λ via a direct solution of the differential Equation (1) is feasible.
For the general SLP of Equations (1) and (2), we supposed that u(a) = 0, u (a) = 0, u(b) = 0, and u (b) = 0. We transformed them into two canonical forms with the Dirichlet and Neumann homogeneous boundary conditions, respectively, and then proposed a BSFM for determining the eigenvalues and finding the corresponding eigenfunctions.

Transforming to Neumann-Type Boundary Conditions
In Equations (1) and (2), suppose that h = 0 and H = 0 such that the coefficients preceding u(a) and u(b) can be set equal to one. When letting µ 1 = −p(a)/h and µ 2 = −p(b)/H, we may need to prove the following result. Theorem 1. For the SLP with the Sturm-Liouville type boundary conditions: where λ is an eigen-parameter. If the following condition is satisfied: then Equations (4) and (5) can be transformed into the Neumann-type SLP: The relation between u(x) and y(x) is given as: Proof. Let: It is apparent that the Neumann boundary conditions in Equation (8) follow from Equation (5) upon using the above equation. It follows from Equation (10) that: If condition (6) holds, Equation (11) reduces to: Thus, y(x) = (x + µ 1 − a)u(x) in Equation (9) is proved. Inserting Equation (9) and its first and second differentials into Equation (4), we can derive: which, when multiplied by 1/(x + µ 1 − a) 2 on both sides, can lead to Equation (7). Motivated by Theorem 1, we generalized it to the Sturm-Liouville type boundary conditions without needing the constraint µ 2 − µ 1 = b − a. For this purpose, we sought a function F(x) in: so that the Sturm-Liouville type boundary conditions in Equation (5) with u (a) = 0 and u (b) = 0 were transformed to: Through some manipulations, we can derive: Theorem 2. The SLP (4) and (5) with u (a) = 0 and u (b) = 0 in the Sturm-Liouville boundary conditions can be transformed to the Neumann-type boundary conditions: where the relation between u(x) and y(x) is given by Equation (14) and F(x) is given by Equations (17) or (18).
Proof. The Neumann boundary conditions in Equation (20) are given by Equations (15) and (16). By inserting Equation (14) and its first differential into Equation (4), we can derive: When Equation (21) is divided by F(x) on both sides and use: we can change it to: The expansion of the second term generates Equation (19).
Note that Niessen and Zettl [23] developed a similar transformation (y(x) = u(x)/ f (x)), where f (x) is a solution of Equation (1) for some λ. The purpose was to transform the singular SLP with a singular nonoscillatory limit circle endpoint into a regular one. However, the present F(x) in Equations (17) and (18) was different from the 1/ f (x) used in [23].

Transforming to Dirichlet-Type Boundary Conditions
This is similar to Theorem 2 in that we can derive the following result:  (4) and (5) can be transformed to the Dirichlet-type SLP: where: Proof. By letting: where p(a) > 0 and p(b) > 0, we transform Equations (4) and (5) into the canonical forms of (23) and (24). It is easy to confirm that y(x) in Equation (26) satisfies Equation (24) when u(x) satisfies Equation (5). For easy operations, y(x) in Equation (26) is written as: where: When taking the derivative of Equation (27) with respect to x, we get: where for simplicity we omit the variable x in the functions, and: is a constant. Inserting Equation (4) for (pu ) into Equation (29) leads to: Using Equations (27) and (31), u and pu in terms of y and y can be solved as follows: where: Further, when taking the derivative of Equation (33) with respect to x and using Equation (4) for (pu ) again, we get: The latter equation can be written as: which, upon using Equations (34) and (32) for u, can be arranged as: Thus, we prove Equation (23).

Transforming to an Initial Value Problem for z(x)
To clearly demonstrate the new BSFM, we began with the SLP with the homogeneous Dirichlet boundary conditions: where p(x), q(x), and ω(x) ∈ L 1 [a, b], with p(x) > 0 and s(x) > 0. In view of Equations (40) and (41), if y(x) is an eigenfunction, then αy(x), α = 0 is also an eigenfunction. For the uniqueness of y(x), we considered a normalization condition as follows: where A 0 is a given nonzero constant. Equation (42) is more straightforward than the usual normalization condition b a y 2 (x)dx = 1. According to the theory of ODE, the solution of y(x) that satisfies Equation (40) with the initial conditions y(a) = 0 and y (a) = A 0 = 0 is unique.
To guarantee that the solution of y(x) could exactly satisfy Equation (41), we introduced a shape function: As a consequence, we could prove the following result: by Equation (43), which becomes: Similarly, inserting x = b into Equation (44) yields: which, in view of Equation (43), becomes: Thus, we complete the proof.
Theorem 4 is crucial in that the newly introduced shape function method guaranteed that the boundary conditions in Equation (41) could be exactly satisfied by y(x) in Equation (44). Corresponding to the shape function s(x), y(x) was called a boundary shape function since it automatically satisfied the boundary conditions (41). Starting from this new idea, we could develop a BSFM to solve Equations (40) and (41) by considering the relation between y(x) and z(x): where: by substituting Equation (43) which can be viewed as an IVP for z(x) in the interval x ∈ [a, b] upon giving the initial values C 1 and C 2 definitely.

Finding the Terminal Value of z(x)
However, z(b) in the function G(x) given by Equation (46) is an unknown constant. For Equations (47) and (48) to be the IVP, we had to ensure that z(b) was not an unknown constant, which was proved as follows: Moreover, we have: where z p and z h are, respectively, the particular and homogeneous solutions of Equation (47).
Proof. Taking the differential of Equation (45) with respect to x and inserting x = a leads to: where G is a constant due to Equation (46); from Equation (48), it follows that: When substituting Equation (52) for G , y (a) = A 0 as shown in Equation (42) and Equation (48) for z (a) = C 2 in Equation (51), we get: hence, Equation (49) is proven. Equation (47) is a nonhomogeneous ODE for z(x). We could decompose z(x) into: where z h (x) and z p (x) are the homogeneous and particular solutions of Equation (47), respectively. The latter means that z p (x) satisfied the nonhomogeneous part of Equation (47): which by observation immediately implied the first result in Equation (50). Inserting Equation (54) into Equation (45) and using z p (x) = G(x) yields: Thus, we end the proof.
The derivation of Equation (49) for z(b) was independent of the governing Equation (4), which meant that the current method of BSFM was applicable to more general eigenvalue problems rather than Equation (4); for instance, nonlinear eigenvalue problems.

Algorithm to Solve Eigenvalues and Eigenfunctions
We have proved the two main Theorems, 4 and 5, which could help us efficiently compute the eigenvalues. The result in Theorem 5 was very important because we could avoid the iteration to determine z(b). Upon letting it follows from Equation (47) that: where the initial condition ξ(a) = p(a)A 0 for ξ(x) is based on Equations (49), (52), and (58): In Equations (58) Because the initial values z(a) = C 1 and ξ(a) = p(a)A 0 in Equation (60) were available, now, from x = a to x = b, we could apply the fourth-order Runge-Kutta method (RK4) to integrate the first-order ODEs (58) and (59) with a given λ to obtain z(b); hence: could be obtained for each λ given. For use in the latter, we called the curve of y(b; λ) versus λ an eigenvalue curve, where y(b; λ) was a target function and λ was an eigen-parameter. We could adjust the value of λ until y(b; λ) satisfied |y(b; λ)| < ε, where ε is a given tolerant error to mismatch the right-boundary condition y(b; λ) = 0. The algorithm based on the BSFM for solving the SLP with the homogeneous Dirichlet boundary conditions can be summarized as follows: initial guess λ (0) . (ii) For k = 1, 2, . . . , applying the RK4 to integrate the ODEs in Equations (58)-(60) with N steps from x = a to x = b. Take: where y(b; λ (k) ) is an implicit function of λ (k) ; then, we employed the fictitious time integration method (FTIM) [24] to compute the next step value: where t k = k∆t is a fictitious time; if λ (k) rendered: then we could find an eigenvalue inside an interval in which we were interested. We could obtain all eigenvalues sequentially by changing the interval and going to (ii). In (i), the initial guess λ (0) could be observed in the data of the eigenvalue curve. A more straightforward method to determine λ is to use a finer tuning technique (FTT).
including the eigenvalues of interest, from which the zero points inside that interval can be detected from the plot of y b; λ (k) versus λ ∈ [c, d]. We could tune the interval to a finer one by gradually reducing the length of the interval [c, d] to locate the minimum of: until the minimum was smaller than ε. At this moment, we could precisely find the eigenvalue inside that interval. We can confirm that the initial values of y(x) and y (x) are: The first one was already proved in Theorem 4. Based on Equations (45), (46), (49), and (52), it followed that: Suppose one wants to solve the corresponding eigenfunction for an eigenvalue λ obtained. In that case, one can employ the initial values in Equation (68) and integrate Equation (40) after inserting the obtained eigenvalue λ to solve the eigenfunction. (59) involves λ as an unknown value, the integrated value y(b; λ) is indeed an implicit function of λ. By meeting the right boundary condition for y(x), an algebraic equation y(b; λ) = 0 can be used to determine the eigenvalue λ, which can be performed by a numerical method based on FTIM [24]. The readers may ask why we did not directly integrate Equation (40) by using the initial values y(a) = 0, y (a) = C 0 , with the guessed C 0 and the guessed eigenvalue λ to match the right boundary condition y(b) = 0 to pick up the correct eigenvalue. However, this approach cannot guarantee that the right boundary condition can be exactly satisfied, which may render an unstable and incorrect solution. We will give an example to demonstrate the failure of this naive approach in the shooting method. We instead integrated z(x) using Equations (58)-(60) in the BSFM, which guaranteed that the right boundary condition y(b; λ) = 0 could be satisfied automatically as proved in Theorem 4. From this aspect, the BSFM was quite different from the shooting method, where the correct value C 0 was not known in advance, and even Equation (40) was integrated by using the correct initial values y(a) = 0, y (a) = C 0 , while the end value y(b) was not guaranteed to satisfy the right-end condition y(b) = 0.

Remark 1. Because Equation
For the general SLPs (4) and (5), the boundary conditions of which were not of the Dirichlet type, we could employ Theorem 3 to transform them into Equations (23) and (24). Then, after using the above BSFM to determine λ, y(x), and y (x), we could determine the eigenfunction u(x) using Equation (32).

Example 1
For the first test example, we considered: where λ k = (k + 1) 2 π 2 , k ∈ N, y k = x sin((k + 1)π lnx), and k = 0, 1, . . .. We chose an adequately large integer N to enhance the accuracy when integrating ODEs using the RK4, the CPU time of which was saved because we merely needed to integrate two first-order ODEs with N steps from x = a to x = b. When a larger eigenvalue was computed, a larger N was demanded. Figure 1a shows a plot of the eigenvalue curve of y(b) versus λ in an interval (0, 300) with M = 601, in which we can see that there were five zero points that corresponded to the first five eigenvalues.

of 22
Suppose one wants to solve the corresponding eigenfunction for an eigenvalue obained. In that case, one can employ the initial values in Equation (68) and integrate Equaion (40) after inserting the obtained eigenvalue to solve the eigenfunction. emark 1. Because Equation (59) involves as an unknown value, the integrated value ( ; ) s indeed an implicit function of . By meeting the right boundary condition for ( ), an algebraic quation ( ; ) = 0 can be used to determine the eigenvalue , which can be performed by a umerical method based on FTIM [24]. The readers may ask why we did not directly integrate quation (40) by using the initial values ( ) = 0, ′( ) = 0 , with the guessed 0 and the uessed eigenvalue to match the right boundary condition ( ) = 0 to pick up the correct eienvalue. However, this approach cannot guarantee that the right boundary condition can be exctly satisfied, which may render an unstable and incorrect solution. We will give an example to emonstrate the failure of this naive approach in the shooting method. We instead integrated ( ) sing Equations (58)-(60) in the BSFM, which guaranteed that the right boundary condition ( ; ) = 0 could be satisfied automatically as proved in Theorem 4. From this aspect, the BSFM as quite different from the shooting method, where the correct value 0 was not known in adance, and even Equation (40) was integrated by using the correct initial values ( ) = 0, ′( ) = 0 , while the end value ( ) was not guaranteed to satisfy the right-end condition ( ) = 0.
For the general SLPs (4) and (5), the boundary conditions of which were not of the irichlet type, we could employ Theorem 3 to transform them into Equations (23) and 24). Then, after using the above BSFM to determine , ( ), and ′( ), we could deterine the eigenfunction ( ) using Equation (32).

.1. Example 1
For the first test example, we considered: here = ( + 1) 2 2 , ∈ ℕ, = sin(( + 1) ln ), and = 0, 1, …. We chose an adequately large integer to enhance the accuracy when integrating DEs using the RK4, the CPU time of which was saved because we merely needed to ntegrate two first-order ODEs with steps from = to = . When a larger eigenalue was computed, a larger was demanded. Figure 1a shows a plot of the eigenvalue urve of ( ) versus in an interval (0, 300) with = 601, in which we can see that here were five zero points that corresponded to the first five eigenvalues. We applied the BSFM to solve this problem with 0 = 1, 1 = 2 = 1, and = 1000. As shown in Figure 2, by tuning to finer intervals using the FTT and applying the BSFM with 0 = ( + 1) (for coinciding with the exact one), 1 = 2 = 1, = 10 −9 , and = 1000, we could sequentially obtain the required eigenvalues. We showed the calculated eigenfunctions for = 0, 2, 9, which coincided with the exact ones. It can be seen in Figure 3 that the numerical errors of the eigenfunctions were very small, on the We applied the BSFM to solve this problem with A 0 = 1, C 1 = C 2 = 1, and N = 1000. As shown in Figure 2, by tuning to finer intervals using the FTT and applying the BSFM with A 0 = (k + 1)π (for coinciding with the exact one), C 1 = C 2 = 1, ε = 10 −9 , and N = 1000, we could sequentially obtain the required eigenvalues. We showed the calculated eigenfunctions for k = 0, 2, 9, which coincided with the exact ones. It can be seen in Figure 3 that the numerical errors of the eigenfunctions were very small, on the order of 10 −10 to 10 −8 . We applied the BSFM to solve this problem with 0 = 1, 1 = 2 = 1, and = 1000. As shown in Figure 2, by tuning to finer intervals using the FTT and applying the BSFM with 0 = ( + 1) (for coinciding with the exact one), 1 = 2 = 1, = 10 −9 , and = 1000, we could sequentially obtain the required eigenvalues. We showed the calculated eigenfunctions for = 0, 2, 9, which coincided with the exact ones. It can be seen in Figure 3 that the numerical errors of the eigenfunctions were very small, on the order of 10 −10 to 10 −8 .  The present problem had three eigenvalues in the range of λ < 100, as shown in Figure 1a by the three intersection points before λ = 100, which are listed in Table 1. If we directly integrated Equation (40) by employing the initial values y(a) = 0 and y (a) = C 0 and picked the zero points in the eigenvalue curve, it led to incorrect results as listed in Table 1. In contrast, we could obtain very accurate eigenvalues by using the FTT, FTIM, and BSFM no matter which value of A 0 = 1 or A 0 = 2 was used; the errors were on the order of 10 −11 to 10 −10 . In the FTIM, we took ∆t = 0.001 and v = −50 for λ 0 , ∆t = 0.0005 and v = 150 for λ 1 , and ∆t = 0.0005 and v = −200 for λ 2 .   Table 1. In contrast, we could obtain very accurate eigenvalues by using the FTT, FTIM, and BSFM no matter which value of 0 = 1 or 0 = 2 was used; the errors were on the order of 10 −11 to 10 −10 . In the FTIM, we took ∆ = 0.001 and = −50 for 0 , ∆ = 0.0005 and = 150 for 1 , and ∆ = 0.0005 and = −200 for 2 . To remedy the shortcoming of the naive approach to the IVP of ( ), one way is to use the shooting method [16,17] to pick up the correct initial value of ′( ) = 0 rather than ( ) = 0 such that the integrated value ( ) can match the true right-end value ( ) = 0. From this point, we can see that the BSFM presented here differed significantly from the shooting method. In contrast to the shooting method, the initial values 1 and 2 for the IVP of ( ) were given definitely and the right boundary condition ( ) = 0 was satisfied automatically and exactly upon solving for ( ) with the correct eigenvalue . Notice the differences between the IVP for ( ) in the shooting method and the BSFM for ( ). In the shooting method, there were two unknown values 0 and to be determined. In contrast, was the only unknown value to be determined in the BSFM.

Remark 2.
To remedy the shortcoming of the naive approach to the IVP of y(x), one way is to use the shooting method [16,17] to pick up the correct initial value of y (x) = 0 rather than y(a) = 0 such that the integrated value y(b) can match the true right-end value y(b) = 0. From this point, we can see that the BSFM presented here differed significantly from the shooting method. In contrast to the shooting method, the initial values C 1 and C 2 for the IVP of z(x) were given definitely and the right boundary condition y(b) = 0 was satisfied automatically and exactly upon solving for z(x) with the correct eigenvalue λ. Notice the differences between the IVP for y(x) in the shooting method and the BSFM for z(x). In the shooting method, there were two unknown values C 0 and λ to be determined. In contrast, λ was the only unknown value to be determined in the BSFM.

Example 2
For this example, we considered the following SLP [8,10,17]: No closed-form solutions of λ and y(x) exist. The BSFM with A 0 = 1, C 1 = C 2 = 1, N = 2000, and M = 3001 was adopted to determine the eigenvalues in the range of 4 < λ < 920 as shown in Figure 1b, in which 30 intersection points could be observed. The results agreed with those obtained by Liu [17]. Table 2 compares the calculated eigenvalues with those obtained by Ghelardoni et al. [8] and Liu [17]. By using the FTT and applying the BSFM with A 0 = 1, C 1 = C 2 = 1, ε = 10 −12 , and N = 2000, we could sequentially obtain the required eigenvalues shown in Table 2. Figure 4 displays the calculated eigenfunctions for k = 4, 9, 29, where the numerical errors of the eigenfunctions to satisfy the right boundary condition were very small, on the order of 10 −15 to 10 −12 .  Table 2 compares the calculated eigenvalues with those obtained by Ghelardoni et al. [8] and Liu [17]. By using the FTT and applying the BSFM with 0 = 1, 1 = 2 = 1, = 10 −12 , and = 2000, we could sequentially obtain the required eigenvalues shown in Table 2. Figure 4 displays the calculated eigenfunctions for = 4, 9, 29, where the numerical errors of the eigenfunctions to satisfy the right boundary condition were very small, on the order of 10 −15 to 10 −12 .

Example 3
We considered the following SLP [25]:

Example 3
We considered the following SLP [25]: We employed the BSFM with A 0 = 1, C 1 = C 2 = 1, N = 2000, and M = 3001 to search for the eigenvalues in the range of 1 < λ < 20. Table 3 compares the calculated eigenvalues through the FTT to adjust the eigenvalues to match the right-boundary condition y(π) = 0 with those obtained by Eggert et al. [25]. The accuracy of the computed eigenvalues by the BSFM could arrive at 10 −10 .

Variable Transformation
In this section, we considered Equation (40) with the following homogeneous Neumann boundary conditions: For the SLP with the homogeneous Neumann boundary conditions, we considered another normalization condition: According to the theory of ODE, the nontrivial solution of y(x) satisfying Equation (40) with the initial conditions y(a) = B 0 = 0 and y (a) = 0 is unique.
We can derive: where s 20 is a given constant, so that: Theorem 6. For the function z(x) ∈ C 2 [a, b]: satisfies Equation (77).
Proof. Taking the differential of Equation (82) with respect to x and inserting x = a leads to: which, when using Equations (79) and (80), becomes: When inserting x = b into the differential of Equation (82), we get: which, when using Equations (79) and (80), becomes: Thus, the proof is completed.
Consider the variable transformation from y(x) to z(x) by: where: Using Equation (83), we can transform Equation (40) into: where C 1 and C 2 are constant initial values.

Algorithm Based on the BSFM
Let: The latter condition could be proved as follows. After taking the differential of Equation (85), inserting x = a, and using Equations (79) and (80), we could obtain z (a) = H (a). Inserting x = a into Equation (91) and using z (a) − H (a) = 0 yielded ξ(a) = 0. Then, we could apply the RK4 to integrate the above ODEs to obtain: The algorithm based on the BSFM for solving the SLP with the homogeneous Neumann boundary conditions can be summarized as follows: (i) Give C 1 , C 2 , B 0 , ∆t, v, the stopping criterion ε, ∆x = (b − a)/N with N given, and initial guess λ (0) . (ii) For k = 1, 2, . . ., applying the RK4 to integrate the ODEs in Equations (91)-(93) with N steps from x = a to x = b. Take: employ the FTIM [24] to obtain: and if λ (k) renders: then we can find an eigenvalue inside an interval; we can then obtain all eigenvalues sequentially by changing the interval and going to (ii). A more straightforward method to determine λ is to use an FTT. Let λ (j) = c + (j − 1)(d − c)/(M − 1), j = 1, . . . , M run in an interval [c, d], including the eigenvalues, from which the zero points inside that interval can be detected from the plot of y b; λ (k) versus λ ∈ [c, d]. We could tune the interval to a finer one to locate the minimum of: until the minimum was smaller than ε. Suppose one wants to solve the corresponding eigenfunction for an eigenvalue λ obtained. In that case, one can employ the initial values y(a) = B 0 and y (a) = 0 and integrate Equation (40) with the obtained eigenvalue λ to solve the corresponding eigenfunction.
For the general SLP (4) and (5), the boundary conditions of which were not the Neumann type, we could employ Theorem 1 or Theorem 2 to transform them into Equations (7) and (8) or Equations (19) and (20). Then, after using the above BSFM to determine λ and y(x), we could determine the eigenfunction u(x) using Equation (9) or Equation (14).

Example 4
We tested the BSFM for the eigenfunction in Equation (76), which had the eigenvalues λ k = (k + 1) 2 π 2 , k = 0, 1, 2, . . .. We applied the BSFM in Section 7.2 to solve this problem with B 0 = 1 or B 0 = 2, s 20 = 1, C 1 = C 2 = 1, N = 1000, and M = 101. By tuning to finer intervals, we compared the computed eigenvalues to the exact ones for k = 0, 1, 2 as shown in Table 4. It can be seen that the numerical errors were minimal (on the order of 10 −11 to 10 −10 ) no matter which value of B 0 = 1 or B 0 = 2 was used.

Numerical Examples of the General Sturm-Liouville Problem
For the general SLP, we could either transform it into the Neumann-type SLP as shown in Theorems 1 and 2 or the Dirichlet-type SLP as shown in Theorem 3. Some examples are given below.

Example 5
The following example was borrowed from Fu et al. [26]: where e 0 is a nonzero parameter. This problem had the eigenvalues λ k = (k + 1) 2 π 2 , k = 0, 1, 2, . . ., but we did not have closed-form solutions for the eigenfunctions.

Example 6
At this moment, we solved the general SLP (1) and (2) by using Theorem 2 and the BSFM in Section 7.2. To test the accuracy, we considered: By using the FTT and the FTIM, we compared the calculated eigenvalues to the exact ones for k = 0, 1, . . . , 4 in Table 5. It can be seen that the numerical errors were minimal-on the order of 10 −11 to 10 −10 when using the FTT and 10 −11 to 10 −7 when using the FTIM. By tuning to finer intervals when using the FTT and applying the BSFM with B 0 = 1 or B 0 = 2, s 20 = 1, C 1 = C 2 = 1, N = 1000, and M = 101, we could sequentially obtain the required eigenvalues; the calculated eigenfunctions for k = 0, 1, 5 are shown in Figure 5b as solid lines. It can be seen that the numerical errors were minimal (on the order of 10 −11 to 10 −8 to satisfy the right boundary condition u (x) + u(1)e 0 /(1 + e 0 ) = 0, where we took e 0 = 1.

Example 6
At this moment, we solved the general SLP (1) and (2) by using Theorem 2 and the BSFM in Section 7.2. To test the accuracy, we considered: where λ k , k = 0, 1, 2, . . . were solved using µ 1 µ 2 λ 2 sin λ + (µ 2 − µ 1 )λ cos λ + sin λ = 0; while the eigenfunctions were given by u k = µ 1 λ k cos(λ k x) − sin(λ k x). For Example 5, we showed that using the FTT was more accurate than using the FTIM. We used the FTT to adjust the eigenvalues for the following examples. We applied the BSFM in Section 7.2 to solve this problem with µ 1 = 1 or µ 2 = −2, s 20 = 1, C 1 = C 2 = 1, N = 1000, and M = 101. By tuning to finer intervals using the FTT, we compared the calculated eigenvalues to the exact ones for k = 0, 2, 4 in Table 6. It can be seen that the numerical errors were very small, on the order of 10 −11 to 10 −10 .

Example 7
We considered Example 5 again; however, we used a different right boundary condition: We employed Theorem 2 and applied the BSFM in Section 7.2 to solve this problem with B 0 = 1, s 20 = 1, C 1 = C 2 = 1, N = 500, and M = 2001. Figure 5a shows a plot of the eigenvalue curve of y (b) versus λ by the dashed line in an interval [1,1000], from which we can see that there were 10 zero points that were close to the first 10 eigenvalues of Example 5. However, there were a few differences as shown in Table 5. The numerical errors of the eigenvalues were adjusted to be less than 10 −6 .
By tuning to finer intervals using the FTT and applying the BSFM in Section 7.2 with B 0 = 1, C 1 = C 2 = 1, ε = 10 −8 , and N = 1000, we could sequentially solve for the required eigenvalues; the calculated eigenfunctions for k = 1, 3, 5 are shown by the dashed lines in Figure 5b. The right boundary condition u(1) − 2u (1) = 0 could be satisfied very well, as shown in Table 7.

Example 8
Now, we solved the general SLP (1) and (2) by using Theorem 3 and the BSFM in Section 5. We tested the following problem given in Example 5: where µ 1 = g 0 and µ 2 = 1 + g 0 . This problem had the eigenvalues λ k = (k + 1) 2 π 2 , k = 0, 1, 2, . . .. We applied the present BSFM based on Theorem 3 to solve this problem with g 0 = 1, A 0 = 1, C 1 = C 2 = 1, N = 1000, and M = 101. By tuning to finer intervals using the FTT, the calculated eigenvalues were compared with the exact ones for λ k = (k + 1) 2 π 2 , k = 0, 1, 2, . . . , 4 in Table 8. As shown in the table, the numerical errors were on the order of 10 −13 to 10 −9 . Remark 3. The Liouville transformation [26]: can transform Equations (1) and (2) to the following canonical form: where: The original SLP became 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) was required.
Upon comparing the transform in Equation (109) to that in Equation (14), we could observe the similarity. However, F(x) was an explicit quadratic or cubic function of x; all the coefficient functions were given explicitly in Equation (19). Moreover, the boundary conditions were more accessible to treat than Equation (112). One major drawback of Equation (111) was that the coefficient function Q(t) was an implicit function of t, which could be known only after inverting the relation in Equation (110) to obtain the function x = x(t) and then inserting it into Equation (113).
Both the transformations in Equations (14) and (109) were congruent, which maintained the linearity in the eigenvalue λ or µ. On the other hand, the transform (25) in Theorem 3 was not a congruent one, where the coefficient function D(x, λ) is a linear function of λ and Q(x, λ) is a nonlinear function of λ. Nevertheless, from the aspect of the numerical method, the Dirichlet boundary conditions in Equation (24) were more uncomplicated and were easier to treat than the boundary conditions (20) and (112).
The SLP is very important in partial differential equations, vibrations in continuum mechanics, and heat conduction problems. We applied the Sturm-Liouville theory to solve the longitudinal wave motion problem when the Young's modulus E was fixed with the cross section area A(x) = 1 + x of a tapered rod, a constant density ρ of mass, and the axial displacement y(x, t): Let y(x, t) = e iωt u(x) with ω as a natural vibration frequency; Equation (116) then reduces to Equation (4) with: We fixed the length of the rod be a unit and took E/ρ = 100. In the design of engineering structures, knowing first few frequencies ω is of utmost importance. Figure 6 shows a plot of u(1) versus λ in a range of λ ∈ [0, 4000], where two intersecting points appeared. Then, by tuning to finer intervals, we could obtain λ 0 = 975.3318666643578 and λ 1 = 3935.600370373466.
We fixed the length of the rod be a unit and took ⁄ = 100. In the design of engineering structures, knowing first few frequencies is of utmost importance. Figure 6 shows a plot of (1) versus in a range of ∈ [0,4000], where two intersecting points appeared. Then, by tuning to finer intervals, we could obtain λ 0 = 975.3318666643578 and λ 1 = 3935.600370373466.

Conclusions
For the SLPs with the Sturm-Liouville boundary conditions, we derived two novel transformations that could transform them into two canonical forms with either the Dirichlet boundary conditions or the Neumann boundary conditions. These two canonical forms appeared in the literature at the same time. Numerically, these two simpler forms were more easily treated than the original Sturm-Liouville-type boundary conditions. The major novelty of the paper was that we employed the shape function to transform the Sturm-Liouville boundary value problem to the corresponding initial value problem for the new variable with terminal values that could be determined explicitly by using the Neumann-type normalization condition for the Dirichlet canonical form and the Dirichlettype normalization condition for the Neumann canonical form. Based on these novel ingredients in the BSFM, we could solve the Sturm-Liouville problems easily and directly. In the BSFM, the only unknown constant was the eigenvalue, which was determined by solving the derived algebraic equation with a finer tuning technique or by using the fictitious time integration method. Eight numerical examples confirmed the high performance of the new BSFM algorithm based on methods with high accuracy for the obtained eigenvalues and eigenfunctions. As compared to the conventional shooting method, we did not need to guess the missing initial value. The advantages of the presented two novel transformations were explored when we compared them to the Liouville transform. Future works can focus on extending the BSFM to solve generalized Sturm-Liouville problems and periodic Sturm-Liouville problems in terms of nonlinear problems.

Conclusions
For the SLPs with the Sturm-Liouville boundary conditions, we derived two novel transformations that could transform them into two canonical forms with either the Dirichlet boundary conditions or the Neumann boundary conditions. These two canonical forms appeared in the literature at the same time. Numerically, these two simpler forms were more easily treated than the original Sturm-Liouville-type boundary conditions. The major novelty of the paper was that we employed the shape function to transform the Sturm-Liouville boundary value problem to the corresponding initial value problem for the new variable with terminal values that could be determined explicitly by using the Neumann-type normalization condition for the Dirichlet canonical form and the Dirichlettype normalization condition for the Neumann canonical form. Based on these novel ingredients in the BSFM, we could solve the Sturm-Liouville problems easily and directly. In the BSFM, the only unknown constant was the eigenvalue, which was determined by solving the derived algebraic equation with a finer tuning technique or by using the fictitious time integration method. Eight numerical examples confirmed the high performance of the new BSFM algorithm based on methods with high accuracy for the obtained eigenvalues and eigenfunctions. As compared to the conventional shooting method, we did not need to guess the missing initial value. The advantages of the presented two novel transformations were explored when we compared them to the Liouville transform. Future works can focus on extending the BSFM to solve generalized Sturm-Liouville problems and periodic Sturm-Liouville problems in terms of nonlinear problems.

Conflicts of Interest:
The authors declare no conflict of interest.