Go to previous page Go up Go to next page

8.4 Spectral integral-iteration algorithms

Nakamura, Kojima, and Oohara [113Jump To The Next Citation Point] developed a functional-iteration spectral algorithm for solving the apparent horizon equation (16View Equation).

This algorithm begins by choosing the usual polar-spherical topology for the angular coordinates (θ,φ ), and rewriting the apparent horizon equation (16View Equation) in the form

-∂θh- ∂-φφh- k Δh ≡ ∂θθh + tan θ + sin2θ = G (∂θφh,∂ φφh,∂θh,∂φh; gij,Kij, Γij), (21 )
where Δ is the flat-space angular Laplacian operator, and G is a complicated nonlinear algebraic function of the arguments shown, which remains regular even at θ=0 and θ= π.

Next we expand h in spherical harmonics (5View Equation). Because the left hand side of Equation (21View Equation) is just the flat-space angular Laplacian of h, which has the Yℓm as orthogonal eigenfunctions, multiplying both sides of Equation (21View Equation) by Y ∗ ℓm (the complex conjugate of Y ℓm) and integrating over all solid angles gives

∫ aℓm = − ---1--- Y ∗G dΩ (22 ) ℓ(ℓ+1 ) ℓm
for each ℓ and m except ℓ = m = 0.

Based on this, Nakamura, Kojima, and Oohara [113Jump To The Next Citation Point] proposed the following functional-iteration algorithm for solving Equation (21View Equation):

  1. Start with some (initial-guess) set of horizon-shape coefficients {aℓm}. These determine a surface shape via Equation (5View Equation).
  2. Interpolate the geometry variables to this surface shape (see Section 7.5).
  3. For each ℓ and m except ℓ = m = 0, evaluate the integral (22View Equation) by numerical quadrature to obtain a next-iteration coefficient aℓm.
  4. Determine a next-iteration coefficient a00 by numerically solving (finding a root of) the equation
    ∫ ∗ Y00G dΩ = 0. (23)
    This can be done using any of the 1-dimensional zero-finding algorithms discussed in Appendix A.
  5. Iterate until all the coefficients {aℓm} converge.

Gundlach [80Jump To The Next Citation Point] observed that the subtraction and inversion of the flat-space angular Laplacian operator in this algorithm is actually a standard technique for solving nonlinear elliptic PDEs by spectral methods. I discuss this observation and its implications further in Section 8.7.4.

Nakamura, Kojima, and Oohara [113] report that their algorithm works well, but Nakao (cited as personal communication in [146Jump To The Next Citation Point]) has argued that it tends to become inefficient (and possibly inaccurate) for large ℓ (high angular resolution) because the Y ℓm fail to be numerically orthogonal due to the finite resolution of the numerical grid. I know of no other published work addressing Nakao’s criticism.

8.4.1 Kemball and Bishop’s modifications of the Nakamura–Kojima–Oohara algorithm

Kemball and Bishop [93Jump To The Next Citation Point] investigated the behavior of the Nakamura–Kojima–Oohara algorithm and found that its (only) major weakness seems to be that the a00-update equation (23View Equation) “may have multiple roots or minima even in the presence of a single marginally outer trapped surface, and all should be tried for convergence”.

Kemball and Bishop [93Jump To The Next Citation Point] suggested and tested several modifications to improve the algorithm’s convergence behavior. They verified that (either in its original form or with their modifications) the algorithm’s rate of convergence (number of iterations to reach a given error level) is roughly independent of the degree ℓmax of spherical-harmonic expansion used. They also give an analysis that the algorithm’s cost is 𝒪 (ℓ4max), and its accuracy ɛ = 𝒪 (1∕ℓmax), i.e. the cost is 𝒪 (1∕ɛ4). This accuracy is surprisingly low: Exponential convergence with ℓmax is typical of spectral algorithms and would be expected here. I do not know of any published work which addresses this discrepancy.

8.4.2 Lin and Novak’s variant of the Nakamura–Kojima–Oohara algorithm

Lin and Novak [104Jump To The Next Citation Point] have developed a variant of the Nakamura–Kojima–Oohara algorithm which avoids the need for a separate search for a 00 at each iteration: Write the apparent horizon equation (16View Equation) in the form

Δh − 2h = λΘ + Δh − 2h, (24 )
where Δ is again the flat-space angular Laplacian operator and where λ is a nonzero scalar function on the horizon surface. Then choose λ as
( )1 ∕3 det-gij [ mn ¯ ¯ ]1∕2 2 λ = det fij g (∇mF )(∇nF ) h , (25 )
where fij is the flat metric of polar spherical coordinates, ∇¯ is the associated 3-covariant derivative operator, and F is the level set function (7View Equation).

Lin and Novak [104Jump To The Next Citation Point] showed that all the spherical-harmonic coefficients aℓm (including a00) can then be found by iteratively solving the equation

∫ ----1------ ∗ a ℓm = − ℓ(ℓ+1) + 2 Y ℓm λ(Θ + Δh − 2h)dΩ. (26 )

Lin and Novak [104Jump To The Next Citation Point] find that this algorithm gives robust convergence and is quite fast, particularly at modest accuracy levels. For example, running on a 2 GHz processor, their implementation takes 3.1, 5.8, 17, 88, and 313 seconds to find the apparent horizon in a test slice to a relative error (measured in the horizon area) of −4 9 × 10, −5 5 × 10, −7 6 × 10, −10 9 × 10, and −12 3 × 10 respectively37. This implementation is freely available as part of the External LinkLorene toolkit for spectral computations in numerical relativity (see Table 2).

8.4.3 Summary of spectral integral-iteration algorithms

Despite what appears to be fairly good numerical behavior and reasonable ease of implementation, the original Nakamura–Kojima–Oohara algorithm has not been widely used apart from later work by its original developers (see, for example, [115114]). Kemball and Bishop [93] have proposed and tested several modifications to the basic Nakamura–Kojima–Oohara algorithm. Lin and Novak [104] have develped a variant of the Nakamura–Kojima–Oohara algorithm which avoids the need for a separate search for the a00 coefficient at each iteration. Their implementation of this variant is freely available as part of the External LinkLorene toolkit for spectral computations in numerical relativity (see Table 2).


  Go to previous page Go up Go to next page