Numerical simulation of the solar system

[phate]

Hi everybody,

I would like to ask your opinion about some results from my solar sistem simulator. First, let me explain the method I used. I started with the classical N-Body problem and wrote the equations of motions in the usual form:

 

[math]\frac{d^{2}\textbf{r}_{i}}{dt^{2}}=G\sum\frac{m_{j}}{r^{3}_{ij}}(\textbf{r}_{j}-\textbf{r}_{i})[/math]

 

where, of course, the sum is valid for j = 1,..,N and j not equal to i. To semplify the implementation I integrate the equations of motion in the inertial reference frame centered on the Sun, being the initial conditions for velocity and position given by the JPL ephemerides software HORIZONS (even if my final target is to perform the integration in the Earth - Moon barycenter..).

To date, my software only takes into account the presence of the Sun, the Earth, The Moon and Jupiter and in order to validate my results I focused my attention on the motion of the Moon around the Earth center. After a short time integration (60 days) I noticed a drift of the Moon orbit of about 5000 km approximately in the direction of the Sun. Now, I am quite surprised about that behavior and I would like to ask you if the gravitational effect of the Sun is strong enough to cause such a perturbation wrt the keplerian unperturbed orbit.

I performed the integration in a SIMULINK environment, using a variable step RK4(5) method with a maximum step size of 10 seconds and a relative tolerance of 1e-13, so my hope is that if there is any problem it will not be related to the solver (which means that some how I can lukily manage it ).

What do you think about it?

Looking forward to hear from you, I thank you for your attention.

Regards

[Sisyphus]

This is just a guess, but are you sure you're thinking properly about eccentricity? Remember, the moon is 40000km farther away from Earth at apogee than perigee, and that cycle corresponds with the anomalistic month, which is different from either the sidereal or the synodic months. In other words, maybe nothing is wrong?

[D H]

Phate, do you know anything about Langrange's planetary equations?

[phate]

This is just a guess, but are you sure you're thinking properly about eccentricity? Remember, the moon is 40000km farther away from Earth at apogee than perigee, and that cycle corresponds with the anomalistic month, which is different from either the sidereal or the synodic months. In other words, maybe nothing is wrong?

 

Sisyphus, I think that using the N-body equations of motion with proper initial conditions should take into consideration the real eccentricities of the bodies. I agree, something is surely wrong...

 

Phate, do you know anything about Langrange's planetary equations?

 

DH,

do you suggest me to use the Lagrange's Planetary equations? As far as I know, the should give the rates for the orbital elements variations, but to be honest I am not so confident with them. I will have a look to the theory behind, in the meantime have you got any kind of suggestions?

 

Thank you everybody for your help.

Regards

[D H]

DH,

do you suggest me to use the Lagrange's Planetary equations? As far as I know, the should give the rates for the orbital elements variations, but to be honest I am not so confident with them. I will have a look to the theory behind, in the meantime have you got any kind of suggestions?

I am not suggesting you use Lagrange's Planetary equations for propagation. I am merely suggesting you use them to understand why the Moon's position is advanced in the n-body model compared to the simple 2-body solution.

 

Qualitatively, the Sun acts as a perturbing force on the 2-body problem. If you just use a first order approximation, this perturbation force more-or-less averages out to zero over one Moon orbit. However, the gravity gradient is a bit higher sunward than anti-sunward. Averaging a second-order expansion over a period of one orbit yields a small sunward acceleration.

[phate]

Qualitatively, the Sun acts as a perturbing force on the 2-body problem. If you just use a first order approximation, this perturbation force more-or-less averages out to zero over one Moon orbit. However, the gravity gradient is a bit higher sunward than anti-sunward. Averaging a second-order expansion over a period of one orbit yields a small sunward acceleration.

 

DH,

so you are somehow saying that the drift I observed in my simulator is correct? I agree with the fact that the Sun acts as a perturbing force on the 2-body problem: my doubt is about the module of this drift, because my feeling is that approx 5000 km (compared to the keplerian orbit) in only 60 days is quite a big number. :confused:

Do you have any number to compare? I thought to take the ephemerides of the Moon for one month and compare the propagated solution with the ephemerides, do you think it is a good check?

 

Regards

[D H]

DH,

so you are somehow saying that the drift I observed in my simulator is correct? I agree with the fact that the Sun acts as a perturbing force on the 2-body problem: my doubt is about the module of this drift, because my feeling is that approx 5000 km (compared to the keplerian orbit) in only 60 days is quite a big number. :confused:

I did a back-of-the-envelope caclulation using Lagrange's planetary equations (I ignored the eccentricity and inclination of the Moon's orbit) and got an advance of 30,000 km via [math]2\pi \,\frac {{\omega_E}^2}{\omega_M}\, t\, r_M[/math]. I think that factor of 2*pi might not be right, in which case that 30,000 km becomes 5,000 km. Either way, 5000 km is not out of line because the simplifications.

 

I thought to take the ephemerides of the Moon for one month and compare the propagated solution with the ephemerides, do you think it is a good check?

Comparing against reality is always a good idea.

[Sisyphus]

Sisyphus, I think that using the N-body equations of motion with proper initial conditions should take into consideration the real eccentricities of the bodies. I agree, something is surely wrong...

 

Yes, the equations would. I was asking if you were. As in, are you sure the moon isn't where it's supposed to be?

[D H]

I did a back-of-the-envelope caclulation using Lagrange's planetary equations (I ignored the eccentricity and inclination of the Moon's orbit) and got an advance of 30,000 km via [math]2\pi \,\frac {{\omega_E}^2}{\omega_M}\, t\, r_M[/math]. I think that factor of 2*pi might not be right, in which case that 30,000 km becomes 5,000 km.

That factor of 2*pi is wrong. Ignoring the eccentricity and inclination of the Moon's orbit, Lagrange's planetary equations dictate that

 

[math]\frac {d\mathcal X}{dt} = -\,\frac 2{\omega_Mr_M}\,\frac{\partial \mathcal R}{d r_M}[/math]

 

where [math]\mathcal X[/math] is the difference between the argument of latitude and that expected from a pure 2-body perspective and [math]\mathcal R[/math] is the perturbative potential energy. The Sun is the driving factor here, so

 

[math]\mathcal R = -GM_S\left(\frac 1{||\mathbf r_E + \mathbf r_M||} - \frac 1{r_E}\right)[/math]

 

Expanding this out and averaging over one lunar orbit yields

 

[math]\frac {d\mathcal X}{dt} \approx \frac {{\omega_E}^2}{\omega_M}[/math]

 

and, yep, that factor of 2*pi doesn't belong. Your 5,000 km is right on the mark.

 

Finally, I reiterate that comparing results against reality is always a good idea.

[phate]

Hi everybody,

in the last days I tried to manage with the problem. This is the situation to date.

I used the equation for N-Body problem ([math]\frac{d^{2}\textbf{r}_{i}}{dt^{2}}=G\sum\frac{m_{j}}{r^{3}_{ij}}(\textbf{r}_{j}-\textbf{r}_{i})

[/math]) to integrate the state of the 9 planets of the solar system with initial condition from the JPL ephemerides tool HORIZONS. Then I compared the numerical data against the ephemerides data (assuming that JPL ephemerides are correct... :D). After a great effort I still have big differences between the predicted data and the real ones even after short integration times: I made some tests, and propagated the Earth "switching off" all the perturbing bodies, that is reducing the model to a simple two body problem, and the results are consistent with this simple model. Now my question is: which kind of perturbing actions JPL put into the ephemerides? I wonder if I am looking for a precision that is out of my possibility...

 

The work is going ahead, I will keep you informed. In the meantime if you all have any kind of suggestion I will greatly appreciate.

 

Regards

[D H]

I have two suspects: Your initial conditions and your state propagation.

 

Initial conditions.

In order to do what you did you need initial position and velocity vectors for the planets. This means your initial positions and velocities must both be very accurate if you wish to have any chance of replicating the "truth". The HORIZONS system you are using in turn uses the JPL DE4xx (probably DE405) ephemeris model. If you delve inside this system (and you can do this; the software is available to the public), you will see that the DE series represents planet position by means of Chebyshev polynomials in time. Planet velocities are not represented directly. Instead, velocities are obtained by computing the time derivative of those Chebyshev polynomials. While JPL claims accuracy of one meter for the position vectors, it makes no claims about the accuracy of the initial velocity vectors.

 

State propagation.

From other posts, I suspect you are using the standard Matlab/Simulink RK45 integrator to propagate state. I have three words: Yech, and double yech. I do use this integrator when I need to do a quick-and-dirty kind of analysis. Even then I will do a sanity check because it the Matlab RK45 integrator oftentimes does a rather lousy job. One solution is to write your own integrator. You can take advantage of the conservative nature of the problem whereas RK45 cannot because it is a general-purpose integrator.

 

An alternative.

Another solution is to abandon the concept of integrating the planetary states. You can instead obtain them from the ephemeris model at each and every each time step. That means replacing your use of HORIZONS with the Chebyshev polynomial coefficients and software to use those coefficients. See http://ssd.jpl.nasa.gov/?planet_eph_export for (some) details. There is a book on the system; I'll see if I can dig up a reference for you later today.

[phate]

I have two suspects: Your initial conditions and your state propagation.

 

Initial conditions.

In order to do what you did you need initial position and velocity vectors for the planets. This means your initial positions and velocities must both be very accurate if you wish to have any chance of replicating the "truth". The HORIZONS system you are using in turn uses the JPL DE4xx (probably DE405) ephemeris model. If you delve inside this system (and you can do this; the software is available to the public), you will see that the DE series represents planet position by means of Chebyshev polynomials in time. Planet velocities are not represented directly. Instead, velocities are obtained by computing the time derivative of those Chebyshev polynomials. While JPL claims accuracy of one meter for the position vectors, it makes no claims about the accuracy of the initial velocity vectors.

 

This is a great surprise to me. I thought JPL was able to know both position and velocity with high accuracy. Even because with such a problem about the initial velocities vector to propagate the state vector of a spacecraft AND the solar system bodies with high precision is not possible at all.

I mada a rude check on the velocities and of course they match with those from the keplerian motion, but I can't check small errors...

 

State propagation.

From other posts, I suspect you are using the standard Matlab/Simulink RK45 integrator to propagate state. I have three words: Yech, and double yech. I do use this integrator when I need to do a quick-and-dirty kind of analysis. Even then I will do a sanity check because it the Matlab RK45 integrator oftentimes does a rather lousy job. One solution is to write your own integrator. You can take advantage of the conservative nature of the problem whereas RK45 cannot because it is a general-purpose integrator.

 

You are right. I use the standard MATLAB/SIMULINK rk4(5) integrator. I studied several other methods, but I have to think about how it could be possibile to "integrate" them in the SIMULINK environment. I don't know much about C++ or JAVA, do you think I should move to another programming language or try to implement another integrator in SIMULINK?

 

An alternative.

Another solution is to abandon the concept of integrating the planetary states. You can instead obtain them from the ephemeris model at each and every each time step. That means replacing your use of HORIZONS with the Chebyshev polynomial coefficients and software to use those coefficients. See http://ssd.jpl.nasa.gov/?planet_eph_export for (some) details. There is a book on the system; I'll see if I can dig up a reference for you later today.

 

I had the same idea before I started (even if I don't now anything about the Chebyshev polynomial coefficients, I thought to get the position of each planet from the ephemerides file in cartesian coordinates) but I moved to a more "dynamical" approach. Following the way you suggest, I should also take into account the possibility to insert a spacecraft in my scenario. In this case, I should get the position of the planets at each instant of time from the ephemerides, then calculate the perturbing force acting on the spacecraft to be added to the keplerian motion of the spacecraft around the primary of motion. Would it be correct? Do you suggest me to change the integrator in this case as well?

 

Again, thank you very much for your patience and your willingness.

[D H]

This is a great surprise to me. I thought JPL was able to know both position and velocity with high accuracy.

What JPL knows is observed right ascension and declination of the planets at various times. These measurements are only 2D and are not particular accurate, at least not in the sense of the meter-level accuracy JPL claims for the planetary positions. They also have a very detailed simulation that models the behavior the planets and minor bodies. I don't know how they are modeling state in this simulation (it might not be Cartesian coordinates), but I do know they go so far as to model relativistic effects. They most likely have the state of the system at epoch time use a batch update to tweak that epoch state to minimize the errors between the propagated state and observations. With this calibrated simulation they can now compute the state at any point in time.

 

There is no dynamics at all in the DExxx ephemeris model. It instead simply determines which set of Chebyshev polynomial coefficients are to be used and computes position and velocity using these coefficients. The dynamics comes into play in generating coefficients. The simulator yields position vectors at various time points picked according to Chebyshev polynomial formalism. The coefficients are rather easy to compute; a dozen or so lines of code does it. I suspect they also compute position vectors at other points. These other points will not be used to compute the coefficients; they are explicitly reserved for testing the validity of the fit.

 

One question: what degree polynomial should be used? Too small and the fit will be lousy. Too big and the fit will be great but the polynomial will be absolutely useless. Think about it this way: You have 99 data points comprising (t,x) values at hand that, when plotted, yield a curve that looks quadratic. You could fit the data points to a 100^th degree polynomial and get an exact fit. However, all this exact-fit polynomial is good for one thing and one thing only: Replicating the data you already have at hand. The polynomial will oscillate wildly between your sampled data points and will quickly diverge outside of the data points. It is much better to fit the data to a quadratic and accept the loss in fit accuracy. The people who generate the coefficients have to make a similar trade. Their goal is to reproduce the simulation results to within one meter. If this requires a high-degree polynomial, so be it. An overfit will be rejected because the evaluation data are distinct from the generation data. However, the evaluation data are still position vectors. This is what JPL guarantees the accuracy of. The velocity vectors are a freebie, and like everything else in life, you get what you pay for.

 

You are right. I use the standard MATLAB/SIMULINK rk4(5) integrator. I studied several other methods, but I have to think about how it could be possibile to "integrate" them in the SIMULINK environment. I don't know much about C++ or JAVA, do you think I should move to another programming language or try to implement another integrator in SIMULINK?

If you want something easy, write it in Matlab. (rk4(5) for example is written in Matlab.) Move to another language if performance is an issue. One thing you will find if you do write your own integrator is that your simulation will speed up immensely. Matlab's rk4(5) makes an inordinate number of calls to your derivative function at the start just to determine what it thinks is the correct time step to use. Another problem with the standard rk4(5) algorithm: It isn't very good with stiff systems. The Moon orbits once per month and Jupiter orbits once per 12 years; the outer planets are even slower. These widely-varying frequencies are the hallmark of a stiff system.

 

I had the same idea before I started (even if I don't now anything about the Chebyshev polynomial coefficients, I thought to get the position of each planet from the ephemerides file in cartesian coordinates) but I moved to a more "dynamical" approach.

Do you know of the work of Cesar Ocampo, University of Texas? The Copernicus system his group at UT is developing in conjunction with NASA has the ability to do both. I suggest you do a literature search. You now have a name.

 

Do you suggest me to change the integrator in this case as well?

I use Matlab's rk4(5) integrator when I just want a quick-and-dirty answer and don't care much for accuracy. Otherwise I use something different.

[h4tt3n]

When it comes to choosing an integrator for planetary simulation, Runge-Kutta schemes aren't preferrable, simply because they're not symplectic, ie. they don't conserve energy. RK4 actually looses energy very slowly.

 

I'd reccomend a symplectic 4th, 6th or even 8th order integrator like these (small C++ code sample):

 

http://www.aoc.nrao.edu/~dwhysong/prog/symp.cpp

 

They handle mass-spring-damper systems very poorly, but they are awe-inspiringly accurate and stable when it comes to gravitational simulations. The theory behind them is described in this paper:

 

http://128.125.4.56/education/phys516/Yoshida-symplectic-PLA00.pdf

 

Cheers,

Michael