<smartass>

And that's why its orbit precesses because of Schwarzschild spacetime ...

</smartass>

?

Actually, the majority of the precession of the perihelion of Mercury occurs because of interaction with the other planets. The curvature of time adds only a small precession.

Correct. ~43 arcseconds/century are on the expense of spacetime curvature.

In any case, the problem was the integrator. Using RK or symplectic integrators did the trick.

2 hours ago, Green_Baron said:

Correct. ~43 arcseconds/century are on the expense of spacetime curvature.

Actually, the kinematic time dilation also affects the precession — it’s not just about the curvature of spacetime. The precession is proportional to 3GM, where it’s one part kinematic time dilafion and two parts gravitational time dilation.

I was really questioning the original item's mention of circular orbits. If a circular orbit was expected, or if that is a suitable approximation for Mercury's orbit in your application, then why is integration considered necessary? In that case the radius and angular velocity would be constant and the angular position would simply be proportional to time. Keep it simple.

Sorry. I wanted to fix the circular orbit before I added eccentricity into the mix. Using either the Runge-Kutta integrator or the symplectic integrator fixed the problem. Note that using Kepler's laws for a parametric ellipse instead of integration is also possible. I used integration, so be it.

P.S. cowcow is a temporary account that I used because I was somehow locked out of this taby account for a while.

http://farside.ph.utexas.edu/teaching/336k/Newtonhtml/node116.html

https://www.mathpages.com/home/kmath527/kmath527.htm

The trouble with integration methods for predicting the path of an object is that eventually errors increase to the point where no faith can be had in the results.
In 1987 Bretagnon and Francou, of the Bureau des Longitudes in Paris took a different approach and produced a method called VSOP87 which enables the heliocentric coordinates of all the major planets to be calculated from long series of terms. As I understand it they did harmonic analysis of historically long series of observations, effectively producing something akin to Fourier series for the calculations.
Jean Meeus' book "Astronomical Algorithms" explains much more about this, giving the algorithm and discussing accuracy. He includes all the necessary data as an appendix to the book. See https://www.willbell.com/math/mc1.htm
I obtained the book in 2010 and implemented the algorithm for all major planets for a certain web site. I originally wrote it in Java for client-side applets but when Oracle ceased to keep Java secure enough in browsers, I converted it all to server-side PHP (yuk!).
The original paper is referenced in Jean's book:
P.Bretagnon, G.Francou, "Planetary theories in rectangular and spherical variables. VSOP87 solutions", Astronomy and Astrophysics, Vol. 202, pp 309-315 (1988).
Searching today I discovered that Jean's book is available online as a PDF (488 pages) from http://www.agopax.it/Libri_astronomia/pdf/Astronomical Algorithms.pdf

Thanks for pointing me to MathPages.com - that looks excellent.