February 11, 2017

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!



7 comments:

  1. Hi there,
    Are you familiar with the Java programming language? I have written the RK4 algorithm for my Solar System simulation but, when I run the program with RK4 the orbits of the planets diverge from their original trajectories very fast as the system appears to gain energy somehow. I can't see my mistake! I have asked a few people with many years of programming experience to look at my code and all of them have said it looks correct! Is there any way you could look at my code?

    ReplyDelete
  2. It is a hard thing to debug if it's not work... Are you using a visual debugger? I would implemented Euler (very simple) and then print or view RK4 and Euler steps for the same body. See if you can pinpoint what component is changing at the wrong rate.

    ReplyDelete
  3. Thanks a lot for this! One quick question though. I tried running your RK4 code for a few solar system planets, and found that they all fly away from the sun. I wrote a Euler n-body integration and did not see this behavior. Do you know what might be going on?

    ReplyDelete
  4. Hey, very good question. Honestly I don't understand what exactly this means, but there was some excellent discussion on reddit around this where a commenter explained RK4 is not symplectic (https://www.reddit.com/r/Physics/comments/5tf61n/nbody_orbit_simulation_with_rungekutta/). I have been meaning to understand the differences better and make this post more accurate and educational!

    ReplyDelete
  5. Euler and RK4 are non-symplectic, which means they do not conserve energy over long periods of time, whatever the time step. As I said in the "Euler" page, you have accidentally created a symplectic (energy-conserving) integrator there ;)

    ReplyDelete
  6. You're not using RK4 to solve the location, you are using RK4 to solve acceleration, and then applying Euler to solve for Velocity and Location

    ReplyDelete
    Replies
    1. And you don't need to use numeric methods of integration to solve for acceleration, since you have the formula from the differential equation

      Delete