Python Help designing an n-body gravity simulator in Python

[Wer900]

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.

 

[Filip Larsen]

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

 

[phyzguy]

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.

 

[jtbell]

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

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 k₁, k₂, 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.

 

[Wer900]

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

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).

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.

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.