"Constrained" numerical solutions of ODEs with conservation laws?
The short answer is yes. Fun fact is that such methods can even be easy to implement, or, in some sense, they look more natural than methods like Heun or Runge-Kutta.
When I was a kid, I programmed a simulation of the motion of a satellite under the influence of a planet's gravity, as part of a game. The scheme was straightforward:
- Compute acceleration $\vec{a}$ as a function of position $\vec{x}$,
- Update velocity $\vec{v}$ by adding $\vec{a}\,\Delta t$,
- Update position $\vec{x}$ by adding $\vec{v}\,\Delta t$,
- Repeat.
That method created a seemingly stable elliptical orbit. With a rather coarse $\Delta t$, the ellipse would rotate noticably, but it would not blow up nor collapse. In fact, I did not notice that integrating motions could generally lead to systematically biased energy conservation errors until I had to study the subject and thus got to know quite some schemes, all of which seemed more complicated in terms of operations per step, yet none seemed appropriate for simulating celestial mechanics because of the energy drift.
Later I learned that I had unknowingly implemented the leapfrog method. I have also heard the method referenced in a talk about the millenium simulation. A general framework is given with so-called symplectic integrators.
While energy conservation is important, specifically in Hamiltonian systems, it is not enough. Rather, it is important to preserve the symplectic structure, which means that your initial conditions should be mapped forward in time using only canonical transformations. That is where the symplectic integrators come in, as mentioned by ccorn.