## Software for Computing the Spheroidal Wave Functions Using Arbitrary Precision Arithmetic

Ross Adelman, Nail A. Gumerov, and Ramani Duraiswami

## Abstract

The spheroidal wave functions, which are the solutions to the Helmholtz equation in spheroidal coordinates, are notoriously difficult to compute. Because of this, practically no programming language comes equipped with the means to compute them. This makes problems that require their use hard to tackle. We have developed computational software for calculating these special functions. Our software is called spheroidal and includes several novel features, such as: using arbitrary precision arithmetic; adaptively choosing the number of expansion coefficients to compute and use; and using the Wronskian to choose from several different methods for computing the spheroidal radial functions to improve their accuracy. There are two types of spheroidal wave functions: the prolate kind when prolate spheroidal coordinates are used; and the oblate kind when oblate spheroidal coordinate are used. In this paper, we describe both, methods for computing them, and our software. We have made our software freely available on our webpage.

## 1 Introduction

Elementary functions, such as powers, sines, and exponentials, are solutions to differential equations arising in various math, science, and engineering problems. Other elementary functions, such as roots, arcsines, and logarithms, are inverses of these. They are called elementary because the differential equations from which they are derived are usually linear, homogenous, and have constant coefficients. This makes them easy to solve and the expressions for their solutions easy to compute. In many problems, however, the differential equations encountered are not so easy to solve, and their solutions can be extremely complicated. The functions arising from these are called special functions. Examples of such functions are the associated Legendre polynomials and the spherical Bessel functions.

In this paper, we explore the spheroidal wave functions. These special functions are the solutions to the differential equations obtained by applying the method of separation of variables to the Helmholtz equation in spheroidal coordinates. There are two cases: the prolate spheroidal wave functions when prolate spheroidal coordinates are used; and the oblate spheroidal wave functions when oblate spheroidal coordinates are used. Unfortunately, there are no simple expressions for computing them. Instead, they must be written as infinite series expansions in terms of various other special functions. For example, the spheroidal

angle functions can be written in terms of the associated Legendre polynomials, and the spheroidal radial functions can be written in terms of the spherical Bessel and Neumann functions. Depending on the method used, there are three or four different sets of expansion coefficients that need to be computed.

The spheroidal wave functions have applications in many disciplines. Our primary motivation for studying them was for computing the solutions to acoustic scattering problems involving prolate spheroids, oblate spheroids, and disks [1]. However, they are also encountered in signal processing [2].

We have developed computational software for calculating the spheroidal wave functions using C++ and MATLAB. Our software is called spheroidal and includes the following features. First, the software uses GNU MPFR, a library for performing arbitrary precision arithmetic [3]. Using arbitrary precision arithmetic provides greater accuracy in many of the computations, especially for higher wavenumbers and modes. Second, the software allows the user to specify the level of precision to use at every computational step. For example, different levels of precision can be used for computing the different sets of expansion coefficients that need to be computed. Third, the software allows the user to specify in two different ways how many expansion coefficients to compute and use. In the first way, the user specifies the number of expansion coefficients exactly (e.g., compute 200 of this type and 300 of that type). In the second way, the user allows the software to choose the number of expansion coefficients adaptively. All of the expansion coefficients decay exponentially in the long run, so in this method, the user specifies the minimum magnitude that the expansion coefficients should reach (e.g., keep computing expansion coefficients until the next one drops below 10 -200 ). Fourth, the expansion coefficients, as well as other special values, are saved to disk. This way, they can be reused later on without having to recompute them from scratch. Fifth, there are several methods for computing the spheroidal radial functions. The actual value of the Wronskian of these functions is easy to compute, so the combination of methods for computing the spheroidal radial functions is chosen so that the computed Wronskian has the smallest error. We have made our software freely available for download on our webpage.

The spheroidal wave functions have been studied for over six decades. Perhaps the most complete description of the spheroidal wave functions is given by Flammer [4]. As pointed out by [5], some of the expressions in [4] are incorrect. Nevertheless, Flammer's book is an invaluable resource. Actually implementing the spheroidal wave functions in code is very involved. In [6], they were implemented in Fortran, and in [7, 8], they were implemented in C. These implementations used double precision. Due to round-off errors, double precision can lead to large errors, especially for higher frequencies and modes. In [9, 10], they were implemented in Fortran using quad precision and expressions that converge faster and more accurately in some cases to obtain better accuracy over a wide range of frequencies, modes, and argument values. These implementations, as well as our own, are only for real frequnecies and integer modes. Other authors have investigated complex frequnecies and non-integer modes [11, 12], and in these particular references, the respective authors used Mathematica, which can work in arbitrary precision. Some have also looked at numerical techniques, such as finite difference approximations and relaxation methods [13, 14]. Many of these authors have released their code free to use.

Figure 1. The prolate (left) and oblate (right) spheroidal coordinate systems. The three colored surfaces are isosurfaces for η = ± 1 / 2 (red), ξ = 3 2 for the prolate case and / ξ = 1 2 / for the oblate case (green), and φ = 0 (blue).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0074v1-figure-7c7affad-7bd3-4457-a5fc-4062fd574051.jpg)

## 2 Spheroidal Coordinates

The prolate spheroidal coordinate system, shown in Figure 1, is related to the Cartesian coordinate system by [4]

$$x = a \left ( 1 - \eta ^ { 2 } \right ) ^ { 1 / 2 } \left ( \xi ^ { 2 } - 1 \right ) ^ { 1 / 2 } \cos \left ( \phi \right ), \quad y = a \left ( 1 - \eta ^ { 2 } \right ) ^ { 1 / 2 } \left ( \xi ^ { 2 } - 1 \right ) ^ { 1 / 2 } \sin \left ( \phi \right ), \quad z = a \eta \xi, \ \ ( 1 )$$

where 2 a is the interfocal distance. The Helmholtz equation, ∇ 2 V + k V 2 = 0, where k is the wavenumber, can be written in prolate spheroidal coordinates as

$$\left ( \frac { \partial } { \partial \eta } \left ( 1 - \eta ^ { 2 } \right ) \frac { \partial } { \partial \eta } + \frac { \partial } { \partial \xi } \left ( \xi ^ { 2 } - 1 \right ) \frac { \partial } { \partial \xi } + \frac { \xi ^ { 2 } - \eta ^ { 2 } } { \left ( 1 - \eta ^ { 2 } \right ) \left ( \xi ^ { 2 } - 1 \right ) } \frac { \partial ^ { 2 } } { \partial \phi ^ { 2 } } + c ^ { 2 } \left ( \xi ^ { 2 } - \eta ^ { 2 } \right ) \right ) V = 0, \quad ( 2 )$$

where c = ka . Applying the method of separation of variables yields three uncoupled ordinary differential equations, one for each coordinate:

$$\frac { \partial } { \partial \eta } \left ( \left ( 1 - \eta ^ { 2 } \right ) \frac { \partial } { \partial \eta } S _ { m n } \left ( c, \eta \right ) \right ) + \left ( \lambda _ { m n } - c ^ { 2 } \eta ^ { 2 } - \frac { m ^ { 2 } } { 1 - \eta ^ { 2 } } \right ) S _ { m n } \left ( c, \eta \right ) = 0,$$

$$\frac { \partial } { \partial \xi } \left ( \left ( \xi ^ { 2 } - 1 \right ) \frac { \partial } { \partial \xi } R _ { m n } \left ( c, \xi \right ) \right ) - \left ( \lambda _ { m n } - c ^ { 2 } \xi ^ { 2 } + \frac { m ^ { 2 } } { \xi ^ { 2 } - 1 } \right ) R _ { m n } \left ( c, \xi \right ) = 0,$$

$$\frac { \partial ^ { 2 } } { \partial \phi ^ { 2 } } \Phi _ { m } \left ( \phi \right ) + m ^ { 2 } \Phi _ { m } \left ( \phi \right ) = 0,$$

where m = 0 1 , , . . . and n = m,m +1 , . . . . While Eq. (5) is easily solved, Eqs. (3) and (4) are much more complicated. The solutions to Eq. (3) are called the prolate spheroidal angle functions, and the solutions to Eq. (4) are called the prolate spheroidal radial functions. Collectively, they are called the prolate spheroidal wave functions. Any solution to Eq. (2) can be written as

$$V = \sum _ { m = 0 } ^ { \infty } \sum _ { n = m } ^ { \infty } S _ { m n } \left ( c, \eta \right ) \left ( A _ { m n } R _ { m n } ^ { \left ( 1 \right ) } \left ( c, \xi \right ) + B _ { m n } R _ { m n } ^ { \left ( 3 \right ) } \left ( c, \xi \right ) \right ) \cos \left ( m \phi \right ),$$

Figure 2. Characteristic and other special values for the prolate spheroidal wave functions for c = 10, m = 0 1 , , . . . , 29, and n = m,m +1 , . . . , m +29.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0074v1-figure-71734666-602e-497f-ad77-d1152a8f97f5.jpg)

where the expansion coefficients, A mn and B mn , depend on the problem being solved.

The expressions arising in the oblate case are very similar to (and sometimes exactly the same as) those arising in the prolate case. In many cases, simply letting c, ξ →-ic, iξ provides a transformation from the prolate case to the oblate case [4]. Indeed, the preceeding paragraphs and equations for the prolate case can be transformed into those for the oblate case by using this transformation. The oblate spheroidal coordinate system is shown in Figure 1.

Figure 3. The prolate spheroidal wave functions and their derivatives for c = 10, m = 10, and n = 10 11 , , . . . , 39.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0074v1-figure-c9c08a83-f342-46cb-be2d-6cb034b64f04.jpg)

## 3 Spheroidal Wave Functions

## 3.1 Prolate Spheroidal Wave Functions

## 3.1.1 Angle Functions

The prolate spheroidal angle functions of the first and second kinds can be written in terms of the associated Legendre polynomials of the first and second kinds, respectively:

$$S _ { m n } ^ { \left ( 1 \right ) } \left ( c, \eta \right ) = \sum _ { r = 0, 1 } ^ { \infty } d _ { r } ^ { m n } \left ( c \right ) P _ { m + r } ^ { m } \left ( \eta \right ),$$

$$S _ { m n } ^ { \left ( 2 \right ) } \left ( c, \eta \right ) = \sum _ { r = - \infty } ^ { \infty } d _ { r } ^ { m n } \left ( c \right ) Q _ { m + r } ^ { m } \left ( \eta \right ),$$

where the prime over the sum means that only the even terms are included when n -m = even and only the odd terms are included when n -m = odd. The angle functions of the first kind are orthogonal over [ -1 1]: ,

$$\int _ { - 1 } ^ { 1 } S _ { m n } ^ { \left ( 1 \right ) } \left ( c, \eta \right ) S _ { m n ^ { \prime } } ^ { \left ( 1 \right ) } \left ( c, \eta \right ) d \eta = \delta _ { n n ^ { \prime } } N _ { m n } \left ( c \right ),$$

where glyph[negationslash]

$$\text{where} \quad \delta _ { n n ^ { \prime } } = \begin{cases} \ 1, & n = n ^ { \prime } \\ 0, & n \neq n ^ { \prime } \end{cases} \, \quad N _ { m n } \left ( c \right ) = 2 \sum _ { r = 0, 1 } ^ { \infty } d _ { r } ^ { m n } \left ( c \right ) ^ { 2 } \frac { ( 2 m + r )! } { ( 2 m + 2 r + 1 ) \, r! }. \quad \ \ ( 1 0 )$$

The angle functions of the first kind can also be written as a power series. This can be done with the help of the hypergeometric function. When n -m = even,

$$S _ { m n } ^ { \left ( 1 \right ) } \left ( c, \eta \right ) = \sum _ { r = 0 } ^ { \infty } d _ { 2 r } ^ { m n } \left ( c \right ) P _ { m + 2 r } ^ { m } \left ( \eta \right ).$$

The associated Legendre polynomials of the first kind can be written as

$$P _ { n } ^ { m } \left ( x \right ) = \frac { \left ( - 1 \right ) ^ { m } \left ( n + m \right )! } { 2 ^ { m } m! \left ( n - m \right )! } \left ( 1 - x ^ { 2 } \right ) ^ { m / 2 } F \left ( m - n, n + m + 1 ; m + 1 ; \frac { 1 - x } { 2 } \right ),$$

where

$$F \left ( \alpha, \beta ; \gamma ; x \right ) = \sum _ { k = 0 } ^ { \infty } \frac { \left ( \alpha \right ) _ { k } \left ( \beta \right ) _ { k } } { k! \left ( \gamma \right ) _ { k } } x ^ { k }$$

is the hypergeometric function and

$$( a ) _ { 0 } = 1, \quad ( a ) _ { k } = a \, ( a + 1 ) \dots ( a + k - 1 )$$

is the rising Pochhammer symbol. Using this, the associated Legendre polynomials of the first kind in Eq. (11) can be written as

$$P _ { m + 2 r } ^ { m } \left ( x \right ) = \frac { \left ( - 1 \right ) ^ { m } \left ( 2 m + 2 r \right )! } { 2 ^ { m } m! \left ( 2 r \right )! } \left ( 1 - x ^ { 2 } \right ) ^ { m / 2 } F \left ( - 2 r, 2 m + 2 r + 1 ; m + 1 ; \frac { 1 - x } { 2 } \right ).$$

Using the following identity for the hypergeometric function, we can rearrange the previous expression slightly:

$$F \left ( \alpha, \beta ; \frac { \alpha + \beta + 1 } { 2 } ; x \right ) = F \left ( \frac { \alpha } { 2 }, \frac { \beta } { 2 } ; \frac { \alpha + \beta + 1 } { 2 } ; 4 x \left ( 1 - x \right ) \right ),$$

$$P _ { m + 2 r } ^ { m } \left ( x \right ) = \frac { \left ( - 1 \right ) ^ { m } \left ( 2 m + 2 r \right )! } { 2 ^ { m } m! \left ( 2 r \right )! } \left ( 1 - x ^ { 2 } \right ) ^ { m / 2 } F \left ( - r, m + r + \frac { 1 } { 2 } ; m + 1 ; 1 - x ^ { 2 } \right ), \quad \left ( 1 7 \right )$$

$$P _ { m + 2 r } ^ { m } \left ( x \right ) = \frac { \left ( - 1 \right ) ^ { m } \left ( 2 m + 2 r \right )! } { 2 ^ { m } m! \left ( 2 r \right )! } \left ( 1 - x ^ { 2 } \right ) ^ { m / 2 } \sum _ { k = 0 } ^ { \infty } \frac { \left ( - r \right ) _ { k } \left ( m + r + \frac { 1 } { 2 } \right ) _ { k } } { k! \left ( m + 1 \right ) _ { k } } \left ( 1 - x ^ { 2 } \right ) ^ { k }.$$

Plugging this into Eq. (11), we have

$$S _ { m n } ^ { ( 1 ) } \left ( c, \eta \right ) & = \sum _ { r = 0 } ^ { \infty } d _ { 2 r } ^ { m n } \left ( c \right ) \frac { \left ( - 1 \right ) ^ { m } \left ( 2 m + 2 r \right )! } { 2 ^ { m } m! \left ( 2 r \right )! } \left ( 1 - \eta ^ { 2 } \right ) ^ { m / 2 } \sum _ { k = 0 } ^ { \infty } \frac { \left ( - r \right ) _ { k } \left ( m + r + \frac { 1 } { 2 } \right ) _ { k } } { k! \left ( m + 1 \right ) _ { k } } \left ( 1 - \eta ^ { 2 } \right ) ^ { k }.$$

Rearranging,

$$\text{gmg}, & S _ { m n } ^ { ( 1 ) } \left ( c, \eta \right ) = \left ( - 1 \right ) ^ { m } \left ( 1 - \eta ^ { 2 } \right ) ^ { m / 2 } \sum _ { k = 0 } ^ { \infty } c _ { 2 k } ^ { m n } \left ( c \right ) \left ( 1 - \eta ^ { 2 } \right ) ^ { k }, & ( 2 0 )$$

where

$$\text{where} \\ c _ { 2 k } ^ { m n } \left ( c \right ) = \frac { 1 } { 2 ^ { m } \left ( m + k \right )! k! } \sum _ { r = 2 k } ^ { \infty } d _ { r } ^ { m n } \left ( c \right ) \frac { ( 2 m +$$

When n -m = odd,

$$S _ { m n } ^ { ( 1 ) } \left ( c, \eta \right ) = \sum _ { r = 0 } ^ { \infty } d _ { 2 r + 1 } ^ { m n } \left ( c \right ) P _ { m + 2 r + 1 } ^ { m } \left ( \eta \right ). \quad \quad \quad$$

Recall the following identity for the associated Legendre polynomials of the first kind:

$$\left ( n - m + 1 \right ) P _ { n + 1 } ^ { m } \left ( x \right ) = \left ( n + 1 \right ) x P _ { n } ^ { m } \left ( x \right ) - \left ( 1 - x ^ { 2 } \right ) P _ { n } ^ { m \prime } \left ( x \right )$$

Setting n = m +2 , plugging in Eq. (18), and rearranging, r

$$P _ { m + 2 r + 1 } ^ { m } \left ( x \right ) = \frac { \left ( - 1 \right ) ^ { m } \left ( 2 m + 2 r + 1 \right )! } { 2 ^ { m } m! \left ( 2 r + 1 \right )! } x \left ( 1 - x ^ { 2 } \right ) ^ {$$

Finally, plugging this into Eq. (22),

$$\dots \dots _ { J }, r \dots \circ \dots \circ \dots \circ \dots \circ \dots \circ \dots \circ \dots \circ \dots \circ \dots \circ \dots \circ \dots \circ \dots \circ \dots \circ \dots \circ \dots \circ \$$

Rearranging,

$$\text{augg}, & S _ { m n } ^ { ( 1 ) } \left ( c, \eta \right ) = \left ( - 1 \right ) ^ { m } \eta \left ( 1 - \eta ^ { 2 } \right ) ^ { m / 2 } \sum _ { k = 0 } ^ { \infty } c _ { 2 k } ^ { m n } \left ( c \right ) \left ( 1 - \eta ^ { 2 } \right ) ^ { k },$$

where

$$c _ { 2 k } ^ { m n } \left ( c \right ) = \frac { 1 } { 2 ^ { m } \left ( m + k \right )! k! } \sum _ { r = 2 k + 1 } ^ { \infty } { ^ { \prime } } d _ { r } ^ { m n } \left ( c \right ) \frac { ( 2 m + r )! } { r! } \left ( - \frac { r - 1 } { 2 } \right ) _ { k } \left ( m + \frac { r } { 2 } + 1 \right ) _ { k }.$$

## 3.1.2 Radial Functions

The prolate spheroidal radial functions of the first and second kinds can be written in terms of the spherical Bessel and Neumann functions, respectively:

$$R _ { m n } ^ { ( 1 ) } \left ( c, \xi \right ) = F _ { m n } \left ( c \right ) ^ { - 1 } \left ( 1 - \frac { 1 } { \xi ^ { 2 } } \right ) ^ { m / 2 } \sum _ { r = 0, 1 } ^ { \infty } \, ^ { \prime } \, ( - 1 ) ^ { ( r - ( n - m ) ) / 2 } \, d _ { r } ^ { m n } \left ( c \right ) \frac { ( 2 m + r )! } { r! } j _ { m + r } \left ( c \xi \right ), \ \ ( 2 8 )$$

$$R _ { m n } ^ { ( 2 ) } \left ( c, \xi \right ) = F _ { m n } \left ( c \right ) ^ { - 1 } \left ( 1 - \frac { 1 } { \xi ^ { 2 } } \right ) ^ { m / 2 } \sum _ { r = 0, 1 } ^ { \infty } \, ^ { \prime } \, ( - 1 ) ^ { ( r - ( n - m ) ) / 2 } \, d _ { r } ^ { m n } \left ( c \right ) \frac { ( 2 m + r )! } { r! } y _ { m + r } \left ( c \xi \right ), \ \ ( 2 9 )$$

where

$$F _ { m n } \left ( c \right ) = \sum _ { r = 0, 1 } ^ { \infty } { ^ { \prime } } d _ { r } ^ { m n } \left ( c \right ) \frac { ( 2 m + r )! } { r! }.$$

The radial functions of the third and fourth kinds are linear combinations of those of the first and second kinds:

$$R _ { m n } ^ { \left ( 3 \right ) } \left ( c, \xi \right ) = R _ { m n } ^ { \left ( 1 \right ) } \left ( c, \xi \right ) + i R _ { m n } ^ { \left ( 2 \right ) } \left ( c, \xi \right ),$$

$$R _ { m n } ^ { \left ( 4 \right ) } \left ( c, \xi \right ) = R _ { m n } ^ { \left ( 1 \right ) } \left ( c, \xi \right ) - i R _ { m n } ^ { \left ( 2 \right ) } \left ( c, \xi \right ).$$

The Wronskian of the radial functions of the first and second kinds is given by

$$W _ { m n } \left ( c, \xi \right ) = R _ { m n } ^ { \left ( 1 \right ) } \left ( c, \xi \right ) \frac { \partial } { \partial \xi } R _ { m n } ^ { \left ( 2 \right ) } \left ( c, \xi \right ) - \frac { \partial } { \partial \xi } R _ { m n } ^ { \left ( 1 \right ) } \left ( c, \xi \right ) R _ { m n } ^ { \left ( 2 \right ) } \left ( c, \xi \right ) = \frac { 1 } { c \left ( \xi ^ { 2 } - 1 \right ) } \quad \left ( 3 3 \right )$$

and is useful for validating computed values of these functions.

The radial functions are related to the angle functions by

$$S _ { m n } ^ { \left ( 1 \right ) } \left ( c, z \right ) = k _ { m n } ^ { \left ( 1 \right ) } \left ( c \right ) R _ { m n } ^ { \left ( 1 \right ) } \left ( c, z \right ),$$

$$S _ { m n } ^ { \left ( 2 \right ) } \left ( c, z \right ) = k _ { m n } ^ { \left ( 2 \right ) } \left ( c \right ) R _ { m n } ^ { \left ( 2 \right ) } \left ( c, z \right ),$$

where k (1) mn ( c ) is given by

$$k _ { m n } ^ { \left ( 1 \right ) } \left ( c \right ) = \frac { \left ( 2 m + 1 \right ) \left ( m + n \right )! F _ { m n } \left ( c \right ) } { 2 ^ { m + n } d _ { 0 } ^ { m n } \left ( c \right ) c ^ { m } m! \left ( \frac { n - m } { 2 } \right )! \left ( \frac { m + n } { 2 } \right )! }, \quad n - m = \text{even},$$

$$k _ { m n } ^ { ( 1 ) } \left ( c \right ) = \frac { \left ( 2 m + 3 \right ) \left ( m + n + 1 \right )! F _ { m n } \left ( c \right ) } { 2 ^ { m + n } d _ { 1 } ^ { m n } \left ( c \right ) c ^ { m + 1 } m! } \left ( \frac { n - m - 1 } { 2 } \right )! \left ( \frac { m + n + 1 } { 2 } \right )! }, \quad n - m = \text{odd}, \quad \left ( 3 7 \right )$$

and k (2) mn ( c ) is given by

$$k _ { m n } ^ { ( 2 ) } \left ( c \right ) = \frac { 2 ^ { n - m } \left ( 2 m \right )! \left ( \frac { n - m } { 2 } \right )! \left ( \frac { m + n } { 2 } \right )! d _ { - 2 m } ^ { m n } \left ( c \right ) F _ { m n } \left ( c \right ) } { \left ( 2 m - 1 \right ) m! \left ( m + n \right )! c ^ { m - 1 } }, \quad n - m = \text{even},$$

$$k _ { m n } ^ { ( 2 ) } \left ( c \right ) = - \frac { 2 ^ { n - m } \left ( 2 m \right )! \left ( \frac { n - m - 1 } { 2 } \right )! \left ( \frac { m + n + 1 } { 2 } \right )! d _ { - 2 m + 1 } ^ { m n } \left ( c \right ) F _ { m n } \left ( c \right ) } { \left ( 2 m - 3 \right ) \left ( 2 m - 1 \right ) m! \left ( m + n + 1 \right )! c ^ { m - 2 } }, \quad n - m = \text{odd}. \\ \text{The envsision for the radial functions of the second kind using} \text{ the spherical Neumann}$$

The expression for the radial functions of the second kind using the spherical Neumann functions converges very slowly for values of ξ near 1 and is, therefore, inaccurate in these cases. While the expression for the radial functions of the first kind is accurate for all values of ξ , having a second method can be used as a check on the first method. Thus, these relationships can be used to construct secondary methods for computing these functions. For the radial functions of the first kind,

$$R _ { m n } ^ { \left ( 1 \right ) } \left ( c, \xi \right ) = k _ { m n } ^ { \left ( 1 \right ) } \left ( c \right ) ^ { - 1 } S _ { m n } ^ { \left ( 1 \right ) } \left ( c, \xi \right ),$$

$$R _ { m n } ^ { \left ( 1 \right ) } \left ( c, \xi \right ) = k _ { m n } ^ { \left ( 1 \right ) } \left ( c \right ) ^ { - 1 } \sum _ { r = 0, 1 } ^ { \infty } \| ^ { \prime } d _ { r } ^ { m n } \left ( c \right ) P _ { m + r } ^ { m } \left ( \xi \right ).$$

Similar to the angle functions of the first kind, this expression can be written as a power series with the help of the hypergeometric function. Because the argument is ξ ≥ 1 as opposed to | η | ≤ 1, the relationship between the associated Legendre polynomials of the first kind and the hypergeometric function is slightly different. In particular, there is no factor of ( -1) m :

$$P _ { n } ^ { m } \left ( x \right ) = \frac { \left ( n + m \right )! } { 2 ^ { m } m! \left ( n - m \right )! } \left ( x ^ { 2 } - 1 \right ) ^ { m / 2 } F \left ( m - n, n + m + 1 ; m + 1 ; \frac { 1 - x } { 2 } \right ).$$

Following a procedure similar to the one followed for the angle functions of the first kind, we have

$$R _ { m n } ^ { ( 1 ) } \left ( c, \xi \right ) = k _ { m n } ^ { ( 1 ) } \left ( c \right ) ^ { - 1 } \left ( \xi ^ { 2 } - 1 \right ) ^ { m / 2 } \sum _ { k = 0 } ^ { \infty } \left ( - 1 \right ) ^ { k } c _ { 2 k } ^ { m n } \left ( c \right ) \left ( \xi ^ { 2 } - 1 \right ) ^ { k }, \quad n - m = \text{even}, \quad \text{(43)}$$

$$R _ { m n } ^ { ( 1 ) } \left ( c, \xi \right ) = k _ { m n } ^ { ( 1 ) } \left ( c \right ) ^ { - 1 } \xi \left ( \xi ^ { 2 } - 1 \right ) ^ { m / 2 } \sum _ { k = 0 } ^ { \infty } \left ( - 1 \right ) ^ { k } c _ { 2 k } ^ { m n } \left ( c \right ) \left ( \xi ^ { 2 } - 1 \right ) ^ { k }, \quad n - m = \text{odd}, \quad \text{(44)}$$

where c mn 2 k ( c ) is the same as before. For the radial functions of the second kind,

$$R _ { m n } ^ { \left ( 2 \right ) } \left ( c, \xi \right ) = k _ { m n } ^ { \left ( 2 \right ) } \left ( c \right ) ^ { - 1 } S _ { m n } ^ { \left ( 2 \right ) } \left ( c, \xi \right ),$$

$$R _ { m n } ^ { \left ( 2 \right ) } \left ( c, \xi \right ) = k _ { m n } ^ { \left ( 2 \right ) } \left ( c \right ) ^ { - 1 } \sum _ { r = - \infty } ^ { \infty } d _ { r } ^ { m n } \left ( c \right ) Q _ { m + r } ^ { m } \left ( \xi \right ).$$

## 3.1.3 Calculating the Characteristic Value and Expansion Coefficients

All of the expressions introduced in the previous sections require, either directly or indirectly, the characteristic value, λ mn ( c ), and expansion coefficients, d mn r ( c ). Below, we derive expressions for computing them.

There are several methods for computing the characterstic value, but here, we use a combination of two: method (1) involves solving for the eigenvalues of a tridiagonal matrix [15]; and method (2) involves solving for the roots of a transcendental equation [4]. To begin, method (1) is used to compute an approximate value for the characteristic value. Then, method (2) is used to compute a more accurate value for the characteristic value using the approximate value computed by method (1) as a starting point. This procedure is similar to the one used in [12]. In our software, we use double precision for method (1) and arbitrary precision for method (2).

Both methods rely on the following recurrence relation, which can be obtained by plugging Eq. (7) into Eq. (3):

$$\alpha _ { r } d _ { r + 2 } ^ { m n } \left ( c \right ) + \left ( \beta _ { r } - \lambda _ { m n } \left ( c \right ) \right ) d _ { r } ^ { m n } \left ( c \right ) + \gamma _ { r } d _ { r - 2 } ^ { m n } \left ( c \right ) = 0,$$

where

$$\alpha _ { r } = \frac { \left ( 2 m + r + 2 \right ) \left ( 2 m + r + 1 \right ) } { \left ( 2 m + 2 r + 5 \right ) \left ( 2 m + 2 r + 3 \right ) } c ^ { 2 },$$

$$\beta _ { r } = \left ( m + r \right ) \left ( m + r + 1 \right ) + \frac { 2 \left ( m + r \right ) \left ( m + r + 1 \right ) - 2 m ^ { 2 } - 1 } { \left ( 2 m + 2 r - 1 \right ) \left ( 2 m + 2 r + 3 \right ) } c ^ { 2 },$$

$$\gamma _ { r } = \frac { r \left ( r - 1 \right ) } { \left ( 2 m + 2 r - 3 \right ) \left ( 2 m + 2 r - 1 \right ) } c ^ { 2 }.$$

For method (1), this recurrence relation is rearranged slightly:

$$\alpha _ { r } d _ { r + 2 } ^ { m n } \left ( c \right ) + \beta _ { r } d _ { r } ^ { m n } \left ( c \right ) + \gamma _ { r } d _ { r - 2 } ^ { m n } \left ( c \right ) = \lambda _ { m n } \left ( c \right ) d _ { r } ^ { m n } \left ( c \right ).$$

In matrix form,

$$\left [ \begin{array} { c c c } \beta _ { 0 } & \alpha _ { 0 } \\ \gamma _ { 2 } & \beta _ { 2 } & \alpha _ { 2 } \\ & \gamma _ { 4 } & \beta _ { 4 } & \alpha _ { 4 } \\ & & \ddots \end{array} \right ] \left [ \begin{array} { c } d _ { 0 } ^ { m n } \left ( c \right ) \\ d _ { 2 } ^ { m n } \left ( c \right ) \\ d _ { 4 } ^ { m n } \left ( c \right ) \\ \vdots \end{array} \right ] = \lambda _ { m n } \left ( c \right ) \left [ \begin{array} { c } d _ { 0 } ^ { m n } \left ( c \right ) \\ d _ { 2 } ^ { m n } \left ( c \right ) \\ d _ { 4 } ^ { m n } \left ( c \right ) \\ \vdots \end{array} \right ], \quad n - m = \text{even,}$$

$$\left [ \begin{array} { c c c } \beta _ { 1 } & \alpha _ { 1 } & & \\ \gamma _ { 3 } & \beta _ { 3 } & \alpha _ { 3 } & \\ & \gamma _ { 5 } & \beta _ { 5 } & \alpha _ { 5 } \\ & & & \ddots \end{array} \right ] \left [ \begin{array} { c } d _ { 1 } ^ { m n } \left ( c \right ) \\ d _ { 3 } ^ { m n } \left ( c \right ) \\ d _ { 5 } ^ { m n } \left ( c \right ) \\ \vdots \end{array} \right ] = \lambda _ { m n } \left ( c \right ) \left [ \begin{array} { c } d _ { 1 } ^ { m n } \left ( c \right ) \\ d _ { 3 } ^ { m n } \left ( c \right ) \\ d _ { 5 } ^ { m n } \left ( c \right ) \\ \vdots \end{array} \right ], \quad n - m = \text{odd.} \quad \quad \quad$$

When n -m = even, the eigenvalues are λ mn ( c ) for n = m,m + 2 , m + 4 , . . . , and when n -m = odd, the eigenvalues are λ mn ( c ) for n = m +1 , m +3 , m +5 , . . . . Thus, we can compute λ mn ( c ) by plugging these tridiagonal matrices into an eigenvalue solver. In our software, we use the eig function in MATLAB.

In method (2), the recurrence relation in Eq. (47) is divided through by d mn r ( c ), which yields

$$\alpha _ { r } \frac { d _ { r + 2 } ^ { m n } \left ( c \right ) } { d _ { r } ^ { m n } \left ( c \right ) } + \beta _ { r } - \lambda _ { m n } \left ( c \right ) + \gamma _ { r } \frac { d _ { r - 2 } ^ { m n } \left ( c \right ) } { d _ { r } ^ { m n } \left ( c \right ) } = 0.$$

Setting

$$N _ { r } ^ { m } = - \alpha _ { r - 2 } \frac { d _ { r } ^ { m n } \left ( c \right ) } { d _ { r - 2 } ^ { m n } \left ( c \right ) }, \quad \gamma _ { r } ^ { m } = \beta _ { r }, \quad \beta _ { r } ^ { m } = \gamma _ { r } \alpha _ { r - 2 }$$

allows us to write Eq. (54) as

$$- \, N _ { r + 2 } ^ { m } + \gamma _ { r } ^ { m } - \lambda _ { m n } \left ( c \right ) - \frac { \beta _ { r } ^ { m } } { N _ { r } ^ { m } } = 0.$$

Rearranging one way leads to a continued fraction in decreasing r :

$$N _ { r } ^ { m } = \gamma _ { r - 2 } ^ { m } - \lambda _ { m n } \left ( c \right ) - \frac { \beta _ { r - 2 } ^ { m } } { \gamma _ { r - 4 } ^ { m } - \lambda _ { m n } \left ( c \right ) - } \frac { \beta _ { r - 4 } ^ { m } } { \gamma _ { r - 6 } ^ { m } - \lambda _ { m n } \left ( c \right ) - } \cdots.$$

Rearranging the other way leads to a continued fraction in increasing r :

$$N _ { r } ^ { m } = \frac { \beta _ { r } ^ { m } } { \gamma _ { r } ^ { m } - \lambda _ { m n } \left ( c \right ) - } \frac { \beta _ { r + 2 } ^ { m } } { \gamma _ { r + 2 } ^ { m } - \lambda _ { m n } \left ( c \right ) - } \cdots.$$

The two expression for N m r should be equal to each other. Setting r = n -m +2 in Eqs. (57) and (58),

$$U _ { 1 } \left ( \lambda _ { m n } \left ( c \right ) \right ) = \gamma _ { n - m } ^ { m } - \lambda _ { m n } \left ( c \right ) - \frac { \beta _ { n - m } ^ { m } } { \gamma _ { n - m - 2 } ^ { m } - \lambda _ { m n } \left ( c \right ) - } \frac { \beta _ { n - m - 2 } ^ { m } } { \gamma _ { n - m - 4 } ^ { m } - \lambda _ { m n } \left ( c \right ) - } \cdots, \quad \text{(59)}$$

$$U _ { 2 } \left ( \lambda _ { m n } \left ( c \right ) \right ) = - \frac { \beta _ { n - m + 2 } ^ { m } } { \gamma _ { n - m + 2 } ^ { m } - \lambda _ { m n } \left ( c \right ) - } \frac { \beta _ { n - m + 4 } ^ { m } } { \gamma _ { n - m + 4 } ^ { m } - \lambda _ { m n } \left ( c \right ) - } \cdots.$$

Adding these together yields a transcendental equation in λ mn ( c ):

$$U \left ( \lambda _ { m n } \left ( c \right ) \right ) = U _ { 1 } \left ( \lambda _ { m n } \left ( c \right ) \right ) + U _ { 2 } \left ( \lambda _ { m n } \left ( c \right ) \right ) = 0.$$

We can compute λ mn ( c ) by solving this transcendental equation. In our software, we use the secant method.

Once λ mn ( c ) is known, Eqs. (55) and (58) can be used to compute d mn r ( c ). The expansion coefficients are unique up to a constant factor, though, so the following normalization scheme is used. When n -m = even,

$$S _ { m n } ^ { \left ( 1 \right ) } \left ( c, 0 \right ) = P _ { n } ^ { m } \left ( 0 \right ),$$

$$\sum _ { r = 0 } ^ { \infty } ^ { \prime } d _ { r } ^ { m n } \left ( c \right ) \frac { \left ( - 1 \right ) ^ { r / 2 } \left ( 2 m + r \right )! } { 2 ^ { r } \left ( \frac { 2 m + r } { 2 } \right )! \left ( \frac { r } { 2 } \right )! } = \frac { \left ( - 1 \right ) ^ { ( n - m ) / 2 } \left ( n + m \right )! } { 2 ^ { n - m } \left ( \frac { n + m } { 2 } \right )! \left ( \frac { n - m } { 2 } \right )! }.$$

When n -m = odd,

$$S _ { m n } ^ { \left ( 1 \right ) \prime } \left ( c, 0 \right ) = P _ { n } ^ { m \prime } \left ( 0 \right ),$$

$$\sum _ { r = 1 } ^ { \infty } ^ { \prime } d _ { r } ^ { m n } \left ( c \right ) \frac { ( - 1 ) ^ { ( r - 1 ) / 2 } \left ( 2 m + r + 1 \right )! } { 2 ^ { r } \left ( \frac { 2 m + r + 1 } { 2 } \right )! \left ( \frac { r - 1 } { 2 } \right )! } = \frac { n \cdots \, \nu \, \nu \, \nu } { ( - 1 ) ^ { ( n - m - 1 ) / 2 } \left ( n + m + 1 \right )! } \cdots \quad \cdots \quad \cdots \quad \cdots \quad \cdots } { 2 ^ { n - m } \left ( \frac { n + m + 1 } { 2 } \right )! \left ( \frac { n - m - 1 } { 2 } \right )! }.$$

To use Eq. (46), d mn r ( c ) must be computed for negative r as well. To begin, Eq. (47) is rewritten as

$$A _ { r + 2 } ^ { m } d _ { r + 2 } ^ { m n } \left ( c \right ) + B _ { r } ^ { m } d _ { r } ^ { m n } \left ( c \right ) + C _ { r - 2 } ^ { m } d _ { r - 2 } ^ { m n } \left ( c \right ) = 0,$$

where

Rearranging,

$$\frac { d _ { r } ^ { m n } \left ( c \right ) } { d _ { r + 2 } ^ { m n } \left ( c \right ) } & = - \frac { A _ { r + 2 } ^ { m } } { B _ { r } ^ { m } + C _ { r - 2 } ^ { m } \frac { d _ { r - 2 } ^ { m n } \left ( c \right ) } { d _ { r } ^ { m n } \left ( c \right ) } }, \\, \,, \quad. \quad, \,, \, e \quad. \quad, \quad.$$

which can be expanded as a continued fraction in decreasing r :

$$\frac { d _ { r } ^ { m n } \left ( c \right ) } { d _ { r + 2 } ^ { m n } \left ( c \right ) } = - \frac { A _ { r + 2 } ^ { m } } { B _ { r } ^ { m } - } \frac { C _ { r - 2 } ^ { m } A _ { r } ^ { m } } { B _ { r - 2 } ^ { m } - } \frac { C _ { r - 4 } ^ { m } A _ { r - 2 } ^ { m } } { B _ { r - 4 } ^ { m } - } \cdots.$$

$$A _ { r } ^ { m } = \alpha _ { r - 2 }, \ \ B _ { r } ^ { m } = \beta _ { r } - \lambda _ { m n } \left ( c \right ), \ \ C _ { r } ^ { m } = \gamma _ { r + 2 }.$$

Because A m r = 0 when r = -2 m or r = -2 m +1, this continued fraction ends:

$$\frac { d _ { r } ^ { m n } \left ( c \right ) } { d _ { r + 2 } ^ { m n } \left ( c \right ) } = - \frac { A _ { r + 2 } ^ { m } } { B _ { r } ^ { m } - } \frac { C _ { r - 2 } ^ { m } A _ { r } ^ { m } } { B _ { r - 2 } ^ { m } - } \frac { C _ { r - 4 } ^ { m } A _ { r - 2 } ^ { m } } { B _ { r - 4 } ^ { m } - } \cdots \frac { A _ { - 2 m + 2 } ^ { m } } { B _ { - 2 m } ^ { m } + C _ { - 2 m - 2 } ^ { m } }$$

when n -m = even, and

$$\frac { d _ { r } ^ { m n } \left ( c \right ) } { d _ { r + 2 } ^ { m n } \left ( c \right ) } = - \frac { A _ { r + 2 } ^ { m } } { B _ { r } ^ { m } - } \frac { C _ { r - 2 } ^ { m } A _ { r } ^ { m } } { B _ { r - 2 } ^ { m } - } \frac { C _ { r - 4 } ^ { m } A _ { r - 2 } ^ { m } } { B _ { r - 4 } ^ { m } - } \cdots \frac { A _ { - 2 m + 3 } ^ { m } } { B _ { - 2 m + 1 } ^ { m } + C _ { - 2 m - 1 } ^ { m } }$$

when n -m = odd. This also means that d mn r ( c ) → 0 when r ≤ -2 m -2 for n -m = even and r ≤ -2 m -1 for n -m = odd. However, Q m m r + ( ξ ) → ∞ in these cases, and d mn r ( c ) Q m m r + ( ξ ) &lt; ∞ :

$$d _ { r } ^ { m n } \left ( c \right ) Q _ { m + r } ^ { m } \left ( \xi \right ) = d _ { r | \epsilon } ^ { m n } \left ( c \right ) P _ { - r - m - 1 } ^ { m } \left ( \xi \right ).$$

The sum in Eq. (46) is, therefore, separated into two pieces:

$$R _ { m n } ^ { ( 2 ) } \left ( c, \xi \right ) & = I _ { m n } ^ { ( 2 ) } \left ( \sum _ { r = - \infty } ^ { - 2 m - 2, - 2 m - 1 } ^ { \prime } d _ { r | c } ^ { m n } \left ( c \right ) P _ { - r - m - 1 } ^ { m } \left ( \xi \right ) + \sum _ { r = - 2 m, - 2 m + 1 } ^ { \infty } d _ { r } ^ { m n } \left ( c \right ) Q _ { m + r } ^ { m } \left ( \xi \right ) \right ), \\ \text{where} \, d ^ { m n } \quad \left ( \vec { \wedge } \, \text{and} \, d ^ { m n } \quad \left ( \vec { \wedge } \, \text{ are committed using} \,.$$

where d mn -2 m - | 2 glyph[epsilon1] ( c ) and d mn -2 m - | 1 glyph[epsilon1] ( c ) are computed using

$$\frac { d _ { - 2 m - 2 | \epsilon } ^ { m n } \left ( c \right ) } { d _ { - 2 m } ^ { m n } \left ( c \right ) } = \frac { c ^ { 2 } } { \left ( 2 m - 1 \right ) \left ( 2 m + 1 \right ) } \frac { 1 } { B _ { - 2 m - 2 } ^ { m } - ^ { - } } \frac { C _ { - 2 m - 4 } ^ { m } A _ { - 2 m - 2 } ^ { m } } { B _ { - 2 m - 4 } ^ { m } - ^ { - } } \frac { C _ { - 2 m - 6 } ^ { m } A _ { - 2 m - 4 } ^ { m } } { B _ { - 2 m - 6 } ^ { m } - ^ { - } } \cdots, \quad \text{(74)}$$

$$\frac { d _ { - 2 m - 1 | \epsilon } ^ { m n } ( c ) } { d _ { - 2 m + 1 } ^ { m n } ( c ) } = - \frac { c ^ { 2 } } { ( 2 m - 1 ) \, ( 2 m - 3 ) } \frac { 1 } { B _ { - 2 m - 1 } ^ { m } - } \frac { C _ { - 2 m - 3 } ^ { m } A _ { - 2 m - 1 } ^ { m } } { B _ { - 2 m - 3 } ^ { m } - } \frac { C _ { - 2 m - 5 } ^ { m } A _ { - 2 m - 3 } ^ { m } } { B _ { - 2 m - 5 } ^ { m } - } \cdots. \quad ( 7 5 )$$

For r &lt; -2 m -2 when n -m = even and r &lt; -2 m -1 when n -m = odd, the remaining expansion coefficients can be computed using

$$\frac { d _ { r | \epsilon } ^ { m n } \left ( c \right ) } { d _ { r + 2 | \epsilon } ^ { m n } \left ( c \right ) } = - \frac { A _ { r + 2 } ^ { m } } { B _ { r } ^ { m } - } \frac { C _ { r - 2 } ^ { m } A _ { r } ^ { m } } { B _ { r - 2 } ^ { m } - } \frac { C _ { r - 4 } ^ { m } A _ { r - 2 } ^ { m } } { B _ { r - 4 } ^ { m } - } \cdots.$$

## 3.2 Oblate Spheroidal Wave Functions

## 3.2.1 Angle Functions

The expressions for the prolate spheroidal angle functions, including the ones written in terms of the associated Legendre polynomials as well as the ones expanded as power series, can be transformed into those for the oblate spheroidal angle functions by letting c →-ic .

Figure 4. Characteristic and other special values for the oblate spheroidal wave functions for c = 10, m = 0 1 , , . . . , 29, and n = m,m +1 , . . . , m +29.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0074v1-figure-c9d4dc30-29c7-470e-81ea-00a8977f0902.jpg)

## 3.2.2 Radial Functions

The oblate spheroidal radial functions of the first and second kinds can be written in terms of the spherical Bessel and Neumann functions, respectively:

$$m / 2$$

$$I _ { m n } ^ { \sim } \left ( - i c, i \xi \right ) & = r _ { m n } \left ( - i c \right ) \quad \left \{ 1 + \frac { 1 } { \xi ^ { 2 } } \right \} \quad \times \\ \sum _ { r = 0, 1 } ^ { \infty } \left ( - 1 \right ) ^ { ( r - ( n - m ) ) / 2 } d _ { r } ^ { m n } \left ( - i c \right ) \frac { ( 2 m + r )! } { r! } j _ { m + r } \left ( c \xi \right ), \\ \\ R _ { m n } ^ { ( 2 ) } \left ( - i c, i \xi \right ) & = F _ { m n } \left ( - i c \right ) ^ { - 1 } \left ( 1 + \frac { 1 } { \xi ^ { 2 } } \right ) ^ { m / 2 } \times \\ \sum _ { r = 0, 1 } ^ { \infty } \left ( - 1 \right ) ^ { ( r - ( n - m ) ) / 2 } d _ { r } ^ { m n } \left ( - i c \right ) \frac { ( 2 m + r )! } {. } y _ { m + r } \left ( c \xi \right ),$$

$$R _ { m n } ^ { ( 1 ) } \left ( - i c, i \xi \right ) = F _ { m n } \left ( - i c \right ) ^ { - 1 } \left ( 1 + \frac { 1 } { \xi ^ { 2 } } \right ) ^ { m / 2 } \times \\ \sum _ { r = 0, 1 } ^ { \infty } ^ { \prime } \left ( - 1 \right ) ^ { \left ( r - \left ( n - m \right ) \right ) / 2 } d _ { r } ^ { m n } \left ( - i c \right ) \frac { \left ( 2 m + r \right )! } { r! } j _ { m + r } \left ( c \xi \right ),$$

$$R _ { m n } ^ { ( 2 ) } \left ( - i c, i \xi \right ) = F _ { m n } \left ( - i c \right ) ^ { - 1 } \left ( 1 + \frac { 1 } { \xi ^ { 2 } } \right ) ^ { m / 2 } \times \\ \sum _ { r = 0, 1 } ^ { \infty } ^ { \prime } \left ( - 1 \right ) ^ { \left ( r - \left ( n - m \right ) \right ) / 2 } d _ { r } ^ { m n } \left ( - i c \right ) \frac { \left ( 2 m + r \right )! } { r! } y _ { m + r } \left ( c \xi \right ),$$

Figure 5. The oblate spheroidal wave functions and their derivatives for c = 10, m = 10, and n = 10 11 , , . . . , 39.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0074v1-figure-94ae1bcc-f13d-44b3-82f8-4f1f79ce1e0f.jpg)

where

$$F _ { m n } \left ( - i c \right ) = \sum _ { r = 0, 1 } ^ { \infty } ^ { \prime } d _ { r } ^ { m n } \left ( - i c \right ) \frac { ( 2 m + r )! } { r! }.$$

The radial functions of the third and fourth kinds are linear combinations of those of the first and second kinds:

$$R _ { m n } ^ { ( 3 ) } \left ( - i c, i \xi \right ) = R _ { m n } ^ { ( 1 ) } \left ( - i c, i \xi \right ) + i R _ { m n } ^ { ( 2 ) } \left ( - i c, i \xi \right ),$$

$$R _ { m n } ^ { ( 4 ) } \left ( - i c, i \xi \right ) = R _ { m n } ^ { ( 1 ) } \left ( - i c, i \xi \right ) - i R _ { m n } ^ { ( 2 ) } \left ( - i c, i \xi \right ).$$

The Wronskian of the radial functions of the first and second kinds is given by

$$W _ { m n } \left ( - i c, i \xi \right ) & = R _ { m n } ^ { \left ( 1 \right ) } \left ( - i c, i \xi \right ) \frac { \partial } { \partial \xi } R _ { m n } ^ { \left ( 2 \right ) } \left ( - i c, i \xi \right ) - \frac { \partial } { \partial \xi } R _ { m n } ^ { \left ( 1 \right ) } \left ( - i c, i \xi \right ) R _ { m n } ^ { \left ( 2 \right ) } \left ( - i c, i \xi \right ) = \frac { 1 } { c \left ( \xi ^ { 2 } + 1 \right ) } \\ \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \end{$$

and is useful for validating computed values of these functions.

The radial functions are related to the angle functions by

$$S _ { m n } ^ { \left ( 1 \right ) } \left ( - i c, i z \right ) = k _ { m n } ^ { \left ( 1 \right ) } \left ( - i c \right ) R _ { m n } ^ { \left ( 1 \right ) } \left ( - i c, i z \right ),$$

$$S _ { m n } ^ { \left ( 2 \right ) } \left ( - i c, i z \right ) = k _ { m n } ^ { \left ( 1 \right ) } \left ( - i c \right ) R _ { m n } ^ { \left ( 2 \right ) } \left ( - i c, i z \right ),$$

where k (1) mn ( -ic ) and k (2) mn ( -ic ) are given by the same expressions as in the prolate case, provided that c →-ic . The expression for the radial functions of the first kind using the

spherical Bessel functions converges and is accurate for all values of ξ , except for ξ = 0, where the expression is undefined due to a divide by zero. The expression for the radial functions of the second kind using the spherical Neumann functions converges very slowly for values of ξ near 0 and is, therefore, inaccurate in these cases. Thus, these relationships can be used to construct secondary methods for computing these functions. The same procedures used in the prolate case for doing so can also be used here, provided that c, ξ → -ic, iξ where necessary.

A tertiary method for computing the radial functions of the second kind can also be constructed using a power series:

$$R _ { m n } ^ { \left ( 2 \right ) } \left ( - i c, i \xi \right ) = Q _ { m n } ^ { * } \left ( - i c \right ) R _ { m n } ^ { \left ( 1 \right ) } \left ( - i c, i \xi \right ) \left ( \arctan \left ( \xi \right ) - \frac { \pi } { 2 } \right ) + g _ { m n } \left ( - i c, i \xi \right ),$$

where

$$Q _ { m n } ^ { * } \left ( - i c \right ) = \frac { \left ( i ^ { - m } k _ { m n } ^ { \left ( 1 \right ) } \left ( - i c \right ) \right ) ^ { 2 } } { c } \sum _ { r = 0 } ^ { m } \alpha _ { r } ^ { m n } \left ( - i c \right ) \frac { \left ( 2 m - 2 r \right )! } { r! \left ( 2 ^ { m - r } \left ( m - r \right )! \right ) ^ { 2 } }, \quad n - m = \text{even}, \quad ( 8 6 )$$

$$Q _ { m n } ^ { * } \left ( - i c \right ) = - \frac { \left ( i ^ { - \left ( m + 1 \right ) } k _ { m n } ^ { \left ( 1 \right ) } \left ( - i c \right ) \right ) ^ { 2 } } { c } \sum _ { r = 0 } ^ { m } \alpha _ { r } ^ { m n } \left ( - i c \right ) \frac { \left ( 2 m - 2 r + 1 \right )! } { r! \left ( 2 ^ { m - r } \left ( m - r \right )! \right ) ^ { 2 } }, \quad n - m = \text{odd, } ( 8 7 )$$

$$\alpha _ { r } ^ { m n } \left ( - i c \right ) = \left [ \frac { d ^ { r } } { d x ^ { r } } \frac { 1 } { \left ( \sum _ { k = 0 } ^ { \infty } c _ { 2 k } ^ { m n } \left ( - i c \right ) x ^ { k } \right ) ^ { 2 } } \right ] _ { x = 0 },$$

$$g _ { m n } \left ( - i c, i \xi \right ) = \xi \left ( \xi ^ { 2 } + 1 \right ) ^ { - m / 2 } \sum _ { r = 0 } ^ { \infty } B _ { 2 r } ^ { m n } \left ( - i c \right ) \xi ^ { 2 r }, \quad n - m = \text{even},$$

$$g _ { m n } \left ( - i c, i \xi \right ) = \left ( \xi ^ { 2 } + 1 \right ) ^ { - m / 2 } \sum _ { r = 0 } ^ { \infty } B _ { 2 r } ^ { m n } \left ( - i c \right ) \xi ^ { 2 r }, \quad n - m = \text{odd}.$$

To compute α mn r ( -ic ), we use the following procedure, which was described in [6] and uses some properties of Cauchy products. Let C k = c mn 2 k ( -ic ), and expand the denominator in Eq. (88) as

$$\left ( \sum _ { k = 0 } ^ { \infty } c _ { 2 k } ^ { m n } \left ( - i c \right ) x ^ { k } \right ) ^ { 2 } = \sum _ { n = 0 } ^ { \infty } \sum _ { k = 0 } ^ { n } C _ { k } x ^ { k } C _ { n - k } x ^ { n - k } = \sum _ { n = 0 } ^ { \infty } B _ { n } x ^ { n },$$

$$B _ { n } = \sum _ { k = 0 } ^ { n } C _ { k } C _ { n - k },$$

$$\frac { 1 } { \left ( \sum _ { k = 0 } ^ { \infty } c _ { 2 k } ^ { m n } \left ( - i c \right ) x ^ { k } \right ) ^ { 2 } } & = \frac { 1 } { \infty } = \sum _ { n = 0 } ^ { \infty } A _ { n } x ^ { n },$$

$$\sum _ { n = 0 } ^ { \infty } A _ { n } x ^ { n } \sum _ { n = 0 } ^ { \infty } B _ { n } x ^ { n } = 1,$$

$$\sum _ { n = 0 } ^ { \infty } \left ( \sum _ { k = 0 } ^ { n } A _ { k } B _ { n - k } \right ) x ^ { n } = 1.$$

In order for this equality to hold,

$$A _ { 0 } B _ { 0 } = 1, \quad \sum _ { k = 0 } ^ { n } A _ { k } B _ { n - k } = 0, \quad n > 0.$$

Rearranging,

$$A _ { 0 } = \frac { 1 } { B _ { 0 } }, \quad A _ { n } = - \frac { 1 } { B _ { 0 } } \sum _ { k = 0 } ^ { n - 1 } A _ { k } B _ { n - k }, \quad n > 0.$$

We can now compute α mn r ( -ic ):

$$\alpha _ { r } ^ { m n } \left ( - i c \right ) = \left [ \frac { d ^ { r } } { d x ^ { r } } \sum _ { n = 0 } ^ { \infty } A _ { n } x ^ { n } \right ] _ { x = 0 } = A _ { r } r!.$$

Plugging Eq. (85) into the oblate version of Eq. (4) yields the following recurrence relation in B mn 2 r ( -ic ):

$$\alpha _ { 2 r } B _ { 2 r + 2 } ^ { m n } \left ( - i c \right ) + \beta _ { 2 r } B _ { 2 r } ^ { m n } \left ( - i c \right ) + \gamma _ { 2 r } B _ { 2 r - 2 } ^ { m n } \left ( - i c \right ) = h _ { 2 r },$$

where

$$\alpha _ { 2 r } & = \left ( 2 r + 2 \right ) \left ( 2 r + 3 \right ), & ( 1 0 0 )$$

$$\beta _ { 2 r } = \left ( 2 r + 1 \right ) \left ( 2 r - 2 m + 2 \right ) + m \left ( m - 1 \right ) - \lambda _ { m n } \left ( - i c \right ), \quad \ \ ( 1 0 1 )$$

$$\gamma _ { 2 r } = c ^ { 2 },$$

$$h _ { 2 r } = - 2 Q _ { m n } ^ { * } \left ( - i c \right ) \left ( i ^ { - m } k _ { m n } ^ { \left ( 1 \right ) } \left ( - i c \right ) \right ) ^ { - 1 } \sum _ { k = r - m + 1 } ^ { \infty } c _ { 2 k } ^ { m n } \left ( - i c \right ) \left ( m + 2 k \right ) \frac { \left ( m + k - 1 \right )! } { \left ( m + k - 1 - r \right )! r! } \ \left ( 1 0 3 \right )$$

when n -m = even, and

$$\alpha _ { 2 r } = \left ( 2 r + 1 \right ) \left ( 2 r + 2 \right ),$$

$$\beta _ { 2 r } & = 2 r \left ( 2 r - 2 m + 1 \right ) + m \left ( m - 1 \right ) - \lambda _ { m n } \left ( - i c \right ), & ( 1 0 5 )$$

$$\gamma _ { 2 r } = c ^ { 2 },$$

$$h _ { 2 r } = - 2 Q _ { m n } ^ { * } \left ( - i c \right ) \left ( i ^ { - \left ( m + 1 \right ) } k _ { m n } ^ { \left ( 1 \right ) } \left ( - i c \right ) \right ) ^ { - 1 } \left ( \sum _ { k = r - m } ^ { \infty } c _ { 2 k } ^ { m n } \left ( - i c \right ) \left ( m + 2 k + 1 \right ) \frac { \left ( m + k \right )! } { \left ( m + k - r \right )! r! }$$

$$- \sum _ { k = r - m + 1 } ^ { \infty } c _ { 2 k } ^ { m n } \left ( - i c \right ) ( m + 2 k ) \, \frac { ( m + k - 1 )! } { ( m + k - 1 - r )! r! } \right ) \\ m = \cap \partial \Delta \, & \, C \sin \, \cos \, \detect \, \max \, \Delta \, & \, \underline { R } ^ { m n } \left ( \Delta \, \partial \, \cap \, \cos \, \cos \, \det \, \max \, \Delta \, \det \, \cos \, \min \, \det \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \sin \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \tan \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \ sin \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \ cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \cos \, \colon$$

when n -m = odd. Given a starting value, B mn 0 ( -ic ), Eq. (99) can be used to compute B mn 2 r ( -ic ). The starting values are

$$B _ { 0 } ^ { m n } \left ( - i c \right ) = \left ( c R _ { m n } ^ { ( 1 ) } \left ( - i c, i 0 \right ) \right ) ^ { - 1 } - Q _ { m n } ^ { * } \left ( - i c \right ) R _ { m n } ^ { ( 1 ) } \left ( - i c, i 0 \right ), \quad n - m = \text{even}, \quad \text{(108)}$$

$$B _ { 0 } ^ { m n } \left ( - i c \right ) = - \left ( c R _ { m n } ^ { \left ( 1 \right ) \prime } \left ( - i c, i 0 \right ) \right ) ^ { - 1 }, \quad n - m = \text{odd}.$$

## 3.2.3 Calculating the Characteristic Value and Expansion Coefficients

The procedures for computing the characteristic value and expansion coefficients for the oblate case are exactly the same as those for the prolate case, provided that c →-ic .

## 4 Software

We implemented our software in C++ and MATLAB. Our code is called spheroidal and is primarily called via the command line. Many programming languages have interfaces to the command line (e.g., the system function in MATLAB), so they can access our code in this way. There are two programs: pro\_sphwv computes the prolate spheroidal wave functions; and obl\_sphwv computes the oblate spheroidal wave functions. Like any piece of software, our code contains many subtleties and nuances. We have commented and documented our code to explain all of these. Below, we describe some of the more important ones.

## 4.1 Using the MPFR Library

Except in a few places, our code uses the GNU MPFR library, which provides interfaces and routines for performing arbitrary precision arithmetic [3]. We created a C++ class, real , which encapsulates many of the features provided by GNU MPFR. These features include: basic arithmetic by overloading the + , -, * , and / operators; comparisons by overloading the &gt; , &gt;= , &lt; , &lt;= , == , and != operators; and some elementary functions, including abs , atan , cos , log , pow , and sin . GNU MPFR allows the programmer to specify the precision to use for these operations by setting the number of bits of precision. As a reference, single and double precision arithmetic found in most programming languages have 24 and 53 bits of precision, respectively. We experimented with several different levels of precision (as low as 24 and as high as 5000 bits of precision). In general, as we increased the precision, the accuracy of the computations increased as well. For lower c , m , and n , single and double precision were good enough. However, for higher values, using such low precision yielded very large errors, and only by increasing the precision were these errors reduced.

## 4.2 Prolate Spheroidal Wave Functions

The code for computing the prolate spheroidal wave functions, as well as the characteristic and other special values and expansion coefficients required for computing them, is called via the command line. Every command uses the same general structure. The program is called pro\_sphwv . Seven arguments are required:

- 1. -max\_memory is the maximum amount of memory, in MB, the program can use before automatically terminating;
- 2. -prec is the number of bits of precision to use;
- 3. -verbose is whether the program should output diagnostic and other information about the computations (accepted values are 'y' and 'n');

- 4. -c is equal to ka , where k is the wavenumber and 2 a is the interfocal distance;
- 5. -m is one modal value;
- 6. -n is the other modal value and is equal to m,m +1 , . . . ; and
- 7. -w is what the program should do.

Following the -w argument are zero or more arguments, the number and type of which depend on the value given for the -w argument. The following sequence of commands computes the characteristic and other special values and expansion coefficients required for computing the prolate spheroidal wave functions:

```
<_C_>
```

- \

The command, ./pro\_sphwv ... -w lambda , uses method (2) from Section 3.1.3 to compute the characteristic value. Method (2) requires an approximate value for the characteristic value, which can be computed by using method (1). We wrote a small MATLAB program for doing so. Once this program has run, the preceeding sequence of commands can be called using one command:

```
       ./pro_sphwv -max_memory 2000 -prec 100 -verbose y -c 10.0 -m 0 -n 0 \
            -w everything -n_dr 10 -dr_min 1.0e-200 -n_dr_neg 10 \
            -dr_neg_min 1.0e-200 -n_c2k 10 -c2k_min 1.0e-200
```

The computed values are stored as ASCII in .txt files in the data directory. This way, they can be reused later on without having to recompute them from scratch. The following two commands compute the prolate spheroidal wave functions over a range of values:

```
       ./pro_sphwv -max_memory 2000 -prec 100 -verbose n -c 10.0 -m 0 -n 0 -w S1 \
            -a -1.0 -b 1.0 -d 0.125 -arg_type eta \
            -p 20 > data/pro_00010000_000_000_S1.txt
       ./pro_sphwv -max_memory 2000 -prec 100 -verbose n -c 10.0 -m 0 -n 0 -w R \
            -a 1.0 -b 9.0 -d 0.125 -arg_type xi -which R1_1,R1_2,R2_1,R2_2 \
            -p 20 > data/pro_00010000_000_000_R.txt
```

- \

where the range is determined by the values passed in for -a (the starting point), -b (the ending point), and -d (the spacing between the points). Normally, the angle functions take η as their argument, and the radial functions take ξ as their argument. The argument, -arg\_type , allows one to use η = cos ( xπ ) for the angle functions, where the range is over x , not η . In this case, -arg\_type should be set to theta/pi . Likewise, for the radial functions, one can use ξ = ( x 2 +1) 1 2 / , where the range is over x , not ξ . In this case, -arg\_type should be set to x . The argument, -which , tells the program which of the different methods for computing the radial functions should be used. The argument, -p , is how many digits of precision to output. There is no guarantee that all of these digits will be accurate. How many digits are actually accurate depends on the level of precision and how many expansion coefficients are computed and used.

## 4.3 Oblate Spheroidal Wave Functions

The calling convention for obl\_sphwv is the same as pro\_sphwv . There are some differences since a couple extra sets of expansion coefficients are computed. The following sequence of commands computes the characteristic and other special values and expansion coefficients required for computing the oblate spheroidal wave functions:

```
                       commands computes the characteristic and other special values and expansion coefficients
                       required for computing the oblate spheroidal wave functions:

                      ./obl_sphwv  -max_memory  2000  -prec  100  -verbose  y  -c  10.0  -m  0  -n  0  \
                           -w  lambda
                       ./obl_sphwv  -max_memory  2000  -prec  100  -verbose  y  -c  10.0  -m  0  -n  0  -w  dr  \
                           -n_dr  10  -dr_min  1.0e-200
                       ./obl_sphwv  -max_memory  2000  -prec  100  -verbose  y  -c  10.0  -m  0  -n  0  \
                           -w  dr_neg  -n_dr_neg  10  -dr_neg_min  1.0e-200
                       ./obl_sphwv  -max_memory  2000  -prec  100  -verbose  y  -c  10.0  -m  0  -w  N
                       ./obl_sphwv  -max_memory  2000  -prec  100  -verbose  y  -c  10.0  -m  0  -w  F
                       ./obl_sphwv  -max_memory  2000  -prec  100  -verbose  y  -c  10.0  -m  0  -w  k1
                       ./obl_sphwv  -max_memory  2000  -prec  100  -verbose  y  -c  10.0  -m  0  -w  k2
                       ./obl_sphwv  -max_memory  2000  -prec  100  -verbose  y  -c  10.0  -m  0  -n  0  \
                            -w  c2k  -n_c2k  10  -c2k_min  1.0e-200
                        ./obl_sphwv  -max_memory  2000  -prec  100  -verbose  y  -c  10.0  -m  0  -w  Q
                        ./obl_sphwv  -max_memory  2000  -prec  100  -verbose  y  -c  10.0  -m  0  -n  0  \
                           -w  B2r  -n_B2r  10  -B2r_min  1.0e-200

                        Like in the prolate case, the command,./pro_sphwv ...  -w  lambda, must be preceeded by
                        running a small MATLAB program to compute an approximate value for the characteristic
```

- \

Like in the prolate case, the command, ./pro\_sphwv ... -w lambda , must be preceeded by running a small MATLAB program to compute an approximate value for the characteristic value. The preceeding sequence of commands can be called using one command:

```
       ./obl_sphwv -max_memory 2000 -prec 100 -verbose y -c 10.0 -m 0 -n 0 \
            -w everything -n_dr 10 -dr_min 1.0e-200 -n_dr_neg 10 \
            -dr_neg_min 1.0e-200 -n_c2k 10 -c2k_min 1.0e-200 -n_B2r 10 \
            -B2r_min 1.0e-200
```

The computed values are stored as ASCII in .txt files in the data directory. This way, they can be reused later on without having to recompute them from scratch. The following two commands compute the oblate spheroidal wave functions over a range of values:

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0074v1-figure-b2177d83-8e86-4fec-83c6-df6b6e321b44.jpg)

Figure 6. The relative error of the computed Wronskian when using different combinations of the methods for computing the prolate spheroidal radial functions for c = 10, m = 10, and n = 39.

```
       ./obl_sphw -max_memory 2000 -prec 100 -verbose n -c 10.0 -m 0 -n 0 -w S1 \
            -a -1.0 -b 1.0 -d 0.125 -arg_type eta \
            -p 20 > data/obl_00010000_000_000_000_S1.txt
       ./obl_sphw -max_memory 2000 -prec 100 -verbose n -c 10.0 -m 0 -n 0 -w R \
            -a 0.0 -b 8.0 -d 0.125 -arg_type xi \
            -which R1_1,R1_2,R2_1,R2_2,R2_31,R2_32 \
            -p 20 > data/obl_00010000_000_000_000_R.txt
```

where the arguments are the same as in the prolate case. The only difference is that, for the radial functions, there is only one valid value for -arg\_type , which is xi .

## 4.4 Using the Wronskian

Two methods were given for computing the prolate spheroidal radial functions of the first kind. One method uses the spherical Bessel functions, and the other method uses a power series. Call these methods R1\_1 and R1\_2 , respectively. Likewise, two methods were given for computing the prolate spheroidal radial functions of the second kind. One method uses the spherical Neumann functions, and the other method uses the associated Legendre polynomials of the second kind. Call these methods R2\_1 and R2\_2 , respectively. In each of these pairs of methods, one is better for certain values of ξ than the other, and vice versa. In particular, R1\_1 is better for larger ξ and R1\_2 is better for smaller ξ . Similarly, R2\_1 is better for larger ξ and R2\_2 is better for smaller ξ . However, when to use which method is not always clear: the exact value of ξ below which one is better and above which the other is better is different for different values of c , m , and n . We solved this dilemma in the following manner. For a given ξ , all four methods are used to compute their respective functions, and the combination that yields the smallest error in the computed Wronskian is used. An example of this can be seen in Figure 6. As expected, the combination, R1\_2, R2\_2 , was superior for smaller ξ and the combination, R1\_1, R2\_1 , was superior for larger ξ .

Figure 7. The relative error of the computed Wronskian when using different combinations of the methods for computing the oblate spheroidal radial functions for c = 10, m = 10, and n = 39.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0074v1-figure-c6db48cc-c118-444b-9a44-fcca81906db3.jpg)

This procedure is also used for the oblate spheroidal radial functions. There are two methods for computing the oblate spheroidal radial functions of the first kind, one that uses the spherical Bessel functions and one that uses a power series. Call these methods R1\_1 and R1\_2 , respectively. There are three methods for computing the oblate spheroidal radial functions of the second kind, one that uses the spherical Neumann functions, one that uses the associated Legendre polynomials of the second kind, and one that uses a power series. Call these methods R2\_1 , R2\_2 , and R2\_3 , respectively. Internally, R2\_3 uses the oblate spheroidal radial functions of the first kind, so when R1\_1 is used, call this method R2\_31 , and when R1\_2 is used, call this method R2\_32 . Thus, there are eight possible combinations, but only six are considered: R1\_1 and R2\_1 ; R1\_1 and R2\_2 ; R1\_1 and R2\_31 ; R1\_2 and R2\_1 ; R1\_2 and R2\_2 ; and R1\_2 and R2\_32 . For a given ξ , of these six, the one that yields the smallest error in the computed Wronskian is used. An example of this can be seen in Figure 7. As expected, the combination, R1\_1, R2\_2 , was superior for smaller ξ and the combination, R1\_1, R2\_1 , was superior for larger ξ .

## 4.5 Solving Forward and Backward Recurrences

Several recurrence relations need to be solved in order to compute the characterstic and other special values, expansion coefficients, and special functions required by many of the expressions for the spheroidal wave functions. These recurrence relations can either have a starting value (e.g., for B mn 2 r ( -ic )) or some kind of normalization scheme (e.g., for d mn r ( ic )). Except for one case, all of the recurrence relations encountered in the previous sections were homogeneous. Depending on whether the solution to a particular recurrence relation grows or decays, the method of solution is different. For solutions that grow (e.g., the spherical Neumann functions), the forward recurrence approach is used (i.e., the recurrence relation is used directly to compute succeeding values). For solutions that decay (e.g., d mn r ( ic )), the forward recurrence approach is numerically unstable. Instead, the continued fraction approach, which is described in [4], is used.

or,

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0074v1-figure-2e74da7f-6302-46f6-b620-ba187467db1e.jpg)

n - m

Figure 8. The order of magnitude of the expansion coefficients, B mn 2 r ( -ic ), for c = 25, m = 49, and n = 49 50 , , . . . , 98.

Nonhomogeneous recurrence relations are more complicated. For solutions that grow, the foward recurrence approach can still be used. However, for solutions that decay, the continued fraction approach no longer works. Instead, we use the tridiagonal matrix method described in [16, 17]: the recurrence relation is written for each index, these are combined into a tridiagonal system of equations, and this system is inverted. The case of computing the expansion coefficients, B mn 2 r ( -ic ), is further complicated by the fact that, for many values of c , m , and n , B mn 2 r ( -ic ) grows for lower r , but decays for higher r . See, for example, Figure 8. Thus, the forward recurrence approach is used when B mn 2 r ( -ic ) is growing, and the tridiagonal matrix method is used when B mn 2 r ( -ic ) is decaying.

## 4.6 Optimizing Factorials and Rising Pochhammer Symbols

Many of the expressions for the characteristic and other special values, expansion coefficients, and spheroidal wave functions involve infinite series expansions. The terms in these series often include factorials and rising Pochhammer symbols. Computing these can be rather expensive: for example, computing a ! and ( a ) k require a -1 and k -1 multiplies, respectively. Because a and k are usually functions of the term number, as the number of terms computed in these series increases, the number of multiplies required increases quadratically. This can be very slow. However, this problem can be solved in the following manner: the products arising from these factorials and rising Pochhammer symbols can be factored across terms, so that only a constant number of multiplies are required for each term.

For example, consider here, again, Eq. (103):

$$h _ { 2 r } = - 2 Q _ { m n } ^ { * } \left ( - i c \right ) \left ( i ^ { - m } k _ { m n } ^ { ( 1 ) } \left ( - i c \right ) \right ) ^ { - 1 } \sum _ { k = r - m + 1 } ^ { \infty } c _ { 2 k } ^ { m n } \left ( - i c \right ) \left ( m + 2 k \right ) \frac { ( m + k - 1 )! } { ( m + k - 1 - r )! r! }, \\ \sim$$

$$\alpha, \\ h _ { 2 r } = - 2 Q _ { m n } ^ { * } \left ( - i c \right ) \left ( i ^ { - m } k _ { m n } ^ { \left ( 1 \right ) } \left ( - i c \right ) \right ) ^ { - 1 } \sum _ { k = r - m + 1 } ^ { \infty } c _ { 2 k } ^ { m n } \left ( - i c \right ) \left ( m + 2 k \right ) a _ { k }, \quad \left ( 1 1 1 \right )$$

where

$$a _ { k } & = \frac { ( m + k - 1 )! } { ( m + k - 1 - r )! r! }.$$

Computing a k requires O m ( + k + ) multiplies. For large r m r , , and k , carrying out these operations for every a k would be time consuming. Fortunately, a k is related to a k -1 by a constant number of operations. Starting with a r -m +1 , every a k can be computed recursively and plugged into Eq. (111). To begin,

$$a _ { r - m + 1 } = 1.$$

For increasing k ,

$$a _ { k } = a _ { k - 1 } \frac { ( m + k - 1 ) } { ( m + k - r - 1 ) }.$$

## 5 Conclusion

The spheroidal wave functions are among the most complicated special functions. There are no simple ways to compute them. However, because the solutions to so many interesting problems require them, software for computing them accurately is needed. We have developed computational software for doing so using C++, MATLAB, and GNU MPFR, a library for performing arbitrary precision arithmetic. In this paper, we have described the prolate and oblate spheroidal coordinate systems and wave functions, methods for deriving analytical expressions for computing them, and our software that implements these expressions. Our software includes many novel features. Some of these features include: using arbitrary precision arithmetic; adaptively choosing the number of expansion coefficients to compute and use; and using the Wronskian to choose from several different methods for computing the spheroidal radial functions to improve their accuracy. We have made our software freely available on our webpage.

## References

- [1] J. J. Bowman, T. B. A. Senior, and P. L. E. Uslenghi, Electromagnetic and acoustic scattering by simple shapes (Hemisphere Publishing Corporation, New York) (1987).
- [2] D. Slepian, 'Some comments on Fourier analysis, uncertainty and modeling', SIAM Review 25(3) (1983).
- [3] L. Fousse, G. Hanrot, V. Lefevre, P. Pelissier, and P. Zimmermann, 'MPFR: A multipleprecision binary floating-point library with correct rounding', ACM Trans. Math. Softw. 33(2) (2007).
- [4] C. Flammer, Spheroidal wave functions (Dover Publications, Mineola) (2005).
- [5] T. Do-Nhat and R. H. MacPhie, 'Accurate values of prolate spheroidal radial functions of the second kind', Can. J. Phys. 75 (1997).

- [6] S. Zhang and J. Jin, Computation of special functions (Wiley-Interscience) (1996).
- [7] W. J. Thompson, Atlas for computing mathematical functions (John Wiley and Sons) (1997).
- [8] W. J. Thompson, 'Spheroidal wave functions', Computing in Science and Engineering 1(3) (1999).
- [9] A. L. V. Buren and J. E. Boisvert, 'Accurate calculation of prolate spheroidal radial functions of the first kind and their first derivatives', Quart. Appl. Math. 60 (2002).
- [10] A. L. V. Buren and J. E. Boisvert, 'Improved calculation of prolate spheroidal radial functions of the second kind and their first derivatives', Quart. Appl. Math. 62 (2004).
- [11] L. Li, M. Leong, T. Yeo, P. Kooi, , and K. Tan, 'Computations of spheroidal harmonics with complex arguments: a review with an algorithm', Phys. Rev. 58 (1998).
- [12] P. E. Falloon, P. C. Abbott, and J. B. Wang, 'Theory and computation of the spheroidal wave functions', J. Physics A 36 (2003).
- [13] J. Caldwell, 'Computation of eigenvalues of spheroidal harmonics using relaxation', J. Phys. A 21 (1988).
- [14] 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', Computer Physics Communications 185 (2014).
- [15] D. B. Hodge, 'Eigenvalues and eigenfunctions of the spheroidal wave equation', J. Math. Phys. 11 (1970).
- [16] F. W. J. Olver, 'Numerical solution of second-order linear difference equations', Journal of Research of the National Bureau of Standards 71B (1967).
- [17] F. W. J. Olver and D. J. Sookne, 'Note on backward recurrence algorithms', Mathematics of Computation 26(120) (1972).
Ross Adelman Nail A. Gumerov Ramani Duraiswami
Software for Computing the Spheroidal Wave Functions Using Arbitrary Precision Arithmetic
The spheroidal wave functions, which are the solutions to the Helmholtz
equation in spheroidal coordinates, are notoriously difficult to compute.
Because of this, practically no programming language comes equipped with the
means to compute them. This makes problems that require their use hard to
tackle. We have developed computational software for calculating these special
functions. Our software is called spheroidal and includes several novel
features, such as: using arbitrary precision arithmetic; adaptively choosing
the number of expansion coefficients to compute and use; and using the
Wronskian to choose from several different methods for computing the spheroidal
radial functions to improve their accuracy. There are two types of spheroidal
wave functions: the prolate kind when prolate spheroidal coordinates are used;
and the oblate kind when oblate spheroidal coordinate are used. In this paper,
we describe both, methods for computing them, and our software. We have made
our software freely available on our webpage.