Elsevier Science Home
Computer Physics Communications Program Library
Full text online from Science Direct
Programs in Physics & Physical Chemistry
CPC Home

[Licence| Download | New Version Template] aeqe_v1_0.tar.gz(156 Kbytes)
Manuscript Title: A finite difference construction of the spheroidal wave functions
Authors: Daniel X. Ogburn, Colin L. Waters, Murray D. Sciffer, Jeff A. Hogan, Paul C. Abbott
Program title: SWF_8thOrder
Catalogue identifier: AEQE_v1_0
Distribution format: tar.gz
Journal reference: Comput. Phys. Commun. 185(2014)244
Programming language: Matlab R2009b.
Computer: Designed to run on any computer capable of running Matlab 2009b with at least 2GB 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 vectorised 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.
Keywords: Slepian function, Spheroidal wave function, Spheroidal harmonics, Spectral concentration problem, Prolate spheroidal wave equation, Oblate spheroidal wave equation, Spheroidal cap harmonics, Basis function expansion, Finite difference method, Sturm-Liouville problem.
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η + (-c2η2 - m2/(1-η2))] Slm = λml(c)Slm

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 "CD2E" and "CD1E", 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 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 discretised into N points, N - 2 of which comprise the 'interior grid' and 2 points are the boundary points. The spheroidal differential operator is discretised which arises from separation of variables of the Laplacian in three dimensional, spheroidal coordinates on the interior grid using 8th order finite difference formulae. 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 discretised eigenvalue problem. The eigenvectors give the discretised 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.

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:
If the spheroidal harmonics are required one must multiply the spheroidal wave functions by a complex exponential expimφ 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.30GHz), with 4.00 GB of RAM, the program takes 149 seconds 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 seconds on the same machine.