1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Trouble with a Solar System simulation

  1. Jan 6, 2013 #1
    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?
  2. jcsd
  3. Jan 9, 2013 #2
    Wow, so I'm the only one in the whole world who's doing this kind of thing? I feel so... special.
  4. Jan 10, 2013 #3


    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    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.)
  5. Jan 10, 2013 #4


    User Avatar
    Science Advisor
    Homework Helper

    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!
  6. Jan 10, 2013 #5


    User Avatar
    2017 Award

    Staff: Mentor

    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?
  7. Jan 19, 2013 #6
    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 [itex]{v^{2} \over c^{2}}[/itex] ([itex]\vec{r}[/itex] . [itex]\vec{v}[/itex]) [itex]\vec{v}[/itex] ? (which I need to do to calculate forces in 3 dimensions)
    I am now interpreting it as 4 [itex]{v^{2} \over c^{2}}[/itex] (rv cos ∅) [x, y, or z of [itex]\vec{v}[/itex]], where
    ∅ is, of course, the angle between [itex]\vec{r}[/itex] and [itex]\vec{v}[/itex]

    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.
    Last edited: Jan 19, 2013
  8. Jan 20, 2013 #7


    User Avatar
    2017 Award

    Staff: Mentor

    $$(\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.
  9. Jan 20, 2013 #8
    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 [itex]\vec{v}^{2}[/itex], which doesn't even make sense. Thanks! I will report back when/if I get this running.
    Last edited: Jan 20, 2013
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook