Trouble with a Solar System simulation

[tfr000]

Hi all.
I have been working, off and on, on a Solar System simulator (in MS Excel) for many months. I am using the Runge-Kutta-Fehlberg 8/9 8th-order integrator and JPL's initial conditions to try to reproduce an accurate model of planetary motions.
This, of course, involves modeling the relativistic effects of gravity. And therein lies my problem... using a couple of different post-Newtonian models of relativity, I seem to be getting exactly twice the effect I should be. In other words, Mercury's perihelion apparently advances by exactly twice the famous 43 arcsecond per century rate.
I have done a lot of checking... Newtonian-only effects are modeled perfectly vs analytical. I ran some test cases of two-body motion and compared the perihelion advance to analytical... exactly 2x.
Obviously, I seem to be applying these models incorrectly. And the fact that a couple of different post-Newtonian models give me exactly the same result points to something I am doing wrong, rather than the models themselves.
I am, however, baffled about what that something might be. Help?

 

[tfr000]

Wow, so I'm the only one in the whole world who's doing this kind of thing? I feel so... special.

 

[haruspex]

Presumably the error is either in the equations or in the implementation of them. Without seeing either it's rather hard to comment. (I'm not saying I would comment if you were to post them... but somebody else might.)

 

[AlephZero]

As Haruspex said you haven't given us anything much to work with.

Assuming there is only one error, either your differential equations are wrong or your integration scheme is wrong. You said you tried two different versions of the equations and got the same wrong answer, so that sugggests it's your integration scheme.

I would try the "standard" 4/5 order RK and see what you get. It's not impossible there was a typo in the source of your data for the 8/9 order scheme!

 

[mfb]

I would be surprised by an integration scheme error which exactly doubles relativistic corrections, but does not influence Newtonian motion. Do you get the correct precession rate from other planets (531+-1 arcsec/century according to Wikipedia)?

The factor 2 reminds me of the bending of light, where a naive calculation gives 1/2 the true value because it neglects the modified time coordinate. Maybe you have a similar effect here?

 

[tfr000]

Oops, never mind... I wrote a new spreadsheet which calculates the Newtonian forces around an orbit analytically, and then plugged in the different relativity approximations for direct comparison with each other. This showed me that I mis-interpreted the notation of the relativity equations, I think. Also, one of the three methods gives very strange results, while the other two agree. I have references for all of the methods, if you'd like to look at them.

So my new question...
How would you resolve 4 {v^{2} \over c^{2}} (\vec{r} . \vec{v}) \vec{v} ? (which I need to do to calculate forces in 3 dimensions)
I am now interpreting it as 4 {v^{2} \over c^{2}} (rv cos ∅) [x, y, or z of \vec{v}], where
∅ is, of course, the angle between \vec{r} and \vec{v}

My errors in the numerical integration must have been an odd coincidence.
BTW, I trust the integrator, because it gives nearly perfect results vs analytical for Newtonian-only forces, even after 100 or more orbits.

 

[mfb]

$$(\vec{r} \cdot \vec{v})\vec{v} = (r_1 v_1 + r_2 v_2 + r_3 v_3)\left( \begin{array}{c} v_1\\ v_2 \\ v_3\end{array}\right) = \left(|r|\, |v| cos(\theta)\right) \,\left( \begin{array}{c} v_1\\ v_2 \\ v_3\end{array}\right) $$

I don't know what you mean with "or", but I think your conversion is correct.

 

[tfr000]

mfb - just how my brain works. What you wrote is what I meant. Not sure what I had read into that notation before - something with \vec{v}^{2}, which doesn't even make sense. Thanks! I will report back when/if I get this running.