Making an N-Body Simulation

Making an n-body gravity simulation without wasting any time. After this 9 minute tutorial, you should have all the knowledge you need to make your own basic n-body simulation in your preferred programming language and most importantly, you might even understand it.

I try to keep my videos relevant, no matter what language you use, as long as there's no language explicitly specified in the title. This means that I don't use any hyper-language-specific features and instead try to stick to things that most widely used languages have. This means no specific libraries exclusive to rust and I instead focus on showing code that does what it looks like it does.

Chapters:
0:00 - Intro
0:18 - Bodies
0:35 - Initialization
1:04 - Interaction
2:27 - Integration
4:42 - Softening
5:15 - Optimization
8:06 - Rendering
8:24 - Outro

[DeadlockCode]

If you're having issues, read this first and if it doesn't solve your problem, leave a comment and I'll respond as soon as I can.

Problem 1:
Description: The bodies either, exist for a while but suddenly disappear, or more likely they exist for only a couple of frames before disappearing.
Reason: Two or more bodies share the exact same floating point position. This results in the distance becoming exactly zero, and we all know division by zero is bad. It's particularly bad here because then the acceleration will become NaN and because NaNs are contagious, the velocity and position will also become NaN. This results in whatever method you have for rendering to stop working, maybe it even crashes.
Solution: Either check if the bodies have the same position, and then just don't compute the acceleration for those bodies, or my preferred approach, do 'mag.max(epsilon)' for the final unclamped magnitude in the acceleration calculation.
Why it works: This solves the problem because there is no longer a division by zero, instead a division by epsilon, which is very tiny depending on implementation. This results in the division becoming huge (but importantly not Inf) and multiplying that with the r vector, which is (0, 0), results in a (0, 0) acceleration towards that body which is like they never tried interacting in the first place. For clarity, epsilon is 1.1920929E-7 in my program, and the max function returns the larger of two floats.

[FredGlt]

4:26 One small improvement you could do for accuracy is to change the equation of the position from (P = Po + V * T) to (P = Po + Vo * T + 0.5 * A * T^2), where Po and Vo are the values computed at the previous step of the simulation.
This is because the accurate way to obtain the speed's and the position's equations is by integrating (a = dv/dt) and then (v = dp/dt) instead of approximating these two equivalences directly.

It might not sound like much, but I’ve seen other content creators with projects such as softbody and gravity simulations where it has greatly helped accuracy as well as stability.

[rorymacdonald]

love it! i've been wanting to get into particle/n-body simulations. gonna take a shot at it over winter!

[axionbuster]

For the Barnes-Hut, using the Z-order curve (by bit interleaving) and then building the tree bottom-up from the lowest layer (where every particle becomes a node) and up (lower level of detail) is a speedy way of accomplishing it. From that point on, you can calculate the approximate view angle (or the tangent of the angle) by a simple division. If the angle is too large (about 10 degrees), you go deeper into the tree, but if it’s small enough, you treat the node as a particle and then move on. Once the tree is built, the depth-first traversal is an embarrassingly parallel problem. I built my tree using the left-child-right-sibling (LCRS) data structure, which is dimensionally independent and memory efficient simultaneously.

[jsbaasi]

Brilliant video, loved the graphics and the soothing music. Looking forward to your implementation of the BH algorithm!

[primalaspid]

Just found your channel and your videos are absolutely amazing! Huge potential there!

[academyofuselessideas]

Fantastic! Many people will say that reinventing the wheel is wrong but those are people often obsessed with efficiency... keep reinventing! what matter is that you enjoy it... It is also great that you chose rust for your simulations. Perhaps, for other videos you can use rust's flamegraph to show what parts of the algorithms are consuming the most resources... Another possible question here is how the simulation is affected by numerical error... one of the issues with the n-body problem is that at some point the numerical errors accumulate fast which forces you to reduce the delta t... another interesting improvement would be to make the delta t adaptive... It is a very interesting problem! keep it up!

[HoSza]

Did you test it with only two bodies? If you fix one object somewhere on the screen, then the other should move in an ellipse like curve and if you don't clear the screen but keep the former states if the other object drawn you might notice it keeps spiraling towards the fixed object.

[chocholatebunny]

The quality of this video has made me a sub, keep going friend, you're making some great content and even better Rust content!

[NeroDefogger]

I also have my simulator, you can check where I explain it a bit, the simulations I do at the end are very scuffed, the simulator works much better now, I'm trying to use it to discover more forces, btw interesting the calculating both objects at the same time to try to halve the complexity

[danielbecker]

Great video, but i believe you can improve the accuracy of the positions, At 4:28 you stated that p_{n+1} = p_n + dt*v_n, but i believe you actually want to update the positon based on the average velocity over the time interval (v_a = 1/2*(v_n + v_{n + 1}), which gives us the formula for position as p_{n+1} = p_n + v_n*dt + 1/2*a*dt^2, (i.e. the displacement formula for uniform acceleration)

[newfilms]

this is so well made, congrats!!

you deserve so many more subs 💪💪

[moritzfechtner]

Continue your good work, looking forward to the next video.

[jonathanbaincosmologyvideo]

Can you perhaps show where this solution is applied?
Try perhaps Mercury's infamous Perihelion Precession !

[uberscott]

I'm really enjoying your content.

[alexandratsankova]

You might be able to use "fast inverse square root" or (i remember it being called) Quake inverse root for the calculations, it might speed up the calculations, but I'm not sure by how much.

Also not sure, but it might be possible to use some type of threading to speed up the calculations. But it might be tedious in order to not get overlaps between the threads (and thus not have to use locks/mutexes)

[unveil]

Wow!!! You are a great teacher!!! ❤ thanks!! This was REALLY helpful 🎉

[Miparwo]

Maybe is faster to calculate the gravitational field, which should not be sensible to distant particles.

[davidpiepgrass]

I look forward to your next video! Have you considered using the faster "Fast-Multipole" method instead of Barnes-Hut? My impression is that you just have to compute interactions between pairs of distant quadtree nodes (or something like that) instead of interactions between individual particles and distant quadtree nodes.

[o]

Is the code available anywhere? I've been trying to replicate this simulation but the particles just get launched away from the center