код для вставки на сайт или в блог

ссылки на документ

Journal of Symbolic Computation 35 (2003) 465?485
www.elsevier.com/locate/jsc
A combined symbolic and numerical algorithm for
the computation of zeros of orthogonal
polynomials and special functions
Amparo Gila,c,?, Javier Segurab,c,1
a Departamento de Matema?ticas, Facultad Ciencias, Universidad Auto?noma de Madrid, 28049-Madrid, Spain
b Departamento de Matema?ticas, Universidad Carlos III de Madrid, 28911-Legane?s, Madrid, Spain
c Universita?t GhK, FB 17 Mathematik-Informatik, 34132-Kassel, Germany
Received 16 November 2001; accepted 16 April 2002
Abstract
A Maple algorithm for the computation of the zeros of orthogonal polynomials (OPs) and
special functions (SFs) in a given interval [x1 , x2 ] is presented. The program combines symbolic
and numerical calculations and it is based on fixed point iterations. The program uses as inputs
the analytic expressions for the coefficients of the three-term recurrence relation and a differencedifferential relation satisfied by the set of OPs or SFs. The performance of the method is illustrated
with several examples: Hermite, Chebyshev, Legendre, Jacobi and Gegenbauer polynomials, Bessel,
Coulomb and Conical functions. Е 2003 Elsevier Science Ltd. All rights reserved.
Keywords: Zeros; Orthogonal polynomials; Special functions; Fixed point iterations
1. Introduction
The accurate and efficient computation of the zeros of orthogonal polynomials (OPs) is
a relevant numerical issue; as is well known, the abscissas of Gaussian quadrature formulas
are the roots of specific OPs. Also, the zeros of special functions (SFs) are important in a
vast number of applications, in particular the zeros of Bessel functions.
The most general methods for the computation of the zeros of OPs are matrix methods,
based on the Golub?Welsch (Golub and Welsch, 1969) algorithm which uses a result
of Wilf (1962). Matrix methods are also available (Ikebe, 1975; Ikebe et al., 1991) for
SFs which are minimal with respect to a 3-term recurrence relation. Methods based on
? Corresponding author.
E-mail addresses: amparo.gil@uam.es (A. Gil), jsegura@math.uc3m.es (J. Segura).
1 Present address: Departamento de Matematicas, Estadistica y Computacion, Universidad de Cantabria,
39005 Santander, Spain.
0014-5793/03/$ - see front matter Е 2003 Elsevier Science Ltd. All rights reserved.
doi:10.1016/S0747-7171(03)00013-0
466
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
specific approximations to the roots (mainly asymptotic expansions) are also available;
for example, in Temme (1979) asymptotic approximations to the zeros of first and second
kind Bessel functions, which are refined by a Newton method, are considered.
In this paper we present a Maple program for the evaluation of the zeros of OPs and SFs
based on a globally convergent fixed point method. The program is available upon request
from the authors.
The method applies for SFs and OPs, yn , which are solutions of linear homogeneous
second order ordinary differential equations (ODEs) and satisfy difference-differential
equations (DDEs) of the type:
yn (x) = an (x)yn (x) + dn (x)yn?1 (x)
yn?1
(x) = bn (x)yn?1 (x) + en (x)yn (x)
where the coefficients an (x), bn (x), dn (x) and en (x) are continuous.
This method proves to be efficient and more general than matrix methods. Fixed point
methods apply for arbitrary solutions of second order linear homogeneous ODEs and
not only to polynomial solutions or minimal solutions of a three-term recurrence relation
(TTRR) (which is, for instance, the case for regular Bessel and Coulomb functions). A further advantage of the fixed point method is that the interval for finding zeros can be arbitrarily chosen, contrary to matrix methods. Furthermore, the symbolic part of the algorithm
provides analytical information concerning bounds on the distance between adjacent zeros.
2. Theory
The method for the computation of zeros in an interval I of solutions of second order
ODEs
yn + Bn (x)yn + An (x)yn = 0
with independent solutions
equation
(1) (2)
{yn , yn }
(1)
is based on the existence of a contrast differential
+ Bn?1 (x)yn?1
+ An?1 (x)yn?1 = 0,
yn?1
(2)
(1)
(2)
, yn?1
} such that there exist first order DDEs:
satisfied by independent solutions {yn?1
yn (x) = an (x)yn (x) + dn (x)yn?1 (x)
yn?1
(x) = bn (x)yn?1 (x) + en (x)yn (x),
(3)
(1)
(1)
with continuous coefficients in I , satisfied both by Y (1) ? {yn , yn?1 } and Y (2) ?
(2) (2)
(1) (2)
{yn , yn?1 }. Such DDEs, satisfied simultaneously by two sets Y (1), Y (2), with {yn , yn }
(1)
(2)
and {yn?1 , yn?1 } independent solutions of two different ODEs, will be called general. The
indices n and n ? 1 will normally denote a parameter of the differential equations, but not
necessarily: n and n ? 1 can be interpreted as labels referring to two different ODEs.
It can be shown that such general DDEs exist and are unique for a pair {Y (1), Y (2) } as
described above. The restriction imposed on the contrast functions yn?1 is the continuity
of the coefficients of the DDEs.
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
467
2.1. Properties of the solutions and the coefficients
Properties of the coefficients of the DDEs can be extracted using well-known properties
of the solutions of ODEs. First, dn and en can never cancel because {yk(1), yk(2) } (k =
n, n ? 1) are independent solutions (see Segura, 2002, Lemma 2.1). Indeed, writing the
first equation in (3) for Y (1) and Y (2) :
(1)
yn
(1)
(1)
= an (x)yn + dn (x)yn?1
(2)
yn(2) = an (x)yn(2) + dn (x)yn?1
and therefore dn (x) = W [yn(1), yn(2) ]/Z n (x), with
(1)
(1) yn
yn?1
Z n (x) = (2)
= det[Y 1 , Y 2 ].
(2) yn
yn?1
(4)
Thus dn (x) = 0 ?x because the Wronskian of two independent solutions of a second
(1)
(2)
order ODE never cancels. Similarly en (x) = ?W [yn?1
, yn?1
]/Z n (x). These arguments,
together with the fact that Z n (x) cannot be identically zero (Marx, 1953, Theorem 1), show
that the DDEs (3) exist and are unique (Marx, 1953).
We also see that the continuity (and differentiability) of the coefficients is equivalent to
the condition Z n (x) = 0 ?x ? I .
On the other hand, if the coefficients are continuous then necessarily the zeros of the
solutions yn and yn?1 are interlaced, and if one (non-trivial) solution of the DDEs has two
or more zeros in I then necessarily en dn < 0 (Segura, 2002, Lemma 2.4).
, we
By differentiating the first DDE and using both DDEs to eliminate yn?1 and yn?1
(1)
(2)
arrive at a second order ODE for yn , satisfied by two independent solutions yn , yn . This
gives the relations:
Bn = ?an ? bn ?
dn
;
dn
An = ?an ? dn en + an
dn
+ an bn
dn
(5)
and similarly
Bn?1 = ?an ? bn ?
en
;
en
An?1 = ?bn ? dn en + bn
en
+ an bn .
en
(6)
In Segura (2002), the case Bn = Bn?1 = B(x) was considered in detail. This is not
a restriction, given
that one can always consider a change of the dependent variables
yk = exp(? 12 Bk dx) y?k , k = n, n ? 1; the functions y?k , with the same zeros as yk ,
satisfy ODEs without the first derivative term (ODEs in normal form).
Considering the first equations in (5) and (6) we see that the condition Bn = Bn?1 is
equivalent to dn /en = constant. Therefore, an alternative way to consider the general case
in which Bn = Bn?1 is to renormalize the solutions yk (x) = ?k (x) y?k (x) in such a way that
the functions y?k satisfy the DDEs (3) with dn /en = constant.
468
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
2.2. Transformation of the DDEs
Let us temporarily change ournotation for the DDEs, by denoting by y our ?problem
functions? and by w the ?contrast functions?:
y (x) = ?(x)y(x) + ?(x)w(x)
w (x) = ?(x)w(x) + ? (x)y(x).
(7)
Next we describe the analytical transformations that are required to build globally
convergent fixed point iterations (FPIs) for the computation of the zeros of the solutions y.
We assume that the coefficients are differentiable, that ? and ? never cancel (Segura, 2002,
Lemma 2.1) and that ?? < 0 (Segura, 2002, Lemma 2.4). As we have discussed, these are
general properties for general DDEs having solutions with at least two zeros.
First, the functions are renormalized:
y(x) = ? y (x) y?(x),
w(x) = ?w (x)w?(x)
(8)
with ? y (x), ?w (x) = 0 ?x ? I ; then the DDEs for the renormalized functions are
y? = ?? y? + ?? w?
w? = ?? w? + ?? y?
with
?y
?? = ? ?
,
?y
?w
?? = ? ,
?y
?w
,
?? = ? ?
?w
?? = ?
?y
.
?w
(9)
The renormalizing functions are chosen in such a way that
?? = ??? ,
that is
?? > 0
?
? y = sign(?) ? ?w .
?
With this selection and the change of variables z(x) =
??(x) dx, we obtain
??
y?? = y? + w?
??
??
w?? = w? ? y?
??
where dots denote derivation with respect to z.
Then, the ratio H (z) = y?/w?, with zeros and singularities interlaced (coinciding, up to
a change of variables, with the zeros of y and w), satisfies:
H? = 1 + H 2 ? 2?H
(10)
? = (?? ? ??)/(2??).
(11)
with
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
469
2.3. Global fixed point iterations
From Eq. (10), it can be shown (Segura, 2002) that the iteration
T (z) = z ? arctan(H (z))
(12)
converges globally to the zeros of y(x(z)) in intervals where ? does not change sign.
It is easy to relate the functions H (z), ?(z) with the original functions y, w and the
coefficients of the DDEs:
? y(x(z))
,
H (z) = sign(?) ?
? w(x(z))
(13)
z(x) =
?? ?,
1 ?
?
1
???+
?
,
?(x) = ?
2 ?
?
2 ???
where ?(z) = ?(x(z)).
Putting all these expressions together, the global iteration reads:
? y(x(z))
T (z) = z ? sign(?) arctan ?
.
? w(x(z))
(14)
As mentioned, T (z) is a globally convergent iteration in intervals where ? does not
change sign. Let us denote the successive zeros of y(x) and w(x) by x ym , x wm ; similarly, we
m
denote the zeros of y(x(z)) and w(x(z)) by z m
y and z w , respectively. The zeros of y(x(z))
and w(x(z)) are interlaced; the indices m enumerate the zeros from the smallest to the
largest ones in the interval I and they are chosen in such a way that:
m
m+1
и и и < z m?1
< zw
< zm
< z m+1
< иии.
y
y < zw
y
Given a value x 0 between two consecutive zeros of w, x wm and x wm+1 , and taking as starting
value z 0 = z(x 0 ) we have Segura (2002)
lim T ( p) (z 0 ) = z m
y,
p??
where x ym = x(z m
y ) is the zero of yn between such two consecutive zeros of w.
m
The iteration T (z) is quadratically convergent because T? (z m
y ) = 0 for any zero z y of
y(x(z)).
Observe that the application of the FPI requires that z(x) is an invertible function (that is,
that z(x) is a change of variables). This is so, because the integrand is positive (Eq. (13)).
Of course, the algorithm will be more efficient if the inversion of z(x) can be obtained
analytically; this is the case for all OPs and SFs we have considered so far.
2.4. Forward and backward iterative schemes
Next we describe the iterative scheme to compute all the zeros inside an interval where ?
does not change sign.
Let us, for instance, assume that ? < 0 in a given interval. It can be shown
m
m+1 ? z m < ?/2 (and therefore
(Segura, 2002) that ?/2 < z m
y ? z w and 0 < z w
y
470
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
z m+1
? zm
y
y > ?/2) for all m such that the zeros are inside this interval. Hence, given that
m+1 < z m + ?/2 < z m+1 , the larger zeros subsequent to z m can be computed iteratively:
zw
y
y
y
m + ?/2) (forward sweep). In contrast, when ? > 0, a backward
=
lim
T
(z
z m+1
p??
y
y
sweep is the way to find the zeros. If ? changes sign, a combination of forward and
backward sweeps is possible (Segura, 2002). As in Segura (2002), we will call expansive
sweep the iterative computation of zeros starting from a zero z m
y and evaluating the zeros
smaller than z m
by
a
backward
sweep
and
the
zeros
greater
than
zm
y
y by a forward sweep;
reciprocally, interchanging forward by backward sweep, we have a contractive sweep.
The current version of the algorithm assumes that only one change of sign in ? can
take place (which is a considerably general condition, see Segura, 2002). More general
situations can be described; however, for the purpose of computing zeros of OPs and SFs,
this is a general enough situation.
2.5. The case of uniparametric families of functions
A further restriction of the algorithm is that it has been implemented for the computation
of zeros of families of functions depending on one parameter and functions which can
be related to this kind of function (like Airy functions, which can be related to Bessel
functions of order 1/3 (Segura, 1998)). This is not an intrinsic limitation of the method, as
we previously discussed.
We consider uniparametric families of ODEs yk + Bk (x)yk + Ak (x)yk = 0, where k is
the parameter, with independent solutions {yk(1)}, {yk(2)} verifying
yk = ak yk + dk yk?1
yk?1
= bk yk?1 + ek yk
with continuous coefficients.
The functions may depend on additional parameters, but only one parameter determines
which is the contrast function appearing together with the problem function yn in the
DDEs.
If allowed by the range of k we can write two systems of DDEs for each problem
function yn , taking k = n in:
yk = ak yk + dk yk?1
yk = bk+1 yk + ek+1 yk+1
(15)
yk?1
= bk yk?1 + ek yk ;
yk+1
= ak+1 yk+1 + dk+1 yk .
This means that a TTRR is verified:
yk+1 = rk yk + sk yk?1 ,
(16)
where
rk =
ak ? bk+1
ek+1
and
sk = dk /ek+1 = 0 ?x ? I.
(17)
If the coefficients of the DDEs are continuous then {yk(1)}, {yk(2)} are necessarily
independent solutions of the TTRR because Z k (x) = 0 ?x, k (Eq. (4)). TTRRs are a useful
tool to compute SFs and OPs.
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
471
Another interesting fact is that we have two alternative FPIs by choosing as contrast
functions w = yn?1 or w = yn+1 . In one case we take (Eqs. (13) and (15)):
y = yn ,
w = yn?1 ,
? = an ,
? = bn
? = dn ,
? = en
(18)
and for the other one
y = yn ,
? = en+1 ,
w = yn+1 ,
? = dn+1 .
? = bn+1 ,
? = an+1
(19)
More explicitly, denoting, as in Segura (2002):
n + 1, i = +1
ni =
n,
i = ?1,
two different FPIs Ti (z) = z ? arctan(Hi (z)) can be considered to compute the zeros of yn ,
where
yn (x(z ni ))
y?n (z)
= i sign(dni )K i
,
y?n+i (z)
yn+i (x(z ni ))
dni i/2
Ki = ?
,
en i
zi =
?dni eni dx,
dn i
1
1 en i
an i ? bn i +
.
?i = i ?
2 en i
dn i
2 ?dni eni
Hi =
(20)
2.6. Improvement of the iteration step
From Eq. (20) one can expect that, generally speaking, when ?+1 > 0 then ??1 < 0
and vice-versa. This means that both forward and backward sweeps (or expansive and
contractive) are generally available to compute the zeros of a function yn .
The symbolic algorithm selects the most appropriate iteration depending on the
monotonicity properties for the second order ODE (in normal form)
е
y?(z)
+ A?(z) y?(z) = 0,
A?(z) = 1 ? ?? + ?2 ,
satisfied by
1
y? = exp ? 2 (?? + ??) dz y?,
where ?? = ??/??, ?? = ??/?? (see Section 2.2). Given the continuity of the coefficients of the
DDEs, the functions y?n (z) have the same zeros as y?.
?
If, for instance, A?(z)
< 0 in a given interval, then the spacing between the zeros of y?(z),
and hence of y(x(z)), increases for increasing z, that is:
m?1
zm
< z m+1
? zm
y ? zy
y
y
472
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
m
m
m?1 < ?/2; then
and if an iteration with ? < 0 has been chosen: z m
y ? z w > ?/2, z w ? z y
it is not difficult to see that
m+1
m
m?1
m+2
zw
< zm
) < z m+1
< zw
.
y + (z y ? z y
y
m
m = (z m ? z m?1 ), provides
This means that the starting value z m
y + z , with z
y
y
m+1 and z m+2 ) and that the step
(which is the zero between z w
convergence to z m+1
y
w
z m improves the default step ?/2 because z m > ?/2; therefore it provides a closer
.
(under)estimation of z m+1
y
This suggests that using the difference of two previously evaluated zeros as iteration
steps instead of ▒?/2 will improve the performance if an iteration (i = +1 or ?1) can be
?
chosen such that A?(z)?
> 0.
For instance, for the case of OPs we observe that A?(z) has a maximum and then we
should choose an expansive iteration.
When possible, the algorithm will consider this improvement of the iteration step.
3. Description of the program
As previously described, the algorithm is based on the existence of two linear first order
DDEs. However, we choose as symbolic input the coefficients of a TTRR and of one of
the DDEs. Of course, both choices are equivalent (Eqs. (16) and (17)). The reason for
giving the TTRR and one DDE is that this option is more often available in books (see
Abramowitz and Stegun, 1972).
The algorithm has several modes of operation. The user can select one of the already
implemented cases for the computation of zeros or introduce other recurrence relations.
The functions which are implemented in the present version of the program are: general
Bessel functions cos ? J? (x) ? sin ?Y? (x) (we denote n = ?), regular Coulomb wave
functions FL (?, ?) (notation: n = L, ? = ?, x = ?), first kind Conical functions
n
P?1/2+i?
(x) (x > 1), where i is the imaginary unit, Hermite polynomials Hn (x), Legendre
(?)
polynomials Pn (x), generalized Laguerre polynomials L n (x), Gegenbauer polynomials
(?,?)
(?)
(x).
Cn (x) and Jacobi polynomials Pn
If one of these functions is considered, we do not need to introduce the coefficients of a
TTRR and a DDE.
The algorithm produces both symbolic and numerical results. The symbolic results can
be useful to implement algorithms for the computation of SFs in compilable languages like
Fortran. This is relevant when the speed of the algorithms is of concern or when the zeros
for different sets of parameters of a same function are needed; a Fortran template will also
be made available in the near future. Also, as described in Segura (2003), this information
can be used to study interlacing properties of the zeros.
The numerical results are the zeros of the function under consideration in the selected
interval. One can select performing only the symbolic task or both the symbolic and
numerical tasks.
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
473
3.1. Inputs
We are prompted by the program to provide the following inputs:
1. Number of digits for the calculation.
2. Selection of one of the implemented functions or introduction of coefficients.
3. If we choose to introduce a function not yet implemented, the program asks for:
(a) Interval of definition of the function.
(b) Number of parameters of the function (the first to be introduced is the parameter
n that is incremented in the TTRR and the DDE).
(c) Name of the parameters and permissible range of the parameters.
(d) Explicit expressions of the coefficients in the TTRR yn+1 = rn yn + sn yn?1 and
the DDE yn = an yn + dn yn?1 . n is one of the parameters that have been defined
before.
4. Values of the parameters: We must fix one particular set of values for the
parameters. Even if only the symbolic computation is required we must give a set
of representative values. The reason for this is that the iteration to be chosen (i = +1
or ?1) may depend on the actual values considered and even if this is not so, the
discussion of the monotonicity properties of A?(z) when free parameters are present
becomes a difficult task.
5. Selection of strictly symbolic or symbolic/numerical computations.
6. If the symbolic/numerical option is chosen:
(a) State the interval (contained in the interval of definition of the functions) where
the zeros are sought.
(b) Choose the method of computation of the ratios yn /yn▒1 . Three possibilities are
given:
(ii.a) Finite continued fraction: This method applies for dominant solutions of
the corresponding TTRR and it corresponds to the forward evaluation of
the recurrence relation. Also, when the TTRR has no dominant solutions,
this option can be generally used when a moderate number of iterations are
needed. In this case, the ratio of starting values y0 /y1 must be given from
which, by forward recursion, the values yn /yn▒1 can be computed. This
method is applied for the OPs already implemented in the algorithm.
(ii.b) Infinite continued fraction: This can be the method of choice when
minimal solutions of the corresponding TTRR are considered. This is the
case, for instance, for the (already implemented) functions: Bessel function
n
(i is
Jn (x), regular Coulomb wavefunction and Conical function P?1/2+i?
the imaginary unit). The continued fraction is computed using the modified
Lenz?Thompson algorithm (Press et al., 1992).
(ii.c) Direct computation: When the functions under consideration are implemented in Maple, Maple can directly compute the ratio yn /yn+1 . The
program then asks us to provide the function yn in the form of the
corresponding Maple command; if the function yn depends on additional
474
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
parameters in addition to n, the values of these extra parameters must be
fixed when the function is provided. The method of direct computation can
be applied to all the already implemented functions in the algorithm except
for the Coulomb functions and Conical functions (not yet implemented in
Maple). However, the algorithms work more efficiently if one of the two
previous approaches can be used. For the already implemented functions,
the option of direct computation is considered for Bessel functions C? (x) =
cos ? J? (x) ? sin ?Y? (x) when cos ? = 0.
(c) Relative precision required, which should never be smaller than 10?Digits (the
number of digits is the first input described above).
(d) Maximum number of iterations allowed.
In the Appendix, we provide an example of a session with the algorithm zeros.mpl.
3.2. Basic tasks
The program zeros.mpl performs the following tasks:
? Calculation of the coefficients bn and en of the difference relation:
yn?1
= bn yn?1 + en yn
?
?
?
?
?
using the coefficients of the TTRR yn+1 = rn yn + sn yn?1 and the differential
relation yn = an yn + dn yn?1 (Eq. (17)).
Check the continuity of the coefficients of the DDEs and the inequality en dn < 0.
Check also the consistency of Eqs. (5) and (6); An and An?1 are computed and it is
verified that they are identical with a shift in the parameter n; the same is done for
Bn . In case the checks fail, an error message appears and the code stops its execution.
Normalization of the solutions; calculation of: change of variable z(x), ?i of the
normalization K i and A?(z), as described in Section 2.5 (initially assuming that the
iteration is i = ?1).
Computation of the extrema of A?(z) and selection of the appropriated iteration
(i = +1, ?1) depending on the monotonicity properties of A?(z) (Section 2.6).
Re-evaluation (if the selected iteration turns out to be i = +1) of the basic functions
of the algorithm (change of variable, ?i , etc).
Selection and initialization of the sweep (forward, backward, contractive or
expansive) depending on the location of the value x ? (the zero of ?i (x)) if it exists
and the sign of ?i . For instance, if ?i is negative in the interval under consideration a
forward sweep is the option.
3.3. Output
The numerical outputs of the code are the zeros computed in the interval and the number
of iterations applied to compute each zero.
The symbolic output is displayed before entering the numerical procedures and it
provides the information needed to implement the fixed point method for the problem
function, namely:
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
475
1. The change of variables z(x) applied.
2. The coefficient A?(z) for the second order ODE in normal form in the z variable (for
brevity, it is expressed in terms of x, which is a function of z).
3. The parameter ?i (x(z)), also expressed in terms of x.
4. The normalization factor K i .
5. Sign of the dni coefficient of the DDE.
6. Type of FPI (index i = 1 or ?1).
Outputs 1, 4, 5 and 6 are all that is needed to implement the numerical method.
The type of sweep (forward/backward, contractive/expansive) depends, as discussed in
more detail in Segura (2002), on the sign of ?i .
4. Examples of symbolic output
Table 1 summarizes the symbolic
For brevity, A?(z)
? outputs for the implemented cases.
?
(n + 1)(n + 1 + ?), Nn,? =
(n + 1 + ?)/(n + 1),
is not shown, where Mn,? =
Rn,? = 2(n + 1) + ? and ?n,? = ((n ? 1/2)2 + ? 2 )?1/2 . Of course the same results
hold for a general solution of the system of DDEs. For instance, the same results apply
to general Bessel functions cos ? Jn (x) ? sin ?Yn (x), for combinations of the regular and
irregular Coulomb functions, and so on.
Table 1
Examples of symbolic outputs provided by the code
yn (x)
z(x)
?
Hn (x)
2(n + 1)x
(Hermite)
Pn (x)
(Legendre)
(?)
L n (x)
(?)
Ki
sign(d)
i
?
2n + 2
+
+1
(n + 1) tanh?1 x
?x
1
+
+1
Mn,? log(x)
2n + 2 + ? ? x
?
2 (n + 1 + ?)(n + 1)
Nn,?
?
+1
+
+1
Nn,?
Nn+?,?
+
+1
(Laguerre)
Cn (x)
?(x)
x
2
?
2 n+1
?
(n + 2?)(n + 1) tanh?1 x
(Gegenbauer)
x(2n + 1 + ?)
??
(n + 2?)(n + 1)
x(Rn,?+? )2 ? ? 2 + ? 2
4Mn,? Mn+?,?
n + 2?
n+1
2Mn,? Mn+?,?
tanh?1 x
Rn,?+?
?
1
cosh?1 x
?n,?
(n ? 1/2)x
?n,?
x2 ? 1
?n,?
?
?1
Jn (x)
x
1
+
?1
(Bessel)
(2n ? 1)
2x
x n2 + ? 2
n
n2 + ? x
x n2 + ? 2
1
+
?1
(?,?)
Pn
(x)
(Jacobi)
n
P?1/2+i?
(x)
(Conical)
Fn (? , x)
(Coulomb)
476
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
All these cases have in common that the normalization factor K i does not depend
on x. As discussed, this reflects the fact that the coefficient Bn (x) of the ODE does not
depend on n and, consequently, the dn and en coefficients are such that dn /en is a constant.
However, the algorithm is able to deal with the more general situation in which Bn (x)
depends on n.
For example, the Kummer function M(a, c; x) satisfies the differential equation
x y + (c ? x)y ? ay = 0
and therefore the differential equations for the family of functions yn (x) = M(a + n, b +
n; x) have coefficients Bn (x) depending on n. In this case, the program gives the following
symbolic output:
z(x) = 2 x(1 ? n ? a);
?
x(1 ? n ? a)
;
Ki =
|b + n ? 1|
?(x) = ?
1 3 ? 2b ? 2n + 2x
?
4
x(1 ? n ? a)
sign(d) = ?1;
i = ?1
and we see that, as expected, K i depends on x.
The type of sweep for each family of functions is easily read from the obtained
expression of ?i , where i has been chosen, as described, to improve the iteration step
by using the monotonicity properties of the function A?(z). For instance, for OPs there is a
transition point (x ? for which ?i (x ? ) = 0) inside their support and ?i is positive on the right
of the interval. This means that expansive sweeps are the option, which is coherent with
the fact that A?(z) has a maximum inside the interval of definition. For Bessel (x > 0) and
Conical functions (x > 1), a forward (backward) sweep would be considered for n < 1/2
(n > 1/2). On the other hand, for Coulomb functions (x > 0) there is the possibility
of backward or expansive sweeps depending on the values of the parameters. In all cases
described with a transition point (x ? ), the expansive sweep will be replaced by a forward or
backward sweep for certain selections of the interval; for instance, if we choose an interval
[x 1 , x 2 ] such that x 1 > x ? with x ? the transition point, then the expansive sweep is replaced
by a forward sweep.
It is interesting to note that the change of variables for OPs maps the support of the
polynomials to (??, +?). This means that the change of variables becomes singular
at the ends of the support, which go to infinity. Therefore, it is expected that the zeros
approaching the ends of the support interval are computationally the most demanding. This
effect becomes more apparent as the order increases and the rest of the parameters approach
singular values (for example, ? ? ?1+ for Laguerre, ? ? ?1/2+ for Gegenbauer and
? ? ?1+ , ? ? ?1+ for Jacobi). To deal with the computation of the extreme zeros of
OPs for extreme values of the parameters, it is convenient to fix the maximum number of
iterations to higher values than the recommended 50 iterations.
On the other hand, the difference equations over the order n are not necessarily the most
efficient option for the computation of zeros of OPs. For instance, by taking into account
the relation of Laguerre polynomials with confluent hypergeometric functions:
n!
L ? = M(?n, ? + 1, x)
(n + ?)n n
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
477
the problem can be translated into finding the zeros of M(a, b, x) for the particular values
of the parameters corresponding to a Laguerre polynomial. We have many possibilities
of defining our sequences of functions yn , an option is yn (x) = M(a + n, b + n, x).
It turns out, as we have already shown, that the associated change of?variables for this
option is not singular as x ? 0+ and that it behaves like z(x) ? x. This selection
behaves much better than the already implemented version for the smallest zero. In a future
publication (Gil et al., in preparation), we will study in detail the use of different sequences
of functions to compute zeros of hypergeometric and confluent hypergeometric functions
depending on the values of the parameters and the range of the zeros.
5. Performance of the code
Let us discuss the performance of the code for the group of polynomials and SFs already
implemented in the algorithm.
The main advantage of Maple algorithms compared to fixed precision programming
languages is that the number of digits in the computation can be chosen at will. This
feature can be used to check the accuracy and efficiency of our algorithm. We have
computed the zeros of OPs with 10?100 relative precision and checked their accuracy by
computing the polynomials at the zeros, p(x z ), and comparing with the estimated value
p(x z (1 + )) p(x z ) + x z p (x z ) = x z p (x z ) (where is the relative error in the
computation of the zero x z ). Typically, the algorithm needs around 10 iterations on average
for each zero, although slower convergence is observed in the extreme cases described in
the previous section. We obtain full agreement within the precision considered in all cases
when enough digits (120?130) are considered for orders up to 100. For even higher orders
and higher precision we expect that the algorithm will also work accurately. The same
check was considered for Bessel functions. Conical functions and Coulomb wavefunctions
are not yet implemented in Maple and this direct check cannot be performed.
Next, we will show the number of required iterations for the computation of the zeros for
a certain selection of parameters. In all cases, the precision demanded for the calculation of
the zeros is 10?12; however, the convergence is so fast that the number of iterations required
increases by few units when much higher precisions are required (for instance 1?2 extra
iterations are required for 10?40 relative precision). By using asymptotic expansions as
initial guesses for the zeros, the performance could be improved in some cases. However,
the algorithm is intended to work under very general circumstances and for arbitrary
solutions of the corresponding ODEs while asymptotic approximations are specific for
particular solutions. For a fast algorithm for the zeros of first and second kind Bessel
functions based on asymptotic expansions see Temme (1979). As we will next show, our
algorithm, being quite general, is also very efficient (see, for instance, Fig. 6).
Figs. 1?8 show the number of iterations needed in order to compute the zeros of Hermite
(Fig. 1), Legendre (Fig. 2), generalized Laguerre (Fig. 3), Gegenbauer (Fig. 4) and Jacobi
(Fig. 5) polynomials, for selected values of the parameters, as a function of the position
n
of the zeros. Figs. 6?8 correspond to Bessel (J? (x)), Conical (P?1/2+i?
(x)) and Coulomb
wavefunctions (Fn (?, x)) respectively. All figures correspond to 10?12 relative precision.
478
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
10
9
Number of iterations
8
7
6
5
4
3
2
1
0
? 10
?5
0
5
10
x
Fig. 1. Number of iterations for computing the zeros of H30 (x).
Figs. 1, 2, 4 and 5 show a similar pattern in the distribution of iterations needed to
compute each zero. In all four cases ? is zero at x = 0 and expansive sweeps are considered
starting from this point. If we had chosen ? = ?, Fig. 5 would not be symmetric; the
transition point would move apart from the origin. The method is expected to converge
better as ? is smaller; in fact, if ? = 0 the ratio Hi is a tangent function in z (see Eqs. (10)
and (20)) as, for instance, happens for Bessel functions of order n = 1/2. In this case the
algorithm gives the correct answer in one iteration. This explains why the smallest zeros
in modulus are the fastest to compute. As we move from x = 0 the number of iterations
increases, and the largest zeros in modulus are the most demanding ones, as we previously
discussed. The fact that close to x = 0 small peaks appear in the figures is an effect of the
improvement of the iteration step, which becomes effective after the zeros corresponding
to these peaks have been evaluated.
For the Laguerre case (Fig. 3) the transition point is at a positive value of x and the
increase in the number of iterations is faster as we move to the left from this transition
point.
The pattern in Fig. 6 (Bessel functions) and 7 (Conical functions) is similar with regard
to the variation in the number of iterations and corresponds to a backward sweep, with a
slight increase in the number of iterations as x becomes smaller. Of course, the distribution
of zeros is quite different in each of these two cases, as could be expected given the
different change of variables z(x) associated to each case. The largest zeros require more
iterations because the improved iteration step starts to operate after these initial zeros have
been computed. Finally, Fig. 8 (Coulomb functions) corresponds again to an expansive
sweep.
Our numerical implementation of the fixed point method proves to be an efficient
algorithm for the evaluation of the zeros of SFs and it appears to be at least a factor 2 faster
than the intrinsic Maple procedures for computing the zeros of Bessel functions. Besides, it
is important to point out that external programs do not have the same privilege as intrinsic
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
479
25
Number of iterations
20
15
10
5
0
?1
? 0.5
0
0.5
1
x
Fig. 2. Number of iterations for computing the zeros of P30 (x).
25
Number of iterations
20
15
10
5
0
0
20
40
60
80
100
120
x
(3/2)
Fig. 3. Number of iterations needed for computing the zeros of L 30
(x).
Maple procedures. Also, the comparison of similar methods (Segura and Gil, 1999)
(implemented in Fortran) with other available codes (Vrahatis et al., 1995), was favourable.
As for OPs, the Maple procedure fsolve is an interesting tool to compute such zeros.
For moderately large orders we observe that our code tends to be faster, in spite of the
fact that, as commented, it is difficult to compare the speed of external algorithms with
intrinsic procedures. For instance, our algorithm is faster than fsolve (using Maple 7) for the
computation of the zeros of Laguerre polynomials of moderately large order and negative
(?1/6)
(x).
parameter, like L 50
The fixed point method converges with certainty even for extreme cases like, for
(?0.99,?0.99)
instance, when computing the zeros of P100
(x). For this polynomial, the
480
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
25
Number of iterations
20
15
10
5
0
?1
? 0.5
0
0.5
1
x
(3/2)
Fig. 4. Number of iterations needed for computing the zeros of C30
(x).
25
Number of iterations
20
15
10
5
0
?1
? 0.5
0
0.5
1
x
(3/2,3/2)
Fig. 5. Number of iterations needed for computing the zeros of P30
(x).
convergence of fsolve is uncertain if a moderate number of digits is chosen (20 digits),
producing zeros outside the interval [?1, 1]. For a larger number of digits (50 digits), this
problem disappears but the convergence rate becomes very slow. This shows that fsolve
(Maple 7) becomes quite unstable for these extreme cases (and for Maple V, fsolve fails to
converge for more than 16 digits), losing too many significant digits. These problems are
not observed for the fixed point method.
Furthermore, we expect to improve the performance for the computation of the extreme
zeros in a forthcoming publication (Gil et al., in preparation). As discussed for Laguerre
polynomials regarding the calculation of the smallest zeros, alternative sequences of
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
481
10
9
Number of Iterations
8
7
6
5
4
3
2
1
0
0
20
40
60
80
100
x
Fig. 6. Number of iterations needed for finding the zeros of J10 (x) in the interval [1, 100].
10
9
Number of Iterations
8
7
6
5
4
3
2
1
0
0
20
40
60
80
100
x
10
Fig. 7. Number of iterations needed for finding the zeros of P?1/2+i10
(x) in the interval [1.1, 100].
functions solve the relatively slow convergence for the extreme zeros when extreme values
of the parameters are considered. The comparison will be even more favourable then.
6. Conclusions
We have described a combined symbolic/numerical algorithm to compute the zeros of
families of functions satisfying a TTRR together with a first order linear DDE. These
conditions are satisfied by a broad family of SFs and OPs. Therefore, the method is
applicable for the computation of the zeros of classical OPss and non-polynomial solutions
of the same ODEs. The algorithm is also successful in computing the zeros of other SFs
482
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
10
9
Number of Iterations
8
7
6
5
4
3
2
1
0
0
20
40
x
Fig. 8. Number of iterations needed for finding the zeros of F20 (?20, x) in the interval [1, 50].
like, for instance, arbitrary solutions of the Bessel equation, the Coulomb wave equation,
Hermite equation (parabolic cylinder functions), Legendre equation with complex degree
(Conical functions) and Airy functions (whose zeros can be computed via the connection
of Airy functions with Bessel functions (Segura, 1998)), among others.
In the literature, there exist several algorithms for the computation of zeros of
the Bessel functions J? and Y? (Piessens, 1990; Segura and Gil, 1999; Temme, 1979;
Vrahatis et al., 1995), but no methods and algorithms are available for other SFs,
although matrix methods have been described to solve the problem for regular Coulomb
wavefunctions (Ikebe, 1975; Ikebe et al., 1991) and a global Newton method was described
in Segura (2001) to compute the zeros of general solutions of the Bessel equation. The
algorithm here described is able to compute all these zeros efficiently.
As a bonus, the method, being based on analytical properties, provides analytical
information on the spacing between the zeros of a given sequence of functions. It also
provides all the information needed to implement the algorithms in compilable languages
like Fortran. A Fortran template will also be provided to build specific programs for the
computation of the zeros of families of functions within this (vast) set of holonomic
functions.
Fixed point methods to compute the turning points of this type of function
(Gil and Segura, submitted) are also available and will be implemented in the next update
of the algorithm.
Acknowledgements
A. Gil acknowledges support from the A. von Humboldt foundation. J. Segura
acknowledges support from DAAD. We thank the editors and the anonymous referees for
valuable suggestions.
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
Appendix
> read ?zeros.mpl?;
This package computes the zeros of Special Functions
and Orthogonal Polynomials
(c) Amparo Gil (Dept. Matematicas, Universidad Autonoma de Madrid, Spain)
and
Javier Segura (Dept. Matematicas, Universidad Carlos III de Madrid, Spain)
Introduce number of digits for the computation
> 20;
Is your function one of the following:
? 1. Bessel functions
c Jn (x) ? 1 ? c2 Yn (x), c in the interval [?1, 1]
2. Conical functions
Pn,?1/2+i,? (x)
3. Coulomb functions
Fn, ? (x)
4. Gegenbauer polynomials
C?,n (x)
5. Generalized Laguerre polynomials
L ?,n (x)
6. Hermite polynomials
Hn (x)
7. Jacobi polynomials
P?,?,n (x)
8. Legendre polynomials
Pn (x)
YES ? ?> 1
NO ? ?> 2
> 1;
Which is your function. Select 1, 2, . . . , 8
> 4;
Values of the parameters of the functions
Parameter n
> 10;
Parameter ?
> .5;
The user must now choose performing only the symbolic tasks
or both symbolic and numerical tasks
0 ? ? > SYMBOLIC
1 ? ? > SYMBOLIC + NUMERIC
> 1;
Provide the following information
a) Interval for the calculation of zeros [xa, xb]
Remember that [xa, xb] must be in the interval ]?1, 1[
483
484
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
> -.9;
> .9;
c) Precision required for the computation of the zeros
10^ (-12);
d) Maximum number of iterations allowed (recommended: 50)
> 50;
First check ok
Second check ok
A has a maximum
1. Change of variables z(x) := (n + 2?)(n + 1) arctanh(x)
2. A := ?
1 ?4n 2 ? 8n? + 4x 2 n 2 + 8x 2 n? ? x 2 + 4x 2? 2 + 2 ? 4?
4
(n + 2?)(n + 1)
1 x(2n + 1 + 2?)
3. Parameter eta(x) := ? ?
?
2 n + 2? n + 1
n + 2?
4. Normalization factor K i :=
n+1
5. Iteration 1
6. Sign of dn i , 1
NUMERICAL OUTPUTS:
1) TYPE OF SWEEP:
Expansive
2) ZEROS:
xcero[
xcero[
xcero[
xcero[
xcero[
xcero[
xcero[
xcero[
1
2
3
4
5
6
7
8
]:=
]:=
]:=
]:=
]:=
]:=
]:=
]:=
-.86506336668898451072
-.67940956829902440623
-.43339539412924719080
-.14887433898163121089
.14887433898163121089
.43339539412924719080
.67940956829902440623
.86506336668898451072
Number of iterations needed for each zero:
it[
it[
it[
it[
it[
it[
it[
it[
1
2
3
4
5
6
7
8
]:=
]:=
]:=
]:=
]:=
]:=
]:=
]:=
8
8
6
4
4
6
8
8
Number of zeros found := 8
A. Gil, J. Segura / Journal of Symbolic Computation 35 (2003) 465?485
485
References
Abramowitz, M., Stegun, I.A. (Eds.), 1972. Handbook of Mathematical Functions With Formulas,
Graphs, and Mathematical Tables. Reprint of the 1972 ed. A Wiley?Interscience Publication.
Selected Government Publications. John Wiley & Sons, Inc, New York; National Bureau of
Standards, Washington, D.C.
Gil, A., Koepf, W., Segura, J., Numerical algorithms for the computation of zeros of hypergeometric
functions (in preparation).
Gil, A., Segura, J., Computing zeros and turning points of special functions from fixed point methods
(submitted for publication).
Golub, G.H., Welsch, J.H., 1969. Calculation of Gauss quadrature rules. Math. Comp. 23, 221?230.
Ikebe, Y., 1975. The zeros of regular Coulomb wave functions and of their derivatives. Math. Comp.
29, 878?887.
Ikebe, Y., Kikuchi, Y., Fujishiro, I., 1991. Computing zeros and orders of Bessel functions. J. Comput.
Appl. Math. 38, 169?184.
Marx, I., 1953. On the structure of recurrence relations II. Michigan Math. J. 2, 99?103.
Piessens, R., 1990. On the computation of zeros and turning points of Bessel functions. Bull. Greek
Math. Soc. 31, 117?122.
Press, W.H., Teukolski, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes in Fortran.
Cambridge University Press, Cambridge.
Segura, J., 1998. A global Newton method for the zeros of cylinder functions. Numer. Algorithms
18, 259?276.
Segura, J., 2001. Bounds on differences of adjacent zeros of Bessel functions and iterative relations
between consecutive zeros. Math. Comp. 70, 1205?1220.
Segura, J., 2002. The zeros of special functions from a fixed point method. SIAM J. Numer. Anal.
40 (1), 114?133.
Segura, J., 2003. On the zeros and turning points of special functions. J. Comput. Appl. Math. (in
press).
Segura, J., Gil, A., 1999. ELF and GNOME: two tiny codes to evaluate the real zeros of the Bessel
functions of the first kind for real orders. Comput. Phys. Comm. 117, 250?262.
Temme, N.M., 1979. An algorithm with ALGOL 60 program for the computation of the zeros of
ordinary Bessel functions and those of their derivatives. J. Comput. Phys. 32, 270?279.
Vrahatis, M.N., Ragos, O., Skiniotis, T., Zafiropoulos, F.A., Grapsa, T.N., 1995. RFSFNS: a portable
package for the numerical determination of the number and the calculation of roots of Bessel
functions. Comput. Phys. Comm. 92, 252?266.
Wilf, H.S., 1962. Mathematics for the Physical Sciences, Wiley, New York (Problem 9).