Why do objects in a gravity simulation experience sudden large accelerations?

Question

I'm trying to create a simple program that simulates gravity. The idea is that I have one central sun and several planets that I can create with a swipe gesture on the screen, and I use the initial swipe to provide the planets with an initial velocity.

After some time the planet will eventually start to move around the sun according to Newton's law.

My current approach is: At any moment I compute the value of the gravitational force that the sun exerts on the planet ($Mp$ is the mass of the planet and $Ms$ is the mass of the sun):

$$F = G \frac{Mp Ms}{r^2}$$

Then I find the acceleration value for the planet with

$$A = \frac{F}{Mp}$$

Then I find the angle, $\theta$, of the line that connects the center of the planet to the center of the sun.

Next, I create an acceleration vector along that line by creating its x and a y component:

$$\begin{align} Ax &= A\cos\theta & Ay &= A\sin\theta \end{align}$$

At the next iteration I use this acceleration vector and the time elapsed since the previous iteration to compute the planet's velocity and then its position. Then the whole thing repeats.

The main problem with this approach is that as soon as the planet reaches a distance close to zero from the sun, it gets a tremendous acceleration and at the next iteration it is simply too far from the sun and just keeps moving along a straight line, right out of the screen, which of course is not what I expect from gravity. Note that I still don't have collision detection so what I'd expect would be the planet to remain kind of still at the center of the sun.

My intuition is shouting that I should use some kind of integration for the acceleration so every acceleration that I miss between one iteration and the next would be taken into account and my planet would stop escaping gravity, so I recovered my math and physics books and tried to figure this out by myself, but no luck.

If I got this right the problem is that my acceleration is a function of distance, so I cannot integrate it in order to get the position for the planet, because that would require an acceleration as a function of time. Am I right? What's the right solution to this?

Answer 1

Brionus has touched on the key - adaptive time steps. When you start getting large accelerations, reduce the size of your time increments. Also, when you are not accelerating much, increase the size.

A fairly standard way to do this is to calculate your position change over one step. Then cut the step in half and, starting from the same starting point, compute 2 successive position steps. Compare the two final position changes, and if they differ by some predetermined factor (let's say 10^-6), replace the original time step with the smaller step, and do the calculation all over again. If the two steps matched closely, try a computation with a time step twice the original.

This takes a lot of extra computation, but it produces a simulation which is neither too precise nor not precise enough. For orbital simulations, the large time steps used during low-gravity trajectories will more than compensate for the extra computation time.

EDIT - In response to the request for "the calculus approach":

The calculation approach which you used will work for a really simple simulator, but it ignores the interactions between the non-sun bodies, as well as assuming that the central sun does not move. To handle this, use a more generalized approach. Store the (x,y,z) (Vx,Vy,Vz) and (Mx,My,Mz) values for your bodies in arrays indexed the same way. Then the gravitational attractions between any two bodies will simply be calculated as (Fx,Fy,Fz), where $$F_n =\frac{\Delta n}{\sqrt{{\Delta x}^2 + {\Delta y}^2 + {\Delta z}^2}} \frac{GM_a M_b}{{\Delta x}^2 + {\Delta y}^2 + {\Delta z}^2}$$ Calculate each component separately, and integrate separately to get new velocities and positions. Also note that you will need to calculate $N(N-1)$ values for N bodies (including the sun) but that body A pulls on body B exactly as hard as body B pulls on body A, so you only need to do the numerical computation half of the apparent total, although you'll need to be careful to reverse the sign to get the other half of the values right. (The reason it's $N(N-1)$ rather than $N^2$ is that you don't calculate the pull of a body on itself.)

For a simple simulator, you can use Euler integration. Given the force on a body and its mass, for a very small $\Delta t$ you can say $$\Delta V = \frac{F \Delta t}{M}$$ and when you know both the original velocity V and the change in velocity, the change in position $\Delta P$ is$$\Delta P = {(V + \Delta V)\Delta t}$$

For shortish time scales and very small $\Delta t$ this will work, but you are fundamentally trying to approximate an irregular curve with straight segments and in the long run you'll see large and growing errors. So for longer periods you'll want to get into more sophisticated algorithms. As zeldridge has commented, Runge-Kutta is a well-known alternative.

SECOND EDIT - Also, when updating your values, do each set of calculations based on the same original conditions for all bodies. That is, if you're calculating Vnew and Pnew for body A, do not calculate the results for body B using the updated Pnew for body A. Calculate an entire new array of V and P values, then replace the old values as a block.

Answer 2

If you really want a general gravitation simulator (i.e. one that will handle more then two bodies), then there are methods for reducing the error involved in the simulation, but there aren't any methods for eliminating the error. Below are a few approaches - none of these approaches are perfect, since there's a balance between physical accuracy, programming & computational speed, and reducing the escaping behavior.

1. Adaptive time steps - when your particle is experiencing large changes in acceleration (i.e. when it's moving fast near a massive body), you can reduce the time step $dt$ so you're "skipping" less time. This will increase the physical accuracy of the simulation, but it will also slow down the simulation (and slow it down unevenly), so it's not a good approach if you want it to look good in real time.

2. Hard spheres - make your masses hard spheres that bounce off each other. This will reduce the buildup of error, since the masses won't get as close to each other. This can increase the physical accuracy of the simulation, unless you object to planets that bounce like rubber balls (which is a fair objection, I suppose - you could also make them inelastic spheres, which is probably more accurate for planets).

3. Speed limit - you could just program in a hard speed limit to limit escaping behavior. Whenever the velocity exceeds the speed limit, resize the velocity magnitude back to the speed limit. This isn't very physically accurate, and could result in some strange looking behavior, but it's easy to program, and will reduce the escaping of your masses.

4. Conservation of energy - at each time step, calculate the total amount of gravitational and kinetic energy all the masses have. At each time step, if the total amount of energy has changed, artificially adjust the velocities of the masses so the total amount of energy stays the same. This is also not perfectly accurate, but it does maintain fidelity with one physical law, and it will reduce the escaping behavior.

If you'd like help understanding the implementation of one of these methods, I can explain in more detail.

Answer 3

If you are not interested in a full n-body simulation and will accept one where the sun is much more massive than all the planets, you can simplify things a lot. We have an analytic solution to the two body problem, so you can apply the sun's gravity that way. For each planet, given its position and velocity, you can compute its angular momentum and energy, then get the orbit. This lets you compute the location of the planet at the end of the time step ignoring the attractions of all the planets. To zeroth order, you are done, because you ignore the planet-planet interactions. If you want to apply the planet-planet interactions, you can compute the force on planet $i$ from all the other planets and apply the acceleration for the entire time step. As the force is much smaller than the sun's gravity, the change in position will be small. This will keep your errors much smaller, roughly as the ratio of the forces on the planet from the other planets to the force from the sun. You can then apply the adaptive time steps and energy conservation that the other responders suggest.

Answer 4

A different programming approach in a step-by-step animation is to assign position and momentum variables to each of the objects. In programming terms we have object name, mass, location x/y/z, and momentum px/py/pz. Between each frame, for each object, first calculate its distance to every other object and calculate the gravitational force on that object based on Newton’s inverse force law; second, add the force to the momentum; finally, add the momentum to the location (px=px+forcex, x=x+px).

You can get very good animations with this method, where the accuracy is only bounded by the accuracy of the calculations, the frequency of the calculations and the starting position/momentum of each of the objects. An example of this method can be found at http://www.animatedphysics.com/planets/Moon_Orbit.htm

WRT your problem "The main problem with this approach is that as soon as the planet reaches a distance close to zero from the sun, it gets a tremendous acceleration and at the next iteration it is simply too far from the sun and just keeps moving along a straight line, right out of the screen". This is a very significant problem. If you allow one object to get arbitrarily close to a second object, the force between the objects will get arbitrarily large (hence the object suddenly jumps out of the field of view as you experienced) no matter how small the timesteps you use.

I think it would be very illuminating if @Brionius could expand on his option 4.

Answer 5

Based on the physical model you're adopting (that only the "sun" exerts gravitational force on the "planets" and that we ignore the gravitational effects of the planets on the sun or on each other), which is certainly a decent enough starting point (not to mention good enough for Kepler w.r.t our own solar system), there are two major points to realize:

1) There ARE legitimate orbits (or, more properly, let's call them trajectories) in which the planet is not bound to the sun! That is, the planet, rather than moving in a path that stays near the sun and revolves periodically around it, can, in complete consistency with Newton's laws, be set on a course which takes it out of the solar system, never to return (and this can occur even after the planet makes an initial close pass near the sun which might lead one to believe it was going to orbit). Nor, given arbitrary initial velocities for the planets, should we expect these trajectories to be rare or unusual, though they may not be the ones you are interested in. Solving Newton's laws of motion together with his law of gravity tells us that the planet's trajectory will be one of three types: elliptical, parabolic, or hyperbolic. Only the first of these, elliptical trajectories, are closed, or bounded. Now, it is quite possible, even likely, that the "runaway" trajectories you observe really are due to the numerical artifacts of a insufficiently precise simulation algorithm, but it would likely behoove you to confirm that this is true as a first step by checking that the initial velocities you're assigning the planets are small enough that they don't put the planet on a parabolic or hyperbolic trajectory. If you want to completely rule out such trajectories, as it seems you do, it would be smart to write your code in such a way that prevents the user from choosing initial velocities that would lead to them. The test is easy; just demand that $v^2 < \sqrt {G M / r}$, where $v$ is the planet's initial total speed and $r$ is the sun-planet distance at the moment when that initial speed is measured. You might want to go step further and demand that the eccentricity of all elliptical orbits is small enough that the planet does not leave the screen (as it might, for example, for a Halley's comet-like object which, though bound to our sun, varies its distance to the sun by a factor of 70 over its 76 year orbital period). By requiring that $a \equiv \frac{1}{\frac{2}{r}-\frac{v^2}{G M}}$, the so-called "semi-major axis" of the orbit, is smaller that the distance from the sun to the edge of the screen, you can be sure that any planets that leave the screen are never coming back, and represent either parabolic/hyperbolic trajectories, or if these have been excluded, numerical artifacts.

2) Perhaps more important, because it has the potential to greatly improve the simplicity, accuracy, and efficiency of your code, is that Newton's laws, under your physical assumptions, are completely solvable, meaning that there is no need to simulate the physics (!). It is called the two-body problem (despite there being multiple planets, we are assuming they don't interact with one another, so each planet's path is completely determined by its interactions with exactly one other object, the sun), and it's essentially the calculation that Kepler did to show that the planets move in elliptical orbits. It's somewhat complicated to write down the $x$ and $y$ positions as a function of time given an arbitrary initial position and velocity for the planet, so I won't do it here, but the Wikipedia article on "Kepler orbits", and links therein provide all the equations one would need to massage to get the desired result (note that there it is assumed that we take the orbit's perihelion/aphelion direction to be aligned with one of the coordinate axes, and one would need to perform a change of basis to recover the more general result). I'm actually an employee at Wolfram Alpha, and I'm surprised that we don't have the parametric results for $x$ and $y$ as a function of time for the Kepler problem available on our website. But I will work to rectify that in the future and update my answer here if/when I do! Suffice it to say, however, that if you don't care about inter-planetary interactions, then simulating Newton's laws is definitely not the right way to approach this problem, and things like adaptive time-steps, while they would work, would be massive overkill.