[...]

One of the biggest issues with these simulations is when particles get close. As a particles distance approaches 0, the force applied approaches infinity. So if a particle passes close by a gravity source, it will be ejected out with a lot more energy that it had going in.

If you give particles a radius, so they represent ghostly spheres that can go through each other, the force between two particles is proportional to the 1/distance^2 if they are farther away than 2 radii, but it's proportional to distance if the spheres intersect. That gets rid of the problem of huge forces when the particles are too close together.

You can also check against a minimum distance in the gravity calculation or, if you're lazy like me, just don't apply any force within a certain radius or impose a speed limit. Wikipedia has a stub on a more formal method that looks good: https://en.wikipedia.org/wiki/Softening

However, these don't really solve the inherent problem: A particle samples the gravity field as it moves. The faster it moves, the fewer samples it can collect. The fewer samples it can collect, the less accurate the result. The only way to get an accurate result is to shorten the timestep. If this results in unacceptable performance for your desired accuracy, then the solution is an adaptive timestep.

I simply cheated: If a particle is too close, it keeps its velocity unchanged. It then goes by the other and regains enough distance to calculate properly. Yeah, before I did that, those things went flying EVERYWHERE!

Not sure I fully understand the advantage of the leapfrog or midpoint? I mean, the implementation seems easy enough, but what do I get out of it?

Asfor the octree, I keep getting more keen on it. I have another experiment running 25,000 particles at okay rate, but at some point, I clearly need to familiarize myself with some octree algos!

However, these don't really solve the inherent problem: A particle samples the gravity field as it moves. The faster it moves, the fewer samples it can collect. The fewer samples it can collect, the less accurate the result. The only way to get an accurate result is to shorten the timestep. If this results in unacceptable performance for your desired accuracy, then the solution is an adaptive timestep.

This is the eternal issue, isn't it. What bothers me is that my brain keeps telling me that there should be a way to compress all the forces into a single equation (possibly one per particle) and find the position of the particle as a function of time. Like if you calculate two bodies influencing each other; you can (to simplify it) calculate a barycenter and simply input the point in time you want to sample their positions and voila, the equation returns the positions. There would be no need for timestep, and CPU pressure is next to nothing (compared tot he current options available, anyway). I have even done some theoretical work on how to move from two to three bodies, hoping that it will let me figure out how to do it with n bodies, but it seems a bit.... challenging. So for now, lesser methods apparently need to be used....

I simply cheated: If a particle is too close, it keeps its velocity unchanged. It then goes by the other and regains enough distance to calculate properly. Yeah, before I did that, those things went flying EVERYWHERE!

Okay, that's one solution indeed! :-) I guess it's important to remind yourself, what are you trying to achieve. Are you trying to model the actual physics or make something that LOOKS real for gaming purposes? If the latter, then tricks like this or softening which was mentioned earlier (which is used in many 'real' science/astro simulations btw), are perfectly valid. Regarding the suggestion about adaptive timesteps, this is again something you would definitely do to follow the physics more accurately but then you have unpredictable rates since your simulation slows down whenever two stars get close and then speeds up again, all very jerky and unstable. So a fixed timestep works fine if your not bothered about getting things right with close interactions. But if you want to model the physics accurately (out of some scientific curiosity for instance) then definitely consider it.

Not sure I fully understand the advantage of the leapfrog or midpoint? I mean, the implementation seems easy enough, but what do I get out of it?

The advantage is you get more accurate integration of the orbits (actually important for both scientifically accurate simulations and for making it look better) since leapfrog/mid-point gives second-order errors compared to Euler integration which is first-order. A leapfrog integrator is also a special kind of integrator (called a symplectic scheme) which has wonderful conservation properties which amongst others preserves circular and elliptical orbits (e.g. planets or binaries) very well. The thing is, all these schemes are as expensive as each other, just one force computation per timestep so in that case you would use the better one right? :-)

This is the eternal issue, isn't it. What bothers me is that my brain keeps telling me that there should be a way to compress all the forces into a single equation (possibly one per particle) and find the position of the particle as a function of time. Like if you calculate two bodies influencing each other; you can (to simplify it) calculate a barycenter and simply input the point in time you want to sample their positions and voila, the equation returns the positions. There would be no need for timestep, and CPU pressure is next to nothing (compared tot he current options available, anyway). I have even done some theoretical work on how to move from two to three bodies, hoping that it will let me figure out how to do it with n bodies, but it seems a bit.... challenging. So for now, lesser methods apparently need to be used....

I'm going to have to disappoint you on this one. You are of course correct for two bodies. This problem was solved by Sir Isaac Newton some 300+ years ago :-) . And in the intermediate centuries, everybody and their cat has tried to add just one more body to solve the so-called 3-body problem. But there is no general solution I'm afraid that can be used like you hope for. There are various special so-called 'restricted' 3-body problems where the motion of 1 body is limited in some way so it becomes a 2-body problem again or has some special solution for some configuration. But in general, and especially when this becomes 4, 5 ... N bodies, you have to do the numerical integration I'm afraid! This is one of the step-ups from linear/soluable equations to non-linear equations. All the 'geniuses' of yesteryear solved the 'easy' linear stuff :-D and left all the difficult non-linear stuff to our generation! :-/ Oh well ...

Regarding the oct-tree, I'm still a little unsure if you really need it for gravity integration (for the reasons discussed earlier about the ~1000 body limit). However, if you want to model some kind of larger Universe with large-scale structure and then 'zoom' into galaxies and further, then having a tree structure to organise this data hierarchically is obviously a good thing. If you can combine it with helping to compute the gravity, then great but I'm still unsure if you'll get that much benefit while maintaining 30/60fps.

However, these don't really solve the inherent problem: A particle samples the gravity field as it moves. The faster it moves, the fewer samples it can collect. The fewer samples it can collect, the less accurate the result. The only way to get an accurate result is to shorten the timestep. If this results in unacceptable performance for your desired accuracy, then the solution is an adaptive timestep.

This is the eternal issue, isn't it. What bothers me is that my brain keeps telling me that there should be a way to compress all the forces into a single equation (possibly one per particle) and find the position of the particle as a function of time. Like if you calculate two bodies influencing each other; you can (to simplify it) calculate a barycenter and simply input the point in time you want to sample their positions and voila, the equation returns the positions. There would be no need for timestep, and CPU pressure is next to nothing (compared tot he current options available, anyway). I have even done some theoretical work on how to move from two to three bodies, hoping that it will let me figure out how to do it with n bodies, but it seems a bit.... challenging. So for now, lesser methods apparently need to be used....

I really can't tell if you are trolling, at this point. I'll feed you a Wikipedia page about the subject: https://en.wikipedia.org/wiki/N-body_problem . No, you will not find an analytical solution to anything beyond n=3.

I see octrees already mentioned, but in case you don't know about it, there's a name for this technique used for this specific purpose: Barnes-Hut. You can trade accuracy for efficiency.

Some links that explain more:

https://en.wikipedia.org/wiki/Barnes-Hut

http://www.cs.princeton.edu/courses/archive/fall04/cos126/assignments/barnes-hut.html

Okay, that's one solution indeed! :-) I guess it's important to remind yourself, what are you trying to achieve. Are you trying to model the actual physics or make something that LOOKS real for gaming purposes? If the latter, then tricks like this or softening which was mentioned earlier (which is used in many 'real' science/astro simulations btw), are perfectly valid. Regarding the suggestion about adaptive timesteps, this is again something you would definitely do to follow the physics more accurately but then you have unpredictable rates since your simulation slows down whenever two stars get close and then speeds up again, all very jerky and unstable. So a fixed timestep works fine if your not bothered about getting things right with close interactions. But if you want to model the physics accurately (out of some scientific curiosity for instance) then definitely consider it.

Definitely going for "close enough"! I will no doubt go for something more advanced later on, but for now, something that is scientifically 'close', but not 'accuracy over usability'.

The advantage is you get more accurate integration of the orbits (actually important for both scientifically accurate simulations and for making it look better) since leapfrog/mid-point gives second-order errors compared to Euler integration which is first-order. A leapfrog integrator is also a special kind of integrator (called a symplectic scheme) which has wonderful conservation properties which amongst others preserves circular and elliptical orbits (e.g. planets or binaries) very well. The thing is, all these schemes are as expensive as each other, just one force computation per timestep so in that case you would use the better one right? :-)

Decisions, decisions. But you do make it sound pretty good, although I will need to study it much closer! I am starting to think I need a mixture...

I'm going to have to disappoint you on this one. You are of course correct for two bodies. This problem was solved by Sir Isaac Newton some 300+ years ago :-) . And in the intermediate centuries, everybody and their cat has tried to add just one more body to solve the so-called 3-body problem. But there is no general solution I'm afraid that can be used like you hope for. There are various special so-called 'restricted' 3-body problems where the motion of 1 body is limited in some way so it becomes a 2-body problem again or has some special solution for some configuration. But in general, and especially when this becomes 4, 5 ... N bodies, you have to do the numerical integration I'm afraid! This is one of the step-ups from linear/soluable equations to non-linear equations. All the 'geniuses' of yesteryear solved the 'easy' linear stuff :-D and left all the difficult non-linear stuff to our generation! :-/ Oh well ...

I am starting to accept that. I am still holding out hope for something similar that is simply 'science-adjacent'. But that is a future thing, I am trying to keep my eyes on the ball for now...

Regarding the oct-tree, I'm still a little unsure if you really need it for gravity integration (for the reasons discussed earlier about the ~1000 body limit). However, if you want to model some kind of larger Universe with large-scale structure and then 'zoom' into galaxies and further, then having a tree structure to organise this data hierarchically is obviously a good thing. If you can combine it with helping to compute the gravity, then great but I'm still unsure if you'll get that much benefit while maintaining 30/60fps.

I am convinced that I will have to do a lot of weighing systems against each other, maybe mix and match. Who would have thought this would be difficult ;-p

I really can't tell if you are trolling, at this point. I'll feed you a Wikipedia page about the subject: https://en.wikipedia.org/wiki/N-body_problem . No, you will not find an analytical solution to anything beyond n=3.

No, not trolling, why would I? I am just hoping that there is something I am ignorant of, but hope is failing me on that particular account... :(

I see octrees already mentioned, but in case you don't know about it, there's a name for this technique used for this specific purpose: Barnes-Hut. You can trade accuracy for efficiency.

Some links that explain more:

https://en.wikipedia.org/wiki/Barnes-Hut

http://www.cs.princeton.edu/courses/archive/fall04/cos126/assignments/barnes-hut.html

Thanks! I am shopping around for good stuff, didn't know the Princeton one!

This is the eternal issue, isn't it. What bothers me is that my brain keeps telling me that there should be a way to compress all the forces into a single equation (possibly one per particle) and find the position of the particle as a function of time.

I'm going to have to disappoint you on this one. You are of course correct for two bodies. This problem was solved by Sir Isaac Newton some 300+ years ago :-) .

Kepler doesn't get the respect he deserves :P
But yes, you are correct. This is an inherently hard problem and there is no analytical solution. Only numerical ones. And lots of trickery and faking-it, as is always the case in game programming.

Not sure I fully understand the advantage of the leapfrog or midpoint? I mean, the implementation seems easy enough, but what do I get out of it?

Better stability in general. Honestly, though, it's probably not worth implementing until you identify a need for it.

https://codepen.io/kemick/pen/YQpepN

Wrote a quick comparison. The particles fall (using a constant gravity force) and then bounce back up using a spring force. The blue lines indicate the starting point and the floor. The red lines indicate the maximum bounds that the particle has moved.

The leftmost box integrates with a bad version of euler. It also explodes very quickly at the given timestep.

 (newPos, newVel) += (oldVel, oldForce)

The middle box integrates with symplectic euler integration. It gains and loses some energy. At the very least, you should be using this form.

(newPos, newVel) += (oldVel, newForce) 

The right box integrates using leapfrog. See the source for details since it's a bit longer.

So as you can see, for the given timestep the leapfrog method conserves its energy very well. However, it may not be necessary depending on your accuracy requirements.

This is the eternal issue, isn't it. What bothers me is that my brain keeps telling me that there should be a way to compress all the forces into a single equation (possibly one per particle) and find the position of the particle as a function of time.

I'm going to have to disappoint you on this one. You are of course correct for two bodies. This problem was solved by Sir Isaac Newton some 300+ years ago :-) .

Kepler doesn't get the respect he deserves :P
But yes, you are correct. This is an inherently hard problem and there is no analytical solution. Only numerical ones. And lots of trickery and faking-it, as is always the case in game programming.

Not sure I fully understand the advantage of the leapfrog or midpoint? I mean, the implementation seems easy enough, but what do I get out of it?

Better stability in general. Honestly, though, it's probably not worth implementing until you identify a need for it.
https://codepen.io/kemick/pen/YQpepN
Wrote a quick comparison. The particles fall (using a constant gravity force) and then bounce back up using a spring force. The blue lines indicate the starting point and the floor. The red lines indicate the maximum bounds that the particle has moved.
The leftmost box integrates with a bad version of euler. It also explodes very quickly at the given timestep.

 (newPos, newVel) += (oldVel, oldForce)

The middle box integrates with symplectic euler integration. It gains and loses some energy. At the very least, you should be using this form.

(newPos, newVel) += (oldVel, newForce) 

The right box integrates using leapfrog. See the source for details since it's a bit longer.
So as you can see, for the given timestep the leapfrog method conserves its energy very well. However, it may not be necessary depending on your accuracy requirements.

Thanks, I needed that kind of more 'tangible' thing. I think you're right, leapfrog may be a bit better in some ways, but it would take an immense effort to rewrtie for it (including learning an dunderstanding it), so that will be for a possible future version.
It still feels utterly bizarre that there is no conversion-to-static-equation solution. But something tells me others have said the same as me, and then lost decades to finding an answer. The bold should beware, I guess :-o

It still feels utterly bizarre that there is no conversion-to-static-equation solution.

The 2-body solution utilizes conic sections (ellipses, in the case of a stable orbit) to approximate the movement of a body. If you wanted, you could model 8 planets orbiting the sun with reasonable accuracy because, with a few exceptions, our planets have stable elliptical orbits. So the issue is not with the number of bodies, but instead how hard their movements are to model.

This becomes more obvious if you look at it from the opposite direction. Instead of looking at the interaction, look at the result.

A single body's trajectory (with no forces applied) can be described by a line. The line only ever needs position/velocity. The velocity is constant.

The trajectories of two bodies can be modeled using conic sections. In the case of stable orbits, that would be an ellipse. An ellipse only ever needs major and minor axes and position/velocity. The velocity changes, but is cyclical and so it repeats.

Now, how would you model the trajectories of three arbitrary bodies doing whatever? In the previous two cases, there was a single equation that you found the parameters to. Here, the velocity can do whatever and so the equation will be some wacky polynomial whose form changes depending on the state of the system. There can be no general equation that you just plug numbers into. And any specific equation that you could generate would grow larger and more wacky as the system progressed, whereas the previous two examples remain constant.

As I've said, this is an inherently hard problem. There are no clever answers to be found, only clever approximations.