# "Constrained" numerical solutions of ODEs with conservation laws?

I know little about numerical methods and I was considering the following problem that possibly has standard solution in the literature.

Suppose you have an ODE for wich we already know that it must satisfy a conservation law for some "energy" (for example hamiltonian dynamical systems such as $\frac 1 2 m \dot{x}^2+U(x)=\mbox{const}$). Then if we try to solve it with standard Runge Kutta methods we are likely to have errors for the energy at each step, which possibly will sum up. But what if I want to sacrifice a little accuracy for the "position" in order to have a greater accuracy for the energy? (This could be relevant if for example I am interested in having good approximation of long time asimptotic behaviors which should not be affected by errors in short time dynamics). Are there numerical methods which tries to give the best approximation given the constrain of keeping the energy "strictly" fixed?

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:

1. Compute acceleration $\vec{a}$ as a function of position $\vec{x}$,
2. Update velocity $\vec{v}$ by adding $\vec{a}\,\Delta t$,
3. Update position $\vec{x}$ by adding $\vec{v}\,\Delta t$,
4. 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.

• Also see here for more detailed answers. Aug 21, 2015 at 12:17

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.