d02kef

d02kef © Numerical Algorithms Group, 2002.

Purpose

D02KEF 2nd order Sturm-Liouville problem, regular/singular system, finite/infinite range, eigenvalue and eigenfunction, user-specified break-points

Synopsis

[elam,hmax,delam,match,maxit,ifail] = d02kef(xpoint,coeffn,bdyval,k,tol,...
elam,monit,report<,hmax,delam,match,maxfun,maxit,ifail>)

Description

 
                                     ~~~~~~~~
 D02KEF finds a specified eigenvalue (lambda) of a Sturm-
 Liouville system defined by a self-adjoint differential equation 
 of the second-order
 
                 (p(x)y')'+q(x;(lambda))y=0, a<x<b
 
 together with the appropriate boundary conditions at the two 
 (finite or infinite) end-points a and b. The functions p and q, 
 which are real-valued, must be defined by a subroutine COEFFN. 
 The boundary conditions must be defined by a subroutine 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 (lambda)=(lambda) of the eigenvalue has 
 been found, the routine integrates the differential equation once
 more with that value of (lambda), and with initial conditions 
 chosen so that the integral
 
                   b                
                   /    2    ddq    
                S= |y(x)  ----------(x;(lambda))dx
                   /      dd(lambda)
                   a                
 
 is approximately one. When q(x;(lambda)) is of the form 
 (lambda)w(x)+q(x), which is the most common case, S represents 
 the square of the norm of y induced by the inner product
 
                             b
                             /
                      <f,g>= |f(x)g(x)w(x)dx,
                             /
                             a
 
 with respect to which the eigenfunctions are mutually orthogonal.
 This normalisation 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 routine supplied by the 
 user is called at each integration mesh point x. Sufficient 
 information is returned to permit the user to compute y(x) and 
 y'(x) for printing or plotting. D02KEF passes across to REPORT, 
 not y and y', but the Pruefer variables (beta), (phi) and (rho) 
 on which the numerical method is based. Their relationship to y 
 and y' is given by the equations
 
 
          ______   ( (rho))   ( (phi)) 
 p(x)y'=\/(beta)exp( -----)cos( -----);
                   (   2  )   (   2  ) 
                                       
            1       ( (rho))   ( (phi))
      y= --------exp( -----)sin( -----).
           ______   (   2  )   (   2  )
         \/(beta)                      
 
 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 non-zero and of one sign throughout the 
       interval (a,b); and,
 
           ddq                                                 
 (b)    ---------- must be of one sign throughout (a,b) for all 
        dd(lambda)                                             
       relevant values of (lambda), and must not be identically 
       zero as x varies, for any (lambda).
 
 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.
 

Parameters

d02kef

Required Input Arguments:

xpoint (:)                            real
coeffn                                function (User-Supplied)
bdyval                                function (User-Supplied)
k                                     integer
tol                                   real
elam                                  real
monit                                 function (User-Supplied)
report                                function (User-Supplied)

Optional Input Arguments:                       <Default>

hmax (2,:)                            real     zeros(2,length(xpoint))
delam                                 real     0.25*max(1,abs(elam))
match                                 integer  0
maxfun                                integer  0
maxit                                 integer  0
ifail                                 integer  -1

Output Arguments:

elam                                  real
hmax (2,:)                            real
delam                                 real
match                                 integer
maxit                                 integer
ifail                                 integer