Many physical systems of current interest are chaotic, which means that numerical errors in their simulation are exponentially magnified with the passage of time. This could mean that a numerical solution of a chaotic system is the result of nothing but magnified noise, which calls into question the value of such simulations. Although this fact has been well known for a long time, its impact on the validity of simulations is not well understood. The study of shadowing may provide an answer. A shadow is an exact trajectory of a chaotic map or ordinary differential equation that remains close to an approximate solution for a nontrivial duration of time. If it can be shown that a numerical solution has a shadow, then the validity of the solution is strong, in the sense that it can be viewed as an experimental observation of the shadow, which is an exact solution. We present a discussion of shadowing, including an algorithm to find shadows, using the gravitational N-body problem as an example.

1.
This quotation is frequently attributed to Mark Twain, but in the Chronicle of Higher Education (4 September 1991, p. A8), the group at the University of California, Berkeley (〈http://bancroft.berkeley.edu/MTP/〉), which is editing all of Twain’s papers, put this one into the list of quotations that Twain did not say.
2.
D. A. Clarke and M. J. West, in The 12th Kingston Meeting: Computational Astrophysics (ASP Conference Series, Astron. Soc. of the Pacific, 1997), Vol. 123.
3.
E.
Bertschinger
, “
Simulations of structure formation in the universe
,”
Annu. Rev. Astron. Astrophys.
36
,
599
654
(
1998
).
4.
J.
Goodman
,
D. C.
Heggie
, and
P.
Hut
, “
On the exponential instability of N-body systems
,”
Astrophys. J.
415
,
715
733
(
1993
).
5.
R. H.
Miller
, “
Irreversibility in small stellar dynamical systems
,”
Astrophys. J.
140
,
250
256
(
1964
).
6.
M.
Lecar
, “
A comparison of eleven numerical integrations of the same gravitational 25-body problem
,”
Bull. Astron.
3
,
91
(
1968
).
7.
H. E.
Kandrup
and
H.
Smith
, “
On the sensitivity of the N-body problem to small changes in initial conditions
,”
Astrophys. J.
374
,
255
265
(
1991
).
8.
H. E.
Kandrup
,
H.
Smith
, and
D.
Willmes
, “
On the sensitivity of the N-body problem to small changes in initial conditions. III
,”
Astrophys. J.
399
,
627
633
(
1992
).
9.
L.
Hernquist
,
P.
Hut
, and
J.
Makino
, “
Discreteness noise versus force errors in N-body simulations
,”
Astrophys. J. Lett.
402
,
L85
L88
(
1993
).
10.
J. E.
Barnes
and
P.
Hut
, “
Error analysis of a tree code
,”
Astrophys. J., Suppl. Ser.
70
,
389
417
(
1989
).
11.
J. A.
Sellwood
, “
The art of N-body building
,”
Annu. Rev. Astron. Astrophys.
25
,
151
-
86
(
1987
).
12.
J. P.
Singh
,
J. L.
Hennessy
, and
A.
Gupta
, “
Implications of hierarchical N-body methods for multiprocessor architecture,” Technical Report
CSL-TR-92-506, Computer Systems Lab, Stanford University,
Stanford, CA
1992
.
13.
D.
Pfenniger
and
D.
Friedli
, “
Computational issues connected with 3D N-body simulations
,”
Astron. Astrophys.
270
,
561
572
(
1993
).
14.
J. G.
Jernigan
and
D. H.
Porter
, “
A tree code with logarithmic reduction of force terms, hierarchical regularization of all variables, and explicit accuracy controls
,”
Astrophys. J., Suppl. Ser.
71
,
871
893
(
1989
).
15.
J.
Barnes
and
P.
Hut
, “
A hierarchical O(N log N) force-calculation algorithm
,”
Nature (London)
324
,
446
449
(
1986
).
16.
J. Binney and S. Tremaine, Galactic Dynamics (Princeton University Press, Princeton, NJ, 1987).
17.
L.
Greengard
and
V.
Rokhlin
, “
A fast algorithm for particle simulation
,”
J. Comput. Phys.
73
,
325
348
(
1987
).
18.
P. J.
Channell
and
C.
Scovel
, “
Symplectic integration of Hamiltonian Systems
,”
Nonlinearity
3
,
231
259
(
1990
).
19.
J. M. Sanz-Serna, “Symplectic Integrators for Hamiltonian Problems: An Overview,” Acta Numerica 1992, edited by A. Iserles (Cambridge University Press, Cambridge, 1992), pp. 243–286.
20.
G. D.
Quinlan
and
S.
Tremaine
, “
On the reliability of gravitational N-body integrations
,”
Mon. Not. R. Astron. Soc.
259
,
505
518
(
1992
).
21.
This is in stark contrast to N-body simulations of our solar system, in which case we are interested in the precise evolution of the particular solar system we live in, not a hypothetical one very similar to it.
22.
Speech at Bristol on Declining the Poll Vol. ii, p. 420. (1780) from Macmillan Book of Proverbs, Maxims, and Famous Phrases (Macmillan, New York, 1965), p. 2080 (quotes about shadows).
23.
C.
Grebogi
,
S. M.
Hammel
,
J. A.
Yorke
, and
T.
Sauer
, “
Shadowing of physical trajectories in chaotic dynamics: Containment and refinement
,”
Phys. Rev. Lett.
65
,
1527
1530
(
1990
).
24.
D. V.
Anosov
, “
Geodesic flows and closed Riemannian manifolds with negative curvature
,”
Proc. Steklov Inst. Math.
90
,
1
235
(
1967
).
25.
R.
Bowen
, “
ω-limit sets for axiom A diffeomorphisms
,”
J. Diff. Eqns.
18
,
333
339
(
1975
).
26.
W. B.
Hayes
and
K. R.
Jackson
, “
Rigorous shadowing of numerical solutions of ordinary differential equations by containment
,”
SIAM (Soc. Ind. Appl. Math.) J. Numer. Anal.
41
,
1948
1973
(
2003
).
27.
W. B. Hayes, “Rigorous shadowing of numerical solutions of ordinary differential equations by containment,” Ph.D. thesis, Department of Computer Science, University of Toronto, 2001. Available at 〈http://www.cs.toronto.edu/NA/reports.html#hayes-01-phd〉.
28.
In other words, let y=f(y(t)) be a first-order ordinary differential equation. Note that yi+1=φ(yi) is the solution of Eq. (3) using yi as the initial condition and integrating f to time ti+1. The Jacobian Df(yi) measures how y changes if y is changed by a small amount. The resolventR(ti+1,ti) is the integral of Df(y) along the path y(t), and describes how a small perturbation δy of yi at time ti is mapped to a perturbation of yi+1 at time ti+1.R(ti+1,ti) is the solution of the variational equation,
∂R∂t=Df(y(t))R(t,ti), R(ti,ti)=I,
where I is the identity matrix. The reason the arguments to R seem to be reversed is for notational convenience: they satisfy the identity R(t2,t0)=R(t2,t1)R(t1,t0), and so a perturbation δy at time t0 gets mapped to a perturbation at time t2 by the matrix–matrix and matrix–vector multiplication R2δy=R1R0δy.29 Finally, the linear map in the GHYS refinement procedure, if φ is the time-h solution operator for Eq. (3), is Dφ(yi)=R(ti+1,ti).
29.
E. Hairer, S. P. Nørsett, and Gerhard Wanner, Solving Ordinary Differential Equations, 2nd ed. (Springer-Verlag, Berlin, 1993), Vols. 1 and 2.
30.
B. W. Char, First Leaves: A Tutorial Introduction to Maple V (Springer-Verlag, Berlin, 1993).
31.
J. D.
Farmer
and
J. J.
Sidorowich
, “
Optimal shadowing and noise reduction
,”
Physica D
47
,
373
392
(
1991
).
32.
S. T.
Fryska
and
M. A.
Zohdy
, “
Computer dynamics and the shadowing of chaotic orbits
,”
Phys. Lett. A
166
,
340
346
(
1992
).
33.
R. M.
Corless
, “
What good are numerical simulations of chaotic dynamical systems?
,”
Comput. Math. Appl.
28
,
107
121
(
1994
).
34.
S. Guazzo, Civile Conversation, Book ii (1574), p. 231.
35.
W. Hayes, “Efficient shadowing of high dimensional chaotic systems with the large astrophysical n-body problem as an example,” Master’s thesis, Department of Computer Science, University of Toronto, January 1995.
36.
W.
Hayes
, “
Shadowing high-dimensional Hamiltonian systems: The gravitational n-body problem
,”
Phys. Rev. Lett.
90
,
054104
-
1
054104
-
4
(
2003
).
37.
J. A. Rice, Mathematical Statistics and Data Analysis (Wadsworth & Brooks/Cole, 1988).
38.
W.
Hayes
, “
Shadowing-based reliability decay in softened n-body simulations
,”
Astrophys. J. Lett.
587
,
L59
L62
(
2003
).
39.
La Fontaine, Fables, Book vi, Fable 17, 1668.
40.
K. T. Alligood, T. D. Sauer, and J. A. Yorke, Chaos: An Introduction to Dynamical Systems (Springer-Verlag, Berlin, 1996).
This content is only available via PDF.
AAPT members receive access to the American Journal of Physics and The Physics Teacher as a member benefit. To learn more about this member benefit and becoming an AAPT member, visit the Joining AAPT page.