D02KEF : NAG Library, Mark 25

Note: before using this routine, please read the Users' Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.

D02KEF finds a specified eigenvalue of a regular or singular second-order Sturm–Liouville system on a finite or infinite interval, using a Pruefer transformation and a shooting method. It also reports values of the eigenfunction and its derivatives. Provision is made for discontinuities in the coefficient functions or their derivatives.

D02KEF has essentially the same purpose as D02KDF with minor modifications to enable values of the eigenfunction to be obtained after convergence to the eigenvalue has been achieved.

It first finds a specified eigenvalue

\tilde{λ}

of a Sturm–Liouville system defined by a self-adjoint differential equation of the second-order

{(p (x) y^{'})}^{'} + q (x; λ) y = 0, a < x < b,

together with appropriate boundary conditions at the two, finite or infinite, end points

a

and

b

. The functions

p

and

q

, which are real-valued, are defined by COEFFN. The boundary conditions must be defined by BDYVAL, and, in the case of a singularity at

a

or

b

, take the form of an asymptotic formula for the solution near the relevant end point.

When the final estimate

λ = \tilde{λ}

of the eigenvalue has been found, the routine integrates the differential equation once more with that value of

λ

, and with initial conditions chosen so that the integral

S = \int_{a}^{b} y {(x)}^{2} \frac{\partial q}{\partial λ} (x; λ) d x

is approximately one. When

q (x; λ)

is of the form

λ w (x) + q (x)

, which is the most common case,

S

represents the square of the norm of

y

induced by the inner product

〈f, g〉 = \int_{a}^{b} f (x) g (x) w (x) d x,

with respect to which the eigenfunctions are mutually orthogonal. This normalization of

y

is only approximate, but experience shows that

S

generally differs from unity by only one or two per cent.

During this final integration the REPORT is called at each integration mesh point

x

. Sufficient information is returned to permit you to compute

y (x)

and

y^{'} (x)

for printing or plotting. For reasons described in Section 9.2, D02KEF passes across to REPORT, not

y

and

y^{'}

, but the Pruefer variables

β

,

ϕ

and

ρ

on which the numerical method is based. Their relationship to

y

and

y^{'}

is given by the equations

p (x) y^{'} = \sqrt{β} \exp (\frac{ρ}{2}) \cos (\frac{ϕ}{2}), y = \frac{1}{\sqrt{β}} \exp (\frac{ρ}{2}) \sin (\frac{ϕ}{2}) .

A specimen REPORT is given in Section 10 below.

For the theoretical basis of the numerical method to be valid, the following conditions should hold on the coefficient functions:

(a)	$p (x)$ must be nonzero and must not change sign throughout the interval $(a, b)$ ; and,
(b)	$\frac{\partial q}{\partial λ}$ must not change sign throughout the interval $(a, b)$ for all relevant values of $λ$ , and must not be identically zero as $x$ varies, for any $λ$ .

Points of discontinuity in the functions

p

and

q

or their derivatives are allowed, and should be included as ‘break-points’ in the array XPOINT.

A good account of the theory of Sturm–Liouville systems, with some description of Pruefer transformations, is given in Chapter X of Birkhoff and Rota (1962). An introduction to the use of Pruefer transformations for the numerical solution of eigenvalue problems arising from physics and chemistry is given in Bailey (1966).

The scaled Pruefer method is described in a short note by Pryce and Hargrave (1977) and in some detail in the technical report by Pryce (1981).

Abramowitz M and Stegun I A (1972) Handbook of Mathematical Functions (3rd Edition) Dover Publications

Bailey P B (1966) Sturm–Liouville eigenvalues via a phase function SIAM J. Appl. Math. 14 242–249

Banks D O and Kurowski I (1968) Computation of eigenvalues of singular Sturm–Liouville Systems Math. Comput. 22 304–310

Birkhoff G and Rota G C (1962) Ordinary Differential Equations Ginn & Co., Boston and New York

Pryce J D (1981) Two codes for Sturm–Liouville problems Technical Report CS-81-01 Department of Computer Science, Bristol University

Pryce J D and Hargrave B A (1977) The scaled Prüfer method for one-parameter and multi-parameter eigenvalue problems in ODEs IMA Numerical Analysis Newsletter 1(3)

If on entry

IFAIL = 0

or

- 1

, explanatory error messages are output on the current error message unit (as defined by X04AAF).

Errors or warnings detected by the routine:

Note: error exits with

IFAIL = 2

,

3

,

6

,

7

,

8

or

11

are caused by being unable to set up or solve the differential equation at some iteration and will be immediately preceded by a call of MONIT giving diagnostic information. For other errors, diagnostic information is contained in

HMAX (2, j)

, for

j = 1, 2, \dots, m

, where appropriate.

See the discussion in Section 9.2.

Not applicable.

The time taken by D02KEF depends on the complexity of the coefficient functions, whether they or their derivatives are rapidly changing, the tolerance demanded, and how many iterations are needed to obtain convergence. The amount of work per iteration is roughly doubled when TOL is divided by

16

. To make the most economical use of the routine, one should try to obtain good initial values for ELAM and DELAM, and, where appropriate, good asymptotic formulae. Also the boundary matching points should not be set unnecessarily close to singular points. The extra time needed to compute the eigenfunction is principally the cost of one additional integration once the eigenvalue has been found.

A shooting method, for differential equation problems containing unknown parameters, relies on the construction of a ‘miss-distance function’, which for given trial values of the parameters measures how far the conditions of the problem are from being met. The problem is then reduced to one of finding the values of the parameters for which the miss-distance function is zero, that is to a root-finding process. Shooting methods differ mainly in how the miss-distance is defined.

D02KEF defines a miss-distance

f (λ)

based on the rotation about the origin of the point

P (x) = (p (x) y^{'} (x), y (x))

in the Phase Plane as the solution proceeds from

a

to

b

. The boundary conditions define the ray (i.e., two-sided line through the origin) on which

p (x)

should start, and the ray on which it should finish. The eigenvalue

k

defines the total number of half-turns it should make. Numerical solution is actually done by ‘shooting forward’ from

x = a

and ‘shooting backward’ from

x = b

to a matching point

x = c

. Then

f (λ)

is taken as the angle between the rays to the two resulting points

P_{a} (c)

and

P_{b} (c)

. A relative scaling of the

p y^{'}

and

y

axes, based on the behaviour of the coefficient functions

p

and

q

, is used to improve the numerical behaviour.

Figure 1

The resulting function

f (λ)

is monotonic over

- \infty < λ < \infty

, increasing if

\frac{\partial q}{\partial λ} > 0

and decreasing if

\frac{\partial q}{\partial λ} < 0

, with a unique zero at the desired eigenvalue

\tilde{λ}

. The routine measures

f (λ)

in units of a half-turn. This means that as

λ

increases,

f (λ)

varies by about

1

as each eigenvalue is passed. (This feature implies that the values of

f (λ)

at successive iterations – especially in the early stages of the iterative process – can be used with suitable extrapolation or interpolation to help the choice of initial estimates for eigenvalues near to the one currently being found.)

The routine actually computes a value for

f (λ)

with errors, arising from the local errors of the differential equation code and from the asymptotic formulae provided by you if singular points are involved. However, the error estimate output in DELAM is usually fairly realistic, in that the actual error

|\tilde{λ} - ELAM|

is within an order of magnitude of DELAM.

We pass the values of

β

,

ϕ

,

ρ

across through REPORT rather than converting them to values of

y

,

y^{'}

inside D02KEF, for the following reasons. First, there may be cases where auxiliary quantities can be more accurately computed from the Pruefer variables than from

y

and

y^{'}

. Second, in singular problems on an infinite interval

y

and

y^{'}

may underflow towards the end of the range, whereas the Pruefer variables remain well-behaved. Third, with high-order eigenvalues (and therefore highly oscillatory eigenfunctions) the eigenfunction may have a complete oscillation (or more than one oscillation) between two mesh points, so that values of

y

and

y^{'}

at mesh points give a very poor representation of the curve. The probable behaviour of the Pruefer variables in this case is that

β

and

ρ

vary slowly whilst

ϕ

increases quickly: for all three Pruefer variables linear interpolation between the values at adjacent mesh points is probably sufficiently accurate to yield acceptable intermediate values of

β

,

ϕ

,

ρ

(and hence of

y, y^{'}

) for graphical purposes.

Similar considerations apply to the exponentially decaying ‘tails’ of the eigenfunctions that often occur in singular problems. Here

ϕ

has approximately constant value whilst

ρ

increases rapidly in the direction of integration, though the step length is generally fairly small over such a range.

If the solution is output through REPORT at

x

values which are too widely spaced, the step length can be controlled by choosing HMAX suitably, or, preferably, by reducing TOL. Both these choices will lead to more accurate eigenvalues and eigenfunctions but at some computational cost.

This point is always one of the values

x_{i}

in array XPOINT. It may be specified using the parameter MATCH. The default value is chosen to be the value of that

x_{i}

,

2 \leq i \leq m - 1

, that lies closest to the middle of the interval

[x_{2}, x_{m - 1}]

. If there is a tie, the rightmost candidate is chosen. In particular if there are no break-points, then

c = x_{m - 1}

(

= x_{3}

); that is, the shooting is from left to right in this case. A break-point may be inserted purely to move

c

to an interior point of the interval, even though the form of the equations does not require it. This often speeds up convergence especially with singular problems.

Note that the shooting method used by the code integrates first from the left-hand end

x_{l}

, then from the right-hand end

x_{r}

, to meet at the matching point

c

in the middle. This will of course be reflected in printed or graphical output. The diagram shows a possible sequence of nine mesh points

τ_{1}

through

τ_{9}

in the order in which they appear, assuming there are just two sub-intervals (so

m = 5

).

Figure 2

Since the shooting method usually fails to match up the two ‘legs’ of the curve exactly, there is bound to be a jump in

y

, or in

p (x) y^{'}

or both, at the matching point

c

. The code in fact ‘shares’ the discrepancy out so that both

y

and

p (x) y^{'}

have a jump. A large jump does not imply an inaccurate eigenvalue, but implies either

(a)	a badly chosen matching point: if $q (x; λ)$ has a ‘humped’ shape, $c$ should be chosen near the maximum value of $q$ , especially if $q$ is negative at the ends of the interval;
(b)	an inherently ill-conditioned problem, typically one where another eigenvalue is pathologically close to the one being sought. In this case it is extremely difficult to obtain an accurate eigenfunction.

In Section 10, we find the 11th eigenvalue and corresponding eigenfunction of the equation

y^{''} + (λ - x - 2 / x^{2}) y = 0 on 0 < x < \infty,

the boundary conditions being that

y

should remain bounded as

x

tends to

0

and

x

tends to

\infty

. The coding of this problem is discussed in detail in Section 9.5.

The choice of matching point

c

is open. If we choose

c = 30.0

as in D02KDF example program we find that the exponentially increasing component of the solution dominates and we get extremely inaccurate values for the eigenfunction (though the eigenvalue is determined accurately). The values of the eigenfunction calculated with

c = 30.0

are given schematically in Figure 3.

Figure 3

If we choose

c

as the maximum of the hump in

q (x; λ)

(see item (a) above) we instead obtain the accurate results given in Figure 4

Figure 4

Coding COEFFN is straightforward except when break-points are needed. The examples below show:

(a)	a simple case,
(b)	a case in which discontinuities in the coefficient functions or their derivatives necessitate break-points, and
(c)	a case where break-points together with the HMAX parameter are an efficient way to deal with a coefficient function that is well-behaved except over one short interval.

(Some of these cases are among the examples in Section 10.)

Example A

The modified Bessel equation

x {(x y^{'})}^{'} + (λ x^{2} - ν^{2}) y = 0 .

Assuming the interval of solution does not contain the origin and dividing through by

x

, we have

p (x) = x

and

q (x; λ) = λ x - ν^{2} / x

. The code could be

SUBROUTINE COEFFN (P,Q,DQDL,X,ELAM,JINT)
...
P = X
Q = ELAM*X - NU*NU/X
DQDL = X
RETURN
END

where NU (standing for

ν

) is a real variable that might be defined in a DATA statement, or might be in user-declared COMMON so that its value could be set in the main program.

Example B

The Schroedinger equation

y^{''} + (λ + q (x)) y = 0,

where

q (x) = \{\begin{cases} x^{2} - 10 & (|x| \leq 4), \\ \frac{6}{|x|} & (|x| > 4), \end{cases}

over some interval ‘approximating to

(- \infty, \infty)

’, say

[- 20, 20]

. Here we need break-points at

\pm 4

, forming three sub-intervals

i_{1} = [- 20, - 4]

,

i_{2} = [- 4, 4]

,

i_{3} = [4, 20]

. The code could be

SUBROUTINE COEFFN(P,Q,DQDL,X,ELAM,JINT)
...
IF (JINT.EQ.2) THEN
  Q = ELAM + X*X - 10.0E0
ELSE
  Q = ELAM + 6.0E0/ABS(X)
ENDIF
P = 1.0E0
DQDL = 1.0E0
RETURN
END

The array XPOINT would contain the values

x_{1}

,

- 20.0

,

- 4.0

,

+ 4.0

,

+ 20.0

,

x_{6}

and

m

would be

6

. The choice of appropriate values for

x_{1}

and

x_{6}

depends on the form of the asymptotic formula computed by BDYVAL and the technique is discussed in Section 9.5.

Example C

y^{''} + λ (1 - 2 e^{- 100 x^{2}}) y = 0, - 10 \leq x \leq 10 .

Here

q (x; λ)

is nearly constant over the range except for a sharp inverted spike over approximately

- 0.1 \leq x \leq 0.1

. There is a danger that the routine will build up to a large step size and ‘step over’ the spike without noticing it. By using break-points – say at

\pm 0.5

– one can restrict the step size near the spike without impairing the efficiency elsewhere.

The code for COEFFN could be

SUBROUTINE COEFFN(P,Q,DQDL,X,ELAM,JINT)
...
P = 1.0E0
DQDL = 1.0E0 - 2.0E0*EXP(-100.0E0*X*X)
Q = ELAM*DQDL
RETURN
END

XPOINT might contain

- 10.0

,

- 10.0

,

- 0.5

,

0.5

,

10.0

,

10.0

(assuming

\pm 10

are regular points) and

m

would be

6

.

HMAX (1, j)

, for

j = 1, 2, 3

, might contain

0.0

,

0.1

and

0.0

.

Quoting from page 243 of Bailey (1966): ‘Usually ... the differential equation has two essentially different types of solution near a singular point, and the boundary condition there merely serves to distinguish one kind from the other. This is the case in all the standard examples of mathematical physics.’

In most cases the behaviour of the ratio

p (x) y^{'} / y

near the point is quite different for the two types of solution. Essentially what you provide through the BDYVAL is an approximation to this ratio, valid as

x

tends to the singular point (SP).

You must decide (a) how accurate to make this approximation or asymptotic formula, for example how many terms of a series to use, and (b) where to place the boundary matching point (BMP) at which the numerical solution of the differential equation takes over from the asymptotic formula. Taking the BMP closer to the SP will generally improve the accuracy of the asymptotic formula, but will make the computation more expensive as the Pruefer differential equations generally become progressively more ill-behaved as the SP is approached. You are strongly recommended to experiment with placing the BMPs. In many singular problems quite crude asymptotic formulae will do. To help you avoid needlessly accurate formulae, D02KEF outputs two ‘sensitivity coefficients’

σ_{l}, σ_{r}

which estimate how much the errors at the BMPs affect the computed eigenvalue. They are described in detail in Section 9.6.

Example of coding BDYVAL:

The example below illustrates typical situations:

y^{''} + (λ - x - \frac{2}{x^{2}}) y = 0, 0 < x < \infty

the boundary conditions being that

y

should remain bounded as

x

tends to

0

and

x

tends to

\infty

.

At the end

x = 0

there is one solution that behaves like

x^{2}

and another that behaves like

x^{- 1}

. For the first of these solutions

p (x) y^{'} / y

is asymptotically

2 / x

while for the second it is asymptotically

- 1 / x

. Thus the desired ratio is specified by setting

YL (1) = x and YL (2) = 2.0 .

At the end

x = \infty

the equation behaves like Airy's equation shifted through

λ

, i.e., like

y^{''} - t y = 0

where

t = x - λ

, so again there are two types of solution. The solution we require behaves as

\exp (- \frac{2}{3} t^{\frac{3}{2}}) / \sqrt[4]{t}

and the other as

\exp (+ \frac{2}{3} t^{\frac{3}{2}}) / \sqrt[4]{t} .

Hence, the desired solution has

p (x) y^{'} / y \sim - \sqrt{t}

so that we could set

YR (1) = 1.0

and

YR (2) = - \sqrt{x - λ}

. The complete subroutine might thus be

SUBROUTINE BDYVAL (XL,XR,ELAM,YL,YR)
real XL, XR, ELAM, YL(3), YR(3)
YL(1) = XL
YL(2) = 2.0E0
YR(1) = 1.0E0
YR(2) = -SQRT(XR-ELAM)
RETURN
END

Clearly for this problem it is essential that any value given by D02KEF to XR is well to the right of the value of ELAM, so that you must vary the right-hand BMP with the eigenvalue index

k

. One would expect

λ_{k}

to be near the

k

th zero of the Airy function

Ai (x)

, so there is no problem estimating ELAM.

More accurate asymptotic formulae are easily found: near

x = 0

by the standard Frobenius method, and near

x = \infty

by using standard asymptotics for

Ai (x)

,

{Ai}^{'} (x)

, (see page 448 of Abramowitz and Stegun (1972)).

For example, by the Frobenius method the solution near

x = 0

has the expansion

y = x^{2} (c_{0} + c_{1} x + c_{2} x^{2} + \dots)

with

c_{0} = 1, c_{1} = 0, c_{2} = \frac{- λ}{10}, c_{3} = \frac{1}{18}, \dots, c_{n} = \frac{c_{n - 3} - λ c_{n - 2}}{n (n + 3)} .

This yields

\frac{p (x) y^{'}}{y} = \frac{2 - \frac{2}{5} λ x^{2} + \dots}{x (1 - \frac{λ}{10} x^{2} + \dots)} .

The sensitivity parameters

σ_{l}

,

σ_{r}

(held in

HMAX (1, m - 1)

and

HMAX (1, m)

on output) estimate the effect of errors in the boundary conditions. For sufficiently small errors

Δ y

,

Δ p y^{'}

in

y

and

p y^{'}

respectively, the relations

\begin{array}{l} Δ λ ≃ {(y . Δ p y^{'} - p y^{'} . Δ y)}_{l} σ_{l} \\ Δ λ ≃ {(y . Δ p y^{'} - p y^{'} . Δ y)}_{r} σ_{r} \end{array}

are satisfied, where the subscripts

l

,

r

denote errors committed at the left- and right-hand BMPs respectively, and

Δ λ

denotes the consequent error in the computed eigenvalue.

This is a pitfall to beware of at a singular point. If the BMP is chosen so far from the SP that a zero of the desired eigenfunction lies in between them, then the routine will fail to ‘notice’ this zero. Since the index of

k

of an eigenvalue is the number of zeros of its eigenfunction, the result will be that

(a)	the wrong eigenvalue will be computed for the given index $k$ – in fact some $λ_{k + k^{'}}$ will be found where $k^{'} \geq 1$ ;
(b)	the same index $k$ can cause convergence to any of several eigenvalues depending on the initial values of ELAM and DELAM.

It is up to you to take suitable precautions – for instance by varying the position of the BMPs in the light of knowledge of the asymptotic behaviour of the eigenfunction at different eigenvalues.

This example finds the 11th eigenvalue and eigenfunction of the example of Section 9.5, using the simple asymptotic formulae for the boundary conditions.

Comparison of the results from this example program with the corresponding results from D02KDF example program shows that similar output is produced from MONIT, followed by the eigenfunction values from REPORT, and then a further line of information from MONIT (corresponding to the integration to find the eigenfunction). Final information is printed within the example program exactly as with D02KDF.

Note the discrepancy at the matching point

c

(

= \sqrt[3]{4}

, the maximum of

q (x; λ)

, in this case) between the solutions obtained by integrations from left- and right-hand end points.

$HMAX (2, 1) = 1$ :	$M < 4$ ;
$HMAX (2, 1) = 2$ :	$K < 0$ ;
$HMAX (2, 1) = 3$ :	$TOL \leq 0.0$ ;
$HMAX (2, 1) = 4$ :	$XPOINT (1)$ to $XPOINT (m)$ are not in ascending order. $HMAX (2, 2)$ gives the position $i$ in XPOINT where this was detected.

NAG Library Routine Document

D02KEF

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

9.1 Timing

9.2 General Description of the Algorithm

9.3 The Position of the Shooting Matching Point $c$

9.4 Examples of Coding the COEFFN

9.5 Examples of Boundary Conditions at Singular Points

9.6 The Sensitivity Parameters $σ_{l}$ and $σ_{r}$

9.7 ‘Missed Zeros’

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG Library Routine DocumentD02KEF

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

9.1 Timing

9.2 General Description of the Algorithm

9.3 The Position of the Shooting Matching Point c

9.4 Examples of Coding the COEFFN

9.5 Examples of Boundary Conditions at Singular Points

9.6 The Sensitivity Parameters σl and σr

9.7 ‘Missed Zeros’

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG Library Routine Document

D02KEF

9.3 The Position of the Shooting Matching Point $c$

9.6 The Sensitivity Parameters $σ_{l}$ and $σ_{r}$