N-Body Orbit Simulation with Runge-Kutta

In a previous post I introduced a simple orbital simulation program written in python. In that post we simulated orbits by simply taking the location, and velocities of a set of masses, and computed the force on each body. Then we calculated where all the bodies would end up under that force a small time step into the future. This process was repeated over and over again, and was used to simulate gravitational systems like our solar system giving outputs like you see below.

This technique is called the Euler method. If you're not familiar with using numerical methods to simulate orbits, I'd recommend taking a look there first!

Part 1: Python N Body Simulation

Orbit paths from the previous example

In this post I will be adding a more advanced time stepping technique called the Fourth Order Runge-Kutta method. Kids these days just call it RK4. I'll walk through the logic behind RK4, and share a python implementation. I will also link to a C++ implementation, and do a brief performance comparison.

Fourth Order Runge-Kutta

So we are already familiar with the Euler method from the previous post. This method is lacking because it assumes constant acceleration through your whole timestep, which is usually not true. It still works though because the steps are small enough it doesn't matter.

RK4 suffers from the same flaw, but improves on the method by taking a few acceleration values and averaging them. As we are fourth order, we're taking four approximations. 

The four locations are called k values and are as follows:

k1 = the acceleration at the begging
k2 = the acceleration at the middle
k3 = the acceleration at the middle if the acceleration at the beginning was k2
k4 = the acceleration at the end if the acceleration at the beginning was k3

In other words:

k1 = acceleration based on your bodies location at the starting point.
k2 = acceleration 0.5 timesteps (t = 0.5) in the future (other bodies stay in their old locations)
k3 = acceleration at t = 0.5, under the acceleration of k2
k4 = acceleration at t = 1 under acceleration k3

Then a weighted average is taken of these values, and you find your final acceleration value as:

acceleration =  (k1 + 2 * k2 + 2 * k3  + k4) / 6

These weightings are determined by a butcher table, which I don't totally understand, so I'll just link to wikipedia. The result of this average is that now our value contains information about how the acceleration will change throughout your timestep, and corrects accordingly.

It's similar to the Euler technique if the timestep were four times smaller, but since we don't update the location of all the other bodies for the k calculations, it requires less computation. Less computation is good!

Python RK4

As we are introducing a new integration method, I decided to make individual classes for the RK4 and Euler integrators. These are created by passing a list of gravitational bodies, and a time step duration. Integrators must also declare a compute_gravity_step function.

In RK4, the location and velocity updates remain the same as in the Euler case, which is to say simply the time step multiplied by the acceleration to get updated velocity, and timestep multiplied by velocity to get updated location.

The interesting bit is inside our new calculate_single_body_acceleration function. In this function, we now add the four k values as described in the previous section.

That looks like this:

The full run-able python script can be found here.

Results

In the following orbit graphs, you can see the benefit of RK4. In these examples, I created an artificial satellite object and placed it near the sun with some random velocity.

Satellite starting info: location (m) (0, 1.0e10, 0), mass (kg) = 1e23, velocity (m/s) (3e4, 0, 0)

When using the Euler integrator, it appears to enter a stable orbit. When using RK4 however, the forces while passing the sun reveal a decaying and unstable orbit.

Euler orbit of a artificial satellite object with a stable orbit

RK4 path of artificial satellite object displaying decaying and unstable orbit

C++ Orbit Simulation

Our python orbit simulation is pretty neat, but if you add all the planets into the mix and it gets pretty slow. While there are a lot of ways to make python run fast, like simply using numpy, most high performance python is really just calling C or C++ code underneath. When it comes to performance, I prefer to skip the middle man and go straight to C++.

(If you feel like telling me about Fortran here, double check your modern performance numbers, and show me a well organized Fortran project with more than 10 files)

I also plan to add some parallelism in future posts, so I decided now was a good time to go to C++.

So I implemented this simulation in C++ including RK4. The program layout is very similar to python, and can be found on github here. There is nothing parallel or performance oriented here yet though, besides the properties of the languages. Checkout the README for build instructions. The main deviation from python is simply syntax, I overloaded the point basic math operations, so now the component wise (a.x += b.x) automatically occurs without explicitly writing it.

To compare the performance of the two, I was going to just show some bar charts. The differences are so huge though, the C++ bar doesn't even appear visible, so I'll go with a table instead.

Description	Python Runtime (s)	C++ Runtime (s)
10 Body Solar System, 1000 iterations, report interval 100	12.2	0.11
10 Body Solar System, 10000 iterations, report interval 1000	116.2	0.63
5 Body Inner Solar System, 1000 iterations, report interval 100	3.48	0.04
5 Body Inner Solar System, 10000 iterations, report interval 1000	29.5	0.19

So we see C++ providing over 100x performance simply by changing languages!

What a rush.

If you'd like to take a crack at getting some python performance numbers, I'd love the discussion. In future posts I'll be discussing threading and maybe even GPU acceleration if I can build a system large enough to require it.

Please let me know if you have any good material regarding the  Fast Multipole Method!