Writing N-body gravity simulations code in Python(alvinng4.github.io) |
Writing N-body gravity simulations code in Python(alvinng4.github.io) |
https://alvinng4.github.io/grav_sim/5_steps_to_n_body_simula...
Naturally, this difference becomes big when two objects are close together, since the acceleration will induce a large change of velocity, and a lower time step is needed to keep the error under control.
But in our simulation (reality, or whatever), uncertainty increases the closer we look.
Does that suggest we don’t live in a simulation? Or is “looking closer increases uncertainty” something that would emerge from a nested simulation?
laughs in 32kb of RAM
Sounds like a great project, though. There's a lot of fundamental concepts involved in something like this, so it's a great learning exercise (and fun as well!) :)
For a few-body system (e.g., the Solar System) this probably won't provide any speedup. But once you get to ~100 bodies you should start to see substantial speedups by running the simulator on a GPU.
For the programmer, yes it is easy enough.
But there is a lot of complexity hidden behind that change. If you care about how your tools work that might be a problem (I'm not judging).
Obviously if you want maximum performance you'll probably still have to roll up your sleeves and write CUDA yourself, but there are a lot of situations where you can still get most of the benefit of using a GPU and never have to worry yourself over how to write CUDA.
I'm trying to evaluate if I can get away with f32 for GPU use for my molecular docking software. Might be OK, but I've hit cases broadly where f64 is fine, but f32 is not.
I suppose this is because the dominant uses games and AI/ML use f32, or for AI even less‽
> Large-scale simulation: So far, we have only focused on systems with a few objects. What about large-scale systems with thousands or millions of objects? Turns out it is not so easy because the computation of gravity scales as . Have a look at Barnes-Hut algorithm to see how to speed up the simulation. In fact, we have documentations about it on this website as well. You may try to implement it in some low-level language like C or C++.
Once particles accelerate they'll just 'jump by' each other of course rather than collide, if you have no repulsive forces. I realized you had to take smaller and smaller time-slices to get things to slow down and keep running when the singularities "got near". I had a theory even back then that this is what nature is doing with General Relativity making time slow down near masses. Slowing down to avoid singularities. It's a legitimate theory. Maybe this difficulty is why "God" (if there's a creator) decided to make everything actually "waves" as much as particles, because waves don't have the problem of singularities.
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.
a = G * m1 * m2 / |r|^2 * r_unit
r_unit = r / |r|
a = G * m1 * m2 / |r|^3 * r F = G * m1 * m2 / |r|^2 * r_unit
F = G * m1 * m2 / |r|^3 * r
a = G * m / |r|^2 * r_unit
a = G * m / |r|^3 * rhttps://benchmarksgame-team.pages.debian.net/benchmarksgame/...
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_...
I vividly remember the Pythagorean three-body problem example, and how it required special treatment for the close interactions.
Which made me very pleased to see that configuration used as the example here.
ParallelNBodyPerformance Public
Now you're making me worry I'm not "physicsing" hard enough
And yes that's a dramatic over-simplification of uncertainty principle, but it conveys the concept perfectly.
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
r_dir = r_vec / r_mag
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_vecThis 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!
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_vecYour solution sounds interesting, but isn't it only practical when you have a small number of bodies?
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") 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.
Optimistic fixed time steps are just going to work well almost always and almost everywhere, accumulating errors behind your back at every episode of close approach.
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.