Conservertion of Energy

Required Reading


As we integrate ODEs which describe physical systems, we expect to conserve energy. As integration progresses, the total energy in the system may change. When this happens, it's likely that the integration scheme contains errors.


To demonstrate, take the simple case of a satellite orbiting a planet. Newton's law of gravity describes the force experienced by both the satellite and the planet:

\begin{equation} \vec{F} = -\frac{m_p m_s G}{|\vec{r}|^3}\vec{r} \label{newtongravity} \end{equation}.

Newton's second law gives the force felt by mass as a function of its acceleration as

\begin{equation} \eqalign{ \vec{F} & = m\vec{a} \cr & = m\frac{d^2\vec{r}}{dt^2} } \label{newton2} \end{equation}.

This gives us, for the satellite:

\begin{equation} -\frac{m_p G}{|\vec{r}|^3}\vec{r} = \frac{d^2\vec{r}}{dt^2} \label{gravity} \end{equation}.

We assume $m_p >> m_s$ such that the force exerted on the planet by the satellite is negligible.

Equation \eqref{gravity} is a second order differential equation. In order to obtain a system of first order differential equations, introduce a new function

\begin{equation} \vec{v} = \frac{d\vec{r}}{dt} \label{velocity} \end{equation},

which you may recognise as the velocity.

Equation \eqref{gravity} can now be rewritten as a system of first order equations:

\begin{equation} \eqalign{ \frac{d\vec{v}}{dt} &= -\frac{m_p G}{|\vec{r}|^3}\vec{r} \cr \frac{d\vec{r}}{dt} &= \vec{v} } \label{gravity2} \end{equation}.

Note that $\frac{d\vec{v}}{dt}$ only depends upon $\vec{r}$ and $\frac{d\vec{r}}{dt}$ only depends upon $\vec{v}$. To solve these equations of motion for the satellite, integrate according to the steps described by the Leapfrog method.

First advance $r$ by half a step to get $r_{1/2}$:

\begin{equation} \vec{r}_{1/2} = \vec{r}_0 + \frac{\Delta t}{2} \vec{v}_0 \label{leapfrog_t0} \end{equation}.

Then, for $i \in [0, 1, 2 \dots]$:

\begin{equation} \eqalign{ \vec{v}_{i+1} &= \vec{v}_i - \Delta t \frac{m_p G}{|\vec{r}_{i + 1/2}|^3} \vec{r}_{i+1/2} \cr \vec{r}_{i+3/2} &= \vec{r}_{i+1/2} + \Delta t \cdot \vec{v}_i } \label{leapfrog_t1} \end{equation}.

Conservation of Energy

An animation demonstrating the conservation properties of the Leapfrog integrator.

To show that the Leapfrog integrator conserves energy, consider what happens when the integration is reversed, i.e. from $\vec{r}_{i+3/2}, \vec{v}_{i+1}$ to $\vec{r}_{i+1/2}, \vec{v}_{i}$,

\begin{equation} \eqalign{ \vec{r}_{i+1/2} &= \vec{r}_{i+3/2} - \Delta t \cdot \vec{v}_{i+1} \cr \vec{v}_i &= \vec{v}_{i+1} + \Delta t \frac{m_p G}{|\vec{r}_{i + 1/2}|^3}\vec{r}_{i+1/2} } \label{reverse} \end{equation}.

Note that this is the same process as in \eqref{leapfrog_t1}, only with the signs of $v$, $r$ and $\Delta t$ reversed. It follows then that by integrating from $t_0$ to $t_e$, and then by integrating in reverse from $t_e$ back to $t'_0 = t_0$

\begin{equation} \eqalign{ \vec{r}'_{1/2} &= \vec{r}_{1/2} \cr \vec{v}'_{0} &= \vec{v}_{0} } \label{tprime} \end{equation}.

Suppose now that the initial integration from $t_0$ to $t_e$ results in an energy error of $\varepsilon$, and the error from $t_e$ to $t'_0$ is $\varepsilon'$. If $\vec{v}'$ and $\vec{r}'$ are exactly equal to $\vec{v}$ and $\vec{r}$ respectively, then $\varepsilon + \varepsilon' = 0$. However, because $\vec{r}$ and $\vec{v}$ do not depend on $t$, the reverse integration is equivalent to the forward integration of the same magnitude; therefore $\varepsilon' = \varepsilon$. It follows that if $2\varepsilon = 0$, then $\varepsilon = 0$. The Leapfrog integrator conserves energy.