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:
Δ ϕ = σ, (7 )
where Δ is the Laplace operator (93View Equation) expressed in spherical coordinates (r,𝜃,φ ) (see also Section 3.2). We want to solve Equation (7View Equation) 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 [97Jump To The Next Citation Point]: the linear part of these equations can be cast in form (7View Equation), 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):

ℓm∑ax m∑= ℓ m σ(r,𝜃,φ ) ≃ sℓm (r )Y ℓ (𝜃,φ), (8 ) ℓ=0m= −ℓ
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:
d2fℓm- 2dfℓm- ℓ(ℓ-+-1)fℓm- ∀(ℓ,m ), dr2 + r dr − r2 = sℓm (r). (9 )
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 (9View Equation) being, thus, a matrix Lℓm acting on the vector:
L ℓmA ℓm = Sℓm, (12 )
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 (7View Equation) is ill posed. One must indeed specify boundary conditions at r = Rmin and r = Rmax. For simplicity, let us suppose

∀(𝜃,φ ), ϕ(r = Rmin,𝜃, φ) = ϕ(r = Rmax, 𝜃,φ ) = 0. (13 )
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:
∑N ∑N ∀(ℓ,m ), (− 1)iaiℓm = aiℓm = 0. (14 ) i=0 i=0
Equations (14View Equation) are equivalent to boundary conditions (13View Equation), 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:

  1. 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).
  2. Get the values of the source σ on these grid points.
  3. Perform a spherical-harmonics transform (for example, using some available library [151Jump To The Next Citation Point]), followed by the Chebyshev transform (using a Fast Fourier Transform (FFT), or a Gauss–Lobatto quadrature) of the source σ.
  4. 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.
  5. 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. [109Jump To The Next Citation Point], 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.

  Go to previous page Go up Go to next page