The basis function n-body force solver is optimal for studying the global response of galaxies to perturbations or stability (Earn & Sellwood 1995). This technique was developed for astrophysical problems by Clutton-Brock (1972, 1973), Kalnajs (1976), Fridman & Polyachenko (1984) and more recently by Hernquist & Ostriker (1992) who dubbed it the self-consistent field (SCF) method. Orthogonal function expansions are attractive Poisson equation solvers for two reasons: 1) the expansions can be chosen to filter the structure over an interesting range of scales and simultaneously suppress small-scale noise; and 2) the algorithm is computationally efficient, scaling linearly with the number of particles. Mathematically, this entire class of algorithms relies on the general properties of the Sturm-Liouville equation (SLE) of which the Poisson equation is a particular case. This same approach is common in perturbation theories and so facilitates direct comparison between n-body simulation and linear perturbation theory. In addition, this approach is straightforward to parallelize (e.g. Hernquist, Sigurdsson & Bryan 1995); we find the algorithm scales linearly with the number of processors with low overhead. If the basis set resembles the equilibrium galaxy, most of the computational work is concentrated on resolving the perturbation rather than the equilibrium.
This last point is also a disadvantage of this technique in applications to date. If the equilibrium does not look like the basis set, the technique becomes less efficient and noisy because the expansion series must be sufficiently long to represent the equilibrium even without the perturbation. This paper describes a general method based on a numerical construction of orthogonal bases which remedies this situation. Solutions to the fundamental equation, the Sturm-Liouville equation, are well-understood and well-behaved. A number of recently published algorithms take advantage of the special properties of this differential equation to yield high-accuracy solutions with low computational work. Harnessing these developments to our needs leads to an algorithm for computing orthogonal bases whose lowest-order function matches any given any regular equilibrium; spherical and three-dimensional cylindrical solutions are described in detail here. The basic algorithm will be described in §2.
For the spherical case, the proposed algorithm is competitive in
performance with evaluation by recursion relation used for the
published bases cited above and has reproduced them with high accuracy
as a check. The cylindrical basis is a bit more cumbersome: one may
rely on the same numerical solution to tailor the basis in the radial
or vertical direction but not both simultaneously. Here, I choose to
derive the radial basis numerically. The lowest-order radial basis
functions then take the form 
.
These may then
be adapted to the background by principal component analysis.
Although more cumbersome to implement and more time consuming to
execute than the spherical case, it is still fast relative to
non-expansion-based solvers. The details of the cylindrical basis are
given in §3.2.1.