The Sturm-Liouville 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 STURM-LIOUVILLE EIGENVALUE PROBLEM - A NUMERICAL SOLUTION USING THE CONTROL VOLUME METHOD
Jarosław Siedlecki1, Mariusz Ciesielski2, Tomasz Błaszczyk1
1Institute
of Mathematics, Czestochowa University of Technology
Częstochowa, Poland
2Institute 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 Sturm-Liouville 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: Sturm-Liouville 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 Sturm-Liouville 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 Sturm-Liouville equations.
The Sturm-Liouville problem appears directly as the eigenvalue problem in a one-dimensional 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 Sturm-Liouville theory, we refer the reader to [1-5].
The Sturm-Liouville 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, 6-8] 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 Sturm-Liouville 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: l1 < l2 < l3 < … and limk®¥lk = ¥.
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 Sturm-Liouville 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 Wi for i = 1,…,N with the central nodes xi. The mesh is presented in Figure 1. The auxiliary nodes xi = 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: xi = a + (i - 0.5)Dx, i = 1,…,N.
Fig. 1. Mesh of control volumes
The integration of Eq. (1) with respect to volume Wi leads to
(3) |
or written in the form (assuming that Wi: [xi-1, xi])
(4) |
All the components in Eq. (4) can be approximated as follows:
(5) |
(6) |
(7) |
The values of ya = y(a) and yb = 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. a1 = 1, a2 = 0) and/or x = b (i.e. b1 = 1, b2 = 0), the coefficients ga and gb are equal to 2, and in the case of the Neumann boundary conditions at x = a (a1 = 0, a2 = 1) and/or x = b (b1 = 0, b2 = 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 W1:
(12) |
· for Wi, for i = 2,…,N - 1:
(13) |
· for WN:
(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 li(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 i-th eigenvalue li 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: li(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 li, i = 1,…,∞ and discrete one for li(N), i = 1,…,N )
|
The Dirichlet B.C. at x = b gb = 2 (for b1 = 1, b2 = 0) |
The Neumann B.C. at x = b gb = 0 (for b1 = 0, b2 = 1) |
The Dirichlet B.C. at x = a
ga = 2 (for a1 = 1, a2 = 0) |
||
The Neumann B.C. at x = a
ga = 0 (for a1 = 0, a2 = 1) |
For more complicated cases such as non-constant 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 Sturm-Liouville problem (1)-(2) there exists a constant C such that the i-th exact eigenvalue li and the approximate eigenvalue li(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 sin2(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 Dirichlet-Dirichlet 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 i-th 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 li 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, a1 = 1, a2 = 0, b1 = 1, b2 = 0.
· Example 2: p(x) = (x + 1)2, q(x) = x2 - 2, w(x) = exp(x), a1 = 1, a2 = 0, b1 = 0, b2 = 1.
· Example 3: p(x) = 2 + sin(2px), q(x) = -10, w(x) = 1 + sqrt(x), a1 = 1, a2 = 0, b1 = 5, b2 = 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 ERCr ((25) or (27)) for all examples are presented, respectively. In addition, in Table 3, the calculated values of ERCs (26) for Example 1 are shown.
Table 2
Eigenvalues and ERCr for Example 1
N |
l1(N) |
ERCr |
l2(N) |
ERCr |
l3(N) |
ERCr |
l4(N) |
ERCr |
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 |
l5(N) |
ERCr |
l6(N) |
ERCr |
l7(N) |
ERCr |
l8(N) |
ERCr |
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
ERCs for Example 1
N |
ERCs(N,1) |
ERCs(N,2) |
ERCs(N,3) |
ERCs(N,4) |
ERCs(N,5) |
ERCs(N,6) |
ERCs(N,7) |
ERCs(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 ERCr for Example 2
N |
l1(N) |
ERCr |
l2(N) |
ERCr |
l3(N) |
ERCr |
l4(N) |
ERCr |
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 |
l5(N) |
ERCr |
l6(N) |
ERCr |
l7(N) |
ERCr |
l8(N) |
ERCr |
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 ERCr for Example 3
N |
l1(N) |
ERCr |
l2(N) |
ERCr |
l3(N) |
ERCr |
l4(N) |
ERCr |
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 |
l5(N) |
ERCr |
l6(N) |
ERCr |
l7(N) |
ERCr |
l8(N) |
ERCr |
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 2-5 indicates that the rate r (ERCr) is close to 2, while the rate s (ERCs) 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 Sturm-Liouville 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 Sturm-Liouville problem and can be applied to the fractional Sturm-Liouville 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 Sturm-Liouville 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 Sturm-Liouville problems by approximating the differential equation, SIAM J. Numer. Anal. 1973, 10, 55-68.
[6] Aceto L., Ghelardoni P., Magherini C., Boundary value methods as an extension of Numerov’s method for Sturm-Liouville eigenvalue estimates, Applied Numerical Mathematics 2009, 59 (7), 1644-1656.
[7] Amodio P., Settanni G., A matrix method for the solution of Sturm-Liouville problems, Journal of Numerical Analysis, Industrial and Applied Mathematics (JNAIAM) 2011, 6 (1-2), 1-13.
[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, Wiley-Interscience, 1969.
[12] Akulenko L.D., Nesterov S.V., High-Precision 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 non-homogenous fractional oscillator equation in integral form, Journal of Theoretical and Applied Mechanics 2015, 53 (4), 959-968.
[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.