Create an Orbital Simulator

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 web site 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

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

tony873004

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.


This next image compares the the precession of Mercury's perihelion between the Newtonian-predicted value of 531 arcseconds per century, Euler's method, and a 2nd-order method:

Euler's method actually out-performs the 2nd order method. It requires time steps of less than 64 seconds before the 2nd order method begins to catch up. But notice that every time the time step is halved, the 2nd order method's error is cut by 4x, and Euler's method is only cut 2x. So there will reach a point where 2nd order is more accurate, but it will be when the time step is very slow.

These next two images show a massless object on a chaotic trajectory in the Earth/Moon system. The simulation is run for a few months, then the color of the object is changed from magenta to green, and time is reversed. If everything is perfect, it will retrace its same path. The time step used is 16 seconds.



You'll notice that the 2nd order method is nearly perfect, while Euler is not. Is this a victory for the 2nd order method? I'm not sure. Some routines have the feature that they're perfectly time-reversable. So even if there are errors, the errors will be mirrored exactly when time is reversed. I'm not sure if my 2nd order routine has this feature. I know that Euler does not.

Just for fun, here's my failed attempt at the Leapfrog method. I obviously did something wrong


Here's a graph created with Gravity Simulator (Euler) showing Earth's changing eccentricity over the course of 1 million years. It clearly shows the small 100,000 year cycles and the larger 400,000 year cycle predicted by analytical methods. Earth's orbit is chaotic enough that I would imagine that the graph would be quickly out of phase regardless of the method used, although the 100,000 year and 400,000 cycles would remain evident despite the phase error.

D H

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


Euler-Cromer method


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

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

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

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 wikipedia article on RK4 is compact and accurate. Other articles on Runge-Kutta include http://www.nrbook.com/a/bookcpdf/c16-1.pdf and http://numericalmethods.eng.usf.edu/...8ode_ebook.pdf.
• Chin and Forest-Ruth. You can find Chin's papers at arXiv.org. 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

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

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/gravit...onbuilder.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 definately 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.


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 wont get an ejection no matter how long I simulate. It takes a pretty big time step to get millions of years into the future.