<p>I like to classify solutions of the problem of the time evolution of the
complete initial state of a set of objects at some epoch time, where the
objects are subject to Newtonian gravitation into two main groups.
One approach is to use orbital elements of some sort.
The other is to use a numerical initial value problem solver,
aka a numerical integrator. The latter is the primary subject of this answer.</p>

<p>Note well: This classification isn't quite perfect as hybrid approaches are
also possible, wherein one uses a numerical integrator to integrate
time-varying orbital elements.</p>

<p><br>
Keplerian orbital elements work quite nicely in the case of two point
masses, or more generally, two objects with a spherical mass distribution.
The anomaly is the only Keplerian orbital element that changes over time.
Keplerian elements can be used in situations where the underlying
assumptions are approximately correct by developing a model of how
those supposedly unchanging elements vary with time. One way of doing
this is to use Lagrange's planetary equations. (There are other related
approaches such as Gauss' planetary equations, Delaunay's planetary
equations, etc.)
Lagrange's planetary equations yield expressions for how Keplerian orbital
elements vary over time given a set of perturbing forces.</p>

<p>Another approach is to use something akin to those Keplerian elements
such as Delaunay elements), coupled with planetary equations for those
alternative elements. Yet another approach is to use orbital elements
(e.g., Brouwer-Lyddane elements, SGP4 elements) in which the planetary
equations are embedded in the orbital element to Cartesian state
transformation algorithm. This final approach is used to this day to
describe vehicles in Earth orbit.</p>

<p><br>
The other approach is to use numerical integration.</p>

<p>I'll start with a discussion of how to solve for the value of a scalar
function $x(t)$ at some time $t_1$ given an initial value $x(t_0) = x_0$
and some well-behaved (continuous and bounded) derivative function
$f(x(t),t) = dx(t)/dt$ that describes the time evolution of $x$.
This falls in the very broad category of initial value problems.</p>

<p>Suppose the ordinary differential equation cannot be solved analytically
and cannot be expressed terms of a useful power series. This doesn't mean
nothing can be done. There are a number of techniques for solving this
problem numerically.</p>

<p>Note the dependence of the derivative function on the dependent variable $x$.
This becomes the much simpler problem of numerical quadrature if the
derivative function can be expressed independent of $x$.
The discussion that follows assumes that the derivative function $f$ does
indeed depend on $x$.
Note well: Newtonian gravitation falls in this category.</p>

<p>It also assumes the derivative function is well-behaved.
Numerically integrating across a discontinuity is a bad idea.</p>

<p>The foundation of the integration-based techniques for a scalar function
is the mean value theorem, which says that at some time $t_c$ between
$t_0$ and $t_1$, the value at $t_1$ is exactly
$x(t_1) = x(t_0) + (t_1-t_0)\,f(x(t_c),t_c)$. If only we could find that
magical $t_c$ and the derivative at that point. There's a chicken and
egg problem here: that magical point in time is not known. Even if it was,
the derivative function depends on state, and that too isn't known.</p>

<p>A very simple approach around this problem is to assume that this magical
point is the initial point:
$$x(t_1) = x(t_0) + (t_1-t_0)f(x(t_0),t0)$$
This works quite nicely for values of $t_1$ that are very close to $t_0$.
It doesn't work very well at all where $|t_1 - t_0|$ is not small.
This suggests splitting the interval $(t_0, t_1)$ into a number of smaller
intervals. This results in Euler's method: Apply the above to advance
state to time $t_0+\Delta t$, then to $t_0+2\Delta t$, and so on, eventually
reaching the desired time. Euler's method is rather lousy, even for a simple
first order scalar ODE. We can do much better than this. The key reason for
discussing Euler's method is that it is the basis for many other integration
techniques. Learn how it works, then toss it.</p>

<p>One approach to approving on Euler's method is to somehow correct the
result from Euler's method. For example, take an Euler step and compute
the derivative at the end point. Then use the average of those two
derivative values (the original value used to make the Euler step, and
the other from the end of the Euler step) to recompute the step from
$t$ to $t+\Delta t$. THis is Heun's method.</p>

<p>Another approach is to guess that the magical point $t_c$ lies
somewhere between $t$ and $t+\Delta t$. Perhaps the middle? We can use
Euler's method to advance state to the midpoint, and then use the
derivative at that point to advance state from $t$ to $t+\Delta t$.
This is the midpoint method.</p>

<p>Both Heun's method and the midpoint method appear to be steps backwards,
computationally. While Euler's method requires but one evaluation of the
derivative function per time step, these improved methods require two.
However, the error growth is in general so much smaller with either
Heun's method or the midpoint method compared to that from Euler's method.
This means that those "improvements" most definitely are improvements.
The expense of calling the derivative function twice per step is more than
offset by the fact that imprpovements enables take steps that are orders
of magnitude larger than one can make with Euler's method.</p>

<p>Both Heun's method and the midpoint method are simple improvements.
This problem has been studied in many guises. There are many more advanced
techniques. One is the class of Runge-Kutta integration techniques.
Both Heun's method and the midpoint method fall into this class.
The most popular of these, classical Runge-Kutta 4, is a significant
improvement on those two methods. There are even higher order Runge-Kutta
integrators than RK4.
Heun's method also falls into the broad class of predictor-correctors,
wherein one method (the predictor) advances state to the end of the interval
and another method (the corrector) uses the derivative at this approximate
endpoint to correct the guess made by the predictor.</p>

<p><br></p>

<p>The above focused solely on first order ODEs involving a scalar function.
What if the problem is multidimensional or involves higher order derivatives?
The mathematics of the techniques described above can easily handle
multidimensional data: Simply use the vector-valued time derivative.
Since a higher order ODE can be converted to a first order ODE via an
augmented vector-valued state, the same approaches used to address
multidimensional data can also be employed to address higher order ODEs.</p>

<p>There's a problem with doing this: It throws out geometry. For example,
consider the rather simple first order ODE $\dot x = -y, \dot y = x$.
The solution to this multivariate ODE is uniform circular motion.
Applying Euler's method to this results in
$$\begin{aligned}
x(t+\Delta t) &amp;= x(t) - \Delta t\,y(t) \\
y(t+\Delta t) &amp;= y(t) + \Delta t\,x(t)
\end{aligned}$$
The square magnitude of this new vector is
$(x(t)^2+y(t)^2)(1+\Delta t^2)$, which is always greater than magnitude
of the vector at the start of the step. This is not uniform circular motion.
The solution obtained via Euler's method instead spirals out. Other techniques
spiral inward. A geometric integrator on the other hand will somehow
maintain the constraint that $x^2+y^2$ is a constant of motion.</p>

<p>The above example showed why we don't want to toss geometry in a very
simple problem. The geometry of Newtonian gravitation, along with much
classical mechanics in general, is symplectic geometry.
This is why symplectic integrators are of great concern. A simple example
again, with Euler's method: Suppose the second derivative of position
is given by some function $\ddot {\vec x}(t) = \vec f(x(t),t)$.
Applying the basic Euler method against 3+3 dimensional phase space
dictates that
$$\begin{aligned}
\vec x(t+\Delta t) &amp;= \vec x(t) + \Delta t \vec v(t) \\
\vec v(t+\Delta t) &amp;= \vec v(t) + \Delta t \vec f(\vec x(t),t)
\end{aligned}$$
A simple change makes this symplectic:
$$\begin{aligned}
\vec v(t+\Delta t) &amp;= \vec v(t) + \Delta t \vec f(\vec x(t),t) \\
\vec x(t+\Delta t) &amp;= \vec x(t) + \Delta t \vec v(t+\Delta t)
\end{aligned}$$
As is the case with the scalar techniques discussed at the start,
Euler's method is a starting point rather than the end with
regard to symplectic integration techniques. Symplectic Euler's method
is rather lousy. But at least orbits don't spiral outward.</p>

<p>In making an N-body gravitation simulation, the size of N (the number of
bodies) is a key concern. Simulating a galaxy is a very different concern
from simulating a star system. The techniques used in simulating the
formation of a galaxy are very different from those used to develop a
solar system ephemeris. Galactic scale simulations cannot afford to
calculate all of the N<sup>2</sup> gravitational interactions amongst all
the particles, and because N is so large, it cannot afford anything
more complex than very simple integrators. A star system model that does
not calculate all N<sup>2</sup> of the gravitational interactions
or that uses very a simple integrator will by viewed in disdain.</p>

CC BY-SA 3.0
2017-05-19T19:50:12.943Z
2017-05-22T10:56:08.637Z
2017-05-22T10:56:08.637Z
David Hammen
<p>I am making a small orbit simulator. I figured out Kepler's laws and know how to work with them to "update" an orbit based on time. So now I want to know how I can describe orbits with newton laws. Given a position vector and a velocity vector of a satellite, What I currently do is this:</p>

<ol>
<li>Calculate gravitational acceleration: (g*mass)/distance^2</li>
<li>Add this number to the velocity vector.</li>
<li>Multiply the new velocity vector with the DeltaTime</li>
<li>Add it to the position vector</li>
<li>Go back to 1.</li>
</ol>

<p>However when converting this to Keplerian Orbit Elements doesn't really make them constant.(Calculations from here: <a href="http://ccar.colorado.edu/asen5070/handouts/cart2kep2002.pdf" rel="nofollow noreferrer">Cartesian To Kepler</a>)</p>

<p>Also the orbit isn't really stable, and also wobbles around. So can anyone help me fix this? Or am I doing something wrong?</p>

CC BY-SA 3.0
2017-05-18T17:56:02.243Z
This might be too broad. The question on "how to integrate a second order ode" could easily fill several books. However, with small step size there shouldn't be too much wobbling, so there may be bugs in your code. If not then investigate verlet integration, and Runge-Kutta methods. There are several existing implementations that can be searched for and compared with yours.
James K
CC BY-SA 3.0
2017-05-18T18:01:47.527Z
You first want to multiply your gravity vector with dt and then add this number to the velocity vector. Otherwise you're not solving Newton's ODE.
AtmosphericPrisonEscape
CC BY-SA 3.0
2017-05-18T18:05:00.410Z
But that leaves me with a huge number as the velocity vector wich I add to position
GR00G0
CC BY-SA 3.0
2017-05-18T18:45:00.943Z
We might solve your problem easier if you'd just post which discretization you tell the computer to do.
AtmosphericPrisonEscape
CC BY-SA 3.0
2017-05-18T18:49:20.273Z
I am not sure what you mean
GR00G0
CC BY-SA 3.0
2017-05-18T19:01:30.353Z
Write down the formula you hack into the computer.
AtmosphericPrisonEscape
CC BY-SA 3.0
2017-05-19T15:54:22.940Z
The process you describe is called Euler integration. It's the most basic numerical  integration method, and notoriously bad. Even if you make the time step tiny it's performance is poor, especially if the orbit is eccentric.
PM 2Ring
CC BY-SA 3.0
2017-05-19T16:06:31.580Z
So you need a better integrator, preferably one that's symplectiic (conserves energy) like leapfrog or verlet.  Lots of people use and recommend Runge-Kutta methods, but they aren't so suitable for pure gravity sims because they aren't symplectic.
PM 2Ring
CC BY-SA 3.0
2017-05-19T16:13:56.150Z
If you're getting unexpectedly huge numbers then you may have a mistake somewhere with your units, eg using seconds where you should have days or years.  And with Euler integration you do need tiny time steps.  IIRC, doing Mars with Euler the orbit won't even close up unless delta t is around a day or smaller.
PM 2Ring
CC BY-SA 3.0
2017-05-19T20:01:12.010Z
Search stackoverflow for "numerically solving kepler's equation". One of the solutions I found helpful for this is: http://stackoverflow.com/questions/5287814/solving-keplers-equation-computationally. The js program at the bottom is a very nice step through to solving the equations using an iterative method.
scrappedcola
CC BY-SA 3.0
2017-05-19T20:29:06.293Z
So what I could do is during physics time(normal delta time where everything happens at realtime, I can use semi-implict euler, and as soon as I switch to faster time, i convert position and velocity to kepler elements and calculate it with kepler laws? It makes sense that larger deltatime values result in less accurate calculations, but could that be the reason that the orbit slowly falls down?
GR00G0
CC BY-SA 3.0
2017-05-19T20:39:17.197Z
Also, Ido know how to solve in keplerian elements, i just dont want to use it while in physics mode, so I can simulate effects of a rocket burn better, that by converting between kepler and cartesian
GR00G0
<orbit><orbital-elements><orbital-mechanics><newtonian-gravity>
Orbits using Newtons laws
GR00G0
# Orbits using Newtons laws

(Asked by: username_0 on 2017-05-18)

I am making a small orbit simulator. I figured out Kepler's laws and know how to work with them to "update" an orbit based on time. So now I want to know how I can describe orbits with newton laws. Given a position vector and a velocity vector of a satellite, What I currently do is this:

1. Calculate gravitational acceleration: (g\*mass)/distance^2
2. Add this number to the velocity vector.
3. Multiply the new velocity vector with the DeltaTime
4. Add it to the position vector
5. Go back to 1.

However when converting this to Keplerian Orbit Elements doesn't really make them constant.(Calculations from here: [Cartesian To Kepler](http://ccar.colorado.edu/asen5070/handouts/cart2kep2002.pdf))

Also the orbit isn't really stable, and also wobbles around. So can anyone help me fix this? Or am I doing something wrong?

**Tags:** orbit, orbital-elements, orbital-mechanics, newtonian-gravity

> **Comments:**
> **username_1**: This might be too broad. The question on "how to integrate a second order ode" could easily fill several books. However, with small step size there shouldn't be too much wobbling, so there may be bugs in your code. If not then investigate verlet integration, and Runge-Kutta methods. There are several existing implementations that can be searched for and compared with yours.  
> **username_2**: You first want to multiply your gravity vector with dt and then add this number to the velocity vector. Otherwise you're not solving Newton's ODE.  
> **username_0**: But that leaves me with a huge number as the velocity vector wich I add to position  
> **username_2**: We might solve your problem easier if you'd just post which discretization you tell the computer to do.  
> **username_0**: I am not sure what you mean  
> **username_2**: Write down the formula you hack into the computer.  
> **username_3**: The process you describe is called Euler integration. It's the most basic numerical  integration method, and notoriously bad. Even if you make the time step tiny it's performance is poor, especially if the orbit is eccentric.  
> **username_3**: So you need a better integrator, preferably one that's symplectiic (conserves energy) like leapfrog or verlet.  Lots of people use and recommend Runge-Kutta methods, but they aren't so suitable for pure gravity sims because they aren't symplectic.  
> **username_3**: If you're getting unexpectedly huge numbers then you may have a mistake somewhere with your units, eg using seconds where you should have days or years.  And with Euler integration you do need tiny time steps.  IIRC, doing Mars with Euler the orbit won't even close up unless delta t is around a day or smaller.  
> **username_4**: Search stackoverflow for "numerically solving kepler's equation". One of the solutions I found helpful for this is: http://stackoverflow.com/questions/5287814/solving-keplers-equation-computationally. The js program at the bottom is a very nice step through to solving the equations using an iterative method.  
> **username_0**: So what I could do is during physics time(normal delta time where everything happens at realtime, I can use semi-implict euler, and as soon as I switch to faster time, i convert position and velocity to kepler elements and calculate it with kepler laws? It makes sense that larger deltatime values result in less accurate calculations, but could that be the reason that the orbit slowly falls down?  
> **username_0**: Also, Ido know how to solve in keplerian elements, i just dont want to use it while in physics mode, so I can simulate effects of a rocket burn better, that by converting between kepler and cartesian  

---

## 1 answers

### Answer by username_5 on 2017-05-19. Score: 7 | Accepted answer

I like to classify solutions of the problem of the time evolution of the
complete initial state of a set of objects at some epoch time, where the
objects are subject to Newtonian gravitation into two main groups.
One approach is to use orbital elements of some sort.
The other is to use a numerical initial value problem solver,
aka a numerical integrator. The latter is the primary subject of this answer.

Note well: This classification isn't quite perfect as hybrid approaches are
also possible, wherein one uses a numerical integrator to integrate
time-varying orbital elements.

Keplerian orbital elements work quite nicely in the case of two point
masses, or more generally, two objects with a spherical mass distribution.
The anomaly is the only Keplerian orbital element that changes over time.
Keplerian elements can be used in situations where the underlying
assumptions are approximately correct by developing a model of how
those supposedly unchanging elements vary with time. One way of doing
this is to use Lagrange's planetary equations. (There are other related
approaches such as Gauss' planetary equations, Delaunay's planetary
equations, etc.)
Lagrange's planetary equations yield expressions for how Keplerian orbital
elements vary over time given a set of perturbing forces.

Another approach is to use something akin to those Keplerian elements
such as Delaunay elements), coupled with planetary equations for those
alternative elements. Yet another approach is to use orbital elements
(e.g., Brouwer-Lyddane elements, SGP4 elements) in which the planetary
equations are embedded in the orbital element to Cartesian state
transformation algorithm. This final approach is used to this day to
describe vehicles in Earth orbit.

The other approach is to use numerical integration.

I'll start with a discussion of how to solve for the value of a scalar
function $x(t)$ at some time $t\_1$ given an initial value $x(t\_0) = x\_0$
and some well-behaved (continuous and bounded) derivative function
$f(x(t),t) = dx(t)/dt$ that describes the time evolution of $x$.
This falls in the very broad category of initial value problems.

Suppose the ordinary differential equation cannot be solved analytically
and cannot be expressed terms of a useful power series. This doesn't mean
nothing can be done. There are a number of techniques for solving this
problem numerically.

Note the dependence of the derivative function on the dependent variable $x$.
This becomes the much simpler problem of numerical quadrature if the
derivative function can be expressed independent of $x$.
The discussion that follows assumes that the derivative function $f$ does
indeed depend on $x$.
Note well: Newtonian gravitation falls in this category.

It also assumes the derivative function is well-behaved.
Numerically integrating across a discontinuity is a bad idea.

The foundation of the integration-based techniques for a scalar function
is the mean value theorem, which says that at some time $t\_c$ between
$t\_0$ and $t\_1$, the value at $t\_1$ is exactly
$x(t\_1) = x(t\_0) + (t\_1-t\_0)\,f(x(t\_c),t\_c)$. If only we could find that
magical $t\_c$ and the derivative at that point. There's a chicken and
egg problem here: that magical point in time is not known. Even if it was,
the derivative function depends on state, and that too isn't known.

A very simple approach around this problem is to assume that this magical
point is the initial point:
$$x(t\_1) = x(t\_0) + (t\_1-t\_0)f(x(t\_0),t0)$$
This works quite nicely for values of $t\_1$ that are very close to $t\_0$.
It doesn't work very well at all where $|t\_1 - t\_0|$ is not small.
This suggests splitting the interval $(t\_0, t\_1)$ into a number of smaller
intervals. This results in Euler's method: Apply the above to advance
state to time $t\_0+\Delta t$, then to $t\_0+2\Delta t$, and so on, eventually
reaching the desired time. Euler's method is rather lousy, even for a simple
first order scalar ODE. We can do much better than this. The key reason for
discussing Euler's method is that it is the basis for many other integration
techniques. Learn how it works, then toss it.

One approach to approving on Euler's method is to somehow correct the
result from Euler's method. For example, take an Euler step and compute
the derivative at the end point. Then use the average of those two
derivative values (the original value used to make the Euler step, and
the other from the end of the Euler step) to recompute the step from
$t$ to $t+\Delta t$. THis is Heun's method.

Another approach is to guess that the magical point $t\_c$ lies
somewhere between $t$ and $t+\Delta t$. Perhaps the middle? We can use
Euler's method to advance state to the midpoint, and then use the
derivative at that point to advance state from $t$ to $t+\Delta t$.
This is the midpoint method.

Both Heun's method and the midpoint method appear to be steps backwards,
computationally. While Euler's method requires but one evaluation of the
derivative function per time step, these improved methods require two.
However, the error growth is in general so much smaller with either
Heun's method or the midpoint method compared to that from Euler's method.
This means that those "improvements" most definitely are improvements.
The expense of calling the derivative function twice per step is more than
offset by the fact that imprpovements enables take steps that are orders
of magnitude larger than one can make with Euler's method.

Both Heun's method and the midpoint method are simple improvements.
This problem has been studied in many guises. There are many more advanced
techniques. One is the class of Runge-Kutta integration techniques.
Both Heun's method and the midpoint method fall into this class.
The most popular of these, classical Runge-Kutta 4, is a significant
improvement on those two methods. There are even higher order Runge-Kutta
integrators than RK4.
Heun's method also falls into the broad class of predictor-correctors,
wherein one method (the predictor) advances state to the end of the interval
and another method (the corrector) uses the derivative at this approximate
endpoint to correct the guess made by the predictor.

The above focused solely on first order ODEs involving a scalar function.
What if the problem is multidimensional or involves higher order derivatives?
The mathematics of the techniques described above can easily handle
multidimensional data: Simply use the vector-valued time derivative.
Since a higher order ODE can be converted to a first order ODE via an
augmented vector-valued state, the same approaches used to address
multidimensional data can also be employed to address higher order ODEs.

There's a problem with doing this: It throws out geometry. For example,
consider the rather simple first order ODE $\dot x = -y, \dot y = x$.
The solution to this multivariate ODE is uniform circular motion.
Applying Euler's method to this results in
$$\begin{aligned}
x(t+\Delta t) &= x(t) - \Delta t\,y(t) \\
y(t+\Delta t) &= y(t) + \Delta t\,x(t)
\end{aligned}$$
The square magnitude of this new vector is
$(x(t)^2+y(t)^2)(1+\Delta t^2)$, which is always greater than magnitude
of the vector at the start of the step. This is not uniform circular motion.
The solution obtained via Euler's method instead spirals out. Other techniques
spiral inward. A geometric integrator on the other hand will somehow
maintain the constraint that $x^2+y^2$ is a constant of motion.

The above example showed why we don't want to toss geometry in a very
simple problem. The geometry of Newtonian gravitation, along with much
classical mechanics in general, is symplectic geometry.
This is why symplectic integrators are of great concern. A simple example
again, with Euler's method: Suppose the second derivative of position
is given by some function $\ddot {\vec x}(t) = \vec f(x(t),t)$.
Applying the basic Euler method against 3+3 dimensional phase space
dictates that
$$\begin{aligned}
\vec x(t+\Delta t) &= \vec x(t) + \Delta t \vec v(t) \\
\vec v(t+\Delta t) &= \vec v(t) + \Delta t \vec f(\vec x(t),t)
\end{aligned}$$
A simple change makes this symplectic:
$$\begin{aligned}
\vec v(t+\Delta t) &= \vec v(t) + \Delta t \vec f(\vec x(t),t) \\
\vec x(t+\Delta t) &= \vec x(t) + \Delta t \vec v(t+\Delta t)
\end{aligned}$$
As is the case with the scalar techniques discussed at the start,
Euler's method is a starting point rather than the end with
regard to symplectic integration techniques. Symplectic Euler's method
is rather lousy. But at least orbits don't spiral outward.

In making an N-body gravitation simulation, the size of N (the number of
bodies) is a key concern. Simulating a galaxy is a very different concern
from simulating a star system. The techniques used in simulating the
formation of a galaxy are very different from those used to develop a
solar system ephemeris. Galactic scale simulations cannot afford to
calculate all of the N2 gravitational interactions amongst all
the particles, and because N is so large, it cannot afford anything
more complex than very simple integrators. A star system model that does
not calculate all N2 of the gravitational interactions
or that uses very a simple integrator will by viewed in disdain.


---

