  We describe an approach for finding upper bounds on an ODE dynamical system's
maximal Lyapunov exponent among all trajectories in a specified set. A
minimization problem is formulated whose infimum is equal to the maximal
Lyapunov exponent, provided that trajectories of interest remain in a compact
set. The minimization is over auxiliary functions that are defined on the state
space and subject to a pointwise inequality. In the polynomial case -- i.e.,
when the ODE's right-hand side is polynomial, the set of interest can be
specified by polynomial inequalities or equalities, and auxiliary functions are
sought among polynomials -- the minimization can be relaxed into a
computationally tractable polynomial optimization problem subject to
sum-of-squares constraints. Enlarging the spaces of polynomials over which
auxiliary functions are sought yields optimization problems of increasing
computational cost whose infima converge from above to the maximal Lyapunov
exponent, at least when the set of interest is compact. For illustration, we
carry out such polynomial optimization computations for two chaotic examples:
the Lorenz system and the H\'enon-Heiles system. The computed upper bounds
converge as polynomial degrees are raised, and in each example we obtain a
bound that is sharp to at least five digits. This sharpness is confirmed by
finding trajectories whose leading Lyapunov exponents approximately equal the
upper bounds.
