Numerical Model of The Solar System - Integration

Integration

In the previous section it was assumed that acceleration remains constant over a small timestep Δt so that the calculation reduces to simply the addition of V × Δt to R and so forth. In reality this is not the case, except when one takes Δt so small that the number of steps to be taken would be prohibitively high. Because while at any time the position is changed by the acceleration, the value of the acceleration is determined by the instantaneous position. Evidently a full integration is needed.

Several methods are available. First notice the needed equations:

This equation describes the acceleration all bodies i running from 1 to N exercise on a particular body j. It is a vector equation, so it is to be split in 3 equations for each of the X, Y, Z components, yielding:

with the additional relationships

,

likewise for Y and Z.

The former equation (gravitation) may look foreboding, but its calculation is no problem. The latter equations (motion laws) seems simpler, but yet it cannot be calculated. Computers cannot integrate, they cannot work with infinitesimal values, so instead of dt we use Δt and bringing the resulting variable to the left:

, and:

Remember that a is still a function of time. The simplest way to solve these is just the Euler algorithm, which in essence is the linear addition described above. Limiting ourselves to 1 dimension only in some general computer language:

a.old = gravitationfunction(x.old) x.new = x.old + v.old * dt v.new = v.old + a.old * dt

As in essence the acceleration used for the whole duration of the timestep, is the one as it was in the beginning of the timestep, this simple method has no high accuracy. Much better results are achieved by taking a mean acceleration, the average between the beginning value and the expected (unperturbed) end value:

a.old = gravitationfunction(x.old) x.expect = x.old + v.old * dt a.expect = gravitationfunction(x.expect) v.new = v.old + (a.old + a.expect) * 0.5 * dt x.new = x.old + (v.new + v.old) * 0.5 * dt

Of course still better results can be expected by taking intermediate values. This is what happens when using the Runge-Kutta method, especially the one of grade 4 or 5 are most useful.

A completely different method is the use of Taylor series. In that case we write:

but rather than developing up to some higher derivative in r only, one can develop in r and v (that is r') by writing and then write out the factors f and g in a series.

All these more advanced methods easily allow for Solar System calculations with a stepsize Δt of 10 days and yet yield satisfactory results.