The FD-method for solving Sturm-Liouville problems with special singular differential operator

In the paper we describe a superexponentially convergent numerical-analytical method for solving the eigenvalue problem for the some class of singular differential operators with boundary conditions. The method (FD-method) was firstly proposed by V. L. Makarov and successfully combines the benefits of using the {\it coefficient approximation methods} (CAM) and the homotopy approach. The sufficient conditions which provides convergence of the proposed method are stated and rigorously substantiated. The algorithm for the software implementation of the proposed method is described too. A lot of numerical examples are included in the paper. The examples confirm the theoretical conclusions. We also have made the comparison between the results obtained by FD-method and results obtained by the powerful software package for solving Sturm-Liouville problems --- SLEIGN2.


Introduction
There is a great scope of the numerical methods for solving eigenvalue problems for differential operators of the second order (see [9]). The known numerical techniques for eigenvalue problems can be divided into two groups: methods based on the direct approximation of the solution of differential equation and methods based on approximation of the coefficients of differential equation (see [7]).
The methods of the first group such as finite difference, finite elements and spectral are extensively treated through both theoretical investigations and well developed soft. The main idea of this approach is replacement of the eigenfunctions by piecewise polynomial functions, which results in the approximation of a differential equation by a system of linear algebraic equations. Because of the different nature of the original differential operator and approximating algebraic operator the numerical solution comes close to the exact one only for the low-indexed eigenvalues and corresponding eigenfunctions, that is defined by the size of the mesh step. Hence, such numerical techniques are effective to find the low-indexed eigenvalues, but not effective to find the high-indexed eigenvalues [7].
The methods of the second group are based on the general idea of the coefficient approximation methods (hereinafter referred to as CAM), also known as the Pruess methods, -to replace the coefficient functions of the problem by simpler approximation and then to solve the approximating problem (see [9]). Firstly this idea was proposed by Kryloff and Bogolioubov (1926) in [4] for the case of piecewise constant approximation. S. Pruess provided a thorough convergence and error analysis of such methods when piecewise polynomial approximation is used. Except for Gordon, who uses linear functions, the earlier references all confine themselves to piecewise constant approximations, and this is the most practical choice: otherwise the approximating problem may be no easier to solve numerically than is the original. It is well known (see [9]) that if the piecewise constant approximation with the maximum step size h is used, the error of the CAM is estimated by O(h). Hence, when we using the CAM, a crucial problem arise: to halve the error we have to double the computational time. This problem in many cases can be solved by using the functional-discrete method (hereinafter referred to as the FD-method) which is the essential generalization of the CAM for solving Sturm-Liouville problems of different types. The algorithm of FD-method was firstly described and justified for linear Sturm-Liouville problem in [6]. Generally speaking, it consists of two main steps: the first one is to find the initial approximation by applying the CAM and the second one is to calculate a sufficient number of corrections to achieve needed accuracy. The second step is substantiated by the homotopy method (also known as perturbation method, continuation method or successive loading method, see, for example, [2]). From that point of view, FD-method is the synthesis of the CAM and homotopy method.
In the present paper we consider the following singular Sturm-Liouville problem Since the problem (1.1),(1.2) has singularities at the points ±1, we can't directly apply the algorithm of the FD-method, proposed in [6]. However, in [5] the modified FDmethod technique for the problems with such singularities have been presented. Let us state briefly the algorithm of the FD-method described in [5]. According to the FDmethod we approximate the eigenvalues and eigenfunctions of the problem (1.1), (1.2) by the following truncated series where the unknown summands λ j n and u the following recurrence problems: (1.4) x ∈ (−1, 1) with the boundary conditions 5) and the orthogonality conditions The initial data λ n (x) which is needed to begin the recurrence process can be found as the solution to the following eigenvalue problem, called the basic problem: The solutions to problem (1.7), (1.8) can be expressed in the following form where P n (x) denotes the Legendre polynomial of the first kind of order n ∈ N (see [3]).
The following statement about the convergence of the FD-method, described above, was proved in [5].
Theorem 1.1 was proved in [5] for the case when q(x) = q(−x) but, using the same technique, it can be proved for the case of an arbitrary function from Q 0 [−1, 1] without considerable difficulties.
If we try to apply the FD-method (1.4) -(1.9) to approximate eigensolution λ n , u n (x) with n = 0, 1, . . . , n 0 − 1 we might find it divergent. The reason of that is the lack of approximation for q(x) in equation (1.7) of the basic problem. To overcome such difficulties we have to use the general scheme of the FD-method which uses the approximation of the function q(x) on [−1, 1] by the piecewise constant functionq(x).
2 The general algorithm of the FD-method for solving the Sturm-Liouville problem with Legendre differential operator:

theoretical justification
To make the general algorithm of the FD-method more easy for understanding, let us consider the following general problem To define the functionq(x), we have to fix some mesh on the interval [−1, 1] : We also require that all points of discontinuity of the function q(x) on the interval [−1, 1] are contained in the setω. There are many possibilities to defineq(x), for example, It is easy to see that if we set t = 1, then problem (2.1), (2.2) would transform to the problem (1.1), (1.2). Assume that the solution to the problem (2.1), (2.2) can be expressed in the series form and (2.6) The unknown summands of the series (2.5) can be found as the solutions to the recurrence problems (j = 0, 1, . . . , m): with the boundary conditions and additional requirement Here square brackets denotes the jump of the function at the point x = x i . Among the problems (2.7) -(2.9) it is worth to emphasize the first one (j = 0): As it was mentioned above, the problem (2.11) -(2.13) is called the basic problem. It is easy to make sure that the differential operator is self-adjoint in the Sobolev space This fact implies that the eigenfunctions u (0) n (x), n = 0, 1, . . . of the problem (2.11) -(2.13) form the complete orthogonal system in W 2,1 (−1, 1) and the corresponding eigenvalues j , when i < j, i, j = 0, 1, 2 . . .) are simple. Suppose that the basic problem is solved and a pare λ n (x) denotes its arbitrary eigensolution.
It is well known that the problems (2.7) -(2.9) for j = 1, 2, . . . , m can be solved if and only if the following equality holds Equality (2.15) together with (2.7) yield us the formula for computing λ (j) n : Taking into account (2.10) we can simplify formula (2.16) in the following way The function u (j) n (x) can be represented as a Fourier series To find the Fourier coefficients let us multiply both sides of the equation (2.7) by u Taking into account boundary conditions (2.8) and using the integration by parts (two times), from (2.19) we obtain And (2.20) immediately lead us to the equality Using (2.21) we can express the formula (2.19) in the following form and Multiplying both sides of the inequality (2.24) by M n q(x) −q(x) ∞,[−1,1] −j and using estimate (2.23), we obtain (2.26) we can represent inequality (2.26) in the following form It is easy to see that the sequence {V i } , defined by the following recurrence formula, Let us consider the generating function f (z) defined by the formula From (2.29) we obtain the equation with respect to f (z) where ρ denotes the convergence radius for the power series (2.31). The solution to equation (2.32), which satisfies the condition f (0) = 1 is And inequality (2.35) yields the following error estimations of the FD-method's convergence rate (2.37) The infinite sum in the right sides of inequalities (2.36) and (2.37) can be easily estimated We have proved the following theorem. Before passing to the next section, let us investigate the question about the asymptotical behavior of parameter M n (2.25) as n tends to +∞. It is easy to see that the basic problem (2.11) -(2.13) is equivalent to the following auxiliary eigenvalue problem when τ = 1. On the other hand, as it follows from the proof of Theorem 1.1 (see [5]), for sufficiently large number n, namely Hence, from equality (2.39) it is easy to obtain x ∈ (−1, 1), τ ∈ [0, 1]. After that, applying to the both sides of equality (2.44) the integral operator Taking into account (2.43) and equality (2.45), it is easy to obtain λ (0) n = n(n + 1) Finally, equality (2.46) yields the following estimation Hence, for the n sufficiently large (see (2.42)) we have the following estimation which describes the asymptotical behavior of the M n : 3 General algorithm of the FD-method: software implementation.
The solution to the basic problem (2.11) -(2.13) can be represented in the following form where and Q ν i (x) denotes the Legendre function of the second kind (see [3]). It is well known that (see, for example, [3]) (2 The sing "+" or "−" can be chosen optionally. To satisfy boundary condition (2.12) we have to set The other constants A i , B i have to be found from the matching conditions (2.13) as the solutions to the following recurrence systems of linear algebraic equations: Let us consider the function where constants A 1 , B 1 are defined by the following recurrence formulas (see (3.5)) And for the computation of the functions u (j) n (x), j = 1, 2, . . . , m it is convenient to use the formula and w n (x) denotes a function, defined by the formula where A i , B i , i = 1, 2, . . . , N can be computed recursively by formulas (3.6) with the initial data Taking into account the properties of the Legendre functions (see [3]), it is easy to verify that ∂K(x, ξ) ∂x and K(x, ξ) is the Cauchy operator for the nonhomogeneous equations (2.7). In general case, the integrals in formulas (3.9) and (3.10) could not be calculated analytically. Moreover, the integrands in this formulas are unbounded at the points ±1. To compute this integrals it is convenient to use numerical methods, for example, the tanh rule (see [8]) and Stenger's formula (see [10] (e −lh sinc /2 + e lh sinc /2 ) 2 (3 , (3.14) (3 The function f (x) is required to be sufficiently smooth on (a, b), see [10].
where z k = a + be h sinc k 1 + e h sinc k , k = −K . . . , K, δ Henceforth we require that the function q(x) is analytical on each interval (x i−1 , x i ), i = 1, 2, . . . , N. Before passing to the algorithm description, we need to introduce the following auxiliary notation n , and r -the depth of the FD-method. output: The arrays of values u input : The arrays of values u (0) n (z i,j ), w n (z i,j ) and F n (ξ)dξ by the Stenger's formula (3.14) n (ξ)dξ by the Stenger's formula (3.14) The described algorithm does not clime to be the most efficient. It is evident that the highest accuracy, that can be achieved by increasing the rank r of the FD-method, is restricted by the accuracy of the quadrature formulas (3.13), (3.14), used in this algorithm. The question about optimal choice of the parameters r and K is still a pressing issue. In the next section the numerical results obtained by the Algorithm 1 are stated. In the first example we apply quadrature formulas (3.13), (3.14) with K = 500 and in the other examples we used K = 350.
4 Numerical examples Example 1. Let us consider the following Sturm-Liouville problem Using the FD-method, described in [5] and the computer algebra system Maple, it is easy to obtain To control the accuracy of the results obtained by the FD-method we use the following functionals Also we compare the FD-method results with eigenvalues obtained by SLAIGN2 (see. [1]) In all examples, presented in this section, we apply the FD-method with uniform mesh on the interval [−1, 1] and N denotes the number of subintervals. The results obtained     by the FD-method are presented in tables 2-4 and figure 1 below. The result obtained by the SLEIGN2 are presented in table 1.
Analyzing the data in tables 2, 3 and 4, we can conclude that the FD-method with N = 1 converges much more slower then the FD-method with N = 3. This fact is in good agreement with the results of Theorem 2.1. Furthermore, as it follows from table 4, the convergence rate of the FD-method increases as the index n of the trial eigenvalue increases. Figure 1 illustrates the exponential nature of the FD-method's convergence.
Example 2. As the second example, let us consider problem (4.1) with The results obtained with SLEIGN2 are presented in table 5 below.  It is worth to emphasize that the problem under consideration do not satisfy the conditions of theorem 2.1. However, from the results presented in table 6 and figure 2 it follows that the FD-method have successfully handled this problem as opposed to SLEIGN2, which gives results with essential errors (see the leftmost colon in the table 6).  We find out that SLAIGN2 does not handle this problem. But the FD-method does. The results obtained by the FD-method are presented in table 7 and figure 3. As before, figure  3 confirms the exponential convergence rate of the method.

Conclusions
In the present paper we construct and theoretically justified the generalized algorithm of the FD-method for solving the Sturm-Liouville problem for differential equation of the second order (1.1), (1.2) with piecewise continuous functional coefficient q(x). As it follows from Theorem 2.1, the generalized FD-method, which uses the piecewise constant approximation of the function q(x), can be applied for the approximation of eigenvalues and eigenfunctions with any nonnegative index n. The convergence rate of the method can be increased by decreasing the value q(x) −q(x) ∞, [−1,1] . In the case whenq(x) ≡ 0 (this case was considered in [5]) the FD-method (if converges) allows us to calculate the approximation to the eigensolution analytically. But in general case, whenq(x) = 0, the analytical calculations are almost always impossible and it is necessary to use numerical integration methods, such as sinc quadratures and Stenger's formula (see [8], [10]).
The problems considered in examples 2 and 3 do not satisfy the conditions of Theorem 1.1 because of the unboundedness of the function q(x) on [−1, 1]. However, as it was shown in the mentioned examples, the method successfully converges whereas the well known in the mathematical world package SLAIGN2 either gives not more then three correct numbers after decimal point (example 2) or can not handle the problem at all (example 3). This examples indicate that the FD-method has a considerable potential which are to be investigated in further mathematical works.