That is, as described by Hughes et al. [88], we numerically integrate the null geodesic equation
(where is an affine parameter) forwards in time from a set of starting events and check which events have “escaping” geodesics. For analytical or semianalytical studies like that of Bishop [31], this is an excellent algorithm.For numerical work it is straightforward to rewrite the null geodesic equation (8) as a coupled system of two firstorder equations, giving the time evolution of photon positions and 3momenta in terms of the 3 + 1 geometry variables , , , and their spatial derivatives. These can then be timeintegrated by standard numerical algorithms^{8}. However, in practice several factors complicate this algorithm.
We typically only know the 3 + 1 geometry variables on a discrete lattice of spacetime grid points, and we only know the 3 + 1 geometry variables themselves, not their spatial derivatives. Therefore we must numerically differentiate the field variables, then interpolate the field variables and their spacetime derivatives to each integration point along each null geodesic. This is straightforward to implement^{9}, but the numerical differentiation tends to amplify any numerical noise that may be present in the field variables.
Another complicating factor is that the numerical computations generally only span a finite region of spacetime, so it is not entirely obvious whether or not a given geodesic will eventually reach . However, if the final numericallygenerated slice contains an apparent horizon, we can use this as an approximation: Any geodesic which is inside this apparent horizon will definitely not reach , while any other geodesic may be assumed to eventually reach if its momentum is directed away from the apparent horizon. If the final slice (or at least its strongfield region) is approximately stationary, the error from this approximation should be small. I discuss this stationarity assumption further in Section 5.3.1.
In spherical symmetry this algorithm works well and has been used by a number of researchers. For example, Shapiro and Teukolsky [141, 142, 143, 144] used it to study event horizons in a variety of dynamical evolutions of spherically symmetric collapse systems. Figure 2 shows an example of the event and apparent horizons in one of these simulations.

In a nonsphericallysymmetric spacetime, several factors make this algorithm very inefficient:
Because of these limitations, for nonsphericallysymmetric spacetimes the “integrate null geodesics forwards” algorithm has generally been supplanted by the more efficient algorithms I describe in the following.
http://www.livingreviews.org/lrr20073 
This work is licensed under a Creative Commons AttributionNoncommercialNo Derivative Works 2.0 Germany License. Problems/comments to 