As the creator of the Unity asset Gravity Engine I was curious about how well Unity modeled orbits without all my extra work. In investigating the way that Unity is doing physics I came across a very dramatic difference in performance based on the order of two lines of code and a numerical method I should have known about!

A Unity project with these demonstrations can be found on github.

The “natural approach” is to use Unity do the physics  with a rigidbody AddForce() to provide the gravitational interaction. This requires a script with a FixedUpdate() that determines the force of gravity due to some other object based on that objects mass and distance. For simplicity we’ll consider just one “heavy” object in the scene at the origin. The result can be seen below. The trail of the object does not form a closed orbit (as it should for Newtonian gravity) and the orbit axis shifts over time. Does this provide a clue about the way Unity is doing physics?

 

 

In order to investigate how Unity is evolving the force (i.e. what numerical approach it is using) we can now create some mini-physics engines of our own and see how they compare to the result from AddForce().

Creating a simple physics integrator for a single body is straight-forward. In a FixedUpdate() method we can update position by a time step, calculate acceleration due to gravity and use it to advance velocity by a time step.

The simplest algorithm (Euler)  updates position and velocity each time step (r, v and a can be Vector3 or this can be taken as pseudo-code for array operations):

// Euler update
r = r + v * dt; 
v = v + a * dt;

In my initial coding, I did this from memory and typed:

// NOT Euler (velocity is updated first)
v = v + a * dt; 
r = r + v * dt;

this is NOT Euler and this make a very big difference!

In addition we’ll compare the simplest of the integrators in Gravity Engine, named LeapFrog (or Verlet depending where you look). It is:

// Leapfrog
v = v + a * 0.5 * dt;

r = r + v * dt;

v = v + a * 0.5 * dt;

Taking a look at these algorithms:

Wow. Euler is terrible; the orbit is not even closed. The “Not Euler” method is pretty good and is quite close to what Unity does. Leapfrog is also pretty good and looks to be about the same as “Not Euler”.

The “Not Euler” method is more correctly called the Semi-Implicit Euler method. The key is that the velocity is advanced to the next time step before it is used. This seems like a weird mix – using the velocity from a time step ahead to update the position in this time step. What surprised me was given how well this works, why it does not come up in the scientific introductions to N-body integration (e.g. Moving Stars Around). I assume they go direct to Leapfrog because they are headed to more complex algorithms and Leapfrog is a useful segue. The Semi-Implicit Euler method shares an attribute with the Leapfrog integrator. They are both good for systems in which energy is conserved, like gravity. Technically they are called symplectic integrators due to the way the preserve area in a phase space of position and velocity. Semi-Implicit Euler is the simplest symplectic integrator, Leapfrog/Verlet is a second order symplectic integrator.

What is AddForce() doing?

It appears that Unity AddForce() is keeping up well with Semi-Implicit Euler and Leapfrog.

Can we figure out which one it is using?

Let’s take a look at how well the integrators are conserving energy. It’s straight forward to determine the kinetic and potential energy at each time step and with a second camera in the Unity scene we can do a very basic plot of energy vs time by using energy and time to position three new objects in the view of the second camera. Here’s the result with the default time step for FixedUpdate:

The energy errors for AddForce() and Semi-Implicit Euler are overlaid on each other, while the Leapfrog method has a more stable energy. All of them maintain the initial energy except during the close approach to the central mass. Unity appears to be using Semi-Implicit Euler.

Getting Better Orbits

What steps do we need to take to make the orbits show a single ellipse that does not precess? [pun intended, changing the number of steps is exactly what we’ll need to do]. The Unity sample scene allows the numerical integrators to choose a smaller dt and then interate more times per fixed update to stay on the same time scale as the AddForce() approach.

Let’s try a factor of 10 reduction:

Those orbits look good. Notice that the energy error for SI Euler is much reduced and the error for Leapfrog is barely detectable.

The specific time step choice will depend on the maximum velocity of the orbit.

There is one detail about symplectic integrators worth knowing: the time step must stay constant otherwise they do not conserve energy. There are other higher order integrators that allow the time step to change based on the maximum velocity in the system. Gravity engine provides one such integrator. More sophisticated scientific code maintains a heirarchy of time steps and allows e.g. fast moving inner planets to use a smaller time step than much slower outer planets. This is on the roadmap for Gravity Engine.

Other Considerations

The numerical code in the sample project is useful to do some simple investigations but does not really scale. A scene with 10 bodies would require each of them to know about the other 9 to compute the total gravitational force. There also be wasted work, since the force of 1 on 2 is equal and opposite to the force of 2 on 1 and both would be calculated. If you plan to do many-body gravity then having a central engine to do it is more effective.

Gravity Engine provides this and MANY more features.