I just tried it myself, using my fixed snippet.
Probably i see the same problem as you do: As planets come close to each other, forces become so huge the integration seemingly fails even with tiny timesteps. If the planets come very close, which they do because gravity pulls them together, they separate very quickly after that. There is a explosive raise in energy, and the simulation clearly isn't right. Maybe i still have bugs.

No wonder we have black holes and i had to invent my own laws of physics to do better :D

I don't have the physics engine running, so i can't quickly compare with RK4. But i doubt this alone would fix it.
It would help to model collisions, so planets can not get too close, but separate before forces become too discontinuous to integrate.
It would also help to model some realistic solar system, initializing all velocities so they represent natural orbits. Then they would not collide either.

		static bool testSimplePlanets = 1; ImGui::Checkbox("testSimplePlanets", &testSimplePlanets);
		if (testSimplePlanets)
		{
			constexpr int count = 10;
			struct Planet
			{
				sVec3 com = sVec3(0);
				sVec3 vel = sVec3(0);
				sVec3 force = sVec3(0);
				float mass = 1;
			};
			
			static Planet planets[count];

			// generate random planets (this code runs only once)
			static bool init = true;
			if (init || ImGui::Button("Reset Simulation"))
			{
				int k = 0;
				auto rnd = [&](float vMin, float vMax)
				{
					float f = PRN::randF(k++);
					return (f * vMax-vMin) + vMin;
				};

				for (int i=0; i<count; i++)
				{
					planets[i].com = sVec3(rnd(-1,1), rnd(-1,1), rnd(-1,1));
					planets[i].vel = sVec3(0);//sVec3(rnd(-1,1), rnd(-1,1), rnd(-1,1)) * 0.001f;//
					planets[i].force = sVec3(0);
					planets[i].mass = rnd(0.001f, 0.001f);
				}

				init = false;
			}

			// function to update forces of a pair of planets
			auto GravityPair = [&](Planet &bodyI, Planet &bodyJ)
			{
				const float G = 1.f; // some number to relate your mass metric to gravity. i just assume we use a mass unit so G becoems 1.

				sVec3 diff = bodyJ.com - bodyI.com;

				float dist = diff.Length();
				float force = G * (bodyI.mass * bodyJ.mass) / (dist * dist);

				sVec3 direction = diff / dist;
				bodyI.force += direction * force;
				bodyJ.force -= direction * force;
			};

			static float timestep = 0.001f; 
			ImGui::SliderFloat("timestep", &timestep, 0.0002f, 0.005f, "%.5f");


			// update (this code runs once every frame)

			// accumulate forces
			for (int i=0; i<count-1; i++)
			for (int j=i+1; j<count; j++)
				GravityPair(planets[i], planets[j]); 

			// integrate 
			for (int i=0; i<count; i++)
			{
				Planet &p = planets[i];

				sVec3 acc = p.force / p.mass;
				p.vel += acc * timestep;
				p.com += p.vel * timestep;

				p.force = sVec3(0); // reset accumulation
			}

			// render
			for (int i=0; i<count; i++)
			{
				const Planet &p = planets[i];
				float vol = p.mass;
				float radius = pow (vol * 3.f / (4.f * float(PI)), 1.f/3.f);
				RenderCircle (radius, p.com, sVec3(0), 1,1,1);
			}
		}

hmmm… it raises some interesting philosophical questions:

Does space time exist so the universe can reduce it's timestep at spots where mass is large, e.g. nearby a black hole? Could we form a universal theory around computational power? Physicians do so regarding ‘information’, so why not regarding ‘required calculations’?

Can we eventually build a quantum Perpetuum Mobile by cheating laws of physics to raise energy somehow because the universe can't do some infinite small timestep to keep up at some point?

But also: Why do i always assume knowledge from looking at printed equations, but then, after i use them in computer simulations, i feel more like still knowing even less than nothing about anything?

hehe, at least enough nonsense for a background science fiction story… : )

Yes JoeJ; very observant of you.

it is all about processes. As you fall into a black hole, your process is interrupted more and more, until you reach the event horizon, where you become fully assimilated.

Time dilation is the interruption of processes.

Time is not necessarily discrete.

The standard RK4 integrator is great, compared to Euler integration, but it doesn’t necessarily conserve energy. For that, you can rely on symplectic integrators.

jonny-b said:
I think Rk4 would improve the accuracy but ultimately if the last frame on one side and the first frame on the other aren't identical velocity, acceleration, and position will always be off.

I don't understand what kind of experiment you try to do. What do you mean with ‘side’?
In earlier posts it sounded you have a static body, e.g. a sun, which has mass but can not move?
Physics simulators often use 1/mass instead mass, and set this to zero for static bodies. This way all the acceleration goes just to the dynamic body. But for gravitational force calculations this would not work because force would become zero for both bodies.
So i initially i would not support static bodies until it works for only dynamic bodies. Sounds easier to get started.

taby said:
The standard RK4 integrator is great, compared to Euler integration, but it doesn’t necessarily conserve energy. For that, you can rely on symplectic integrators.

I did a whole physics engine around the idea of verlet integration when i was young. This was the time where i realized my missing math background will cause me some problems in the future. : )

But i doubt it could guarantee to conserve energy here either. Maybe it tends to remove energy instead adding it, which might be good in this case.
The main problem to me seems the function or signal (not sure how to call it) we try to approximate here. Any integrator we use treats such function as piece wise linear, but this will always fail to hit those steep spikes of forces local maxima in time. Midpoint / subdivision methods will reduce the error, but it will be still big.

There are some fluid simulators using variable timesteps, so some regions in space have more accuracy than others where not much goes on. Maybe that's required to simulate galaxies in practice.

Another problem i see is the simple model of using point masses or rigid bodies. That's fine at astronomical distances, but once planets come close we might need to represent them with increasing numbers of particles to get a more accurate representation of gravitational force.
Not sure, but i have a similar problem in my global illumination solver. It represents the surface with surfels, so oriented disks. The form factor calculation to get energy transfer between two disks is based on distance, area, and angle of the disks, and you can find the formula in many books. But is not accurate. If the disks are close to each other, the standard formula misses the effect of perspective, so the disk from the perspective of another is no longer an ellipse, but a distorted ellipse, and there is no analytical way to calculate visible area precisely. This causes some error, and i guess the same is true for gravity between spheres.

Seems a surprisingly difficult problem. But the universe is sparse, so the problem i see can probably addressed well enough with some hack.

I just tried to add the planets radii to distance, and this already solves the explosion problem for me:

			auto CalcPlanetRadius = [](const Planet &p)
			{
				float vol = p.mass;
				float radius = pow (vol * 3.f / (4.f * float(PI)), 1.f/3.f);
				return radius;
			};

			// function to update forces of a pair of planets
			auto GravityPair = [&](Planet &bodyI, Planet &bodyJ)
			{
				const float G = 1.f; // some number to relate your mass metric to gravity. i just assume we use a mass unit so G becoems 1.

				sVec3 diff = bodyJ.com - bodyI.com;

				float dist = diff.Length();
				float hackedDist = dist + (CalcPlanetRadius(bodyI) + CalcPlanetRadius(bodyJ)) * .5f;
				float force = G * (bodyI.mass * bodyJ.mass) / (hackedDist * hackedDist);

				sVec3 direction = diff / dist;
				bodyI.force += direction * force;
				bodyJ.force -= direction * force;
			}; 

It's wrong, but it won't affect astronomical distances much.

jonny-b said:
I started out with a simple example to teach my self the basics and get the simulation working. I figured a static body that doesn't move and a dynamic body that moves from “left” to “right” (this is what I meant by that) according to the laws of motion.

The simpler example would be to make both bodies dynamic. Because static bodies do not exist in the real world, you add the complexity to define faked laws of physics handling impossible cases, and so proofing for correctness also becomes harder.

jonny-b said:
Here is an example of what I'd like to build.

This really looks like the results i get from my code example, so you could use that for reference.

Looking at your snippets, maybe your'e not used yet to work with 2d / 3d vectors? That's quite easy, but just ask if you have any questions about details.
There are also libraries for JavaScript which have that, together with matrices and all the basic math stuff, e.g. Three.js (not sure if this also supports 2D).
I guess they also help with performance, which will quickly become a problem. But ofc. you can do it all yourself at least initially. Anything where i use sVec3 would just become 3 lines of code instead one (or 2 for 2D).

jonny-b said:
My ultimate goal is to create a bunch of stars on a 2d plane that all interact with each other.

This gives us a O(N^2) performance problem, so that's not real time for a larger number of stars.
Once the simulation works, this becomes the really interesting problem.
The obvious solution would be to build a tree over the stars. The internal nodes store averaged position and summed mass of all its child nodes. Then, for distant sections of the universe, you can interact with just one internal node, instead of iterating all of it's 1000 child stars. It's an approximation, but you have good control over error, and i'm sure any scientific simulation does this.
But we can still go much faster than that, which becomes pretty advanced then: You could interact with the distant galaxy not for each star, but instead you do this interaction only at some parent node of that current star. This way the interaction work and result is shared for all of the child nodes, which may be again 1000 stars.
Combining those two, we just got a speedup of 1000*1000 for interacting our current galaxy with a distant galaxy. That's nice, so we will accept the small error it causes. The error can be further reduced if we use some spational filter, e.g. a bilinear texture lookup if our tree would be a regular multi level grid. (Just, because the universe is so sparse, a simple multi level grid is not really ideal and a waste of memory and processing)

Barnes-Hut exists to reduce the complexity to below O(n^2).