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)

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

  • 1
    $\begingroup$ 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. $\endgroup$
    – James K
    May 18, 2017 at 17:56
  • $\begingroup$ 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. $\endgroup$ May 18, 2017 at 18:01
  • $\begingroup$ But that leaves me with a huge number as the velocity vector wich I add to position $\endgroup$
    – GR00G0
    May 18, 2017 at 18:05
  • 1
    $\begingroup$ Everytime I've tried this, I've ended up using too large a DeltaTime. Try a smaller value and see if things get better. $\endgroup$
    – user21
    May 19, 2017 at 0:28
  • 2
    $\begingroup$ 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. $\endgroup$
    – PM 2Ring
    May 19, 2017 at 15:54

1 Answer 1


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.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .