  In this paper we consider discrete gradient methods for approximating the
solution and preserving a first integral (also called a constant of motion) of
autonomous ordinary differential equations. We prove under mild conditions for
a large class of discrete gradient methods that the numerical solution exists
and is locally unique, and that for arbitrary $p \in \mathbb{N}$ we may
construct a method that is of order $p$. In the proofs of these results we also
show that the constants in the time step constraint and the error bounds may be
chosen independently from the distance to critical points of the first
integral. In the case when the first integral is quadratic, for arbitrary $p
\in \mathbb{N}$, we have devised a new method that is linearly implicit at each
time step and of order $p$. This new method has significant advantages in terms
of efficiency. We illustrate our theory with a numerical example.
