A simple example

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 $(r,𝜃,φ )$ (see also Section 3.2). We want to solve Equation (7

) in the domain where $0 < Rmin ≤ r ≤ Rmax, 𝜃 ∈ [0,π ],φ ∈ [0,2 π)$ . This Poisson equation naturally arises in numerical relativity when, for example, solving for initial conditions or the Hamiltonian constraint in the 3+1 formalism [97

]: 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 $m {Yℓ }$ (see Section 3.2.2):

with a similar formula relating $ϕ$ to the radial functions $fℓm(r)$ . 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 $(ℓ,m )$ , in terms of the coordinate $r$ :

We then map

[Rmin, Rmax ] → [− 1,1 ] 2r − R − R r ↦→ ξ = -------max----min-, (10 ) Rmax − Rmin

and decompose each field in a (finite) basis of Chebyshev polynomials ${Ti}i=0...N$ (see Section 2.4.3):

∑N sℓm(ξ) = ciℓmTi(ξ), i=0 N ∑ fℓm(ξ) = aiℓmTi (ξ). (11 ) i=0

Each function $fℓm (r)$ can be regarded as a column-vector $A ℓm$ of its $N + 1$ coefficients $aiℓm$ in this basis; the linear differential operator on the left-hand side of Equation (9

) being, thus, a matrix $Lℓm$ acting on the vector:

with $S ℓm$ being the vector of the $N + 1$ coefficients $c iℓm$ of $s (r) ℓm$ . This matrix can be computed from the recurrence relations fulfilled by the Chebyshev polynomials and their derivatives (see Section 2.4.3 for details).

The matrix $L$ is singular because problem (7 ) is ill posed. One must indeed specify boundary conditions at $r = Rmin$ and $r = Rmax$ . For simplicity, let us suppose

To impose these boundary conditions, we adopt the tau methods (see Section 2.5.2): we build the matrix $L¯$ , taking $L$ 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 $¯L$ , which can now be inverted and give the coefficients of the solution $ϕ$ .

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 [152 ]), followed by the Chebyshev transform (using a Fast Fourier Transform (FFT), or a Gauss–Lobatto quadrature) of the source $σ$ .
For each couple of values $(ℓ,m )$ , build the corresponding matrix $¯ L$ 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 coefficients.

A numerical implementation of this algorithm has been reported by Grandclément et al. [109 ], who have observed that the error decayed as $−ℓmax −N e ⋅ e$ , provided that the source $σ$ was smooth. Machine round-off accuracy can be reached with $ℓmax ∼ N ∼ 30$ , 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 sections.