Hacker News Board Hacker News, browsed like a classic message board
Live data from Hacker News

Writing N-body gravity simulations code in Python

alvinng4.github.io

Post reply on HN
Print view View original

41–50 of 51 posts

Re: Writing N-body gravity simulations code in Python

#41

Post by markstock »

hermitcrab wrote: post #39
markstock wrote: post #34

Earlier quoted context omitted.

Yes, the author uses a globally-adaptive time stepper, which is only efficient for very small N. There are adaptive time step methods that are local, and those are used for large systems. If you see bodies flung out after close passes, three solutions are available: reduce the time step, use a higher order time integrator, and (the most common method) add regularization. Regularization (often called "softening") remo…

>Regularization (often called "softening") removes the singularity by adding a constant to the squared distance. So 1 over zero becomes one over a small-ish and finite number. IIRC that is what I did in the end. It is fudge, but it works.

It is a fudge if you really are trying to simulate true point masses. Mathematically, it's solving for the force between fuzzy blobs of mass.

Re: Writing N-body gravity simulations code in Python

#42

Post by itishappy »

itishappy wrote: post #17
munchler wrote: post #5

In Step 2 (Gravity), why are we summing over the cube of the distance between the bodies in the denominator? Edit: To answer myself, I think this is because one of the factors is to normalize the vector between the two bodies to length 1, and the other two factors are the standard inverse square relationship.

You got it. The familiar inverse square formula uses a unit vector: a = G * m1 * m2 / |r|^2 * r_unit r_unit = r / |r| a = G * m1 * m2 / |r|^3 * r

Oops, pretty big mistake on my part. I gave the formula for force and labeled it acceleration. Way too late to edit, but the correctly labeled formulas should have been:

    F = G * m1 * m2 / |r|^2 * r_unit
    F = G * m1 * m2 / |r|^3 * r
Or as acceleration, by dividing out an `m` using `F=ma`:

    a = G * m / |r|^2 * r_unit
    a = G * m / |r|^3 * r

Re: Writing N-body gravity simulations code in Python

#44

Post by itishappy »

quantadev wrote: post #21
itishappy wrote: post #17

Earlier quoted context omitted.

You got it. The familiar inverse square formula uses a unit vector: a = G * m1 * m2 / |r|^2 * r_unit r_unit = r / |r| a = G * m1 * m2 / |r|^3 * r

The force itself is `G * m1 * m2 / (r^2)`. That's a pure magnitude. The direction of the force is just the unit vector going from m1 to m2. You need it to be a unit vector or else you're multiplying up to something higher than that force. However, I don't get why you'd ever cube the 'r'. Never seen that. I don't think it's right, tbh.

> I don't get why you'd ever cube the 'r'.

It's pulled out of the unit vector. Might be more clear if I notated the vector bits a bit:

    old    : new
    r      : r_vec
    |r|    : r_mag
    r_unit : r_dir
As you know, a vector is a magnitude and direction:

    r_dir = r_vec / r_mag
So the formulas from before become (also correctly labeled as `F` per my other comment):

    F = G * m1 * m2 / r_mag^2 * r_dir
    F = G * m1 * m2 / r_mag^2 * r_vec / r_mag
    F = G * m1 * m2 / r_mag^3 * r_vec

Re: Writing N-body gravity simulations code in Python

#45

Post by mkoubaa »

markstock wrote: post #41
hermitcrab wrote: post #39

Earlier quoted context omitted.

>Regularization (often called "softening") removes the singularity by adding a constant to the squared distance. So 1 over zero becomes one over a small-ish and finite number. IIRC that is what I did in the end. It is fudge, but it works.

It is a fudge if you really are trying to simulate true point masses. Mathematically, it's solving for the force between fuzzy blobs of mass.

You are never simulating pure anything. All computational models are wrong. Some are useful

Re: Writing N-body gravity simulations code in Python

#46

Post by quantadev »

itishappy wrote: post #44
quantadev wrote: post #21

Earlier quoted context omitted.

The force itself is `G * m1 * m2 / (r^2)`. That's a pure magnitude. The direction of the force is just the unit vector going from m1 to m2. You need it to be a unit vector or else you're multiplying up to something higher than that force. However, I don't get why you'd ever cube the 'r'. Never seen that. I don't think it's right, tbh.

> I don't get why you'd ever cube the 'r'. It's pulled out of the unit vector. Might be more clear if I notated the vector bits a bit: old : new r : r_vec |r| : r_mag r_unit : r_dir As you know, a vector is a magnitude and direction: r_dir = r_vec / r_mag So the formulas from before become (also correctly labeled as `F` per my other comment): F = G * m1 * m2 / r_mag^2 * r_dir F = G * m1 * m2 / r_mag^2 * r_vec / r_mag…

Ok, I see what you're doing. Your multiplying the force vector by a non unit-vector, and then dividing back out the linear amount to correct for it. You never see this in a physics book because it's a computational hack, probably because it saves you the CPU cost of not having to do the 3 division operations it takes to get each component (X,Y,Z) of the unit vector.

This makes sense to do in computer code also because if you were going to raise r_mag to a power, you might as well raise it to 3 instead of 2, because it's not extra cost, but you do avoid the three divisions, by never calculating a unit vector. Back when I was doing this work, was decades ago and I had no idea about cost of floating points. Thanks for explaining!

Re: Writing N-body gravity simulations code in Python

#47

Post by itishappy »

quantadev wrote: post #46
itishappy wrote: post #44

Earlier quoted context omitted.

> I don't get why you'd ever cube the 'r'. It's pulled out of the unit vector. Might be more clear if I notated the vector bits a bit: old : new r : r_vec |r| : r_mag r_unit : r_dir As you know, a vector is a magnitude and direction: r_dir = r_vec / r_mag So the formulas from before become (also correctly labeled as `F` per my other comment): F = G * m1 * m2 / r_mag^2 * r_dir F = G * m1 * m2 / r_mag^2 * r_vec / r_mag…

Ok, I see what you're doing. Your multiplying the force vector by a non unit-vector, and then dividing back out the linear amount to correct for it. You never see this in a physics book because it's a computational hack, probably because it saves you the CPU cost of not having to do the 3 division operations it takes to get each component (X,Y,Z) of the unit vector. This makes sense to do in computer code also becaus…

Glad I could help!

Also fun is that taking the magnitude involves a square root that can sometimes be avoided, but that doesn't really help us here because of the power of three. If the denominator were squared we could just use `r_mag^2 = r_x^2 + r_y^2`, but we still need the root to get the direction. It is kinda interesting though that in 2d it expands to a power of `3/2`:

    F_vec = G * m1 * m2 / (r_x^2 + r_y^2) ^ (3/2) * r_vec

Re: Writing N-body gravity simulations code in Python

#48

Post by mkoubaa »

kbelder wrote: post #37
JKCalhoun wrote: post #14

Earlier quoted context omitted.

Decades ago ... I think it was Computer Recreations column in "Scientific American" ... the strategy, since computers were less-abled then, was to advance all the bodies by some amount of time — call it a "tick" — when bodies got close, the "tick" got smaller: therefore the calculations more nuanced, precise. Further apart and you could run the solar system on generalities.

Cut the time interval shorter and shorter, and then BAM you've independently discovered calculus.

You need to invent calculus to get partial differential equations. You need invent PDEs to get approximate solutions to PDEs. You need to invent approximate solutions to PDEs to get computational methods. You need to invent computational methods to explore time integration techniques.

Re: Writing N-body gravity simulations code in Python

#49

Post by quantadev »

itishappy wrote: post #47
quantadev wrote: post #46

Earlier quoted context omitted.

Ok, I see what you're doing. Your multiplying the force vector by a non unit-vector, and then dividing back out the linear amount to correct for it. You never see this in a physics book because it's a computational hack, probably because it saves you the CPU cost of not having to do the 3 division operations it takes to get each component (X,Y,Z) of the unit vector. This makes sense to do in computer code also becaus…

Glad I could help! Also fun is that taking the magnitude involves a square root that can sometimes be avoided, but that doesn't really help us here because of the power of three. If the denominator were squared we could just use `r_mag^2 = r_x^2 + r_y^2`, but we still need the root to get the direction. It is kinda interesting though that in 2d it expands to a power of `3/2`: F_vec = G * m1 * m2 / (r_x^2 + r_y^2) ^ (…

Yeah, on paper (or mathematical symbolics) it comes down to what's more clear and representing reality. That's why I initially said I know there's no cubic relations in the physics of this, which was correct.

But that doesn't mean that therefore there's no correct physics equations (for gravity) involving the cube of a distance, even when there's only squares in these "laws" of physics.

In both cases the power of 2, as well as 3/2, is there merely to "cancel out" the fact that you didn't use a unit vector (in the numerator) and therefore need to divide that out in the denominator, to end up scaling the force magnitude against a unit vector.

Re: Writing N-body gravity simulations code in Python

#50

Post by [deleted] »

igouy wrote: post #38
bgwalter wrote: post #33

For comparison, the famous Debian language benchmark competition: https://benchmarksgame-team.pages.debian.net/benchmarksgame/...

And for those who have access 2022 "N-Body Performance With a kD-Tree: Comparing Rust to Other Languages" https://ieeexplore.ieee.org/document/10216574 2025 "Parallel N-Body Performance Comparison: Julia, Rust, and More" https://link.springer.com/chapter/10.1007/978-3-031-85638-9_...

[deleted]
Post reply on HN