SLEIG

Determines eigenvalues, eigenfunctions and/or spectral density functions for Sturm-Liouville problems in the form

 

with boundary conditions (at regular points)

 

Required Arguments

CONS — Array of size eight containing

a1, a1′, a2, a2′, b1, b2, a, and b

in locations CONS(1) through CONS(8), respectively. (Input)

COEFFN — User-supplied subroutine to evaluate the coefficient functions. The usage is CALL COEFFN (X, PX, QX, RX)

X — Independent variable. (Input)

PX — The value of p(x) at X. (Output)

QX — The value of q(x) at X. (Output)

RX — The value of r(x) at X. (Output)

COEFFN must be declared EXTERNAL in the calling program.

ENDFIN — Logical array of size two. ENDFIN(1) = .true. if the endpoint a is finite, ENDFIN(2) = .true. if endpoint b is finite. (Input)

INDEX — Vector of size NUMEIG containing the indices of the desired eigenvalues. (Input)

EVAL — Array of length NUMEIG containing the computed approximations to the eigenvalues whose indices are specified in INDEX. (Output)

Optional Arguments

NUMEIG — The number of eigenvalues desired. (Input)
Default: NUMEIG = size (INDEX,1).

TEVLAB — Absolute error tolerance for eigenvalues. (Input)
Default: TEVLAB = 10 * machine precision.

TEVLRL — Relative error tolerance for eigenvalues. (Input)
Default: TEVLRL = SQRT(machine precision).

FORTRAN 90 Interface

Generic: CALL SLEIG (CONS, COEFFN, ENDFIN, INDEX, EVAL [, …])

Specific: The specific interface names are S_SLEIG and D_SLEIG.

FORTRAN 77 Interface

Single: CALL SLEIG (CONS, COEFFN, ENDFIN, NUMEIG, INDEX, TEVLAB, TEVLRL, EVAL)

Double: The double precision name is DSLEIG.

Description

This subroutine is designed for the calculation of eigenvalues, eigenfunctions and/or spectral density functions for Sturm-Liouville problems in the form

 

with boundary conditions (at regular points)

 

We assume that

 

when aʹ1 ≠ 0 and aʹ2 ≠ 0. The problem is considered regular if and only if

• a and b are finite,

• p(x) and r(x) are positive in (a, b),

• 1/p(x), q(x) and r(x) are locally integrable near the endpoints.

Otherwise the problem is called singular. The theory assumes that p, pʹ, q, and r are at least continuous on (a, b), though a finite number of jump discontinuities can be handled by suitably defining an input mesh.

For regular problems, there are an infinite number of eigenvalues

λ0 < λ1 < … < λk, k→∞

Each eigenvalue has an associated eigenfunction which is unique up to a constant. For singular problems, there is a wide range in the behavior of the eigenvalues.

As presented in Pruess and Fulton (1993) the approach is to replace (1) by a new problem

 

with analogous boundary conditions

 

where

 

are step function approximations to p, q, and r, respectively. Given the mesh a = x1 < x2 < … < xN+1 = b, the usual choice for the step functions uses midpoint interpolation, i. e.,

 

for x in (xn, xn+1) and similarly for the other coefficient functions. This choice works well for regular problems. Some singular problems require a more sophisticated technique to capture the asymptotic behavior. For the midpoint interpolants, the differential equation (2) has the known closed form solution in (xn, xn+1)

 

with

 

where

 

and

 

Starting with,

 

consistent with the boundary condition,

 

an algorithm is to compute for n = 1, 2, …, N,

 

which is a shooting method. For a fixed mesh we can iterate on the approximate eigenvalue until the boundary condition at b is satisfied. This will yield an O(h2) approximation

 

to some λk.

The problem (2) has a step spectral function given by

 

where the sum is taken over k such that

 

and

 

Comments

1. Workspace may be explicitly provided, if desired, by use of S2EIG/DS2EIG. The reference is:

CALL S2EIG (CONS, COEFFN, ENDFIN, NUMEIG, INDEX, TEVLAB, TEVLRL, EVAL, JOB, IPRINT, TOLS, NUMX, XEF, NRHO, T, TYPE, EF, PDEF, RHO, IFLAG, WORK, IWORK)

The additional arguments are as follows:

JOB — Logical array of length five. (Input)

JOB(1) = .true. if a set of eigenvalues are to be computed but not their eigenfunctions.

JOB(2) = .true. if a set of eigenvalue and eigenfunction pairs are to be computed.

JOB(3) = .true. if the spectral function is to be computed
over some subinterval of the essential spectrum.

JOB(4) = .true. if the normal automatic classification is overridden. If JOB(4) = .true. then TYPE(*,*) must be entered correctly. Most users will not want to override the classification process, but it might be appropriate for users experimenting with problems for which the coefficient functions do not have power-like behavior near the singular endpoints. The classification is considered sufficiently important for spectral density function calculations that JOB(4) is ignored with JOB(3) = .true..

JOB(5) = .true. if mesh distribution is chosen by SLEIG. If JOB(5) = .true. and NUMX is zero, the number of mesh points are also chosen by SLEIG. If NUMX > 0 then NUMX mesh points will be used. If JOB(5) = .false., the number NUMX and distribution XEF(*) must be input by the user.

IPRINT — Control levels of internal printing. (Input)
No printing is performed if IPRINT = 0. If either JOB(1) or JOB(2) is true:

 

IPRINT	Printed Output
1	Initial mesh (the first 51 or fewer points), eigenvalue estimate at each level.
4	The above and at each level matching point for eigenfunction shooting, X(*), EF(*) and PDEF(*) values.
5	The above and at each level the brackets for the eigenvalue search, intermediate shooting information for the eigenfunction and eigenfunction norm.

If JOB(3) = .true.

 

IPRINT	Printed Output
1	The actual (a, b) used at each iteration and the total number of eigenvalues computed.
2	The above and switchover points to the asymptotic formulas, and some intermediate (t) approximations.
3	The above and initial meshes for each iteration, the index of the largest eigenvalue which may be computed, and various eigenvalue and RN values.
4	The above and values at each level.
5	The above and RN add eigenvalues below the switchover point

If JOB(4) = .false.

 

IPRINT	Printed Output
2	Output a description of the spectrum.
3	The above and the constants for the Friedrichs' boundary condition(s).
5	The above and and intermediate details of the classification calculation.

TOLS — Array of length 4 containing tolerances. (Input)

TOLS(1) — absolute error tolerance for eigenfunctions

TOLS(2) — relative error tolerance for eigenfunctions

TOLS(3) — absolute error tolerance for eigenfunction derivatives

TOLS(4) — relative error tolerance for eigenfunction derivatives

The absolute tolerances must be positive.The relative tolerances must be at least 100 * amach(4)

NUMX — Integer whose value is the number of output points where each eigenfunction is to be evaluated (the number of entries in XEF(*)) when JOB(2) = .true.. If JOB(5)= .false. and NUMX is greater than zero, then NUMX is the number of points in the initial mesh used. If JOB(5) = .false., the points in XEF(*) should be chosen with a reasonable distribution. Since the endpoints a and b must be part of any mesh, NUMX cannot be one in this case. If JOB(5) = .false. and JOB(3) = .true., then NUMX must be positive. On output, NUMX is set to the number of points for eigenfunctions when input NUMX = 0, and JOB(2) or JOB(5) = .true.. (Input/Output)

XEF — Array of points on input where eigenfunction estimates are desired, if JOB(2) = .true.. Otherwise, if JOB(5) = .false. and NUMX is greater than zero, the user’s initial mesh is entered. The entries must be ordered so that a = XEF(1) < XEF(2) < < XEF(NUMX) = b. If either endpoint is infinite, the corresponding XEF(1) or XEF(NUMX) is ignored. However, it is required that XEF(2) be negative when ENDFIN(1) = .false., and that XEF(NUMX‑1) be positive when ENDFIN(2) = .false.. On output, XEF(*) is changed only if JOB(2) and JOB(5) are true. If JOB(2) = .false., this vector is not referenced. If JOB(2) = .true. and NUMX is greater than zero on input, XEF(*) should be dimensioned at least NUMX + 16. If JOB(2) is true and NUMX is zero on input, XEF(*) should be dimensioned at least 31.

NRHO — The number of output values desired for the array RHO(*). NRHO is not used if JOB(3) = .false.. (Input)

T — Real vector of size NRHO containing values where the spectral function RHO(*) is desired. The entries must be sorted in increasing order. The existence and location of a continuous spectrum can be determined by calling SLEIG with the first four entries of JOB set to false and IPRINT set to 1. T(*) is not used if JOB(3) = .false.. (Input)

TYPE — 4 by 2 logical matrix. Column 1 contains information about endpoint a and column 2 refers to endpoint b.
TYPE(1,*) = .true. if and only if the endpoint is regular
TYPE(2,*) = .true. if and only if the endpoint is limit circle
TYPE(3,*) = .true. if and only if the endpoint is nonoscillatory for all eigenvalues
TYPE(4,*) = .true. if and only if the endpoint is oscillatory for all eigenvalues
Note: all of these values must be correctly input if JOB(4) = .true..
Otherwise, TYPE(*,*) is output. (Input/Output)

EF — Array of eigenfunction values. EF((k‑1)*NUMX + i) is the estimate of u(XEF(i)) corresponding to the eigenvalue in EV(k). If JOB(2) = .false. then this vector is not referenced. If JOB(2) = .true. and NUMX is greater than zero on entry, then EF(*) should be dimensioned at least NUMX * NUMEIG. If JOB(2) = .true. and NUMX is zero on input, then EF(*) should be dimensioned 31 * NUMEIG. (Output)

PDEF — Array of eigenfunction derivative values. PDEF((k‑1)*NUMX + i) is the estimate of (puʹ) (XEF(i)) corresponding to the eigenvalue in EV(k). If JOB(2) = .false. this vector is not referenced. If JOB(2) = .true., it must be dimensioned the same as EF(*). (Output)

RHO — Array of size NRHO containing values for the spectral density function ρ(t), RHO(I) = ρ(T(I)). This vector is not referenced if JOB(3) is false. (Output)

IFLAG — Array of size max(1, numeig) containing information about the output. IFLAG(K) refers to the K-th eigenvalue, when JOB(1) or JOB(2) = .true.. Otherwise, only IFLAG(1) is used. Negative values are associated with fatal errors, and the calculations are ceased. Positive values indicate a warning. (Output)
These are the possible values for IFLAG(K)

IFLAG(K)	Description
-1	Too many levels needed for the eigenvalue calculation; problem seems too difficult at this tolerance. Are the coefficient functions nonsmooth?
-2	Too many levels needed for the eigenfunction calculation; problem seems too difficult at this tolerance. Are the eigenfunctions ill-conditioned?
-3	Too many levels needed for the spectral density calculation; problem seems too difficult at this tolerance.
-4	The user has requested the spectral density function for a problem which has no continuous spectrum.
-5	The user has requested the spectral density function for a problem with both endpoints generating essential spectrum, i.e. both endpoints either OSC or O-NO.
-6	The user has requested the spectral density function for a problem in spectral category 2 for which a proper normalization of the solution at the NONOSC endpoint is not known; for example, problems with an irregular singular point or infinite endpoint at one end and continuous spectrum generated at the other.
-7	Problems were encountered in obtaining a bracket.
-8	Too small a step was used in the integration. The TOLS(*) values may be too small for this problem.
-9	Too small a step was used in the spectral density function calculation for which the continuous spectrum is generated by a finite endpoint.
-10	An argument to the circular trig functions is too large. Try running the problem again with a finer initial mesh or, for singular problems, use interval truncation.
-15	p(x) and r(x) are not positive in the interval (a, b).
-20	Eigenvalues and/or eigenfunctions were requested for a problem with an OSC singular endpoint. Interval truncation must be used on such problems.
1	Failure in the bracketing procedure probably due to a cluster of eigenvalues which the code cannot separate. Calculations have continued but any eigenfunction results are suspect. Try running the problem again with tighter input tolerances to separate the cluster.
2	There is uncertainty in the classification for this problem. Because of the limitations of floating point arithmetic, and the nature of the finite sampling, the routine cannot be certain about the classification information at the requested tolerance.
3	There may be some eigenvalues embedded in the essential spectrum. Use of IPRINT greater than zero will provide additional output giving the location of the approximating eigenvalues for the step function problem. These could be extrapolated to estimate the actual eigenvalue embedded in the essential spectrum.
4	A change of variables was made to avoid potentially slow convergence. However, the global error estimates may not be as reliable. Some experimentation using different tolerances is recommended.
6	There were problems with eigenfunction convergence in a spectral density calculation. The output ρ(t) may not be accurate.

WORK — Array of size MAX(1000, NUMEIG + 22) used for workspace.

IWORK — Integer array of size NUMEIG + 3 used for workspace.

Examples

Example 1

This example computes the first ten eigenvalues of the problem from Titchmarsh (1962) given by

p(x) = r(x) = 1

q(x) = x

[a, b] = [0, ∞]

u(a) = u(b) = 0

The eigenvalues are known to be the zeros of

 

For each eigenvalue λk, the program prints k, λk and f(λk).

 

USE SLEIG_INT

USE CBJS_INT

 

IMPLICIT NONE

! SPECIFICATIONS FOR LOCAL VARIABLES

INTEGER I, INDEX(10), NUMEIG, NOUT

REAL CONS(8), EVAL(10), LAMBDA, TEVLAB,&

TEVLRL, XNU

COMPLEX CBS1(1), CBS2(1), Z

LOGICAL ENDFIN(2)

! SPECIFICATIONS FOR INTRINSICS

INTRINSIC CMPLX, SQRT

REAL SQRT

COMPLEX CMPLX

! SPECIFICATIONS FOR SUBROUTINES

! SPECIFICATIONS FOR FUNCTIONS

EXTERNAL COEFF

!

CALL UMACH (2, NOUT)

! Define boundary conditions

CONS(1) = 1.0

CONS(2) = 0.0

CONS(3) = 0.0

CONS(4) = 0.0

CONS(5) = 1.0

CONS(6) = 0.0

CONS(7) = 0.0

CONS(8) = 0.0

!

ENDFIN(1) = .TRUE.

ENDFIN(2) = .FALSE.

! Compute the first 10 eigenvalues

NUMEIG = 10

DO 10 I=1, NUMEIG

INDEX(I) = I - 1

10 CONTINUE

! Set absolute and relative tolerance

!

CALL SLEIG (CONS, COEFF, ENDFIN, INDEX, EVAL)

!

XNU = -1.0/3.0

WRITE(NOUT,99998)

DO 20 I=1, NUMEIG

LAMBDA = EVAL(I)

Z = CMPLX(2.0/3.0*LAMBDA*SQRT(LAMBDA),0.0)

CALL CBJS (XNU, Z, 1, CBS1)

CALL CBJS (-XNU, Z, 1, CBS2)

WRITE (NOUT,99999) I-1, LAMBDA, REAL(CBS1(1) + CBS2(1))

20 CONTINUE

!

99998 FORMAT(/, 2X, 'index', 5X, 'lambda', 5X, 'f(lambda)',/)

99999 FORMAT(I5, F13.4, E15.4)

END

!

SUBROUTINE COEFF (X, PX, QX, RX)

! SPECIFICATIONS FOR ARGUMENTS

REAL X, PX, QX, RX

!

PX = 1.0

QX = X

RX = 1.0

RETURN

END

Output

 

index     lambda     f(lambda)

 

    0       2.3381    -0.8285E-05

    1       4.0879    -0.1651E-04

    2       5.5205     0.6843E-04

    3       6.7867    -0.4523E-05

    4       7.9440     0.8952E-04

    5       9.0227     0.1123E-04

    6      10.0401     0.1031E-03

    7      11.0084    -0.7913E-04

    8      11.9361    -0.5095E-04

    9      12.8293     0.4645E-03

Example 2

In this problem from Scott, Shampine and Wing (1969),

p(x) = r(x) = 1

q(x) = x2 + x4

[a, b] = [-∞,∞]

u(a) = u(b) = 0

the first eigenvalue and associated eigenfunction, evaluated at selected points, are computed. As a rough check of the correctness of the results, the magnitude of the residual

 

is printed. We compute a spline interpolant to uʹ and use the function CSDER to estimate the quantity ‑(p(x)uʹ)ʹ.

 

USE S2EIG_INT

USE CSDER_INT

USE UMACH_INT

USE CSAKM_INT

 

IMPLICIT NONE

! SPECIFICATIONS FOR LOCAL VARIABLES

 

INTEGER I, IFLAG(1), INDEX(1), IWORK(100), NINTV, NOUT, NRHO, &

NUMEIG, NUMX

REAL BRKUP(61), CONS(8), CSCFUP(4,61), EF(61), EVAL(1), &

LAMBDA, PDEF(61), PX, QX, RESIDUAL, RHO(1), RX, T(1), &

TEVLAB, TEVLRL, TOLS(4), WORK(3000), X, XEF(61)

LOGICAL ENDFIN(2), JOB(5), TYPE(4,2)

! SPECIFICATIONS FOR INTRINSICS

INTRINSIC ABS, REAL

REAL ABS, REAL

! SPECIFICATIONS FOR SUBROUTINES

EXTERNAL COEFF

! Define boundary conditions

CONS(1) = 1.0

CONS(2) = 0.0

CONS(3) = 0.0

CONS(4) = 0.0

CONS(5) = 1.0

CONS(6) = 0.0

CONS(7) = 0.0

CONS(8) = 0.0

! Compute eigenvalue and eigenfunctions

JOB(1) = .FALSE.

JOB(2) = .TRUE.

JOB(3) = .FALSE.

JOB(4) = .FALSE.

JOB(5) = .FALSE.

!

ENDFIN(1) = .FALSE.

ENDFIN(2) = .FALSE.

! Compute eigenvalue with index 0

NUMEIG = 1

INDEX(1) = 0

!

TEVLAB = 1.0E-3

TEVLRL = 1.0E-3

TOLS(1) = TEVLAB

TOLS(2) = TEVLRL

TOLS(3) = TEVLAB

TOLS(4) = TEVLRL

NRHO = 0

! Set up mesh, points at which u and

! u' will be computed

NUMX = 61

DO 10 I=1, NUMX

XEF(I) = 0.05*REAL(I-31)

10 CONTINUE

!

CALL S2EIG (CONS, COEFF, ENDFIN, NUMEIG, INDEX, TEVLAB, TEVLRL, &

EVAL, JOB, 0, TOLS, NUMX, XEF, NRHO, T, TYPE, EF, &

PDEF, RHO, IFLAG, WORK, IWORK)

!

LAMBDA = EVAL(1)

20 CONTINUE

! Compute spline interpolant to u'

!

CALL CSAKM (XEF, PDEF, BRKUP, CSCFUP)

NINTV = NUMX - 1

!

CALL UMACH (2, NOUT)

WRITE (NOUT,99997) ' lambda = ', LAMBDA

WRITE (NOUT,99999)

! At a subset of points from the

! input mesh, compute residual =

! abs( -(u')' + q(x)u - lambda*u ).

! We know p(x) = 1 and r(x) = 1.

DO 30 I=1, 41, 2

X = XEF(I+10)

CALL COEFF (X, PX, QX, RX)

!

! Use the spline fit to u' to

! estimate u'' with CSDER

!

RESIDUAL = ABS(-CSDER(1,X,BRKUP,CSCFUP)+QX*EF(I+10)- &

LAMBDA*EF(I+10))

WRITE (NOUT,99998) X, EF(I+10), PDEF(I+10), RESIDUAL

30 CONTINUE

!

99997 FORMAT (/, A14, F10.5, /)

99998 FORMAT (5X, F4.1, 3F15.5)

99999 FORMAT (7X, 'x', 11X, 'u(x)', 10X, 'u''(x)', 9X, 'residual', /)

END

!

SUBROUTINE COEFF (X, PX, QX, RX)

! SPECIFICATIONS FOR ARGUMENTS

REAL X, PX, QX, RX

!

PX = 1.0

QX = X*X + X*X*X*X

RX = 1.0

RETURN

END

Output

 

lambda = 1.39247

x u(x) u'(x) residual

-1.0 0.38632 0.65019 0.00189

-0.9 0.45218 0.66372 0.00081

-0.8 0.51837 0.65653 0.00023

-0.7 0.58278 0.62827 0.00113

-0.6 0.64334 0.57977 0.00183

-0.5 0.69812 0.51283 0.00230

-0.4 0.74537 0.42990 0.00273

-0.3 0.78366 0.33393 0.00265

-0.2 0.81183 0.22811 0.00273

-0.1 0.82906 0.11570 0.00278

0.0 0.83473 0.00000 0.00136

0.1 0.82893 -0.11568 0.00273

0.2 0.81170 -0.22807 0.00273

0.3 0.78353 -0.33388 0.00267

0.4 0.74525 -0.42983 0.00265

0.5 0.69800 -0.51274 0.00230

0.6 0.64324 -0.57967 0.00182

0.7 0.58269 -0.62816 0.00113

0.8 0.51828 -0.65641 0.00023

0.9 0.45211 -0.66361 0.00081

1.0 0.38626 -0.65008 0.00189