As described in the previous chapter, Grebogi, Hammel, Yorke, and Sauer (GHYS) , invented a two-dimensional shadowing procedure consisting of two parts:
Quinlan and Tremaine (QT) , wanted to test for the existence of shadows in large N-body systems. They generalized the refinement portion of the GHYS algorithm to work on arbitrary-dimensional Hamiltonian systems. QT showed that, using 100 decimal digits of precision, if the noise level of a four-dimensional orbit could be reduced to about , then the noise could also be reduced all the way down to , with the number of correct digits in the solution increasing geometrically. They conclude that it is reasonable to assume that further refinement would continue to reduce the noise indefinitely, leading in the limit to a true orbit. This seems less rigourous that containment, but in practice it is probably almost as reliable.
QT then applied their refinement procedure to a simplified N-body system in which a single particle moves, unsoftened, among 100 fixed particles that are distributed in a standard Plummer distribution. For their ``cheap'' integrator, they used a 5th-order predictor-corrector routine that is commonly used by astronomers [Appendix B of ], arranged to have 1-step errors of about . For their expensive integrator, they used a Bulirsch-Stoer method  with tolerance . They found they could shadow the particle for several crossing times. Unfortunately, their refinement procedure was too inefficient to work on significantly higher dimensional systems. They tried two and three moving particles, but found they were pushing the bounds of feasibility that their shadowing procedure could handle.
It is perhaps important to note before continuing that it is not clear that refinement alone, even if optimized as introduced in this chapter, and even if it is as reliable as containment, is more efficient than using containment. But if refinement is to be used at all, it needs to be made more efficient than the direct, naïve method of GHYS and QT.
It is not at all clear that being able to shadow one particle moving ``pinball'' fashion amongst 100 fixed ones implies that shadows of large N-body systems exist. To establish the shadowability of large N-body simulations requires shadowing tests to be applied to more realistic systems. At present it is far beyond feasibility to attempt shadowing systems of or particles that are commonly being used today by astronomers. Nevertheless, we can learn much about how shadowing results change as the dimension increases, using systems with just tens of moving particles. In such systems it may be reasonable to eliminate the fixed particles entirely, but to extend QT's results, most of the experiments performed in this thesis (and reported in the next chapter) were done on a system with a total of 100 particles, M of which move, for M ranging from 1 to 25.
It is important to emphasize the main scientific question that we will be studying once we have a reasonably fast shadowing implementation. The question is not only a question concerning N-body systems, but a general question for chaotic systems: how do shadowing results change as the number of dimensions increases?
For the remainder of this thesis, the following definitions will be in effect:
As has already been mentioned, most of the expense of refinement of an N-body orbit is in computing the linear maps . See a previous footnote for an explanation of , and the Appendix for the derivation of the Jacobian of an N-body system. Manually counting floating-point operations in the Jacobian gives an operation count of roughly . The total operation count to compute the right hand side of the variational equation, including the matrix-matrix multiply and the computation of the Jacobian, is about . The term is probably negligible, and the MN term is only relevant if N >> M. Thus, if there are S shadow steps, and the function computing the right hand side is invoked an average of k times per shadow step, the total cost of computing all the resolvents is about . I have found empirically that, using an Adams's integrator, k seems to be independent of M, and roughly equal to about 500 for the shadow steps I use in my N-body simulations. Figure 2.0 shows some experiments comparing the time to compute one resolvent to M, demonstrating the relationship.
Figure 2.0: Time in seconds to compute one resolvent as a function of the number of moving particles M. For comparison, is also plotted.