#############################################
Note* : I have not made changes to this code since the end of my PhD (December 2016). There are no guarantees.
#############################################
The purpose of this Julia package is to provide an efficient fast general purpose differential eigenvalue solver package, supporting the canonical interval, and semi-infinite and infinite domains for Sturm-Liouville problems. The following algorithm utilizes the Double Exponential Sinc-Collocation method. This package allows the user to consider other domains by declaring a new instance of the SincFun type Domain.
The primary function of this module computes the eigenvalues of a general Sturm-Liouville problem of the form:
P1 : (-D^2 + q(x) ) u(x) = λ ρ(x) u(x), a < x < b, with u(a) = u(b) = 0,
where -∞ ≦ a < b ≦ ∞ .
In the problem P1,
1. D is the differential operator, an
2. q(x) is a continuous bounded function defined on (a,b)
3. ρ(x) is a continuous positive function defined on (a,b)
4. u(x) are the eigenfunction corresponding to the eigenvalues λ.In this algorithm, we use the transformation developed by Eggert et al. in reference [1]. This transformation results in a symmetric generalized eigenvalue problem defined as:
P2 : (-D^2 + qtilde(x))v(x) = λ ρtilde(x) v(x), -∞ < x < ∞, with v(±∞) = 0
In the problem P2,
1. D is the differential operator, an
2. qtilde(x) is the resulting transformed function defined on (-∞,∞)
3. ρtilde(x) is the resulting transformed function defined on (-∞,∞)
4. v(x) are the eigenfunction corresponding to the eigenvalues λ.
v(x) now have double expoenential decay at both infinities.See reference [2] for more details of the form of qtilde(x) and ρtilde(x).
The type Domain is used to select from the conformal maps depending on the domain of the problem and the decay rate of the solution of P1.
For S-L problems on a finite domain I=(a,b) with algebraic decay at the endpoints:
Domain = Finite(a,b).
For S-L problems on a infinite domain I=(-∞,∞) with algebraic decay at the endpoints:
Domain = Infinite1{Float64}().
For S-L problems on a infinite domain I=(-∞,∞) with single-exponential decay at the endpoints:
Domain = Infinite2{Float64}().
For S-L problems on a semi-infinite domain I=(0,∞) with single-exponential decay at infinity and algebraic decay at 0:
Domain = SemiInfinite1{Float64}().
To use this package, once simply writes:
using SincFun, DESincEigThe main function of this package is SincEigen. This function computes the eigenvalues of the the Sturm-Liouville (P1) using the double exponential Sinc-Collocation method.
For optimal results, it is worth performing an asymptotic analysis of the solution of P2.
The resulting analysis will lead to the following bounds for constants βL,βR,γL,γR>0:
|v(x)| < AL exp(-βL exp ( γL |x| ) ) , on (-∞,0)
|v(x)| < AR exp(-βR exp ( γR |x| ) ) , on (0, ∞)The parameters βL, βR, γL, γR are used in the calculation of the mesh size, h, for the DESCM.
Another important parameter involved in the calculations of the mesh size for the DESCM is the width, d, of the strip D_{d}.
Let S denote the set of complex singularities of the functions qtilde(x) and ρtilde(x):
S = { z ∈ C : qtilde(z) or ρtilde(z) does not exist.}
let s denote the positive imaginary part of nearest singularity to the real axis:
s = min | Im{S} |,
then d = min{ π/2max{γL,γR} , s }.
Input:
Necessary parameters
1. q(x):: Function, The function in P1
2. ρ(x):: Function, The function in P1
3. domain:: Domain, Finite(a,b), Infinite1{Float64}(), Infinite2{Float64}() or SemiInfinite1{Float64}()
4. βopt:: Vector{T}, [βL,βR]
5. γopt:: Vector{T}, [γL,γR]
6. d:: Number, min{ π/2max{γL,γR} , s }The result from using the main function SincEigen is a DESincEig.SincResults object containing five different components.
1. RESULTS : Matrix containing all the converged eigenvalue
First column contains the eigenvalue number
Second column contains the matrix size needed to achieve desired acurracy: tol
Third column contains the value of the conveged eigenvalue
Fourth column contains an approximation to the absolute error.
2. N : Vector containing the values of N used in the DESCM
3. hoptimal : Vector containing the optimal meshsizes used for the values of N
4. MatrixSizes : Vector(s) of Matrix sizes at every iteration. if Centro = true is use, two vectors are returned, matrix sizes for the even eigenvalues and matrix size for the odd eigenvalues.
5. All_Abs_Error_Approx : Matrix(ces) of approximations to the absolute errors at every iteration. if Centro = true is use, two Matrices are returned, one for the even eigenvalues and one for the odd eigenvalues.Suppose we are interested in computing the eigenvalues of the Laguerre equation:
(-D^2 + (35/4)/x^2 - 2 + x^2 /16 ) u(x) = λ u(x), with u(0) = u(∞) = 0we use the package function SincEigen to calculate the eigenvalues:
S = SincEigen( x -> (35/4)./x.^2 .- 2 .+ x.^2 /16 , ones , SemiInfinite1{Float64}() , [1.5,0.03125] , [1.0,2.0] , pi/4 )
S.RESULTSSuppose we are interested in computing the energy eigenvalues E of Schrödinger equation:
(-D^2 + V(x) ) ψ(x) = E ψ(x), with ψ(±∞) = 0for the quartic anharmonic oscillator potential:
V(x) = x.^2 + x.^4we use the package function SincEigen to calculate the eigenvalues:
S = SincEigen(V, ones, Infinite2{Float64}(), [0.125,0.125] , [2.0,2.0] , pi/4 )
S.RESULTS