Numerical model of the Solar System: Difference between revisions

Revision as of 08:00, 9 April 2009

Latest revision as of 07:06, 14 August 2025

(45 intermediate revisions by 33 users not shown)

Line 1:

{{Short description|Equations to predict the position of planets}}

{{No footnotes|article|date=April 2009}}A '''numerical model of the [[solarSolar system]]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.

==Older efforts==

The simulations can be done in either [[Cartesian coordinate system|Cartesian]] or in [[Spherical coordinate system|spherical]] coordinates. The former are easier, but extremely calculation intensive, and only doablepractical on an electronic computer. As such only the latter was used in former times. Strictly speaking, the latter was not much less calculation intensive, but it was possible to start with some simple approximations and then to add [[Perturbation (astronomy)|perturbations]], as much as needed to reach the wanted accuracy.

In essence this mathematical simulation of the solar[[Solar systemSystem]] is a form of the ''[[N-body problem]]''. The symbol '''''N''''' represents the number of bodies, which can grow quite large if one includes 1the sunSun, 8 planets, dozens of moons, and countless planetoids, comets and so forth. However the influence of the sunSun 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. (SmallThese are small compared to a full planetary orbit,. someSome 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 (''v_x'', ''v_y'', ''v_z'') for each of the bodies involved. When also the mass of each body is known, the acceleration (''a_x'', ''a_y'', ''a_z'') can be calculated from [[Newton's Lawlaw of Gravitationgravitation]]. Each body attracts each other body, the total acceleration being the sum of all these attractions. Next one chooses a small timesteptime-step Δ''t'' and applies [[Newton's Secondsecond Lawlaw of Motionmotion]]. 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 timesteptime-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 timestepstime-steps are to be calculated too.

==Integration==

Line 22 ⟶ 23:

<math>\vec{a}_j = \sum_{i \neq j}^n G \frac{M_i}{|\vec{r}_i - \vec{r}_j|^3} (\vec{r}_i - \vec{r}_j)</math>

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:

<math>(a_j)_x = \sum_{i \neq j}^n G \frac{M_i}{( (x_i - x_j)^2 + (y_i - y_j)^2 + (z_i - z_j)^2 )^{3/2}} (x_i - x_j)</math>

Line 32 ⟶ 33:

likewise for Y and Z.

The former equation (gravitation) may look forbiddingforeboding, but its calculation is no problem. The latter equations (motion laws) seemsseem simpler, but yet itthey 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:

<math>\Delta v_x = a_{x} \Delta t </math>, and: <math>\Delta x = v_{x} \Delta t </math>

Line 49 ⟶ 50:

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. The most common method used is the [[leapfrog method]] due to its good long term energy conservation.

A completecompletely different method is the use of [[Taylor series]]. In that case we write: <math>r = r_0 + r'_0 t + r''_0 \frac{t^2}{2!} + (etc)... </math>

but rather than developing up to some higher derivative in r only, one can develop in r and v (that is r') by writing <math>r = f r_0 + g r'_0</math>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==

Line 62 ⟶ 61:

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 solarSolar systemSystem the [[astronomical unit]] is most straightforward. If this is not done one is almost certain to see a simulation abortedabandoned in the middle of a calculation on a [[floating point]] [[arithmetic overflow|overflow]] or [[arithmetic underflow|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 solarSolar systemSystem 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

<math>a = \frac{G M}{r^2 + e}</math>

==Complications==

InIf casethe highest possible accuracy is needed, the thingscalculations getbecome much more complex. In the case of comets, non gravitationalnongravitational forces, such as (radiation pressure and gas drags)drag, are tomust be taken into account. In the case of Mercury, and other planets for long term calculations, relativityrelativistic effects cannot be ignored. Then also the total energy is no longer a constant (because the four vector energy with linear momentum is). The finalfinite speed of light also makes it important to allow for light-time effects, both classical and relativityrelativistic. Planets can no longer be considered as particles, but their shape and density must also be considered too. For example, the flattening of the Earth causes precession, which causes the axial tilt to change, which has its effects onaffects the longtermlong-term movements of all planets.

Long term models, going beyond a few tens of millions of years, are not possible due to the lack of [[stability of the Solar System]].

==See also==

*[[Ephemeris]]

*[[Astronomical algorithm]]

*[[VSOP (planets)]]

*[[Variations séculaires des orbites planétaires]] (VSOP)

==References==

*{{citeCite book |first=Dan L. |last=Boulet |title=Methods of orbit determination for the microcomputer |publisher=Willmann-Bell, Inc |___location=[[Richmond, Virginia]] |year=1991 |pages= |isbn=978-0-943396-34-47 |oclc=23287041}}{{pnPage needed|date=September 2010}}

{{Solar System models}}

{{DEFAULTSORT:Numerical Model Of Solar System}}

[[Category:Numerical analysis]]

[[Category:Computational physics]]

[[Category:Dynamical systems]]

[[Category:Dynamics of the Solar System]]

[[Category:Solar System models]]