  We discuss several numerical methods for calculating Lyapunov exponents (a
quantitative measure of chaos) in systems of ordinary differential equations.
We pay particular attention to constrained systems, and we introduce a variety
of techniques to address the complications introduced by constraints. For all
cases considered, we develop both deviation vector methods, which follow the
time-evolution of the difference between two nearby trajectories, and Jacobian
methods, which use the Jacobian matrix to determine the true local behavior of
the system. We also assess the merits of the various methods, and discuss
assorted subtleties and potential sources of error.
