Improving the accuracy of a Java simulation of n-orbital bodies

1. Feb 20, 2010

CyJackX

1. The problem statement, all variables and given/known data

So, I'm creating a 2d graphical simulation of orbiting bodies.
Right now I'm just working with 2 bodies, one stationary, but I would hope to extend this program to be able to include any number of moving bodies and have them interact accordingly.

After combing wikipedia's articles on orbiting bodies and this very helpful site,
http://www.braeunig.us/space/orbmech.htm" [Broken]
I have not been able to come to a conclusion of what I need to do.

Very shortly after beginning my puzzling, I learned that the n-body problem is a terrifying ordeal, at least, to me. But, my quest is not for a purely accurate representation, but something that will be, ideally, computationally less intense.

The first thing that came to my mind was to calculate the Force of Gravity acceleration between the two bodies every time interval, then adjust their velocities and positions accordingly using the classic accelerated motion equations. This works, but naturally, because this calculates linear acceleration, essentially moving my object in a bunch of tiny straight lines that form an oval, it is very inaccurate unless dt starts approaching 0. Under this formula, I cannot hold an orbit for very long; it only takes several periods before the satellite falls into the stationary planet.

Ideally, what I had in mind BEFORE seeing the highly inaccurate results, I wanted the n-bodies to calculate their accelerations based on their gravitational pull to each other, then adjust their positions and velocity.

This, I believe, would give an acceleration to each particle that doesn't necessarily point to a specific body, and each instant this point would be changing as all the bodies realigned themselves again.

So, after all this backstory, my question is this:

How does one determine the next position of a satellite if they are given only the satellite's position, velocity, acceleration, and change in time? Assuming that the direction of acceleration will remain the same. (It won't, in actuality, because of all the moving bodies, but I believe that finding a formula that gives mini-curved paths instead of mini-linear paths will be more accurate.)

2. Relevant equations

F = GM1M2/r2
I'm just using the basic gravitational equation to determine accelerations.

I've been using standard uniform acceleration equations for the movements.

3. The attempt at a solution

See above. I'm sorry my explanation kind of eschewed the template, but I don't think it's the kind of question that most commonly gets asked here.

Last edited by a moderator: May 4, 2017
2. Feb 20, 2010

CyJackX

I suddenly thought of something that might help me. I hadn't thought of it yet, but finding the center of gravity of all the other planets would give me a relative position to work use! So for n-bodies, treat one satellite as one body and all the other bodies as a cumulative point mass!

But then I still don't know how to get the orbit from that, although I'm sure it's much easier when you have both points and the satellite's acceleration and velocity...

3. Feb 21, 2010

Ben Niehoff

Finding the center of gravity of all the other bodies doesn't help extremely much, because when you add up the forces from each body individually, you are doing the exact same thing.

The naive method to approach this kind of problem is called the Euler method: basically, you calculate accelerations, and using the current velocities and positions, you apply the standard kinematics equation:

$$x(t + \Delta t) = x(t) + v(t) \Delta t + \frac12 a(t) \Delta t^2$$

But as you have observed, the errors quickly build up this way, and you are likely to run into instabilities.

A better algorithm is called the Velocity Verlet algorithm. This is not necessarily more accurate, but it is more stable (i.e., the errors tend to oscillate between positive and negative, thus preventing runaway solutions). In this algorithm, each time step is divided into halves, and slightly different operations are done at each half step. Look it up to learn more about it. Velocity Verlet is one of what are known as "symplectic integration" algorithms; that is, the algorithm is designed to explicitly conserve energy, so that even if you don't get a correct solution, at least the energy is constant and the solution will not run away.

4. Feb 21, 2010

CyJackX

Yes, Euler method. I knew that system of approximation had a name.

I will try to implement the velocity verlet algorithm, but I am still wondering, is it possible to find exactly the next position of a satellite in relation to a stationary body if you know both masses, positions, the velocity, and the acceleration?

5. Feb 21, 2010

Ben Niehoff

If you only have two bodies, then you have an explicit solution already. But you're looking at n-body dynamics, so no, the only option in general is approximation.

Think of it as a Taylor series. If you look at my Euler's method formula above, it's just the Taylor series for the next position, expanded up to second order in terms of t. But as you know, the Taylor series is in general infinite. In order to know exactly the next position, you need to expand the Taylor series to ALL orders; i.e., you must know not just the velocity and acceleration, but also all higher-order derivatives of the position vector as well.

There are more accurate methods than velocity Verlet, such as Runge-Kutta. But the more accuracy you want, the more computational overhead needed for each time step. If your only goal is to produce pretty pictures, I've done great-looking n-body simulations with velocity-verlet with up to 500 bodies.

6. Feb 21, 2010

D H

Staff Emeritus
Welcome to PhysicForums, CyJackX!

Ben has made some good suggestions. A few more:
1. Don't have any of the bodies stationary. That approximation is valid only in the extreme case where the central body is many orders of magnitude more massive than any of the other objects in the system. Consider our solar system. This approximation is only fair. Jupiter's mass is about 1/1000 that of the Sun. A decent solar system simulation will have all of the bodies in motion.
2. Go to three dimensions.

The approach you take to numerically solve the equations of motion will depend on several factors. These include
1. How many bodies are involve in the simulation.
Thousands of bodies or more will require some very advanced techniques. With a smallish number of bodies, say a few dozen or less, you can chase after accuracy and precision. With an intermediate number, such as Ben's 500 body simulation, you can use simple technique such as Verlet, but the accuracy will suffer.

2. How long of a span of time (simulated time) you want the sim to cover.
Some numerical integration techniques more-or-less conserve energy, others definitely do not. Consider two variations on the simplest of the integration techniques, the Euler method. The basic Euler method involves propagating position and velocity via $\mathbf x(t+\Delta t) = \mathbf x(t) + \mathbf v(t)\Delta t$ and $\mathbf v(t+\Delta t) = \mathbf v(t) + \mathbf a(t)\Delta t$. Suppose you do this the other way around -- update velocity and then use the updated velocity to update position: $\mathbf v(t+\Delta t) = \mathbf v(t) + \mathbf a(t)\Delta t$ and $\mathbf x(t+\Delta t) = \mathbf x(t) + \mathbf v(t+\Delta t)\Delta t$. This simple change makes a *huge* difference. The first approach, basic Euler, does not conserve energy. The planets will spiral out of the solar system. The second approach, symplectic Euler, comes very close to conserving energy. Orbits propagated with the symplectic Euler approach will be stable for a long time. While the propagation is stable with this approach, the accuracy is rather lousy. The planet's won't spiral out, but they won't be in the right place.

3. How accurate you want your simulation to be.
Ideally your simulation will be accuracy and stable. That's usually a pipe dream. Ben mentioned Runge Kutta techniques as an example of high accuracy. While RK integrators are very good, they typically do not conserve energy. The same spiraling problem that occurs with the basic Euler approach also happens with the Runge Kutta techniques. Even more advanced techniques will do a better job of balancing accuracy and stability.

4. Whether you can play games with time.
Every numerical technique has some optimal time step. If you have done some reading, you will have found statements such as Euler having $O(\Delta t)$ accuracy, velocity Verlet $O(\Delta t^2)$ accuracy, and so on. Making the time step smaller will improve the accuracy -- up to a point. Make the time step extremely small and the planets will stop moving in your simulation. All of the integration techniques involve a basic step along the lines of $x(t+\Delta t) = x(t) + \Delta x$. There is a basic problem with the way computers do floating point arithmetic. Make that $\Delta x$ too small and x won't change one bit.

What this means is that every numerical integration technique has some optimal time step. For example, consider a satellite orbiting the Earth in low Earth orbit. The best accuracy with fourth order Runge Kutta is achieved with a time step of about a second or so, or about 5400 steps per orbit. Making the simulation time step much larger than this optimal value and the errors inherent in the RK4 technique will start to creep in. Making the time step much smaller than this optimal value and numerical computation errors start to creep in.

This becomes a problem with propagating our solar system. Mercury's orbital rate is about 1000 times that of Pluto. The optimal time step for propagating Mercury's orbit is suboptimal (too fast) when it comes to propagating Pluto's orbit, while the optimal time step for propagating Pluto's orbit is suboptimal (too slow) for propagating Mercury's orbit. One way around this problem is to use different time steps for the different bodies. Doing this right can be rather tricky.

As far as numerical techniques go, the first step is to use understand how/why the basic Euler works because almost all integration techniques can be viewed as refinements to this basic approach. The very next step is to throw this basic Euler approach out with the trash. While understanding it is important, avoiding the use of it is also very important. Symplectic Euler is just as simple to set up and is a lot more accurate.

With regard to more advanced techniques, Ben mentioned velocity Verlet. I don't particularly like that one because of the (my opinion) goofy half-time steps for velocity. An even better that is essentially just as simple as velocity verlet is Beeman's algorithm. I particularly like the predictor-corrector form of Beeman's algorithm.

The simulations I work with tend to use fourth-order Runge Kutta. This is mostly a domain-specific choice. When simulating a spacecraft's flight computer is a part of the simulation, the rate at which the spacecraft operates pretty much dictates the simulation time step. If that dang flight computer weren't in the simulation I would use much better techniques. Some techniques such as Encke-Nystrom techniques and the Bond-Gottlieb 14-element method can accurately and stably propagate an orbit with only a few steps per orbit.

Last edited: Feb 21, 2010
7. Feb 21, 2010

CyJackX

Alright, I will try RK4, but if it does not preserve energy, should I use it? What if my system will end up changing a lot?

Also, where can I find more instructional documentation on BG14 and Encke-Nystrom?

8. Jan 28, 2011

jaiii

Could you please send me this simulator as PM or attachment, or by e-mail.

Thank you very much.