  We describe methods for proving bounds on infinite-time averages in
differential dynamical systems. The methods rely on the construction of
nonnegative polynomials with certain properties, similarly to the way nonlinear
stability can be proved using Lyapunov functions. Nonnegativity is enforced by
requiring the polynomials to be sums of squares, a condition which is then
formulated as a semidefinite program (SDP) that can be solved computationally.
Although such computations are subject to numerical error, we demonstrate two
ways to obtain rigorous results: using interval arithmetic to control the error
of an approximate SDP solution, and finding exact analytical solutions to
relatively small SDPs. Previous formulations are extended to allow for bounds
depending analytically on parametric variables. These methods are illustrated
using the Lorenz equations, a system with three state variables $(x,y,z)$ and
three parameters $(\beta,\sigma,r)$. Bounds are reported for infinite-time
averages of all eighteen moments $x^ly^mz^n$ up to quartic degree that are
symmetric under $(x,y)\mapsto(-x,-y)$. These bounds apply to all solutions
regardless of stability, including chaotic trajectories, periodic orbits, and
equilibrium points. The analytical approach yields two novel bounds that are
sharp: the mean of $z^3$ can be no larger than its value of $(r-1)^3$ at the
nonzero equilibria, and the mean of $xy^3$ must be nonnegative. The interval
arithmetic approach is applied at the standard chaotic parameters to bound
eleven average moments that all appear to be maximized on the shortest periodic
orbit. Our best upper bound on each such average exceeds its value on the
maximizing orbit by less than 1%. Many bounds reported here are much tighter
than would be possible without computer assistance.
