  We design a novel, exactly energy-conserving implicit non-symplectic
integration method for an eight-dimensional Hamiltonian system with four
degrees of freedom. In our algorithm, each partial derivative of the
Hamiltonian with respect to one of phase-space variables is discretized by the
average of eight Hamiltonian difference terms. Such a discretization form is a
second-order approximation to the Hamiltonian gradient. It is shown numerically
via simulations of an FPU-$\beta$ system and a post-Newtonian conservative
system of compact binaries with one body spinning that the newly proposed
method has extremely good energy-conserving performance, compared to the
Runge-Kutta, implicit midpoint symplectic method and extended phase-space
explicit symplectic-like integrators. The new method is advantageous over very
long times and for large time steps compared to state-of-the-art Runge-Kutta
method in the accuracy of numerical solutions. Although such an
energy-conserving integrator exhibits a higher computational cost than any one
of the other three algorithms, the superior results justify its use for
satisfying some specific purposes on the preservation of energies in numerical
simulations with much longer times, e.g. obtaining a high enough accuracy of
the semimajor axis in a Keplerian problem in the Solar System or accurately
grasping the frequency of a gravitational wave from a circular orbit in a
post-Newtonian system of compact binaries. The new integrator will be
potentially applied to model time-varying external electromagnetic fields or
time-dependent spacetimes.
