1.3 A simple example
Before entering the details of spectral methods in Sections 2, 3 and 4, let us give here their spirit with
the simple example of the Poisson equation in a spherical shell:
where is the Laplace operator (93) expressed in spherical coordinates (see also Section 3.2).
We want to solve Equation (7) in the domain where . This
Poisson equation naturally arises in numerical relativity when, for example, solving for initial conditions or
the Hamiltonian constraint in the 3+1 formalism : the linear part of these equations can be
cast in form (7), and the nonlinearities put into the source , with an iterative scheme on
First, the angular parts of both fields are decomposed into a (finite) set of spherical harmonics
(see Section 3.2.2):
with a similar formula relating to the radial functions . Because spherical harmonics are
eigenfunctions of the angular part of the Laplace operator, the Poisson equation can be equivalently solved
as a set of ordinary differential equations for each couple , in terms of the coordinate :
We then map
and decompose each field in a (finite) basis of Chebyshev polynomials (see Section 2.4.3):
Each function can be regarded as a column-vector of its coefficients in this
basis; the linear differential operator on the left-hand side of Equation (9) being, thus, a matrix
acting on the vector:
with being the vector of the coefficients of . This matrix can be computed
from the recurrence relations fulfilled by the Chebyshev polynomials and their derivatives (see Section 2.4.3
The matrix is singular because problem (7) is ill posed. One must indeed specify boundary
conditions at and . For simplicity, let us suppose
To impose these boundary conditions, we adopt the tau methods (see Section 2.5.2): we build the matrix
, taking and replacing the last two lines by the boundary conditions, expressed in terms of the
coefficients from the properties of Chebyshev polynomials:
Equations (14) are equivalent to boundary conditions (13), within the considered spectral approximation,
and they represent the last two lines of , which can now be inverted and give the coefficients of the
If one summarizes the steps:
- Setup an adapted grid for the computation of spectral coefficients (e.g., equidistant in the
angular directions and Chebyshev–Gauss–Lobatto collocation points; see Section 2.4.3).
- Get the values of the source on these grid points.
- Perform a spherical-harmonics transform (for example, using some available library ),
followed by the Chebyshev transform (using a Fast Fourier Transform (FFT), or a
Gauss–Lobatto quadrature) of the source .
- For each couple of values , build the corresponding matrix with the boundary
conditions, and invert the system (using any available linear-algebra package) with the
coefficients of as the right-hand side.
- Perform the inverse spectral transform to get the values of on the grid points from its
A numerical implementation of this algorithm has been reported by Grandclément et al. , who
have observed that the error decayed as , provided that the source was smooth. Machine
round-off accuracy can be reached with , which makes the matrix inversions of
step 4 very cheap in terms of CPU and the whole method affordable in terms of memory usage.
These are the main advantages of using spectral methods, as shall be shown in the following