Numerical solution of the Sturm-Liouville eigenvalue problem — matslise (MATLAB) and pyslise (python) solvers
In our previous article, we introduced some basic ideas about direct and inverse eigenvalue problems. Here, we will discuss in more detail about the Sturm-Liouville eigenvalue problem and the basic properties of its spectrum. But since we are also interested in the numerical solution, we will present two useful eigenvalue solvers in MATLAB and python, suitable for our boundary value problem.
The Sturm-Liouville eigenvalue problem — properties
We consider the two-point boundary value problem, defined in the unit interval, of the following general form:
together with the boundary conditions:
where we assume that the functions p, p’, q and r are continuous on the interval 0 ≤ x ≤ 1 and that p(x) > 0 and r(x) > 0 at all points in 0 ≤ x ≤ 1. We also assume that not both alphas and betas are equal to zero. In this case the boundary value problem is said to be regular. We can equivalently write the above as a differential operator eigenvalue problem for D:
Each pair of (λ,y) is the corresponding eigenpair of the operator D.
As we explained in the previous article, the direct eigenvalue consists on finding the eigenpairs, for known functions p, q and r. But DO these eigenpairs EXIST? How many are they? Are they real, complex, positive or negative? The following result answers the above questions:
1. Let y and ψ are two eigenfunctions of the Sturm–Liouville problem corresponding to eigenvalues λ and μ, respectively. If λ is not equal with μ, then the eigenfunctions are orthogonal. That is:
2. The eigenvalues form an infinite real sequence and can be ordered according to increasing magnitude so that:
and furthermore
The proof of the above is a consequence of the Hilbert-Schmidt theory for compact and self-adjoint operators and is out of the scope of this article.
Enough with theory (for the time being). Lets see how to use two practical solvers to calculate these eigenpairs!
Numerical solution of the Sturm-Liouville eigenvalue problem
Since boundary value problems and especially Sturm-Liouville eigenvalue problem is of high interest in applications, a lot of numerical methods have been developed to solve these problems.
1. Using matslise in MATLAB
Matslise 2.0 is a MATLAB package which can be used to compute the eigenvalues and eigenfunctions of regular and singular self-adjoint Sturm-Liouville boundary value problems. It was created by Veerle Ledoux and Marnix van Daele.
As an example, we consider the following Sturm-Liouville eigenvalue problem:
with Dirichlet boundary conditions y(0)=y(π/2)=0. The potential q(x) is shown in the following figure:
The code using the matslise built-in slp class:
function TestSturmLiouville
%illustration of the use of the MATSLISE-functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Sturm-Liouvile problem, with functions p,q,w and dirichet boundary conditions
%%y(0)=y(pi/2)=0
o=slp(@p,@q,@w,0,pi/2,1,0,1,0);
%display classification information:
disp(classificationInfo(o))
%compute eigenvalues with indices between 0 and 10:
[E,meshData] = o.computeEigenvalues(0,9,1e-12,true);
%first eigenfunction evaluated in the mesh points:
[x,y,yp] = o.computeEigenfunction(E.eigenvalues(1),meshData);
figure
plot(x,y,x,yp)
xlabel('x')
legend('y','yp')
%the first 3 eigenfunctions evaluated in a range of user-specified points:
[x,y,yp] = o.computeEigenfunction(E.eigenvalues(1),meshData,linspace(0,pi/2,500));
[x,y2,yp] = o.computeEigenfunction(E.eigenvalues(2),meshData,linspace(0,pi/2,500));
[x,y3,yp] = o.computeEigenfunction(E.eigenvalues(3),meshData,linspace(0,pi/2,500));
figure
%plot(x,y,x,yp)
plot(x,y,x,y2,x,y3)
legend('y_{1}','y_{2}','y_{3}')
xlabel('x')
axis tight
%legend('y','yp')
plotMesh(o,meshData);
disp('The first 10 eigenvalues:')
line='-----------------------------------------';
disp(line)
disp(' k E_k estimated error')
disp(line)
for i=1:length(E.eigenvalues)
fprintf('%4.0f %16.12f %+5.2e\n', E.indices(i),E.eigenvalues(i), E.errors(i))
end
disp(line)
disp('---------- orhtogonality test ---------- ')
int1 = trapz(y.*y2);
int2 = trapz(y.*y3);
int3 = trapz(y2.*y3);
fprintf('integral of the y1 and y2 %s \n', int1);
fprintf('integral of the y1 and y3 %s \n', int2);
fprintf('integral of the y2 and y3 %s \n', int3);
disp(line)
end
function r=p(x)
r=1;
end
function r=q(x)
r=(x-1)^2 + exp(2*x);
end
function r=w(x)
r=1;
In the output, we calculate the first 10 eigenvalues. We verify that they are an increasing sequence of positive numbers. We also verify the orthogonality between the first 3 eigenfunctions.
>> TestSturmLiouville
The spectrum is simple, purely discrete and bounded below.
There is no continuous spectrum.
At endpoint A
The problem is regular.
At endpoint B
The problem is regular.
The first 10 eigenvalues:
-----------------------------------------
k E_k estimated error
-----------------------------------------
0 9.125618302908 -8.08e-13
1 22.845873279784 -1.39e-12
2 43.207474751000 -7.25e-13
3 71.278107445561 +1.38e-12
4 107.294741135566 +4.65e-12
5 151.299367990309 +3.44e-12
6 203.300719604859 +2.84e-14
7 263.301043330762 -3.41e-13
8 331.301022502680 -4.77e-12
9 407.300889778476 -2.67e-12
-----------------------------------------
---------- orhtogonality test ----------
integral of the y1 and y2 7.565334e-13
integral of the y1 and y3 4.835973e-12
integral of the y2 and y3 1.979751e-11
-----------------------------------------
Finally, matslise includes a nice GUI for those that don’t like coding:
2. Using pyslise in Python
Based on matslise, Toon Baeyens and Marnix Van Daele created a similar python package, called pyslice. It is a package (written in C++) to solve Schrödinger equations of the form:
i.e. Sturm Liouville equations with p(x)=r(x)=1.
We install pyslise
pip install pysliseand test the above example to get the first 10 eigenvalues:
from pyslise import Pyslise
from math import pi, exp
import numpy as np
import math
import matplotlib.pyplot as plt
problem = Pyslise(lambda x: (x-1)**2 + exp(2*x), 0, pi/2, tolerance=1e-12)
left = (0, 1)
right = (0, 1)
eigenvalues = problem.eigenvaluesByIndex(0, 10, left, right)
xs = np.linspace(0, pi/2, 300)
print('------------- the first 10 eigenvalues -----------')
print('index eigenvalue error')
for index, E in eigenvalues:
error = problem.eigenvalueError(E, left, right)
print(f'{index:>5} {E:>11.5f} {error:>9.1e}')and the results follow:
------------- the first 10 eigenvalues -----------
index eigenvalue error
0 9.12562 7.5e-13
1 22.84587 1.9e-12
2 43.20747 4.3e-12
3 71.27811 1.1e-12
4 107.29474 7.4e-12
5 151.29937 1.3e-11
6 203.30072 1.1e-11
7 263.30104 3.2e-12
8 331.30102 5.7e-12
9 407.30089 1.2e-11
Process finished with exit code 0By comparing the outputs, we verify that both solvers give very close results.
Pyslise also includes a simple GUI.
So, taking advantage of the nice properties of the Sturm-Liouville eigenvalue problems, both matslise and pyslise are very useful in solving the direct problems.
References
- Boyce W.E. and Diprima R.C. (2001), Elementary Differential Equations and Boundary Value Problems, 7th edn, John Wiley & Sons Inc.
- Renardy M. and Rogers R. (2004), An Introduction to Partial Differential Equations, Springer.
- Ledoux V. and Van Daele M. (2016), Matslise 2.0: A Matlab toolbox for Sturm-Liouville computations, ACM Trans. Math. Softw.
- https://sourceforge.net/projects/matslise
