Let's add an approximations to take into account some of the [General Relativity][1] (GR) effects — at least for bodies orbiting close to the massive Sun — and start to look at $J_2$ the lowest order multipole term for a body's gravitational potential beyond the monopole term $-GM/r$. **Newton:** The acceleration of a body in the gravitation field of another body of [standard gravitational parameter][2] $GM$ can be written: $$\mathbf{a_{Newton}} = -GM \frac{\mathbf{r}}{|r|^3},$$ where $r$ is the vector from the body $M$ to the body who's acceleration is being calculated. Remember that the acceleration of each body depends only on the mass of *the other body*, even though the force depends on both masses, because the first mass cancels out by $a=F/m$. **General Relativity (approximate):** Although I'm not familliar with GR, I'm going to recommend an equation that seems to work well and seems to be supported by several links. It is an approximate relativistic correction to Newtonian gravity that is used in orbital mechanics simulations. You will see various forms in the following links, most with additional terms not shown here: - https://physics.stackexchange.com/q/313146/83380 - Eq. 1 in https://www.lpi.usra.edu/books/CometsII/7009.pdf - Eq. 27 in https://ipnpr.jpl.nasa.gov/progress_report/42-196/196C.pdf - Eq. 4-26 in https://descanso.jpl.nasa.gov/monograph/series2/Descanso2_all.pdf - also see David Hammen's discussion/advice in [this answer][3]. - Eq. 3.11 in http://adsabs.harvard.edu/full/1994AJ....107.1885S (click "print") see [this answer][4] - Eq. A.8 in G. F. Rubilar1 & A. Eckart 2001 [Periastron shifts of stellar orbits near the Galactic Center](https://www.aanda.org/articles/aa/abs/2001/28/aah2613/aah2613.html), except that the $\mathbf{r}/r^3$ Newtonian is included rather than treated separately. - Eq. 2 in M. Parsa et al. 2017 [Investigating the Relativistic Motion of the Stars Near the Supermassive Black Hole in the Galactic Center](https://www.eso.org/public/archives/announcements/pdf/ann17051a.pdf), except that the $\mathbf{r}/r^3$ Newtonian is included rather than treated separately. The following approximation should be *added to* the Newtonian term: $$\mathbf{a_{GR}} = GM \frac{1}{c^2 |r|^3}\left(4 GM \frac{\mathbf{r}}{|r|} - (\mathbf{v} \cdot \mathbf{v}) \mathbf{r} + 4 (\mathbf{r} \cdot \mathbf{v}) \mathbf{v} \right),$$ **Oblateness ($J_2$ only, approximate):** I'm just using the math from [Wikipedia's article on the Geopotential Model][5] with a very important-to-remember approximation; I am assuming that the oblateness is in the plane of the ecliptic — that the rotational axis of the oblate body is in the $\mathbf{\hat{z}}$ direction, perpendicular to the ecliptic. **Don't forget that this is an approximation!** The full vector equations are messier than this, I'll try to come back and update this once I'm sure I've got it correct. In the mean time, here is an approximation: $$\mathbf{r} = x \mathbf{\hat{x}} + y \mathbf{\hat{y}} + z \mathbf{\hat{z}} $$ $$a_x = J_2 \frac{x}{|r|^7} (6z^2 - 1.5(x^2+y^2)) $$ $$a_y = J_2 \frac{y}{|r|^7} (6z^2 - 1.5(x^2+y^2)) $$ $$a_z = J_2 \frac{z}{|r|^7} (3z^2 - 4.5(x^2+y^2)) $$ The following should be *added to* the Newtonian term: $$\mathbf{a_{J2}} = a_x \mathbf{\hat{x}} + a_y \mathbf{\hat{y}} + a_z \mathbf{\hat{z}} $$ **Tidal Forces:** It gets even more complicated when looking at terms that involve non-spherical mass distributions in both bodies at the same time, whether they are static or not. At this point it's probably necessary to hit the books. ---- Here's a test run. I've compared to downloaded data from [JPL's Horizons][6]. For the outer planets I use the Horizons data for each planet's barycenter, which makes sure it's the velocity for the center of mass of the planet plus all of its moons. I haven't added the correction to the planet's masses, but that's a much smaller effect since it only affects the movement of other, distant bodies. For the Earth data, make sure to download the Earth's geocenter and the Moon separately (not the Earth-Moon barycenter). For the outer planets remember to download the barycenters. [![Screenshot of JPL Horizons - Earth][7]][7] [![Screenshot of JPL Horizons - Jupiter][8]][8] I've integrated for a year, and everything is within about one kilometer of the Horizons data except for Earth's Moon. That's not a surprise considering all the intimate tidal effects between these two. Adding Earth's $J_2$ to the potential felt by the Moon only helps partially, it's really not the right way to do it, especially considering that the Earth's axis (and therefore oblateness) is so far away from the ecliptic. So this is overall an illustration that the more physics you put in, the closer you can get to a really serious JPL simulation! This is the absolute distance between the simulated positions here and the Horizons output from `2017-01-01 00:00` to `2018-01-01 00:00`. Following that is the Python script I used to generate it. [![Screenshot of Python Output][9]][9] ---- Based on discussion of stiffness in comments below, here's a quick plot of step size versus time. I think the stiffness is coming from the Earth-Moon system, these frequent excursions look a bit like the Earth and Moon's error excursions. I think I am going to try to rescale the problem to dimensionless units, right now the numerical tolerance applies to all velocities and positions equally, not a good idea. [![enter image description here][10]][10] def deriv_Newton_Only(X, t): x, v = X.reshape(2, -1) xs, vs = x.reshape(-1, 3), v.reshape(-1, 3) things = zip(bodies, xs, vs) accs, vels = [], [] for a, xa, va in things: acc_a = np.zeros(3) for b, xb, vb in things: if b != a: r = xa - xb acc_a += -b.GM * r * ((r**2).sum())**-1.5 accs.append(acc_a) vels.append(va) return np.hstack((np.hstack(vels), np.hstack(accs))) def deriv_sunGRJ2EarthJ2(X, t): x, v = X.reshape(2, -1) xs, vs = x.reshape(-1, 3), v.reshape(-1, 3) things = zip(bodies, xs, vs) accs, vels = [], [] for a, xa, va in things: acc_a = np.zeros(3) for b, xb, vb in things: if b != a: r = xa - xb acc_a += -b.GM * r * ((r**2).sum())**-1.5 if a.do_SunGR and not a.name == 'Sun': a.flag0 = True r = xa - xs[0] v = va - vs[0] rsq = (r**2).sum() rm3 = rsq**-1.5 rm1 = rsq**-0.5 # https://physics.stackexchange.com/q/313146/83380 # Eq. 1 in https://www.lpi.usra.edu/books/CometsII/7009.pdf # Eq. 27 in https://ipnpr.jpl.nasa.gov/progress_report/42-196/196C.pdf # Eq. 4-26 in https://descanso.jpl.nasa.gov/monograph/series2/Descanso2_all.pdf # Eq. 3.11 in http://adsabs.harvard.edu/full/1994AJ....107.1885S term_0 = Sun.GM / (clight**2) * rm3 term_1 = 4.*Sun.GM * r * rm1 term_2 = -np.dot(v, v) * r term_3 = 4.*np.dot(r, v) * v accGR = term_0*(term_1 + term_2 + term_3) acc_a += accGR if a.do_SunJ2 and not a.name == 'Sun': a.flag1 = True r = xa - xs[0] # position relative to Sun x, y, z = r xsq, ysq, zsq = r**2 rsq = (r**2).sum() rm7 = rsq**-3.5 # https://en.wikipedia.org/wiki/Geopotential_model#The_deviations_of_Earth.27s_gravitational_field_from_that_of_a_homogeneous_sphere accJ2x = x * rm7 * (6*zsq - 1.5*(xsq + ysq)) accJ2y = y * rm7 * (6*zsq - 1.5*(xsq + ysq)) accJ2z = z * rm7 * (3*zsq - 4.5*(xsq + ysq)) accJ2 = J2s * np.hstack((accJ2x, accJ2y, accJ2z)) acc_a += accJ2 if a.do_EarthJ2 and not a.name == 'Earth': a.flag2 = True r = xa - xs[3] # position relative to Earth x, y, z = r xsq, ysq, zsq = r**2 rsq = (r**2).sum() rm7 = rsq**-3.5 # https://en.wikipedia.org/wiki/Geopotential_model#The_deviations_of_Earth.27s_gravitational_field_from_that_of_a_homogeneous_sphere accJ2x = x * rm7 * (6*zsq - 1.5*(xsq + ysq)) accJ2y = y * rm7 * (6*zsq - 1.5*(xsq + ysq)) accJ2z = z * rm7 * (3*zsq - 4.5*(xsq + ysq)) accJ2 = J2e * np.hstack((accJ2x, accJ2y, accJ2z)) acc_a += accJ2 accs.append(acc_a) vels.append(va) return np.hstack((np.hstack(vels), np.hstack(accs))) import numpy as np import matplotlib.pyplot as plt from scipy.integrate import odeint as ODEint names = ['Sun', 'Mercury', 'Venus', 'Earth', 'Moon', 'Mars', 'Ceres', 'Pallas', 'Vesta', 'Jupiter', 'Saturn', 'Uranus', 'Neptune'] GMsDE430 = [1.32712440040944E+20, 2.203178E+13, 3.24858592E+14, 3.98600435436E+14, 4.902800066E+12, 4.2828375214E+13, 6.28093938E+10, 1.3923011E+10, 1.7288009E+10, 1.267127648E+17, 3.79405852E+16, 5.7945486E+15, 6.83652719958E+15 ] # https://ipnpr.jpl.nasa.gov/progress_report/42-178/178C.pdf # for masses also see ftp://ssd.jpl.nasa.gov/pub/xfr/gm_Horizons.pck # https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/ # https://naif.jpl.nasa.gov/pub/naif/JUNO/kernels/spk/de436s.bsp.lbl # https://astronomy.stackexchange.com/questions/13488/where-can-i-find-visualize-planets-stars-moons-etc-positions # https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/jup310.cmt # ftp://ssd.jpl.nasa.gov/pub/xfr/gm_Horizons.pck GMs = GMsDE430 clight = 2.9979E+08 # m/s halfpi, pi, twopi = [f*np.pi for f in [0.5, 1, 2]] # J2 values J2_sun = 2.110608853272684E-07 # unitless R_sun = 6.96E+08 # meters J2s = J2_sun * (GMs[0] * R_sun**2) # is there a minus sign? J2_earth = 1.08262545E-03 # unitless R_earth = 6378136.3 # meters J2e = J2_earth * (GMs[3] * R_earth**2) # is there a minus sign? JDs, positions, velocities, linez = [], [], [], [] use_outer_barycenters = True for name in names: fname = name + ' horizons_results.txt' if use_outer_barycenters: if name in ['Jupiter', 'Saturn', 'Uranus', 'Neptune']: fname = name + ' barycenter horizons_results.txt' with open(fname, 'r') as infile: lines = infile.read().splitlines() iSOE = [i for i, line in enumerate(lines) if "$$SOE" in line][0] iEOE = [i for i, line in enumerate(lines) if "$$EOE" in line][0] # print name, iSOE, iEOE, lines[iSOE], lines[iEOE] lines = lines[iSOE+1:iEOE] lines = [line.split(',') for line in lines] JD = np.array([float(line[0]) for line in lines]) pos = np.array([[float(item) for item in line[2:5]] for line in lines]) vel = np.array([[float(item) for item in line[5:8]] for line in lines]) JDs.append(JD) positions.append(pos * 1000.) # km to m velocities.append(vel * 1000.) # km/s to m/s linez.append(lines) JD = JDs[0] # assume they are identical class Body(object): def __init__(self, name): self.name = name bodies = [] for name, GM, pos, vel in zip(names, GMs, positions, velocities): body = Body(name) bodies.append(body) body.GM = GM body.daily_positions = pos body.daily_velocities = vel body.initial_position = pos[0] body.initial_velocity = vel[0] x0s = np.hstack([b.initial_position for b in bodies]) v0s = np.hstack([b.initial_velocity for b in bodies]) X0 = np.hstack((x0s, v0s)) ndays = 365 nperday = 144 time = np.arange(0, ndays*24*3600+1, 24*3600./nperday) days = time[::nperday]/(24*3600.) for body in bodies: body.do_SunGR = False body.do_SunJ2 = False body.do_EarthJ2 = False body.flag0 = False body.flag1 = False body.flag2 = False Sun, Mercury, Venus, Earth, Moon, Mars = bodies[:6] Ceres, Pallas, Vesta = bodies[6:9] Jupiter, Saturn, Uranus, Neptune = bodies[9:] Mercury.do_SunGR = True Venus.do_SunGR = True Earth.do_SunGR = True Moon.do_SunGR = True Mars.do_SunGR = True Ceres.do_SunGR = True Pallas.do_SunGR = True Vesta.do_SunGR = True Mercury.do_SunJ2 = True Moon.do_EarthJ2 = True rtol = 1E-12 # this is important!!! answer, info = ODEint(deriv_sunGRJ2EarthJ2, X0, time, rtol = rtol, full_output=True) print answer.shape nbodies = len(bodies) xs, vs = answer.T.reshape(2, nbodies, 3, -1) for body, x, v in zip(bodies, xs, vs): body.x = x body.v = v body.x_daily = body.x[:, ::nperday] body.v_daily = body.v[:, ::nperday] body.difference = np.sqrt(((body.x_daily - body.daily_positions.T)**2).sum(axis=0)) if True: for body in bodies[:6]: print body.name, body.flag0, body.flag1, body.flag2 if True: plt.figure() for i, body in enumerate(bodies[:12]): # Sorry Neptune!!! plt.subplot(4, 3, i+1) plt.plot(days, 0.001*body.difference) plt.title(body.name, fontsize=14) plt.xlim(0, 365) plt.suptitle("calc vs JPL Horizons (km vs days)", fontsize=16) plt.show() [1]: https://en.wikipedia.org/wiki/General_relativity [2]: https://en.wikipedia.org/wiki/Standard_gravitational_parameter [3]: https://space.stackexchange.com/a/21472/12102 [4]: https://astronomy.stackexchange.com/a/25304/7982 [5]: https://en.wikipedia.org/wiki/Geopotential_model#The_deviations_of_Earth.27s_gravitational_field_from_that_of_a_homogeneous_sphere [6]: https://ssd.jpl.nasa.gov/horizons.cgi [7]: https://i.sstatic.net/d1KPF.png [8]: https://i.sstatic.net/Lx7b0.png [9]: https://i.sstatic.net/l32NY.png [10]: https://i.sstatic.net/FVVS9.png