A finite difference construction of the spheroidal wave functions

Abstract

A fast and simple finite difference algorithm for computing the spheroidal wave functions is described. The resulting eigenvalues and eigenfunctions for real and complex spheroidal bandwidth parameter, c, agree with those in the literature from four to more than eleven significant figures. The validity of this algorithm in the extreme parameter regime, up to c2=1014, is demonstrated. Furthermore, the algorithm generates the spheroidal functions for complex order m. The coefficients of the differential equation can be simply modified so that the algorithm may solve any second-order differential equation in Sturm–Liouville form. The prolate spheroidal functions and the spectral concentration problem in relation to band-limited and time-limited signals is discussed. We review the properties of these eigenfunctions in the context of Sturm–Liouville theory and the implications for a finite difference algorithm. A number of new suggestions for data fitting using prolate spheroidal wave functions with a heuristic for optimally choosing the value of c and the number of basis functions are described. Program summaryProgram title: SWF_8thOrderCatalogue identifier: AEQE_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEQE_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 1081No. of bytes in distributed program, including test data, etc.: 160,312Distribution format: tar.gzProgramming language: Matlab R2009b.Computer: Designed to run on any computer capable of running Matlab 2009b with at least 2 GB of RAM in order to handle moderate grid sizes. With an appropriate change of syntax, the program may be easily adapted to Maple, Mathematica, Fortran, IDL and any other software capable of performing the diagonalization of large matrices.Operating system: Any operating system which will run Matlab, Mathematica, Fortran or any other language capable of performing large matrix diagonalizations.Has the code been vectorized or parallelized?: Tested with dual core and quad core systems. The algorithm will also work with single core systems.RAM: 1372 MB total used by Matlab for a grid with 4001 points, 41 MB used to store eigenfunctions, grid and spectrum arrays (4001 points). More RAM is used for larger grids. For example, with 8000 grid points, 162 MB of RAM is required to store the eigenfunction and eigenvalue arrays.Classification: 4.3.External routines: The program uses Matlab’s internal ‘eig’ routine.Nature of problem:The problem is to construct the angular eigenfunctions of the Laplacian in three dimensional, spheroidal coordinates. These are the prolate, oblate and generalized spheroidal wave functions and to compute the corresponding eigenvalues. Equivalently, the task can be seen as generating the angular functions which arise when solving the Helmholtz wave equation by separation of variables in three dimensional, spheroidal coordinates: [∂η(1−η2)∂η+λlm(c)−c2η2−m21−η2]Slm=0. This task often arises in the solution of problems with axial symmetry although setting c=0 restores spherical symmetry. More generally, the coefficient function handles, CD1E and CD2E, can be redefined by the user to match the coefficients of the 2nd and 1st derivatives of any second order Sturm–Liouville type equation, for example, the harmonics of a tri-axial ellipsoid. Hence this algorithm and manuscript provide a foundation for solving a range of problems that have application beyond the spheroidal problems considered here.Solution method:The method of solution is the ‘finite difference method’. The spatial grid is discretized into N points, N−2 of which comprise the ‘interior grid’ and 2 points are the boundary points. The spheroidal differential operator is discretized which arises from separation of variables of the Laplacian in three dimensional, spheroidal coordinates on the interior grid using 8th order finite difference formulas. The boundary conditions for the spheroidal wave functions are implemented implicitly via finite difference operators which relate the boundary points back to the interior points. This is done using sliding off-centered differences at 8th order. The discretization of the Laplace operator and implicit implementation of the boundary conditions leads to a discretized eigenvalue problem. The eigenvectors give the discretized eigenfunctions of the spheroidal differential operator and the eigenvalues give the spectrum of the differential operator. Points at the boundary are reconstructed using forward and backward difference operators. The eigenfunctions are numerically normalized using a 6th order “Boole’s Rule” integration procedure.Restrictions:The current version of this algorithm implements the option of both Dirichlet and Neumann boundary conditions, which is chosen by the user by the “BC” switch. Mixed boundary conditions can also be implemented by the user, by modifying the boundary condition ‘IF’ statements. When solving with a very large concentration parameter, |c|, one must use a large number of grid points which would result in longer computational times and require more RAM.Unusual features:This program solves for the angular eigenfunctions of the spheroidal wave equation in three dimensions. Due to the finite difference approach, one is able to solve over non-standard domains such as spheroidal caps or spheroidal belts. The program is capable of solving for the generalized spheroidal wavefunctions with complex geometric parameter c. Given a sufficiently large number of grid points, the program can generate spheroidal eigenfunctions and eigenvalues for extreme parameter regimes, for example, for |c|∼108 with a grid of 20,000 points. Furthermore, the program can generate spheroidal wavefunctions for non-integer and complex order parameter, m, which may correspond to some analytic continuation of the spheroidal wave functions. In particular, for real, non-integer values of m, this corresponds to an axially symmetric ellipsoid with a defect angle where the 2π periodicity symmetry about the rotation axis becomes 2πα periodicity, where α is some defect factor.Additional comments:Main User-Input Parameters:The main input parameters are located at the top of the code. The parameter C2 sets the value of the square of the spheroidal concentration parameter, c2. The parameter m is the order parameter which is typically an integer such as for the Legendre functions, but can also be non-integer and complex valued. The switch BC, sets the boundary conditions for the differential equation. A value BC=0 gives Dirichlet boundary conditions and BC=1 gives Neumann boundary conditions. The values of MinTheta and MaxTheta set the minimum and maximum values of the angular coordinate, θ, which specify the domain of the differential equation. The parameter Npts sets the number of grid points, with larger grids resulting in more accurate eigenfunctions and spectra.General comments:To obtain spheroidal harmonics the eigenfunctions need to be multiplied by a complex exponential factor eimφ or sinusoidal functions for real solutions. The eigenfunctions generated by this program may be normalized numerically. Since normalization conventions differ by application, normalization is left for the user to implement. For this purpose, we have encoded a plain, 6th order Boole rule unit normalization which is easy to modify.Finally, we note that if the user modifies the function handles, CD1E and CD2E then the program will solve any 2nd order ordinary differential equation with Dirichlet or Neumann boundary conditions. Non-homogeneous terms can be added by modifying the entries in the finite difference matrix where c2 and m2 appear.Running time:On a laptop with an Intel Core i3-2350M CPU (2.30 GHz), with 4.00 GB of RAM, the program takes 149 s for a grid with 4000 points. Running time increases for larger grid sizes. For example, increasing the grid size to 8000 points increased the run time to 1138 s on the same machine.