Numerical model of solar system
Encyclopedia
A numerical model of the Solar System
is a set of mathematical equations, which, when solved, give the approximate positions of the planets as a function of time. Attempts to create such a model established the more general field of celestial mechanics
. The results of this simulation can be compared with past measurements to check for accuracy and then be used to predict future positions. Its main use therefore is in preparation of almanacs.
or in spherical
coordinates. The former are easier, but extremely calculation intensive, and only practical on an electronic computer. As such only the latter was used in former times. Strictly speaking not much less calculation intensive, but it was possible to start with some simple approximations and then to add perturbations
, as much as needed to reach the wanted accuracy.
In essence this mathematical simulation of the Solar System is a form of the N-body problem
. The symbol N represents the number of bodies, which can grow quite large if one includes 1 sun, 8 planets, dozens of moons and countless planetoids, comets and so forth. However the influence of the sun on any other body is so large, and the influence of all the other bodies on each other so small that the problem can be reduced to the analytically solvable 2-body problem. The result for each planet is an orbit, a simple description of its position as function of time. Once this is solved the influences moons and planets have on each other are added as small corrections. (Small compared to a full planetary orbit, some corrections might be still several degrees large, while measurements can be made to an accuracy of better than 1″).
Although this method is no longer used for simulations, it is still useful to find an approximate ephemeris as one can take the relatively simple main solution, perhaps add a few of the largest perturbations, and arrive without too much effort at the wanted planetary position. The disadvantage is that perturbation theory is very advanced mathematics.
The result is a new value for position and velocity for all bodies. Then, using these new values one starts over the whole calculation for the next time-step Δt. Repeating this procedure often enough, and one ends up with a description of the positions of all bodies over time.
The advantage of this method is that for a computer it is a very easy job to do, and it yields highly accurate results for all bodies at the same time, doing away with the complex and difficult procedures for determining perturbations. The disadvantage is that one must start with highly accurate figures in the first place, or the results will drift away from the reality in time; that one gets x, y, z positions which are often first to be transformed into more practical ecliptical or equatorial coordinates before they can be used; and that it is an all or nothing approach. If one wants to know the position of one planet on one particular time, then all other planets and all intermediate time-steps are to be calculated too.
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.
By far the most important trick is the use of a proper integration method, as already outlined above.
The choice of units is important. Rather than to work in SI units, which would make some values extremely small and some extremely large, all units are to be scaled such that they are in the neighbourhood of 1. For example for distances in the Solar System the astronomical unit
is most straightforward. If this is not done one is almost certain to see a simulation aborted in the middle of a calculation on a floating point
overflow
or underflow
, and if not that bad, still accuracy is likely to get lost due to truncation
errors.
If N is large (not so much in Solar System simulations, but more in galaxy simulations) it is customary to create dynamic groups of bodies. All bodies in a particular direction and on large distance from the reference body, which is being calculated at that moment, are taken together and their gravitational attraction is averaged over the whole group.
The total amount of energy
and angular momentum
of a closed system are conserved quantities. By calculating these amounts after every time step the simulation can be programmed to increase the stepsize Δt if they do not change significantly, and to reduce it if they start to do so. Combining the bodies in groups as in the previous and apply larger and thus less timesteps on the faraway bodies than on the closer ones, is also possible.
To allow for an excessively rapid change of the acceleration when a particular body is close to the reference body, it is customary to introduce a small softness parameter e so that
Solar System
The Solar System consists of the Sun and the astronomical objects gravitationally bound in orbit around it, all of which formed from the collapse of a giant molecular cloud approximately 4.6 billion years ago. The vast majority of the system's mass is in the Sun...
is a set of mathematical equations, which, when solved, give the approximate positions of the planets as a function of time. Attempts to create such a model established the more general field of celestial mechanics
Celestial mechanics
Celestial mechanics is the branch of astronomy that deals with the motions of celestial objects. The field applies principles of physics, historically classical mechanics, to astronomical objects such as stars and planets to produce ephemeris data. Orbital mechanics is a subfield which focuses on...
. The results of this simulation can be compared with past measurements to check for accuracy and then be used to predict future positions. Its main use therefore is in preparation of almanacs.
Older efforts
The simulations can be done in either CartesianCartesian coordinate system
A Cartesian coordinate system specifies each point uniquely in a plane by a pair of numerical coordinates, which are the signed distances from the point to two fixed perpendicular directed lines, measured in the same unit of length...
or in spherical
Spherical coordinate system
In mathematics, a spherical coordinate system is a coordinate system for three-dimensional space where the position of a point is specified by three numbers: the radial distance of that point from a fixed origin, its inclination angle measured from a fixed zenith direction, and the azimuth angle of...
coordinates. The former are easier, but extremely calculation intensive, and only practical on an electronic computer. As such only the latter was used in former times. Strictly speaking not much less calculation intensive, but it was possible to start with some simple approximations and then to add perturbations
Perturbation (astronomy)
Perturbation is a term used in astronomy in connection with descriptions of the complex motion of a massive body which is subject to appreciable gravitational effects from more than one other massive body....
, as much as needed to reach the wanted accuracy.
In essence this mathematical simulation of the Solar System is a form of the N-body problem
N-body problem
The n-body problem is the problem of predicting the motion of a group of celestial objects that interact with each other gravitationally. Solving this problem has been motivated by the need to understand the motion of the Sun, planets and the visible stars...
. The symbol N represents the number of bodies, which can grow quite large if one includes 1 sun, 8 planets, dozens of moons and countless planetoids, comets and so forth. However the influence of the sun on any other body is so large, and the influence of all the other bodies on each other so small that the problem can be reduced to the analytically solvable 2-body problem. The result for each planet is an orbit, a simple description of its position as function of time. Once this is solved the influences moons and planets have on each other are added as small corrections. (Small compared to a full planetary orbit, some corrections might be still several degrees large, while measurements can be made to an accuracy of better than 1″).
Although this method is no longer used for simulations, it is still useful to find an approximate ephemeris as one can take the relatively simple main solution, perhaps add a few of the largest perturbations, and arrive without too much effort at the wanted planetary position. The disadvantage is that perturbation theory is very advanced mathematics.
Modern method
The modern method consists of numerical integration in 3-dimensional space. One starts with a high accuracy value for the position (x, y, z) and the velocity (vx, vy, vz) for each of the bodies involved. When also the mass of each body is known, the acceleration (ax, ay, az) can be calculated from Newton's Law of Gravitation. Each body attracts each other body, the total acceleration being the sum of all these attractions. Next one chooses a small time-step Δt and applies Newton's Second Law of Motion. The acceleration multiplied with Δt gives a correction to the velocity. The velocity multiplied with Δt gives a correction to the position. This procedure is repeated for all other bodies.The result is a new value for position and velocity for all bodies. Then, using these new values one starts over the whole calculation for the next time-step Δt. Repeating this procedure often enough, and one ends up with a description of the positions of all bodies over time.
The advantage of this method is that for a computer it is a very easy job to do, and it yields highly accurate results for all bodies at the same time, doing away with the complex and difficult procedures for determining perturbations. The disadvantage is that one must start with highly accurate figures in the first place, or the results will drift away from the reality in time; that one gets x, y, z positions which are often first to be transformed into more practical ecliptical or equatorial coordinates before they can be used; and that it is an all or nothing approach. If one wants to know the position of one planet on one particular time, then all other planets and all intermediate time-steps are to be calculated too.
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
Taylor series
In mathematics, a Taylor series is a representation of a function as an infinite sum of terms that are calculated from the values of the function's derivatives at a single point....
. 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.
Approximations
To calculate the accelerations the gravitational attraction of each body on each other body is to be taken into account. As a consequence the amount of calculation in the simulation goes up with the square of the number of bodies: Doubling the number of bodies increases the work with a factor four. To increase the accuracy of the simulation not only more decimals are to be taken but also smaller timesteps, again quickly increasing the amount of work. Evidently tricks are to be applied to reduce the amount of work. Some of these tricks are given here.By far the most important trick is the use of a proper integration method, as already outlined above.
The choice of units is important. Rather than to work in SI units, which would make some values extremely small and some extremely large, all units are to be scaled such that they are in the neighbourhood of 1. For example for distances in the Solar System the astronomical unit
Astronomical unit
An astronomical unit is a unit of length equal to about or approximately the mean Earth–Sun distance....
is most straightforward. If this is not done one is almost certain to see a simulation aborted in the middle of a calculation on a floating point
Floating point
In computing, floating point describes a method of representing real numbers in a way that can support a wide range of values. Numbers are, in general, represented approximately to a fixed number of significant digits and scaled using an exponent. The base for the scaling is normally 2, 10 or 16...
overflow
Arithmetic overflow
The term arithmetic overflow or simply overflow has the following meanings.# In a computer, the condition that occurs when a calculation produces a result that is greater in magnitude than that which a given register or storage location can store or represent.# In a computer, the amount by which a...
or underflow
Arithmetic underflow
The term arithmetic underflow is a condition in a computer program that can occur when the true result of afloating point operation is smaller in magnitude...
, and if not that bad, still accuracy is likely to get lost due to truncation
Truncation
In mathematics and computer science, truncation is the term for limiting the number of digits right of the decimal point, by discarding the least significant ones.For example, consider the real numbersThe result would be:- Truncation and floor function :...
errors.
If N is large (not so much in Solar System simulations, but more in galaxy simulations) it is customary to create dynamic groups of bodies. All bodies in a particular direction and on large distance from the reference body, which is being calculated at that moment, are taken together and their gravitational attraction is averaged over the whole group.
The total amount of energy
Energy
In physics, energy is an indirectly observed quantity. It is often understood as the ability a physical system has to do work on other physical systems...
and angular momentum
Angular momentum
In physics, angular momentum, moment of momentum, or rotational momentum is a conserved vector quantity that can be used to describe the overall state of a physical system...
of a closed system are conserved quantities. By calculating these amounts after every time step the simulation can be programmed to increase the stepsize Δt if they do not change significantly, and to reduce it if they start to do so. Combining the bodies in groups as in the previous and apply larger and thus less timesteps on the faraway bodies than on the closer ones, is also possible.
To allow for an excessively rapid change of the acceleration when a particular body is close to the reference body, it is customary to introduce a small softness parameter e so that
Complications
If the highest possible accuracy is needed, things become much more complex. In the case of comets, nongravitational forces (radiation pressure and gas drag) must be taken into account. In the case of Mercury, relativistic effects cannot be ignored. Then also the total energy is no longer a constant (because the four vector energy with linear momentum is). The finite speed of light also makes it important to allow for light-time effects, both classical and relativistic. Planets can no longer be considered as particles, but their shape and density must also be considered. For example, the flattening of the Earth causes precession, which causes the axial tilt to change, which affects the long-term movements of all planets.See also
- EphemerisEphemerisAn ephemeris is a table of values that gives the positions of astronomical objects in the sky at a given time or times. Different kinds of ephemerides are used for astronomy and astrology...
- Astronomical algorithmAstronomical algorithmAstronomical algorithms are the algorithms used to calculate ephemerides, calendars, and positions . Examples of large and complex astronomical algorithms are those used to calculate the position of the Moon...
- Secular variations of the planetary orbitsSecular variations of the planetary orbitsThe secular variations of the planetary orbits is a concept describing long-term changes in the orbits of the planets Mercury to Neptune. If one ignores the gravitational attraction between the planets and only models the attraction between the Sun and the planets, then with some further...