Simulating orbits using laws of physics

问题:

Over the past couple of weeks I've been trying to simulate orbits in a solar system simulation I am making as part of a University module. To cut things short, my simulation is written in C++ using the Ogre3D rendering engine. I have attempted to implement orbits using Newton's law of universal gravitation which made my planet head towards the sun in a straight line, pass through the sun and then come back to its starting position. I also tried the steps from 'Position as a function of time' section of this wikipedia article, but that did not work for me either.

I am driving the simulation with a simple Euler integration method. If anyone has any experience with this kind of simulation, or just generally knows a lot about these physics laws then any help or pointing me in the right direction would be greatly appreciated.

回答1:

You must give the planet initial velocity v = (vx, vy, vz) tangent to the desired orbit. If the position of the sun is s and planet is p, then there is always a force acting between the two: the one on the planet points toward the sun, vector t=(s - p) and vice versa. The magnitude of this force is g Ms Mp / (t dot t), where "dot" is the dot product, g is the standard acceleration due to gravity, and Ms, Mp are the respective masses.

If you are doing a detailed model where all bodies can exert pull on all other bodies, then the algorithm is to accumulate all the pairwise forces to get a single resultant force vector acting on each body (planet or sun). Otherwise you might settle for an approximation where only the sun pulls on planets and other forces are deemed too small to matter.

So the algorithm is:

Choose dt, the time step, a small interval.
Set initial positions and velocities 
    (for planets, velocity is tangent to desired orbit. Sun has velocity zero.)
loop
  Accumulate all forces on all bodies with current positions.
  For each body position=p, velocity=v, net resultant force=f, mass=m, 
    update its velocity: v_new = v + f / m dt
    update position p_new = p + 0.5 * (v + v_new) 
    v = v_new; p = p_new
  Render
end loop

As has been mentioned, Euler is simple but requires a very small time step to get even reasonable accuracy. Sometimes you can introduce just a tiny bit of drag in the system (multiply velocity by a factor just a tad below 1) to keep things stable where otherwise they blow up.

回答2:

Consult the project "Moving stars around", there is an old C and a modernized Ruby version. (And now a C++ version?)

Short advice: Eulers methods are bad for energy conservation. explicit Euler increases energy, implicit Euler reduces energy. Just check it on the phase space picture of the harmonic oscillator y''+y=0.

Use symplectic integrators, the most simple and famous is the Leapfrog or Verlet method that already Newton used to reason about planetary movement.

回答3:

You could try Runge Kutta 4, but there will still be some "drift" in the orbits. Keeping the total energy of the system constant is needed to stabilize the system, but I'm not sure how this is done. A better method for celestial mechanics is symplectic integrator.