The SturmLiouville eigenvalue problem  a numerical solution using the Control Volume Method
Jarosław Siedlecki
,Mariusz Ciesielski
,Tomasz Błaszczyk
Journal of Applied Mathematics and Computational Mechanics 
Download Full Text 
THE STURMLIOUVILLE EIGENVALUE PROBLEM  A NUMERICAL SOLUTION USING THE CONTROL VOLUME METHOD
Jarosław Siedlecki^{1}, Mariusz Ciesielski^{2}, Tomasz Błaszczyk^{1}
^{1}Institute
of Mathematics, Czestochowa University of Technology
Częstochowa, Poland
^{2}Institute of Computer and Information Sciences
Czestochowa University of Technology, Częstochowa, Poland
jaroslaw.siedlecki@im.pcz.pl,
mariusz.ciesielski@icis.pcz.pl, tomasz.blaszczyk@im.pcz.pl
Abstract. The solution of the 1D SturmLiouville problem using the Control Volume Method is discussed. The second order linear differential equation with homogeneous boundary conditions is discretized and converted to the system of linear algebraic equations. The matrix associated with this system is tridiagonal and eigenvalues of this system are an approximation of the real eigenvalues of the boundary value problem. The numerical results of the eigenvalues for various cases and the experimental rate of convergence are presented.
Keywords: SturmLiouville problem, eigenvalues, numerical methods, Control Volume Method
1. Introduction
This paper is concerned with the computation of eigenvalues of regular eigenvalue problems occurring in ordinary differential equations. The SturmLiouville problem arises within in many areas of science, engineering and applied mathematics. It has been studied for more than two decades. Many physical, biological and chemical processes are described using models based on the SturmLiouville equations.
The SturmLiouville problem appears directly as the eigenvalue problem in a onedimensional space. It also arises when linear partial differential equations are separable in a certain coordinate system. For a more detailed study of the integer order SturmLiouville theory, we refer the reader to [15].
The SturmLiouville problem can be solved by using either analytical or numerical methods. One of the most common approaches to a numerical solution of the considered problem is the finite difference method [3, 68] where each derivative is discretized at each grid point with an adequate difference scheme. Apart from the finite difference method, several analytical ones, such as the variational or decomposition methods, are proposed to find an approximate solution.
2. Statement of the problem
We consider the following problem defined on the bounded interval x [a, b]
(1) 
(2) 
where p(x) > 0, dp/dx and q(x) are continuous, w(x) > 0 on [a, b], and . The above problem is called the regular SturmLiouville Problem (SLP). The solution to SLP consists of a pair l and y, where l is a constant  called an eigenvalue, while y is a nontrivial (nonzero) function  called an eigenfunction, and together they satisfy the given SLP. For each SLP, the eigenvalues form an infinite increasing sequence: l_{1} < l_{2} < l_{3} < … and lim_{k}_{®}_{¥}l_{k} = ¥.
For arbitrary choices of the functions p(x), q(x) and w(x) in Eq. (1), the computation of the exact values of the eigenvalues for which SLP (1)(2) has a nontrivial eigensolution y(x) which is very complicated or it is practically impossible to determine. Numerical methods should, therefore, be used for computing the approximate values of l. In many cases, the importance of a numerical approximation to the SLP described by a differential eigenvalue problem is to reduce the problem to that of solving the eigenvalue problem of a matrix equation (an algebraic problem).
In this paper, we apply the Control Volume Method (also known as the Finite Volume Method) to compute the eigenvalues of the SturmLiouville problem numerically.
3. Control Volume Method
In this numerical method [9] the considered domain of SLP: x [a, b] is divided into N control volumes W_{i} for i = 1,…,N with the central nodes x_{i}. The mesh is presented in Figure 1. The auxiliary nodes x_{i} = a + i Dx, for i = 0,…,N and Dx = = (b  a) / N have also been introduced. Then, the positions of central nodes are the following: x_{i } = a + (i  0.5)Dx, i = 1,…,N.
Fig. 1. Mesh of control volumes
The integration of Eq. (1) with respect to volume W_{i} leads to
(3) 
or written in the form (assuming that W_{i}: [x_{i}_{1}, x_{i}])
(4) 
All the components in Eq. (4) can be approximated as follows:
(5) 
(6) 
(7) 
The values of y_{a} = y(a) and y_{b} = y(b) are determined on the basis of approximations of the boundary conditions (2)
(8) 
and hence, these values are equal to
(9) 
Substituting (9) into (8), the following approximation of (7) has been obtained
(10) 
where
,  (11) 
For example, in the case of the Dirichlet boundary conditions at boundaries x = a (i.e. a_{1} = 1, a_{2} = 0) and/or x = b (i.e. b_{1} = 1, b_{2} = 0), the coefficients g_{a} and g_{b} are equal to 2, and in the case of the Neumann boundary conditions at x = a (a_{1} = 0, a_{2} = 1) and/or x = b (b_{1} = 0, b_{2} = 1), these coefficients take the values of 0, respectively.
After substitution of (5), (6) and (10) into (4), the following system of the discrete equations for every control volume:
· for W_{1}:
(12) 
· for W_{i}, for i = 2,…,N  1:
(13) 
· for W_{N}:
(14) 
is obtained.
This system can be also written in the matrix form as
(15) 
where
(16) 
(17) 
(18) 
and . Then, SLP (1)(2) is equivalent to the matrix eigenvalue problem (15). The ordered eigenvalues of (15) are denoted by l_{i}^{(N)}, i = 1,…,N.
In order to evaluate eigenvalues of large matrix eigenvalue problem (15), the numerical methods implemented in mathematical software can be used.
3.1. Particular case
Assuming the functions p(x) = 1, q(x) = 0 and w(x) = 1 in the SLP problem, then matrices P, Q and W are reduced to the following forms: Q = O (Zero matrix), W = I (Identity matrix) and
(19) 
In the case of the mixed Dirichlet and/or Neumann boundary conditions at both boundaries, one can find in literature (e.g. [4]) the explicit formulas for the ith eigenvalue l_{i} of the SLP (1)(2) and they are presented in Table 1. For the discrete case, the eigenvalues of the matrix P can be determined in an analytical way. The eigenvalue problems of tridiagonal matrices (in a similar form as the matrix P) are considered in [10, 11]. On the basis of these results, the eigenvalues of the matrix P: l_{i}^{(N)}, i = 1,…,N are adopted to the analysed problem for four mixed boundary conditions and they are also given in Table 1.
Table 1
Eigenvalues of Eq. (1) for functions p(x) = 1, q(x) = 0 and w(x) = 1 and four mixed cases of boundary conditions (explicit formulas for l_{i}, i = 1,…,∞ and discrete one for l_{i}^{(N)}, i = 1,…,N )

The Dirichlet B.C. at x = b g_{b} = 2 (for b_{1} = 1, b_{2} = 0) 
The Neumann B.C. at x = b g_{b} = 0 (for b_{1} = 0, b_{2} = 1) 
The Dirichlet B.C. at x = a
g_{a} = 2 (for a_{1} = 1, a_{2} = 0) 

The Neumann B.C. at x = a
g_{a} = 0 (for a_{1} = 0, a_{2} = 1) 
For more complicated cases such as nonconstant functions p(x), q(x), w(x) occurring in Eq. (1), the eigenvalues of system (15) should be determined in a numerical way (i.e. using mathematical software  here the Maple is used).
3.2. Error of numerical approximation of the eigenvalues
Theorem 1. For the SturmLiouville problem (1)(2) there exists a constant C such that the ith exact eigenvalue l_{i} and the approximate eigenvalue l_{i}^{(N)} satisfy the following relation
(20) 
Proof: For the presented formulas in Table 1 (in these cases, the exact and discrete eigenvalues are known), one can estimate the error of approximation of the discrete eigenvalues. Let us start from the Taylor series expansion for the function sin^{2}(x) at about point x = 0
(21) 
Next, if we take into account the first two terms in the Taylor series (21), then we estimate the following expression as
(22) 
For the case of the DirichletDirichlet boundary conditions, the error is evaluated by using the estimation (22) in the following way:
(23) 
In the other cases (from Table 1), the results are similar.
Another method of error estimation is based on the investigation of the Experimental Rate of Convergence (ERC). The total error in the estimate of the eigenvalues is composed of both the error resulting from the discretization of the equation and the error of the numerical algorithm for finding eigenvalues of system (15). Here, we assume that the error is
(24) 
where the parameters r and s are to be determined experimentally.
If the ith exact eigenvalue l to SLP is known, then the parameters r and s can be determined using the following formulas for the ERC for variable values of N:
(25) 
(26) 
Whereas, if the exact eigenvalue is unknown then, we determine the parameter r from the following formula
(27) 
The parameter s can be estimated using (26) and assuming that l_{i} is numerically determined for sufficiently high value of N.
4. Example of numerical simulations
In tests of verification of the numerical solutions, three cases are taken into account:
· Example 1: p(x) = 1, q(x) = 0, w(x) = (1 + x)^{‒2}, a_{1} = 1, a_{2} = 0, b_{1} = 1, b_{2} = 0.
· Example 2: p(x) = (x + 1)^{2}, q(x) = x^{2}  2, w(x) = exp(x), a_{1} = 1, a_{2} = 0, b_{1} = 0, b_{2} = 1.
· Example 3: p(x) = 2 + sin(2px), q(x) = 10, w(x) = 1 + sqrt(x), a_{1} = 1, a_{2} = 0, b_{1} = 5, b_{2} = 1.
In all examples, the values of a = 0 and b = 1 have been assumed.
In the case of Example 1, the exact eigenvalues are given by [12], while for the remaining cases the exact eigenvalues are unknown.
In Tables 2, 4 and 5, the numerical values of the first 8 eigenvalues for different values of N and the calculated ERC_{r} ((25) or (27)) for all examples are presented, respectively. In addition, in Table 3, the calculated values of ERC_{s} (26) for Example 1 are shown.
Table 2
Eigenvalues and ERC_{r} for Example 1
N 
l_{1}^{(N)} 
ERC_{r} 
l_{2}^{(N)} 
ERC_{r} 
l_{3}^{(N)} 
ERC_{r} 
l_{4}^{(N)} 
ERC_{r} 
100 
20.7904297 
82.3888961 
184.976919 

328.440408 


200 
20.7918238 
2.0001 
82.4115895 
2.0000 
185.092175 
1.9999 
328.805046 
1.9998 
400 
20.7921723 
2.0000 
82.4172627 
2.0000 
185.120991 
2.0000 
328.896222 
1.9999 
800 
20.7922594 
2.0000 
82.4186811 
2.0000 
185.128195 
2.0000 
328.919017 
2.0000 
1600 
20.7922812 
2.0067 
82.4190356 
1.9995 
185.129996 
2.0001 
328.924716 
2.0000 
3200 
20.7922866 
1.9994 
82.4191243 
2.0007 
185.130446 
2.0000 
328.926140 
2.0000 
anal. 
20.7922885 
82.4191538 
185.130596 

328.926615 


N 
l_{5}^{(N)} 
ERC_{r} 
l_{6}^{(N)} 
ERC_{r} 
l_{7}^{(N)} 
ERC_{r} 
l_{8}^{(N)} 
ERC_{r} 
100 
512.619722 
737.309749 
1002.26001 

1307.17479 


200 
513.510262 
1.9996 
739.156478 
1.9994 
1005.68095 
1.9992 
1313.00956 
1.9989 
400 
513.732969 
1.9999 
739.618392 
1.9999 
1006.53680 
1.9998 
1314.46964 
1.9997 
800 
513.788651 
2.0000 
739.733885 
2.0000 
1006.75080 
1.9999 
1314.83475 
1.9999 
1600 
513.802571 
2.0000 
739.762760 
2.0000 
1006.80430 
2.0000 
1314.92603 
2.0000 
3200 
513.806051 
2.0000 
739.769978 
2.0000 
1006.81768 
2.0000 
1314.94885 
2.0000 
anal. 
513.807211 
739.772384 
1006.82213 

1314.95646 

Table 3
ERC_{s} for Example 1
N 
ERC_{s}(N,1) 
ERC_{s}(N,2) 
ERC_{s}(N,3) 
ERC_{s}(N,4) 
ERC_{s}(N,5) 
ERC_{s}(N,6) 
ERC_{s}(N,7) 
ERC_{s}(N,8) 
100 
4.0249 
4.008 
4.0037 
4.0017 
4.0006 
3.9997 
3.9989 
3.9982 
200 
4.0250 
4.0082 
4.0040 
4.0023 
4.0014 
4.0008 
4.0004 
4.0000 
400 
4.0250 
4.0082 
4.0041 
4.0024 
4.0016 
4.0011 
4.0007 
4.0005 
800 
4.0250 
4.0082 
4.0041 
4.0024 
4.0016 
4.0011 
4.0008 
4.0006 
1600 
4.0322 
4.0073 
4.0043 
4.0024 
4.0016 
4.0012 
4.0009 
4.0007 
3200 
4.0309 
4.0086 
4.0042 
4.0024 
4.0017 
4.0011 
4.0008 
4.0007 
Table 4
Eigenvalues and ERC_{r} for Example 2
N 
l_{1}^{(N)} 
ERC_{r} 
l_{2}^{(N)} 
ERC_{r} 
l_{3}^{(N)} 
ERC_{r} 
l_{4}^{(N)} 
ERC_{r} 
100 
1.17035201 
 
26.8573037 
 
78.5053965 
 
155.915777 
 
200 
1.17045757 
2.0000 
26.8618516 
1.9999 
78.5381317 
1.9998 
156.039051 
1.9995 
400 
1.17048396 
2.0000 
26.8629886 
2.0000 
78.5463168 
1.9999 
156.069879 
1.9999 
800 
1.17049056 
2.0006 
26.8632729 
2.0001 
78.5483631 
2.0000 
156.077587 
2.0000 
1600 
1.17049221 
1.9924 
26.8633440 
2.0000 
78.5488747 
1.9999 
156.079514 
2.0000 
3200 
1.17049262 
 
26.8633617 
 
78.5490026 
 
156.079996 
 
N 
l_{5}^{(N)} 
ERC_{r} 
l_{6}^{(N)} 
ERC_{r} 
l_{7}^{(N)} 
ERC_{r} 
l_{8}^{(N)} 
ERC_{r} 
100 
259.010064 
 
387.685285 
 
541.813108 
 
721.239858 
 
200 
259.344073 
1.9992 
388.427244 
1.9988 
543.256336 
1.9982 
723.792747 
1.9976 
400 
259.427622 
1.9998 
388.612893 
1.9997 
543.617583 
1.9996 
724.432011 
1.9994 
800 
259.448513 
1.9999 
388.659316 
1.9999 
543.707922 
1.9999 
724.591893 
1.9999 
1600 
259.453735 
2.0000 
388.670922 
2.0000 
543.730508 
2.0000 
724.631867 
2.0000 
3200 
259.455041 
 
388.673823 
 
543.736155 
 
724.641861 
 
Table 5
Eigenvalues and ERC_{r} for Example 3
N 
l_{1}^{(N)} 
ERC_{r} 
l_{2}^{(N)} 
ERC_{r} 
l_{3}^{(N)} 
ERC_{r} 
l_{4}^{(N)} 
ERC_{r} 
100 
2.89692971 
 
24.0701022 
 
66.8696550 
 
130.023860 
 
200 
2.89762086 
1.9999 
24.0778646 
1.9998 
66.9119064 
1.9997 
130.164588 
1.9996 
400 
2.89779365 
2.0000 
24.0798055 
1.9999 
66.9224715 
1.9999 
130.199780 
1.9999 
800 
2.89783685 
1.9999 
24.0802907 
2.0000 
66.9251129 
2.0000 
130.208579 
2.0000 
1600 
2.89784765 
1.9993 
24.0804120 
1.9999 
66.9257733 
2.0000 
130.210778 
2.0000 
3200 
2.89785036 
 
24.0804424 
 
66.9259384 
 
130.211328 
 
N 
l_{5}^{(N)} 
ERC_{r} 
l_{6}^{(N)} 
ERC_{r} 
l_{7}^{(N)} 
ERC_{r} 
l_{8}^{(N)} 
ERC_{r} 
100 
213.781678 
 
318.196036 
 
443.187196 
 
588.625456 
 
200 
214.140243 
1.9994 
318.965739 
1.9992 
444.653772 
1.9990 
591.185500 
1.9998 
400 
214.229921 
1.9998 
319.158267 
1.9998 
445.020663 
1.9997 
591.826047 
1.9997 
800 
214.252343 
2.0000 
319.206405 
1.9999 
445.112402 
1.9999 
591.986219 
1.9999 
1600 
214.257948 
2.0000 
319.218441 
2.0000 
445.135338 
2.0000 
592.026264 
2.0000 
3200 
214.259350 
 
319.221449 
 
445.141072 
 
592.036275 
 
The analysis of the results presented in Tables 25 indicates that the rate r (ERC_{r}) is close to 2, while the rate s (ERC_{s}) is close to 4. Thus we can confirm that the relationship (20) is satisfied. The errors in the approximation of eigenvalues increase rapidly as the index of the eigenvalue grows.
5. Conclusions
In this paper, the new approach based on the control volume method for finding the eigenvalues of the SturmLiouville problem was discussed. The continuous problem described by the differential equation with the adequate boundary conditions was converted to the corresponding discrete one. The rate of convergence of the proposed numerical scheme is order 2. The presented results of the approximation of eigenvalues are in close agreement with the results obtained in an analytical way or in the mathematical software field. In the future, the presented approach can be extended to apply high order of accuracy difference schemes for approximating eigenvalues of the SturmLiouville problem and can be applied to the fractional SturmLiouville problem which is related to the corresponding fractional Euler‑Lagrangian equation [13].
References
[1] Agarwal R.P., O’Regan D., An Introduction to Ordinary Differential Equations, Springer, New York 2008.
[2] Atkinson F.V., Discrete and Continuous Boundary Value Problems, Academic Press, New York, London 1964.
[3] Pryce J.D., Numerical Solution of SturmLiouville Problems, Oxford Univ. Press, London 1993.
[4] Zaitsev V.F., Polyanin A.D., Handbook of Exact Solutions for Ordinary Differential Equations, CRC Press, New York 1995.
[5] Pruess S., Estimating the eigenvalues of SturmLiouville problems by approximating the differential equation, SIAM J. Numer. Anal. 1973, 10, 5568.
[6] Aceto L., Ghelardoni P., Magherini C., Boundary value methods as an extension of Numerov’s method for SturmLiouville eigenvalue estimates, Applied Numerical Mathematics 2009, 59 (7), 16441656.
[7] Amodio P., Settanni G., A matrix method for the solution of SturmLiouville problems, Journal of Numerical Analysis, Industrial and Applied Mathematics (JNAIAM) 2011, 6 (12), 113.
[8] Ascher U.M., Numerical Methods for Evolutionary Differential Equations, SIAM, 2008.
[9] Ascher U.M., Mattheij R.M.M., Russell R.D., Numerical Solution of Boundary Value Problems for ODEs, Classics in Applied Mathematics 13, SIAM, Philadelphia 1995.
[10] Elliott J.F., The characteristic roots of certain real symmetric matrices, Master’s Thesis, University of Tennessee, 1953.
[11] Gregory R.T., Karney D., A Collection of Matrices for Testing Computational Algorithm, WileyInterscience, 1969.
[12] Akulenko L.D., Nesterov S.V., HighPrecision Methods in Eigenvalue Problems and Their Applications (series: Differential and Integral Equations and Their Applications), Chapman and Hall/CRC, 2004.
[13] Ciesielski M., Blaszczyk T., Numerical solution of nonhomogenous fractional oscillator equation in integral form, Journal of Theoretical and Applied Mechanics 2015, 53 (4), 959968.
[14] Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P., Numerical Recipes: The Art of Scientific Computing (3rd ed.), Cambridge University Press, New York 2007.