What is Orbiter Space Flight Simulator?

[Ulysees]

May I ask, for a given accuracy of 10^-5 say, and for a given orbit lasting 1 year, which integration method is the fastest?

 

[csprof2000]

I would suspect that you can't do much better than the fourth order adaptive stepsize Runge Kutta algorithm.

 

[Ulysees]

Does that mean repeating with a smaller time interval until all 3 steps of 4th order Runge-Kutta agree within an error margin?

 

[csprof2000]

I'm not sure precisely what happens, but I'm pretty sure that a timestep is taken with one RK method and then another is taken with one which is of one higher order, and if the difference is within the error bars, it goes with it. Otherwise, it tells you what the timestep should be to get your accuracy. It has the benefit of remaining as accurate as you tell it (within reason) and not being *too* when the forces are negligibly small (for instance, why use a timestep of h = 1.0e-10 both when the force is large and changing quickly AND when the force is next to nothing and staying put?)

 

[tony873004]

You should start with Euler's method because it is simple and you understand it. Then if you like, you can upgrade your method later.

I wrote Gravity Simulator: www.gravitysimulator.com , which sounds like what you want to do. It uses Euler's method. Using Euler's method I can:

*Predict moon phases at least 50 years into the future
*Propogate Apophis 22 years into the future, and miss by only 200 km its close approach distance to Earth in 2029 (which btw is also the amount the higher-order Picard's method misses by, so the inaccuracy seems to have more to do with forces not modeled, rather than with integrator inaccuracies).
*Correctly predict the Newtonian predicted value of the precession of Mercury's perihelion
*accurately recreate virtually any point-mass Newtonian simulation done by professional astronomers except for the ones that last billions of years.
*and the list goes on. My website has more examples

I'm not an expert on the higher order routines, but most of the bad things I've heard about Euler's method seem to be overrated and seem to disappear with a slow enough time step.

In the real universe, there are lots of other forces besides Newtonian gravity from point-mass sources: radiation pressure, relativity, non-spherical objects, etc. I'm guessing that errors caused by ignoring these overwhelm the errors caused by Euler's method. For example, the toughest thing to compute is the orbits of tight moons such as Phobos and Deimos. After 2 years of simulating, even at a slow time step, they're not where they are supposed to be. The size, shape and orientation of their orbits are still fine, but their position on that orbit is off, by lots for Phobos, and by a little for the more distant Deimos. But this is more likely due to treating Mars as spherical, than with numerical errors. So even a higher order routine will miss this one too.

I recently wrote a 2nd order routine. It was more accurate than Eulers at time steps less than 64 seconds, but less accurate than Euler at time steps higher than 64 seconds. In the 2nd order routine, the planets would spiral out. It was explained to me like this: In Euler (1st order) every time you cut your time step in half, you cut your error in half. In 2nd order routines, every time you cut your time step in half, you cut your error 4x. But the inverse of this is true too. Every time you double your time step, your error grows by 4x for the 2nd order routine and by only 2x for Euler. Since most simulations I like to perform require a time step much high than 64 seconds, I won't find much use for the 2nd order routine. Next I'd like to try RK4. It will be interesting to compare to the current Euler routine. Will I like it better, given that it seems like Euler is very underrated?

Don't be discouraged from using Euler if you don't understand the higher order methods. You'll have lots of fun with it. And don't be scared away from 3d. It's not that much different from 2d. Pothaorean theorem works in 3d too: d=sqrt(x^2+y^2+z^2), a=GM/d^2, r=r+a

 

[Redbelly98]

You should start with Euler's method because it is simple and you understand it. Then if you like, you can upgrade your method later.

Definitely good advice here. Get something working so you can start having fun with it. Then worry about what is the "best" method.

 

[tony873004]

...I've been using Cartesian coordinates for all of my work so far, and for the actual calculations i plan to continue - but i need to be able to use keplerian orbital elements instead of state vectors for initial conditions. I was wondering if anyone out there had experience or direction on how to convert between orbital elements and state vectors besides drowning myself in geometry. Thanks!

Why do you need to use Keplerian orbital elements as your initial conditions? From JPL Horizons, you can get state vectors.

I've got code for converting both directions from state vectors to orbital elements and back again. Send me a PM if you want a copy.

 

[tony873004]

One thing that needs to be asked is why does accuracy matter? Are you trying to predict solar eclipses far into the future? Are you trying to pilot a spacecraft to Mars where missing by just a few kilometers will doom your mission? Or are you simply trying to demonstrate a gravity assist, Kozai mechanism, precession of nodes, etc.? For most of these accuracy isn't really that important. I'll show a few examples of tests I've performed:

This is an animation of the Moon's orbit around the Earth. It was made with Gravity Simulator's 1st order Euler's method. The time step was 16000 seconds. At this time step, the simulation is very innacurate. This would be evident if you look at the Moon's position at a known future date. For example, on August 12, 2045 there will be a total eclipse of the Sun across America. So if I paused the simulation there, I should see the Moon directly between the Earth and Sun. But with the errors introduced by the numerical method, the Moon will likely be elsewhere in its orbit at this time. However, despite the error, many features of the Moon's orbit remain accurate: Its semi-major axis has held constant. Its eccentricity and inclination remained constant as well. And its nodes still precess at their 18 year cycle known as the saros cycle, as evident in this animation. In fact, the saros cycle remains consistent no matter how large of a time step I take as long as it is slow enough that the Moon remains in its orbit.

 

[D H]

In the real universe, there are lots of other forces besides Newtonian gravity from point-mass sources: radiation pressure, relativity, non-spherical objects, etc. I'm guessing that errors caused by ignoring these overwhelm the errors caused by Euler's method.

No, they don't. Moreover, if your orbit simulator maintains the Moon's semimajor axis, then you are not using the basic Euler method. You are most likely using the Euler-Cromer method instead. The two methods:

Euler method
\begin{aligned}<br /> \mathbf x(t+\Delta t) &amp;= \mathbf x(t) + \mathbf v(t)\Delta t \\<br /> \mathbf v(t+\Delta t) &amp;= \mathbf v(t) + \mathbf a(\mathbf x(t))\Delta t<br /> \end{aligned}

Euler-Cromer method
\begin{aligned}<br /> \mathbf v(t+\Delta t) &amp;= \mathbf v(t) + \mathbf a(\mathbf x(t))\Delta t \\<br /> \mathbf x(t+\Delta t) &amp;= \mathbf x(t) + \mathbf v(t+\Delta t)\Delta t<br /> \end{aligned}

The two methods differ only in the calculation of the updated position vector. In the basic Euler method, the pre-update velocity is used while in Euler-Cromer, the position is updated using the post-update velocity. This slight change makes a big difference. Euler-Cromer is a symplectic integrator, meaning that it more-or-less conserves energy in situations with a time-independent Hamiltonian (e.g., gravity). The degree to which it fails to conserve energy is a function of step size and numerical errors inherent in using finite-precision arithmetic. The basic Euler method is not symplectic; orbits spiral out.

 

[tony873004]

Thanks for that, DH. My method updates velocity and then position. I never knew that was called Euler-Cromer. I actually made it up myself. (can we call it Euler-Cromer-Tony? :) ).

I had always heard that Euler should cause the orbits to spiral out, so I wondered why mine never did.

 

[Ulysees]

I should add the more important problem with Euler variants is not accuracy. It is stability. For example a mass attached to a spring without any damping should oscillate forever (that's the analytical solution), but with Euler it oscillates with an increasing amplitude. No matter how small you make the integration step, the amplitude of oscillation is increasing.

4th order Runge Kutta is not so unstable. So its 3 steps are better than 3 normal steps of Euler.

 

[D H]

I should add the more important problem with Euler variants is not accuracy. It is stability. For example a mass attached to a spring without any damping should oscillate forever (that's the analytical solution), but with Euler it oscillates with an increasing amplitude. No matter how small you make the integration step, the amplitude of oscillation is increasing.

That is always a problem with the basic Euler method. It is not so much a problem with the Euler-Cromer method, which like RK4, exhibits conditionally stable (but the stability region is considerably smaller with Euler-Cromer than RK4). The simple act of updating the position with the already-updated velocity makes Euler-Cromer conditionally stable and symplectic. RK4 is considerably more stable and considerably more accurate than Euler-Cromer, but it is not symplectic.

 

[Ulysees]

And what's the benefit of symplectic-ness?

 

[D H]

And what's the benefit of symplectic-ness?

Ignoring errors resulting from the use of finite-precision arithmetic, a symplectic integrator comes very close to conserving energy. How close depends on the order of the integrator and the product of the integrator step size and the largest critical frequency of the system. When applied to a two body system, orbits will remain at more-or-less the size and shape with a symplectic integrator.

This is not the case for non-symplectic integrators such as basic Euler or the RK methods. With Euler the deviation is immediate. Orbits spiral out quickly. With RK4 the degradation is much, much slower, but still occurs. An RK4-based integrator will be closer to truth than the simplest symplectic integrator, Euler-Cromer, when used over just a few orbits. On the other hand, using an RK4 integrator to propagate over hundreds of orbits is a bad idea. At some point the energy will have changed enough to make all results garbage. Euler-Cromer suffers from reduced accuracy but benefits from increased global stability.

Just as RK4 is much better than RK1 (i.e., Euler's method), higher-order symplectic integrators are a lot better than Euler-Cromer. Higher-order methods include Verlet methods, Forest-Ruth, and Chin. Chin's Algorithm C is the new gold standard.

 

[tony873004]

Your replies to my posts and Ulysees posts are very facinating. Do you have any links describing Forest-Ruth, Verlet or Chin?

Something I wanted to work on this summer was a higher order integrator for Gravity Simulator. I thought I wanted to do RK4, but you may have discouraged me from that.

Why would anybody want to use Euler's method for anything, being that if you just solve for v before you solve for x, it becomes much nicer? Didn't Euler realize this?

 

[D H]

Your replies to my posts and Ulysees posts are very facinating. Do you have any links describing Forest-Ruth, Verlet or Chin?

Something I wanted to work on this summer was a higher order integrator for Gravity Simulator. I thought I wanted to do RK4, but you may have discouraged me from that.

Why would anybody want to use Euler's method for anything, being that if you just solve for v before you solve for x, it becomes much nicer? Didn't Euler realize this?

RK4 is a very, very good single-step integrator. (All of these techniques are single-step integrators which means that they don't keep a history. The three internal steps in RK4 are all part of the same step.) RK4 is typically more accurate than all the techniques mentioned but Chin's Algorithm C. That RK4 is not symplectic only becomes a significant factor when the integration is performed over a very long time span.

Ideally your simulator should be constructed so that the propagation technique is completely independent of the rest of the simulation. You should be able to switch propagators without changing any code. This is a completely achievable task and how easily it can be done is a good test of the quality of your code. (Chin's Algorithm C is an exception as it requires the gradient of the acceleration for one internal step.)

Euler died 225 years ago. Give him a break. He predates the concept of a Hamiltonian by 100 years and predates the concept of a vector by even more than that. The simple act of using the post-update velocity to update the position looks like a dumb mistake from a mathematical perspective -- unless one looks at things from the perspective of a Hamiltonian.

As far as references,

• Verlet. Use Google. The first page of hits are all true hits, including the two Wiki articles.
• RK4. The http://en.wikipedia.org/wiki/Runge-Kutta#The_classical_fourth-order_Runge.E2.80.93Kutta_method".
• Chin and Forest-Ruth. You can find Chin's papers at http://arxiv.org/find/grp_physics/1/au:+chin_siu/0/1/0/all/0/1". His papers also discuss FR since his method is an offshoot of that method.

 

[Ulysees]

Have you done any long-term simulations like the formation of a solar system? I was looking for a simulator suitable for testing Titius-Bode law.

Also, there are some who say that if a planet is hit by a large enough object or otherwise explodes, the entire solar system will collapse eventually. Any thoughts on this?

 

[tony873004]

Euler died 225 years ago. Give him a break.

Fair enough. I'll cut his some slack! Thanks for the links.

Have you done any long-term simulations like the formation of a solar system? I was looking for a simulator suitable for testing Titius-Bode law.

The formation of the solar system is too complicated to simulate with point-mass n-body code. When the primordial particles are small, gravity does not play a roll. Once they're large enough for gravity to play a role, certain collisions cause them to stick together, while others break them apart. And there's millions of them. My simulator does not do very well with even hundreds of particles.

However, if you're willing to make some assumptions, there's a few simulations you can run. Check this one out. It's called "moonbuilder": http://www.orbitsimulator.com/gravity/articles/moonbuilder.html . It begins with 100 particles each with 1/50 the mass of the real moon. All collisions cause them to merge together. It only takes a few weeks for the hundred to merge together into 2 moons, one that remains in orbit, and one that gets ejected. And different starting conditions will yield different results, so if you try it again and again, you may get systems with 2 or 3 moons. The more you try it, the more you should see some trends develop.

I've tried the same thing, except placing 100 particles at 1 AU from the Sun to form a planet. Sometimes I get trojans or horseshoe objects. Again, every run is different, but long-term trends can be spotted. So if you're clever, and willing to make some simplifications and assumptions, there are definitely some things you can try.

I don't think it would do to well trying to confirm or reject Titus-Bode. You'd be trying to build multiple planets at the same time, which would require lots of particles. You could play around with "Moonbuilder", perhaps spacing the initial particles further apart, and trying to find the ideal conditions that lead to 2 or 3 stable moons, and see if there's any kind of trend to the spacing. But I'd be careful about the conclusions you draw with such a simplified model.

Also, there are some who say that if a planet is hit by a large enough object or otherwise explodes, the entire solar system will collapse eventually. Any thoughts on this?

I've heard things like Earth's presence supresses a long-term resonance between Venus and Jupiter, and without Earth, Venus would not be stable. Maybe that's what you're referring to? I'm not sure how accurate that statement is. When I heard it, the first thing I tried to do was simulate the solar system without Earth. Venus stayed put, although I only simulated for a few million years. It's possible that it takes longer than this for Venus to get ejected, and its also possible that my time step was high enough that numerical errors dominated, and I won't get an ejection no matter how long I simulate. It takes a pretty big time step to get millions of years into the future.

 

[Ulysees]

Ιsn't it reasonable that Titius-Bode would be close to being satisfied with just 1000 particles? Why do you believe the statistics would only converge to an analytical pattern such as Titius Bode only with zillions of particles? We don't need to poll millions of voters in order to predict the outcome of an election, just a 1000 voters converge to a pattern too.

 

[tony873004]

That was the point I was trying to make: set up what you feel would be a representative system using fewer particles. But you need to understand how your model differs from reality. The newly-formed planets interact with the intermediate belts of trillions of lefover particles. This causes their orbits to circularize and to migrate to their final positions. That would have to be modeled too.

 

[Ulysees]

That would have to be modeled too.

And what's stopping you from imagining a handful of planets made up from 100 particles each, rotating in a cloud of 500 other particles, and this occurring automatically starting from an initial cloud. I only want to test Titius Bode as a mathematical property that should originate from the inverse square law in any scale, atomic or galactic. Not confine it to the actual masses, distances, etc of a solar system.

 

[tony873004]

And what's stopping you...

A lack of interest in Titus-Bode prevents me from trying to verify it through simulations. But if it interests you, that doesn't need to stop you. Anyone with a Windows computer can run my simulator.

...I only want to test Titius Bode as a mathematical property that should originate from the inverse square law in any scale, atomic or galactic...

There are 6 opportunities in our solar system for Titus-Bode to be revealed: the planetary system, and 5 moon systems (>2 moons: Jupiter, Saturn, Uranus, Neptune, Pluto), none of which demonstrate Titus-Bode. The planetary system fails because there's no major planet between Mars and Jupiter, and Neptune is out of place.

 

[D H]

Modeling the formation of planetary system is a vastly more complex task than a simple orbital simulator. Planetary growth is a complex and not yet fully understood process. I assume Tony's simulator treats the planets as point masses. This is a good assumption for a stabilized system with only a small number of orbiting objects. It is not a good assumption for modeling planet growth. Tony's simulator runs on a single off-the-shelf computer. Modeling protoplanetary systems requires a lot of coupled, high powered computers because one needs to model millions of tiny particles. These little particles collide (point masses don't collide), have angular momentum due to rotation (point masses don't have angular momentum due to rotation) and exchange angular momentum with the Sun and each other. Once again, point masses don't do this.

 

[tony873004]

...I assume Tony's simulator treats the planets as point masses...These little particles collide (point masses don't collide),..

They're treated as point masses for the sake of computing acceleration, but they do have size, so they can collide and stick. But every collision simply combines their mass and momentum into a new object if their distance is less than their combined radii. But the points you make are what I was getting at a few posts ago (#48). The real system has lots of stuff going on besides point mass interactions. There are more types of collisions than the simple one my code models. And the real system has gas that behaves differently than solid matter. So I can only agree with you that this type of a simulator is too simple for task. It doesn't mean than you can't make some assumptions and generate simplified models to explore various aspects. But to build a solar system from scratch, which you can draw confident conclusions about, is well beyond its capabilities.

 

[Ulysees]

I assume Tony's simulator treats the planets as point masses. This is a good assumption for a stabilized system with only a small number of orbiting objects. It is not a good assumption for modeling planet growth.

Well don't forget Newton's law of gravity only makes sense for point masses. Everything else is approximated by point masses, for the purposes of this law. Conglomerates of point masses can have angular momentum but there is no need to introduce angular momentum in the model. Just point masses connected with breakable and rejoinable spring-damper links, should be enough to model just about anything we might be interested in gravity-wise, and I think Tony is not far from this. The difficulty is with the estimation of suitable parameters for something that matches the real thing, but then we are not always interested in matching the real thing but just want to explore possibilities.

Regarding Titius-Bode that I am interested in testing, there is no reason to assume it applies exclusively to the scale of our solar system. I would expect it to also apply to much bigger and much smaller systems, because it seems to be a statistical result, like Gaussian distributions that appear in nature.

By the way, Titius-Bode -like equations predict a planet between Mars and Jupiter, and indeed there is the asteroid belt there, ie a planet that failed to form or a planet that formed and was later broken apart somehow, perhaps in a collision with a large object like the one that hit the Earth when the moon was formed.

Modeling protoplanetary systems requires a lot of coupled, high powered computers because one needs to model millions of tiny particles.

There's a simulation of the collision of the Mars sized object with the Earth when the moon was formed. Judging from the number of particles, it seems able to run on todays pc's. It ends up with the moon at the right distance from earth, and the right masses for the Earth and moon.

 

[D H]

Well don't forget Newton's law of gravity only makes sense for point masses.

Wrong. It is an trivial matter to compute the gravitational influence of a non-point mass body that has spherical shape and a uniform density on some object. An ellipsoidal body is not too much harder. We use complex mathematical models such spherical harmonics to describe real body such as the Earth. Some spherical harmonics models of the Earth's gravity field include http://cddis.nasa.gov/926/egm96/egm96.html" and more recent ones based on Lunar Prospector data. JPL has developed spherical harmonics models for Venus, Mars, and even some asteroids.

Regarding Titius-Bode that I am interested in testing, there is no reason to assume it applies exclusively to the scale of our solar system.

One of the cornerstones of the scientific method is the concept of falsification. Neptune falsifies the Titus-Bode law. Most astronomers view the Titus-Bode law as mere numerology. I suggest you post your concepts on this law in the scepticism and debunking forum at PF.

 

[Ulysees]

Newtons law of gravity does not apply to ellipsoid shapes last time I checked, only spherical shells are equivalent to point masses, and only on the outside (shell theorem).

But if you have proof that ellipsoid shapes are equivalent to point masses too, as required by Newton's law, I'd like to see it.

Or if you have proof that ALL spherical harmonics are equivalent to point masses for Newton's law to apply, I'd like to see it.

Cause what I have seen from my numerical calculations (in thread
Numerical integration of gravity at surface of an object), Newton's law of gravity only matches spherical shells. And therefore it can only be used exactly for spherically symmetrical distributions of matter like an onion, everything else is by appoximation.

 

[D H]

I never said ellipsoidal shapes are equivalent to point masses. I said one can compute the gravitational attraction toward an ellipsoidal body. The gravitational attraction at some point \mathbf x outside some body is volume integral

<br /> \mathbf a(\mathbf x) = G\int_V \frac{\rho(\mathbf {\xi})}{||\mathbf {\xi}-\mathbf x||^3}(\mathbf {\xi}-\mathbf x)\,d{\mathbf {\xi}<br />

where the integration is taken over the massive body in question.

 

[tony873004]

Newton's Laws: (from memory)
1. An object at rest remains at rest, or an object in motion stays in linear constant motion unless acted upon by a force.

2. F=ma
3. For every action, there is an equal but opposite reaction.

There's nothing here that excludes ellipsoids. Newton spent a lot of time developing Calculus simply so he could mathamatically justify his assumption that a sphere simplifies into a point mass. He did this by integrating the volume of a sphere. Now we can save a lot of time by using this nice shortcut. But there's no reason you can't integrate an ellipsoid, or any other shape, without going beyond Newton's laws. And if you wanted to do it numerically, if I created an ellipsoid from a collection of 1 million point masses, and then simulated a particle orbiting the ellipsoid, using Newton's laws, I imagine it would work quite nicely, showing the proper orbital precession that gives us "walking" orbits and "sun synchronous" orbits.

 

[Ulysees]

I never said ellipsoidal shapes are equivalent to point masses.

You said "wrong" below this:

Well don't forget Newton's law of gravity only makes sense for point masses.

Wrong.