You all are more than welcome to implement a symplectic or Runge-Kutta integrator. I have codes for those, just not integrated into this project. I must stick with Euler integration, for the sake of transparency and simplicity.

Orbits are stable in the face of perturbation when it's a 2-body problem, or effectively a 1-body problem. Only in 3-body+ does chaos make an appearance.

You'll know it's working when the precession angle is like 0 for Newtonian orbit where alpha = beta = 1.

Aressera said:
Angular velocity refers to the spin of the body

Yeah, but i think he only changes position of the planet according to the simulated angle around the sun, while there is no linear velocity or spin around the planets com. Thus angular velocity actually is the proper term, somehow. Also there is no error from approximating a circular trajectory with linear motion, since a circle is all it can do eventually. But not sure. Too much Einstein stuff which did not make sense to me. : )

Edit: Actually no. Looking at the code again there is regular integration of 3D velocity to the planets position. Sorry about the confusion.

The Einstein stuff is what makes the orbit shape deviate from the Newtonian ellipse.

The space around a black hole contracts, making travel through that space slower than through the space at infinite distance. Basically, Mercury dwells more when it's closer to the Sun, which causes more gravitation than in Newtonian gravity, which means the bend in the path is sharper than in Newtonian gravity, which causes an overall deviation from the ellipse.

Also, the strength of gravitation relies on the speed of the gravitated body. Just a smidgen.

The 4th-order symplectic integrator is like:

void proceed_symplectic4(custom_math::vector_3& pos, custom_math::vector_3& vel, long double G, long double dt)
{
	static double const cr2 = pow(2.0, 1.0 / 3.0);

	static const double c[4] =
	{
		1.0 / (2.0 * (2.0 - cr2)),
		(1.0 - cr2) / (2.0 * (2.0 - cr2)),
		(1.0 - cr2) / (2.0 * (2.0 - cr2)),
		1.0 / (2.0 * (2.0 - cr2))
	};

	static const double d[4] =
	{
		1.0 / (2.0 - cr2),
		-cr2 / (2.0 - cr2),
		1.0 / (2.0 - cr2),
		0.0
	};

	{
		const custom_math::vector_3 grav_dir = sun_pos - pos;
		const double distance = grav_dir.length();
		const double Rs = 2 * grav_constant * sun_mass / (speed_of_light * speed_of_light);

		const double alpha = 2.0 - sqrt(1 - (vel.length() * vel.length()) / (speed_of_light * speed_of_light));

		const double beta = sqrt(1.0 - Rs / distance);
		const double beta_truncated = truncate_normalized_double(beta);

		pos += vel * c[0] * dt * beta_truncated;
		vel += grav_acceleration(pos, vel, G) * d[0] * dt * alpha;
	}

	{
		const custom_math::vector_3 grav_dir = sun_pos - pos;
		const double distance = grav_dir.length();
		const double Rs = 2 * grav_constant * sun_mass / (speed_of_light * speed_of_light);

		const double alpha = 2.0 - sqrt(1 - (vel.length() * vel.length()) / (speed_of_light * speed_of_light));

		const double beta = sqrt(1.0 - Rs / distance);
		const double beta_truncated = truncate_normalized_double(beta);

		pos += vel * c[1] * dt * beta_truncated;
		vel += grav_acceleration(pos, vel, G) * d[1] * dt * alpha;
	}

	{
		const custom_math::vector_3 grav_dir = sun_pos - pos;
		const double distance = grav_dir.length();
		const double Rs = 2 * grav_constant * sun_mass / (speed_of_light * speed_of_light);

		const double alpha = 2.0 - sqrt(1 - (vel.length() * vel.length()) / (speed_of_light * speed_of_light));

		const double beta = sqrt(1.0 - Rs / distance);
		const double beta_truncated = truncate_normalized_double(beta);

		pos += vel * c[2] * dt * beta_truncated;
		vel += grav_acceleration(pos, vel, G) * d[2] * dt * alpha;
	}

	{
		const custom_math::vector_3 grav_dir = sun_pos - pos;
		const double distance = grav_dir.length();
		const double Rs = 2 * grav_constant * sun_mass / (speed_of_light * speed_of_light);

		const double alpha = 2.0 - sqrt(1 - (vel.length() * vel.length()) / (speed_of_light * speed_of_light));

		const double beta = sqrt(1.0 - Rs / distance);
		const double beta_truncated = truncate_normalized_double(beta);

		pos += vel * c[3] * dt * beta_truncated;
		//	vel += grav_acceleration(pos, vel, G) * d[3] * dt * alpha; // last element d[3] is always 0
	}
}

So, I am now using the symplectic integrator, and it allows me to run the simulation where dt = 0.1 instead of dt = 0.01. The result is produced quite a bit faster, thanks to this.

Clang compiler might help as well with perf. Easy to install for VS using its installer app, but maybe a 4GB download. On average it's 25% faster for me. But in one case it was 200 times faster. :O

Btw, what is your result? How do you judge how correct it is?

The desired solution gives a result of around 43 arc seconds per Earth century of precession. I get that with both the Euler and the symplectic integration.

I have thought about it a lot and I believe that it has to do with the quantization. so basically we are snapping the double to a float value. It is this snapping to of values that causes the correct precession.

I tried g++ on Ubuntu, but I’ll also try clang.

I tried clang++ on Ubuntu, but no big difference.

Oh, then there was a big misconception on my side. Still referring to this:

If I comment the casting out, I get an answer of 7.75. If the casting remains, I get an answer of 42.66. The analytical solution gives an answer of 42.94.

I thought you get this big difference within one function call.

But now i see you get it only after many integration steps.

taby said:
I have thought about it a lot and I believe that it has to do with the quantization. so basically we are snapping the double to a float value. It is this snapping to of values that causes the correct precession.

Yeah, that's what i meant with joking: ‘On the other hand, maybe reducing precision is needed to emulate quantized nature? :D’

But i don't believe it. I rather believe the ‘wrong’ result you get if not doing the cast is more correct than the seemingly accurate result from casting.

How can you be so sure that 43 is right for the model you use? Real world measurements i guess?

And still having your code ported, how long would i need to run it to replicate?
I did so only for seconds, and it did not print any numbers in that time. Would if i let it run longer?

Any thoughts on how to achieve snapping-to without casting?