A naive procedure for computing angular spheroidal functions
Abstract
An algorithm for computing eigenvalues and eigenfunctions of the angular spheroidal wave equation, based on a known but scarcely used method, is developed. By requiring the regularity of the wave function, represented by its series expansion, the eigenvalues appear as the zeros of a one variable function easily computable. The iterative extended Newton method is suggested as especially suitable for determining those zeros. The computation of the eigenfunctions is then immediate. The usefulness of the method, applicable also in the case of complex values of the “prolateness” parameter, is illustrated by comparing its results with those of procedures used by other authors.
MSC [2010] 33E10 (Primary) 33F05, 34L16, 65D20 (Secondary)
PACS 02.30.Hq 02.30.Gp 03.65.Ge
1 Introduction
The usefulness of spheroidal functions in many branches of Physics, like Quantum Mechanics, General Relativity, Signal Processing, etc., is well known and it does not need to be stressed. Due to that usefulness, the description of the spheroidal equation and of the main properties of its solutions deserves a chapter in handbooks of special functions like that by Abramowitz and Stegun [2, Chap. 21], the best known one, or the NIST Digital Library of Mathematical Functions [18, Chap. 30], the most recent one.
A review of the procedures used in the past century for obtaining the eigenvalues and eigenfunctions of the spheroidal wave equation can be found in a paper by Li et al. [12], where also an algorithm, implemented with the software package Mathematica
, is provided.
In the present century, articles dealing with the solutions of the angular spheroidal wave equation have continued appearing. Without aiming to be exhaustive, let us mention the papers by Aquino et al. [3], Falloon et al. [7], Boyd [5], Barrowes et al. [4], Walter and Soleski [27], Abramov and Kurochkin [1], Kirby [14], Karoui and Moumni [13], Gosse [9], Tian [26], Rokhlin and Xiao [24], Osipov and Rokhlin [19], Ogburn et al. [16] and Huang et al. [11], and the books by Hogan and Lakey [10], and by Osipov, Rokhlin and Xiao [20].
Different strategies have been used to solve the angular spheroidal wave equation. The classical procedure starts with the angular spheroidal wave function written as a series of solutions of another similar differential equation, commonly the Legendre one, with coefficients obeying a three term recurrence relation. The resulting expansion becomes convergent only when such coefficients constitute a minimal solution of the recurrence relation. The eigenvalue problem encountered in this way is solved either as a transcendental equation involving a continued fraction, or written in a matrix form. Procedures based on the direct solution of the angular spheroidal equation, without having recourse to comparison with other differential equations, have been less frequently used. The relaxation method proposed by Caldwell [6] and reproduced, as a worked example, in the Numerical Recipes [23, Sec. 17.4], and the finite difference algorithm, described in the recently appeared paper by Ogburn et al. [16], deserve to be quoted. Here we suggest to follow a procedure, based also on the direct treatment of the spheroidal equation, which benefits from an idea that can be found in a paper by Skorokhodov and Khristoforov [25] dealing with the singularities of the eigenvalues considered as function of the (complex) prolateness parameter . A shooting method is used. But, instead of imposing the boundary conditions to a numerically obtained solution, algebraic regular solutions around the regular point or around the regular singular point are written. Smooth matching of both solutions, i. e. cancelation of their Wronskian, at any point determines the eigenvalues. In our implementation of the procedure, we choose as matching point.
A discomfort, when dealing with spheroidal wave functions, is the lack of universality of the notation used to represent them. The Digital Library of Mathematical Functions [18, Chap. 30] provides information about the different notations found in the bibliography. Here we adopt, for the eigenvalues and eigenfunctions, the notation of the Handbook of Mathematical Functions [2, Chap. 21]. The same notation is used in Ref. [12], a paper whose results we will try to reproduce, for comparison, with the method here developed.
In the next section, we recall the angular spheroidal equation and write its solutions in the form of power series expansions around the origin and around the singular point . The procedure for computing the eigenvalues is presented in Section 3. The results of its application in some typical cases are also given. Section 4 shows that normalized eigenfunctions can be trivially obtained. Some figures illustrate the procedure. A few final comments are contained in Section 5.
2 The differential equation
The angular spheroidal wave function , defined in the interval , satisfies the differential equation [2, Eq. 21.6.2]
(1) 
stemming from the separation of the wave equation in spheroidal coordinates, with separation constants and . Periodicity of the azimuthal part of the wave restricts the values of to the integers and, given the invariance of the differential equation in the reflection , only nonnegative integer values of need to be considered. The other separation constant, , commonly referred to as eigenvalue, must be such that becomes finite at the singular points . Their different values, for given and , are labeled by the integer . In most applications, the external parameter is real, positive in the case of prolate coordinates and negative for oblate ones. There are, however, interesting cases corresponding to complex values of [1, 4, 12, 16, 17, 25].
Instead of solving directly Eq. (1), it is convenient to introduce the change of function
(2) 
and to solve the differential equation
(3) 
where
(4) 
is considered as the new eigenvalue.
Two independent solutions about the ordinary point , valid in the interval , are
(5) 
with coefficients given by the recurrence relation
(6)  
Obviously, and are respectively even and odd functions of .
Solutions about the regular singular point can also be written. In terms of the variables
(7) 
the differential equation (3) turns into
(8) 
The solution of this equation which makes to be regular at is, except for an arbitrary multiplicative constant,
(9) 
with coefficients given by
(10)  
In terms of the variable , this regular solution, valid for , is
(11) 
3 The eigenvalues
The problem of finding the eigenvalues , for given and , reduces to require the regularity of at . (The regularity at is then implied by the symmetry of .) In the particular case of being , the problem can be solved algebraically. The recurrence relation in (6) reduces in this case to
(12) 
Obviously, the series in the right hand side of (5) is divergent for unless the value of is such that one of the of even subindex, say (with ), becomes zero, in which case turns out to be a polynomial of degree . This happens for
(13) 
By using the habitual notation
(14) 
for the degree of the polynomial, we obtain for the eigenvalues in the case of
(15) 
that is, in view of (4),
(16) 
For , a convenient way of guaranteeing the regularity of at is to require the cancelation of the Wronskian of and ,
(17)  
at an arbitrarily chosen point of the interval . From the computational point of view, an interesting choice of seems to be , in which case
(18) 
The set of coefficients and are solutions of the difference equations (6) and (10), respectively. According to the PerronKreuser theorem on difference equations [15, 21, 22],
(19) 
that is, for any given , constants and can be found such that
(20) 
This makes evident that the series in the right hand side of Eq. (18) converges as fast as the geometric series . We consider, however, that a better choice of the value of in Eq. (17) is . In this case,
(21) 
Obviously, the cancelation of this Wronskian occurs when either or vanish at the origin, as it occurs for respectively even or odd functions of . Needless to say, the right hand side of (21) depends on the variable , introduced in (4), through the coefficients . Therefore, the problem of finding the eigenvalues of the angular spheroidal equations, for given and , reduces to the determination of the zeros of the function
(22) 
where we have indicated the dependence of the on .
Different procedures can be used in the determination of the zeros of . A useful iterative method is the Newton one. Starting with an initial approximate value, , of a certain zero, repeated application of the algorithm
(23) 
allows one to get the value of the zero with the desired accuracy. The extended Newton method, which uses
(24) 
is even more efficient. For the first and second derivatives of with respect to we have the expressions
(25)  
(26) 
The first and second derivatives of the coefficients are easily obtained by means of the recurrence relations
(27)  
(28)  
stemming from (10). From these difference equations, inequalities analogous to the second one in (20) can be deduced. Such inequalities guarantee the convergence, as fast as the geometric series , of the series in the right hand sides of Eqs. (22), (25), and (26).
We have applied the procedure just described for obtaining the behaviour of the lowest eigenvalues of the spheroidal equation when the parameter varies in the interval and for values of 0, 1, and 2. The results are shown in Figs. 1 to 3. (Remember that .) A glance at Fig. 1 suggests a quasiconfluence of trajectories of the eigenvalues and for sufficiently large negative values of . One may conjecture that, for larger negative values of other pairs of trajectories, those of and , present such quasiconfluence. Table 1 shows that this is the case, and that a similar phenomenon occurs for other values of .
In order to compare with results published by other authors, we have applied our method to the computation of for a sample of values of the parameters considered by Li et al. [12] and by Ogburn et al. [16]. The comparison, shown in Table 2, allows one to conclude that, for moderate real values of , the procedure used in Ref. [12] is more reliable than the finite difference algorithm of Ref. [16].
Ref. [12]  Ref. [16]  this work  

1  4  11  131.5600809  131.560080918303  131.56008091940694 
0.1  2  2  6.014266314  6.014266356124070  6.0142663139415926 
1  1  1  2.195548355  2.195612369653500  2.1955483554130039 
1  2  2  6.140948992  6.140948969717170  6.1409489918576905 
1  2  5  30.43614539  30.436145317468500  30.436145388713659 
4  1  1  2.734111026  2.73415086499219  2.7341110256122556 
4  2  2  6.542495274  6.54249530312951  6.5424952743905705 
16  1  1  4.399593067  4.399599760664940  4.3995930671655061 
16  2  5  36.99626750  36.996267483327900  36.996267500847930 
Obviously, the procedure is also applicable in the case of complex . Table 3 shows our results for different values of , , and considered in Ref [12]. As it can be seen, the eigenvalues given by Li et al. are confirmed. Nevertheless, in the neighbourhood of each one of those eigenvalues, we have found another one, reported also in Table 3. This result is not surprising, because the values of considered are the approximations found by Oguchi [17] to what he calls “the branch points of the eigenvalues as functions of ”, that is, in our formalism, values of for which a double zero of exists. According to the results of Skorokhodov and Khristoforov [25], there is a double eigenvalue for . This is a much better approximation to the branch point unveiled by Oguchi. As an illustration of what happens in the vicinity of those values of , we present in Figure 4 a modulusphase plot of the function for and , the first of the cases considered in Table 3. The two eigenvalues reported in the table appear as zeros of . As the value of moves from the approximation to the branch point found by Oguchi towards the more precise value given by Skorokhodov and Khristoforov, the two zeros of shown in Fig. 4 approach to each other and eventually collide at a point in the close neighbourhood of the saddle point of suggested by its modulusphase plot. Similar plots of are obtained for the other cases in Table 3.
Ref. [12]  this work  

0  0  
2  
0  0  
4  
0  2  
4  
0  1  
3  
1  1  
3  
1  2  
4  
2  2  
4 
A comment concerning the values of the label of reported in Table 3 is in order. For real or pure imaginary , i.e. for real , the eigenvalues for given are real and can be ordered by increasing value. The label reflects that order. For complex , instead, the values of become complex and such ordination is no more possible. Nevertheless, a label can be assigned to those complex values of , as done by Skorokhodov and Khristoforov. By keeping constant the real part of and continuously decreasing its imaginary part, describes, in the complex plane, a trajectory which intersects the real axis at a certain for . This label can be attached to the whole trajectory described by as varies in the complex plane. The paper by Skorokhodov and Khristoforov contains a very lucid discussion of those trajectories and shows that the branch points correspond to singular values of such that , with .
4 The eigenfunctions
Once the eigenvalues have been calculated, the corresponding eigenfunctions, in the interval , can be obtained immediately by means of the series expansion
(29) 
the coefficients being given by the recurrence relation (10) with . The normalization constant should be adjusted to the normalization scheme preferred. A discussion of the different normalizations used in the literature can be found in the paper by Kirby [14], where the advantage of the unit normalization
(30) 
is made evident. By choosing this normalization, one has
(31) 
where the asterisk indicates complex conjugation. The constant phase in the right hand side of (29) may be taken at will. In the case of real , it is natural to take . For complex , can be chosen in such a way that becomes real at , or at , or at any other point. Needless to say, , according to the even or odd nature of .
Figures 5 and 6 show two examples of the application of the method to the computation of spheroidal angular wave functions in the case of real . The first one is an even prolate angular wave function of parameters , and , and eigenvalue , a case considered in Table III of Ref. [12] (with the Flammer [8] normalization scheme). The second one is an odd oblate angular wave function corresponding to the first of the cases considered in our Table 2. In both figures, the functions have been normalized according to Eq. (30). We have applied our procedure also in some cases of complex . Figure 7 shows the real and imaginary parts and the squared modulus of the wave function in a case considered by Falloon et al. [7], namely the first one in their Table 2. The parameters are , and , and the eigenvalue, in our notation, is . (We give here the eigenvalue with only 15 decimal digits, but our procedure is able to reproduce the 25 decimal digits given in Ref. [7] and to obtain even more.) Finally, Figures 8, 9 and 10 correspond to the second of the cases considered in Table 2 of Ref. [16], of parameters , , , and eigenvalue . (Notice the discrepancy, in the six last digits of both real and imaginary parts, with the value given in Ref. [16].)
5 Conclusions
We have developed a rarely used method which allows to find the eigenvalues and eigenfunctions of the angular spheroidal equation. Instead of having recourse to comparison with other differential equations, the procedure deals with a direct solution, expressed in the form of a convergent series. Requiring it to be regular gives the eigenvalues, which appear as the zeros of a one variable function, . This function and its derivatives and with respect to the variable can be computed, to the desired precision, by summing rapidly convergent series. This fact makes possible the application of the extended Newton method for the determination of the zeros of , i. e., the eigenvalues of the spheroidal equation. Then, the computation of the corresponding eigenfunctions, conveniently normalized, becomes trivial. For the normalization, one benefits from the fact that the squared modulus of the wave function can be integrated algebraically.
The fact that, for given and , the eigenvalues are the zeros of an easily computable function, , makes possible to get an initial approximate location of all of them by a tabulation or a graphical representation of that function. Repeated application of the extended Newton method allows then to calculate the eigenvalues with the desired accuracy.
We have shown the applicability of the method not only in the cases of prolate (real ) and oblate (imaginary ) spheroidal wave equations, but also when is complex. The procedure provides in all cases a very precise determination of the eigenvalues. This has allowed us to resolve the quasiconfluence of pairs of evenodd eigenvalues for large imaginary values of (Table 1) and of pairs of eveneven or oddodd eigenvalues for complex values of in the neighbourhood of “branch points” (Table 3).
The efficiency of the procedure proposed in this paper is subordinate to the capability of computing with sufficient accuracy. The convergence of the series in the right hand side of Eq. (22) is guaranteed in all cases, since for all larger than a certain , and the series may replaced by a sum up to say . Nevertheless, for large values of and/or , the coefficients increase (in modulus) rapidly with before starting to decrease, and the adequate value of may become very large. Even worse, the values of the terms to be summed may cover so many orders of magnitude that the resulting sum is not reliable, unless many significant digits are carried along the computation. This drawback is not outside other procedures. However, the algorithms proposed by Kirby [14] and by Ogburn et al. [16], and procedures collected in Refs. [20], seem to be able to tackle the issue properly. Asymptotic methods [4, 24] have also been forwarded for the mentioned cases of large values of and/or .
Acknowledgements
The work has been supported by Departamento de Ciencia, Tecnología y Universidad del Gobierno de Aragón (Project 226223/1) and Ministerio de Ciencia e Innovación (Project MTM201564166)
Footnotes
 javier@unizar.es
References
 A. A. Abramov and S. V. Kurochkin Highly accurate calculation of angular spheroidal functions. Comput. Math. Math. Phys. 46 (2006), 10–15.
 M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Dover, New York, 1965.
 N. Aquino, E. Castaño and E. LeyKoo, Spheroidal functions revisited: matrix evaluation and generating functions, Rev. Mex. Fis. 48 (2002), 277–282.
 B. E. Barrowes, K. O’Neill, T. M. Grzegorczyk and J. A. Kong, On the asymptotic expansion of the spheroidal wave function and its eigenvalues for complex size parameter, Stud. Appl. Math. 113 (2004), 271–301.
 J. P. Boyd, Prolate spheroidal wavefunctions as an alternative to Chebysev and Legendre polynomials for spectral element and pseudospectral algorithms, J. Comput. Phys. 199 (2004), 688–716.
 J. Caldwell, Computation of eigenvalues of spheroidal harmonics using relaxation, J. Phys. A: Math. Gen. 21 (1988), 3685–3693.
 P. E. Falloon, P. C. Abbott and J. B. Wang, Theory and computation of spheroidal wavefunctions, J. Phys. A: Math. Gen. 36 (2003), 5477–5495.
 C. Flammer, Spheroidal Wave Functions, Stanford University Press, Stanford, CA, 1957.
 L. Gosse, Effective bandlimited extrapolation relying on Slepian series and regularization, Comput. Math. Appl. 60 (2010), 1259–1279.
 J. A. Hogan and J. D. Lakey, Duration and Bandwith Limiting: Prolate Functions, Sampling and Applications, Birkhäuser, Boston, 2011.
 Z. Huang, J. Xiao, and J. P. Boyd, Adaptive radial basis function and Hermite function pseudospectral methods for computing eigenvalues of the prolate spheroidal wave equation, J. Comput. Phys. 281 (2015), 269–284.
 L.W. Li, M.S. Leong, T.S. Yeo, P.S. Kooi, and K.Y. Tan, Computations of spheroidal harmonics with complex arguments: A review with an algorithm, Phys. Rev. E 58 (1998), 6792–6806.
 A. Karoui and T. Moumni, New efficient methods of computing the prolate spheroidal wave functions and their corresponding eigenvalues, Appl. Comput. Harmon. Anal. 24 (2008), 269–289.
 P. Kirby, Calculation of spheroidal wave functions, Comput. Phys. Comm. 175 (2006), 465–472.
 P. Kreuser, Über das Verhalten der Integrale homogener linearer Differenzengleichungen im Unendlichen, Diss. Tübingen, 1914.
 D. X. Ogburn, C. L. Waters, M. D. Sciffer, J. A. Hogan, and P. C. Abbott, A finite difference construction of the spheroidal wave functions, Comput. Phys. Comm. 185 (2014), 244–253.
 T. Oguchi, Eigenvalues of spheroidal wave functions and their branch points for complex values of propagation constants, Radio Sci. 5 (1970), 1207–1214.
 F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark (Eds.), NIST Handbook of Mathematical Functions, Cambridge Univ. Press, Cambridge, 2010. Available at http://dlmf.nist.gov/.
 A. Osipov and V. Rokhlin, On the evaluation of prolate spheroidal wave functions and associated quadrature rules, arXiv:1301.1707.
 A. Osipov, V. Rokhlin, and H. Xiao, Prolate Spheroidal Wave Functions of Order Zero, Springer, New York, 2013.
 O. Perron, Über lineare Differenzengleichungen, Acta math. Stockh. 34 (1910), 109–137.
 O. Perron, Über lineare Differenzengleichungen und eine Anwendung auf lineare Differentialgleichungen mit Polynomkoeffizienten, Math. Zeitschr. 72 (1959), 16–24.

W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes in
FORTRAN
77: The Art of Scientific Computing, Cambridge Univ. Press, New York, 1992.  V. Rokhlin and H. Xiao, Approximate formulae for certain prolate spheroidal wave functions valid for large values of both order and bandlimit, Appl. Comput. Harmon. Anal. 22 (2007), 105–123.
 S. L. Skorokhodov and D. V. Khristoforov, Calculation of the branch points of the eigenfunctions corresponding to wave spheroidal functions. Comput. Math. Math. Phys. 46 (2006), 1132–1146.
 G. Tian, New investigation on the speroidal wave equations, arXiv:1004.1524.
 G. Walter and T. Soleski, A new friendly method of computing prolate spheroidal wave functions and wavelets, Appl. Comput. Harmon. Anal. 19 (2005), 432–443.