Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Help designing an n-body gravity simulator in Python

  1. Nov 29, 2015 #1
    Hi everyone:

    I wanted to design a 2-dimensional gravity simulator as part of a larger research project. I wrote it in Python using numpy to store positions, velocities, and accelerations, and matplotlib to display the trajectories; code is available upon request.

    In order to solve the problem I've used Euler's method, which I know is numerically unstable over long timescales. In particular, I've been observing that over time many bodies are propelled away, in straight lines, from their original positions. By imposing an artificial minimum limit of 0.3 meters on effective gravitational radius (the masses and gravitational constant I'm using for testing purposes are very different from what I'll actually use) I'm able to reduce this mostly, although this reduces the accuracy of predictions at small radii.

    I've read up on Runge-Kutta integration and other, more exotic methods of calculating orbits, but I have no idea of how to implement them in an orbital simulator. Would anyone here be able to help? Thanks in advance to everyone.
     
  2. jcsd
  3. Nov 29, 2015 #2

    Filip Larsen

    User Avatar
    Gold Member

    If you want a fairly easy improvement over explicit Euler, you may want to try the Leapfrog method [1], which for small enough time steps do not add false energy to the system you are simulating.

    [1] https://en.wikipedia.org/wiki/Leapfrog_integration
     
  4. Nov 29, 2015 #3

    phyzguy

    User Avatar
    Science Advisor

    How do you know that the bodies that are "propelled away in straight lines" are incorrect? Close two body interactions should result in some bodies gaining energy and being ejected. This happens in real clusters and results in the clusters losing energy and becoming more tightly bound over time.
     
  5. Nov 29, 2015 #4

    jtbell

    User Avatar

    Staff: Mentor

    The magic phrase to Google for is "runge kutta for system of differential equations". You'll find pages like this one:

    http://www.physics.buffalo.edu/phy410-505-2009/topic3/lec-3-2.pdf

    Start with Newton's Second Law for each pair of objects. Convert the second-order differential equations into pairs of first-order differential equations. Put the components of all positions and velocities into a single array or vector. The intermediate values k1, k2, etc. in the Runge-Kutta algorithm are now also vectors of the same dimension. You basically re-interpret the scalar R-K equations as vector equations. Write functions that carry out the steps of the calculation, taking vectors as input and returning vectors as output.

    I did something like this about 25 years ago, for a solar-system simulator, but it was in Fortran, not Python, so I can't help with the actual code. Instead of standard 4th-order Runge-Kutta, I used a variation called Runge-Kutta-Fehlberg which varies the step-size dynamically.

    I had the problem that orbits tended to drift outwards because the R-K algorithm tends to not conserve energy over long time spans. Loosely speaking, the errors introduced by the discrete steps tend to reinforce each other instead of cancelling. I understand there are algorithms in which the errors do tend to cancel and energy is at least approximately conserved, but I didn't go any further with my simulator.

    [added] Ah, now I see that Filip also mentioned this "false energy" problem. Maybe a good start would be to vectorize the leapfrog algorithm that he mentioned, and look for higher-order versions if necessary.
     
  6. Nov 29, 2015 #5
    Thanks for the link. I think I tried this before, but maybe reversed the order in which I set velocity and acceleration (causing eventual drifting outward of the system).

    Yes, I understand, but I sometimes saw bodies (representing the majority of the mass of the overall system) propelled outward even in 2 or 3 body systems. This doesn't make a whole lot of sense.

    Also: Thanks jtbell for the info. I'm reading up on it right now, hopefully I can get something working soon.
     
    Last edited: Nov 29, 2015
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Help designing an n-body gravity simulator in Python
  1. N-body problem in Python (Replies: 10)

Loading...