3
$\begingroup$

Not 100% sure where this question belongs, since I'm not sure if the problem is code related or not.

I've written a program to model the solar system. I am now testing its accuracy. I've checked conservation of angular momentum, and it is perfectly conserved, as is energy - I'm using the velocity verlet algorithm. I've also plotted the path of the planets, their orbits are stable over 1000 years and are definitely circular (I'm sticking to modelling circular orbits).

The problem is that the planets aren't where they started after one orbital period. They do get back to where they started, just not at the right time. I have no idea why. And it renders a lot of my checks useless. Why might this happen? The effect is significant, about 0.2*10^11m per year. And yet after 1000 years the orbit is still circular, with roughly the same radius.

|cite|improve this question
$\endgroup$
  • 3
    $\begingroup$ What is your definition of orbital period? Mine is "the moment at which the planet gets back where it started", which would make the question moot. $\endgroup$ – Federico Poloni Jan 3 '16 at 16:07
  • $\begingroup$ @FedericoPoloni I assume the OP means that the orbital period in the simulation is not the same as the expected orbital period (computed analytically). Unfortunately I don't think there's enough information here to figure out what might be going wrong. $\endgroup$ – Doug Lipinski Jan 3 '16 at 17:54
  • 3
    $\begingroup$ @DougLipinski Agree. Hard to tell without more details, but my first guess is that he is not properly taking into account the fact that the sun moves, too (reduced mass and that kind of things). $\endgroup$ – Federico Poloni Jan 3 '16 at 18:04
  • $\begingroup$ @FedericoPoloni My planets move about their common centre of mass. They are modelled as point masses. It is very basic. The problem I'm having is that I used the known orbital periods to calculate what the velocities of the planets should be, using v = 2*pi*r/T, however this velocity doesn't actually result in planets completing a full orbit after one orbital period. $\endgroup$ – user13948 Jan 4 '16 at 15:05
  • $\begingroup$ @DougLipinski My simulation is incredibly basic. Circular orbits, point mass planets, only force acting is gravity. And the problem is exactly as you said, my orbital periods don't match the expected ones. I could post the code, but not sure how helpful that would be - specifically which info would be helpful? $\endgroup$ – user13948 Jan 4 '16 at 15:07
11
$\begingroup$

This is expected behavior with the Verlet algorithm. It is a symplectic integrator, which means that it will preserve quadratic invariants to within roundoff error -- thus the form planetary orbits is maintained quite accurately (I assume you mean to say that the orbits are elliptical, not circular). However, the method is still a numerical approximation and does introduce error. The errors introduced affect only the phase of each orbit. This means that the planets will move in the correct orbits, but not at the correct speeds. In order to reduce the phase errors, you can take smaller time steps -- but there is eventually a limit to this, since it will increase the amount of rounding error in the solution.

To learn more about these issues, I recommend the book Geometric numerical integration.

Studying the long-term stability of the solar system (over millions to billions of years) is a subtle and active research area that makes use of much more accurate numerical methods and even multiple-precision arithmetic.

|cite|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.