Create an Orbital Simulator

[bkelly]

I am toying with the idea of creating a general orbital simulator. Possibly lacking the proper terminology, here is a basic description of my goal, make that wild dream.

Step 1. Working in 2 dimensions, user enters the mass of 1 to N object, their initial position, their initial velocity, and a time calculation interval. For each time interval, the simulator calculates the forces that each body applies to each other body, the summed vector of that force, applies it to the initial velocity, and determine the position at the end of the time interval. It then shows the position of each body. Repeat.

Step 2. Add niceties like lines showing where each body has been for some number of iterations so we can better see the activity.

Some future step, make it three dimensional.

This seems a rather ambitious project so before I start:
Do you know of any simulators of this nature?
Do you have any advice for a novice in celestial mechanics?
Do you think this would be better to attempt on Microsoft Windows or under Linux, I have both.

Thank you.

 

[Nabeshin]

I contemplated doing this but found it to be beyond my current physics (and probably computer programming) abilities. I'd be really interesting to know the method you use to accomplish all this!

 

[D H]

Step 1. Working in 2 dimensions, user enters the mass of 1 to N object, their initial position, their initial velocity, and a time calculation interval. For each time interval, the simulator calculates the forces that each body applies to each other body, the summed vector of that force, applies it to the initial velocity, and determine the position at the end of the time interval. It then shows the position of each body. Repeat.

This approach will work fine with one object. The object will move in a straight line because there no forces act on it. This approach will not work very well with two or more objects. The approach you outlined is called the Euler method for solving an initial value problem (i.e., a system whose complete system state is known at some point and whose evolution is described by a set of ordinary differential equations). The Euler method is very simple to implement but unfortunately does an incredibly lousy job except for the simplest of systems.

Fortunately, a lot of numerical integration techniques exist that do yield good accuracy. They are all more complex than the Euler method, but many are only a bit more complex. A couple of nice balances between complexity and accuracy are the velocity verlet technique, which is widely favored in chemistry simulations, and the fourth-order Runge-Kutta (RK4) technique, which is widely favored in orbital simulations. A quick google on velocity verlet and Runge-Kutta will yield a lot of information.

 

[bkelly]

Hello D H,
We quickly delve into the reasons for my post. I don’t understand why the use of the simple equations f = ma, and f = ( G * M1 * M2 ) / r squared will not work. Comparing two bodies, there is a force between them that will change their velocity. If the masses and velocities are anywhere close to “right” then the two masses should orbit. Assume I iterate fast enough to get sufficient resolution but slow enough to prevent thrashing with no results. I will probably start with a resolution of one hour, 8,760 iterations will simulate a year on the scale of our solar system.

I looked at Runge-Kutta on Wiki and do not understand it at all. I am not trying to predict or smooth our trajectories, just create an instance and calculate the values for the next instance.

My first hurdle will be how to express mass in an equation. I will be looking that up.

Thanks for taking the time to reply.

 

[LoneDragon]

Runge Kutta, Burlish Stoer, and Euclid

Both bkelly and DH are correct. If you use small timesteps (but very very small ones), the basic f=ma for N^2 interactions, will produce a proper simulation. But if two objects in that N body simulation *ever* approach closely enough so that the time-space resolution causes "thrashing" or a large error term in the step wise calculation, the simulation will only be characteristically correct, but will diverge greatly from reality. For the solar system, at a one hour resolution, you should have no problems for the solar system.

Runge-Kutta (RK) methods do not calculate smoothed out orbital calculations, or predictions, per-se. Rather, the multiple differnce terms you see in Runge-Kutta calculations, are for calculaing a quadratically accurate polynomial *approximation* of the orbit. Euler only calculates a first order polynomial *approximation* of the orbit.

Please read about Taylor series, rectangular integration rule, simpson's rule of integration, Euler method, and RK method, in detail and trials in code, to understand how Euler method, is like a first order rectangular integration, and RK on the fourth order method, is somewhat-like simpson's rule that calculates mini-chained parabolic integrals each step of the second order. It all has to do with the infinite Taylor series, and error terms of asymptotic convergent functions like gravitational fields. With RK, you could make the steps 3-5 days long, and be far more accurate with Mercury and close approaching comets, than you would using the Euler method. RK even comes in higher order methods that let you further increase the step size, as each step is even higher order accurate, with an exponentially shrinking residue in the Taylor series like infinite series. Trust me, if you learn RK, you will never regret it. Numerical Recipies in C, by Press, Teukolsky, et al, is a good book with such algoritms available in ready to use C code for either Linux or Microsoft C compiler environments.

Another method you may want to consider, if RK is too complex, is to use Burlish-Stoer (BS) methods. In short, first, you can take Euler method at 16 hours*1 step, 8*2, 4*4, 2*8, and 1 hours*16 step. Then do a least squares fitting of the 5 exponential refined ranges to extrapolate the asymptotic-convergent value using an appropriate asymptotic function, which you can vary it's parameters, until it fits the projected 16, 8, 4, 2, 1 series to find the convergent value at "0 hours". Then take that convergent value and take the 16 hour step at the infinitessimally projected 16 hour time step, using the asymptotic (0 hours*infinite step) projected step with approximated "perfect" precision. You will have to find the right equation for this application for the least squares fitting of the exponential approximations, and make sure the exponential ranged projections are not oscilatory, but asymptotic, so that the curve fitting is always accurate each step.

Incidentally, make sure any Euler, RK, or BS curve fit and planetary position and velocities are all done in double precision, and even then, that your numbers retain large precision in the LSB of the numbers. Study the IEEE (754) format for float and the one for doubles, paying special attention to the mantissa for precision checking. The Euler method, especially, if the time steps are too small, may become susceptible to quantization errors in extreme distance or closeness calculations, without checking precision in the mantissa. RK is the best, and BS can be used intermediately.

 

[Phrak]

Step 1. Working in 2 dimensions, user enters the mass of 1 to N object, their initial position, their initial velocity, and a time calculation interval. For each time interval, the simulator calculates the forces that each body applies to each other body, the summed vector of that force, applies it to the initial velocity, and determine the position at the end of the time interval. It then shows the position of each body. Repeat.

I wrote an 3D, N-body orbital simulator, with displayed output, in basically the same manner you outline. Variables were updated, for the most part, in parallel. The model suffered orbits that would gradually increase in radius over iteration. I repaired it by breaking the parallel-set rule. Near collisions introduced a problem that I never did fix before loosing interest. A close approach would throw the two bodies out towards infinity. I attempted a conservation-of-energy kludge to fix the problem. It was a bust.

 

[bkelly]

Wow, that was a lot to say. This will be more complicated that I imagined, even considering that I knew it would be more complicated than I knew. (Does anyone know of a happy face icon with a recursive grin?)

Simulating a solar system would be nice, but what I really want is the ability to start up a simulation with 3 or more stars in orbit around each other and see what happens. As noted by LoneDragon, that will involve some close approaches, high velocities, and obviously some time-space alterations (as compared to our perceptions of what is normal).

After reading Phrak’s post, I had best expect some anomalies in the behavior.

My thought is to use C++ or C# and create a class for celestial objects. For each body in a simulation, instantiate an object giving it the requisite mass, size, location, and velocity. (and other things I don’t know about yet.) The class would have a method that takes as its input parameter a pointer to another object. That method would then compare the parameters of the two bodies and calculate the resultant force at that point in time. In doing so, it would calculate the position and velocity for some specified time in the future. That time would be variable and should probably decrease based on distance, velocity, and I don’t know what else.

It seems like if I start it simple without close approaches, then when I get the mechanics of being able to calculate the resultant vectors of two bodies, and displaying an image of the bodies, then I can refine the equation within the class to any degree I want. Obvious performance will suffer, but that is the way it goes.

That reminds me. Back in the 1980s I created a program to display views of the Mandelbrot set. (This was when the standard speed was 4 MHz and a fast machine was 8 MHz.) I quickly had to use the 8087 for precision and speed. Using Turbo Pascal, I then dropped into assembly language to program the 8087 for speed and more accuracy. Finally when that was not enough resolution, I wrote some 128 bit fixed point arithmetic routines in assembly language and ran that. I never did find the threads between what appeared to be islands of points within the set.

But I digress.

I have seen images from simulations for galaxies, but I don’t know of any right now. It would be interesting to know how they ran those simulations, on what machines, what language, and how long it tool. A bit of searching is in my immediate future.

Thanks for your replies.

 

[D H]

Both bkelly and DH are correct. If you use small timesteps (but very very small ones), the basic f=ma for N^2 interactions, will produce a proper simulation.

No, it won't. I work in this field and we have very explicit rules barring the use of the Euler method. The basic problem is that the Euler method conserves neither angular momentum nor energy. (RK4 has this problem too, but the error is vastly diminished compared to the Euler method.)

When writing any simulator, it is always a good idea to verify and validate the resultant simulation. One way to do this is to apply simulation to a case where the correct answer can be computed analytically. The case of 2 bodies has an analytical solution in Kepler's laws. The simplest case is two bodies in circular orbits about each other. The problem with the Euler method becomes apparent in this simple simulation: The bodies spiral away from one another. One can reduce the error by taking very small steps. Suppose you place two bodies in circular orbits about each other, simulate for what should be ten orbits, and examine the relative position error. To get a relative error of 1% or smaller, an Euler integrator would have to divide that 10 orbit time interval into well over 3 million subintervals. In comparison, an RK4 integrator would only need to divide that 10 orbit time interval into about 400 subintervals. Achieving a relative error of less than 0.1% requires about 900 subintervals for RK4 and is simply unachievable with Euler integration. The problem is that as one makes the time steps smaller and smaller, at some point the error starts to grow because (1+10^-16)-1 = 0 using the ordinary double-precision arithmetic available on almost all computers now. Make the steps too small and the state freezes.

Edited to add:
RK4 is not the be-all and end-all of numerical integration techniques. It is the most widely-used numerical integration technique because it is easy to implement and offers fairly good accuracy. There are more complex Runge-Kutta techniques that do a much better job than RK4. For example, a Runge-Kutta-Feldberg 7/8 integrator achieves better than 10^-12 relative accuracy for the 10 orbit scenario described above by dividing the 10 orbit interval into about 1000 intervals. In comparison, a RK4 integrator run with the same time step will achieve a relative accuracy of about 5*10^-5. The relative error with the Euler method at this coarse step size is greater than one. There are other more advanced techniques that do much, much better than any Runge-Kutta integrator.

 

[Phrak]

what LD doesn't tell you is that all n-body sumulations are only characteristically correct. how correct you want yours depends upon your intent in doing the simulation.
If you want some initial results you don't need polynomial approximations.
btw, a break-down at close approach is characterised by the displacements obtained between iterations per the radial displacement between objects (rather than the radius you assign the objects).

 

[LoneDragon]

what LD doesn't tell you is that all n-body sumulations are only characteristically correct. how correct you want yours depends upon your intent in doing the simulation.
If you want some initial results you don't need polynomial approximations.
btw, a break-down at close approach is characterised by the displacements obtained between iterations per the radial displacement between objects (rather than the radius you assign the objects).

Actually, I did. All taylor series family like methods are asymptotic approximations, with an error residue. More than once I referred to approximation, error residue, and close approach gross errors.

LD: It all has to do with the infinite Taylor series, and error terms of asymptotic convergent functions " sums all of this the best, and even higher order terms like RK4 or RK6 or RK... predictor corrector methods are all approximations. My one "proper" simulation remark, only refers to a simulation like the solar system where things don't get too close to each other and are near circular orbits, excluding moons. A complex stellar environment, with many close approaches, agreed, does require the more precise methods, like RK4 or RK6 and so forth. All of them do not preserve everything, if you want to be technical, as nearly all natural three-body problems are analytically unsolveable, and N body too.

For bkelly,
For an approximation-near-analytic text, to test a simulation of yours, there are astronomical ephemerides, with heliocentric velocity and position information. JPL has online ephemerides which may also be used to test your method, before doing complex N body stellar simulations, with close approaches, which will require higher order terms, near-accuracy methods.

 

[bkelly]

Wow. This is beginning to sound like an incredibly complex task. I don't know enough to really understand what you guys are saying, but I will continue looking to determine the odds of getting something running that remotely resembles what would happen in reality.

Going to www dot ask dot com (this thing seems to tell me I cannot post a URL so I "said" it instead) and asking what are ephemerides yeilds nothing useful. A google search tells me that they give the location of satellites for each day of the year. Haveing been a systems engineer tracking rocket launches and using tracking antennas, I know what a TLE is. The JPS site appears beyond me, but I will wander around there and see what I can understand.

However, to my knowledge, I don't want to predict anything. I want to run a simulation that takes steps based on the current conditions and calculates movements to the next position or iteration.

However, again, as I alluded earlier, I am vaguely aware of the severe limits of my knowledge in this field. The core question has become:
Concerning the task of creating a basic orbital or stellar simulator, providing the ability to designate a known starting point for each of N bodies such as stars: How difficult will it be to get something that reasonably approximates how two to a dozen stars or planets will orbit each other? This is just the mathematics of it all and ignores the problems of display on a computer screen.

 

[LoneDragon]

Ephemeris and Euler = RK1 for testing character.

Try Wikipedia. An astronomical ephemeris (blue books in the library for an american ephemeris I remember for each year) simply lists the positions and velocities of planets, and many other astronomical tables. For a test...

Euler method will *characteristically* work fine for you if all of the bodies don't have close elliptical approaches with large velocities that go though large angles in one step, which integrate increasingly improperly under those conditions in product. Then, if you have properly designed the input mass-position-velocity matrix, and next-step output mass-position-velocity matrix, and Euler step algorithm, and are happy you have worked out the bugs, you **may** replace the Euler test method step calculator, once you see initial Euler results, with a better algorithm. Modularity being key here.

Your code will work fine-enough, but inaccurately, so long as you are not trying to accurately predict, say, precise galactic arm characters, or precise galactic collision simulations.

Also, you can make a stop-gap measure that reduces the step size of close, high velocity, high angle around other body masses, by reducing the step size for *only* those bodies, leaving all of the other bodies at their normal step size. This is what is called an adaptive step size algorithm.

Incidentally, interestingly for the mathematicians only I guess, is that in the differential equations classifications, Euler's method is a first order RK, a second order modified Euler method is a second order RK (RK2), RK4 is obviously a fourth order RK and is closely related to the conceptually easy to understand Simpson's Rule for integration when done in one dimension, and there are higher order RK with more Taylor series expansion terms. But you can start with RK1 (Euler) and work your way up as you pick up the math.

I hope your simulation works well enough for you. I remember doing a solar system simulation on a commodore 64 in junior high which was a lot of "fun" (please don't calculate my age) *grins*. Do make sure to calculate all of the bodies at once, from input array to output array in ping pong buffers, for if you advance through just one array and place the answers back in the array for the next step, it causes even more errors in a simple characteristic simulation as half the bodies are advanced, say, while the others are distorted in their interactions. To save time from memory swapping, just toggle referencing the ping pong buffers input-output back and forth with pointers, if your using C. C#, yer on yer own!

regards

 

[rizzy]

Man you guys are smart. I too hope to become as smart as you guys one day. As of Right now I will use the excuse of being young. hehe

 

[LoneDragon]

Man you guys are smart. I too hope to become as smart as you guys one day. As of Right now I will use the excuse of being young. hehe

LOL. Just study hard, read *everything* you can, and take whatever classes are possible *while paying attention*. Time will do the rest. But join the club, as I couldn't even multiply numbers until 3rd grade (I know, that sucks bad) ... hmmm, there must have been something in the water, causing the math spurt over the decades since *then*. LOL, then groans at self.

For bkelly:
To reiterate, as an engineer, I can say Euler will "git 'er done", for a quick peek at simulation programming, sufficient for a *basic* gravitation simulation, but with a strong emphasis on the *basic* simulation. So, "caveat emptor", or, "let the buyer beware", if this is a college math, and coding precision test,, versus a high school, personal interest, or introductory simulation level project. Still, it is good to study the other stuff, if you have the time, when regarding a related technical career.

 

[Phrak]

Man you guys are smart. I too hope to become as smart as you guys one day. As of Right now I will use the excuse of being young. hehe

Hey! Don't include me in this category. I'm not BSing my way over someones head to give my gotta-be-smart-ego a plateform. But, hey, I could give it a whirl I 'spose. Let's all meet for tea on friday and we'll try it again.

 

[Redbelly98]

I wrote an 3D, N-body orbital simulator, with displayed output, in basically the same manner you outline. Variables were updated, for the most part, in parallel. The model suffered orbits that would gradually increase in radius over iteration. I repaired it by breaking the parallel-set rule. Near collisions introduced a problem that I never did fix before loosing interest. A close approach would throw the two bodies out towards infinity. I attempted a conservation-of-energy kludge to fix the problem. It was a bust.

Too bad yours didn't work out. I did something similar (in Visual Basic, on a Windows PC) 3 or 4 years ago and was reasonably happy with the results. I remember modifying the force so that it didn't diverge when the bodies got too close to each other, for the reason you saw: all hell breaks loose otherwise. Of course that meant some inaccuracies in the trajectories, but for example I was able to give a cool demonstration of the gravity-assist technique for boosting satellites.

If you're interested I can try to look up exactly how I did the modification to the force.

FYI, I only worked with up to 10 or 20 objects at a time, the computations really bogged down beyond that (Remember, this was Windows and Visual Basic).

 

[LoneDragon]

Hey! Don't include me in this category. I'm not BSing my way over someones head to give my gotta-be-smart-ego a plateform. But, hey, I could give it a whirl I 'spose. Let's all meet for tea on friday and we'll try it again.

Well, I am curious for some enlightenment, as you're speaking above *me*, now. Where are the inaccuracies, in the long list of posts, on this thread? Everyone seems to be on an even keel, as far as I can tell, in my humble engineer's math. I see you are an aerodynamics student, or engineer, from your other posts, so, sir, you can take a whirl if you know *that* level of simulation math, as I, at least, would like to know your own numerical opinions, and knowledge, to correct the errant information on this thread. You know, so that we may all be made higher.

 

[Redbelly98]

Wow, that was a lot to say. This will be more complicated that I imagined, even considering that I knew it would be more complicated than I knew. (Does anyone know of a happy face icon with a recursive grin?)

Maybe the best thing right now is to use some basic algorithm, that you can easily understand, to get something working first. Then worry about making it more accurate or complex as a later project.

 

[bkelly]

I found a simulator here: http :// kornelix . squarespace . com / galaxy /
(it seems I don't have URL priviledge, but google on "stellar simulation")
It runs under Linux so I fired up that computer and installed it. It runs well, but not much in the way of features. I haven't reviewed the code closely, but it seems to just use the F=MA type of approach. If I can get version 1.0 running, then I can look into better accuracy for close approaches.

There are two primary segments of this projects
Calculating the positions
Displaying the field
In this forum I want to discuss the calculations.

What type of coordinate system would be best, polar, Cartesian, or something else I don't know about? It seems to me that Cartesian would be better. If for no other reason, I have never messed with polar coordinates. But if polar would be much better, I would like to know.

If Cartesian, I suspect that it would be easier to deal with if the origin is in one of the corners rather than in the middle. My first thought would be the lower left corner. Thoughts about selecting the range: Astronomical Units, millions of miles, something else?

BTW: This site seems to respond much faster under Linux with Firefox than under XP with IE-7.

 

[Phrak]

Too bad yours didn't work out. I did something similar (in Visual Basic, on a Windows PC) 3 or 4 years ago and was reasonably happy with the results. I remember modifying the force so that it didn't diverge when the bodies got too close to each other, for the reason you saw: all hell breaks loose otherwise. Of course that meant some inaccuracies in the trajectories, but for example I was able to give a cool demonstration of the gravity-assist technique for boosting satellites.

If you're interested I can try to look up exactly how I did the modification to the force.

FYI, I only worked with up to 10 or 20 objects at a time, the computations really bogged down beyond that (Remember, this was Windows and Visual Basic).

bkelley's the one you want to talk to about code. I wrote mine a few moons ago before they has invented rocks. Did you solve or notice the walk-away problem where circular orbits are really outward spirals? You can check if it's happening with a running total-energy value.

I hope I didn't PO lonedragon excessively...

 

[Redbelly98]

bkelly,

In Cartesian coordinates it's easier to convert from the (x,y) calculated positions to displayed pixel locations, as only simple scaling factors are involved. Might also be easier to put the origin at pixel (0,0) in the display in the upper left, but if you don't it's still pretty straightforward to rescale from (x,y) to pixel location.

I remember an algorithm a former physics teacher gave me. Better than Euler, not as accurate as Runge-Kutta, but pretty straightforward in terms of understanding it.

In 1-D it works like this ...

Do a Taylor series expansion to find $$x_{n+1}$$ given $$x_n$$ , where $$x_n = x(n \cdot \Delta t)$$ :

$$x_{n+1} = x_n + v_n \cdot \Delta t + \frac{1}{2} a_n \cdot \Delta t^2 + \frac{1}{6}x'''_n \cdot \Delta t^3 +$$ ...

Similarly,

$$x_{n-1} = x_n - v_n \cdot \Delta t + \frac{1}{2} a_n \cdot \Delta t^2 - \frac{1}{6}x'''_n \cdot \Delta t^3 +$$ ...

where $$v_n, a_n, x'''_n$$ are the velocity, acceleration, and 3rd derivitive of x, respectively, at time $$t= n \cdot \Delta t$$

Add these two equations:

$$x_{n+1} + x_{n-1} = 2 \cdot x_n + a_n \cdot \Delta t^2 +$$ ...

or

$$x_{n+1} = 2 \cdot x_n - x_{n-1} + a_n \cdot \Delta t^2 +$$ ...

Note the cancellation of the odd-order terms, so that the error is actually of order $$\Delta t^4$$, not $$\Delta t^3$$. All you need to calculate $$x_{n+1}$$ are the previous 2 values of x, as well as the acceleration at time $$n \cdot \Delta t$$, which you calculate from the Forces and masses in the system.

In 2-D you'd need to calculate the x- and y-components of the accelerations, then apply the algorithm to find $$x_{n+1}$$ and $$y_{n+1}$$ .

There might be a common name for this algorithm, if anybody knows it feel free to post the info.

 

[Redbelly98]

bkelley's the one you want to talk to about code. I wrote mine a few moons ago before they has invented rocks. Did you solve or notice the walk-away problem where circular orbits are really outward spirals? You can check if it's happening with a running total-energy value.

I didn't check circular orbits, that would have been a good test. What I did notice was that if two objects got very close, they could suddenly go flying off the screen at incredibly high speed.

p.s. I just looked over my code, and found that I used the algorithm I described in message #21. Somehow I thought I later changed this to Runge-Kutta but apparently not.
I can't figure out, from a quick glance at the code, just what I'm doing to the force to make it "behave" at close range. I would have to find some 3-year-old handwritten notes that I have somewhere to figure that out again.

 

[Shukie]

Anyone feel like sharing one of their simulators? Would be interesting.

 

[Phrak]

The guys at NASA, NSA, and whoever, should have devoted quit a few man-years at this. So they may have found a different method than interating t, and updating V_i and X_i. There's the Hamiltonian, if you've heard of it. It's an alternative, equivalent to Newtonian mechanics. What I'm trying to say, is that it may be possible to iterate values in some other mathematical space, rather than velocity and position. Angular momentum and Kinetic energy perhaps; I have no idea really. Just offinging possibilities.

 

[D H]

I wrote an 3D, N-body orbital simulator, with displayed output, in basically the same manner you outline. Variables were updated, for the most part, in parallel. The model suffered orbits that would gradually increase in radius over iteration.

This is a direct consequence of using the Euler method. The Euler method displays this problem even with two bodies and even without any eccentricity.

I remember an algorithm a former physics teacher gave me. Better than Euler, not as accurate as Runge-Kutta, but pretty straightforward in terms of understanding it. ... There might be a common name for this algorithm, if anybody knows it feel free to post the info.

The method you described is the basic verlet method. As you noted, it is better than Euler. An even better approach and still very simple to understand is Heun's method.

I found a simulator here: http:// kornelix.squarespace.com/galaxy/ [Broken]. It runs well, but not much in the way of features. I haven't reviewed the code closely, but it seems to just use the F=MA type of approach.

The cited simulator updates the velocity and then updates the position using the updated velocity. This is not a very good integrator. A simple improvement is to update the position using the average of the velocity before and velocity after the velocity update. This simple improvement is Heun's method, a second-order Runge-Kutta integrator. Heun's method is considerably better than the midpoint method, another second-order Runge-Kutta integrator.

I haven't reviewed the code closely, but it seems to just use the F=MA type of approach.

*All* of the Runge-Kutta techniques as well as the verlet method "just use the F=Ma" type of approach. Some just do it a lot better than others.

I didn't check circular orbits, that would have been a good test.

That is always a good test. A simulator that can't do a very good job reproducing a circular orbit isn't worth very much at all. Highly elliptical orbits and hyperbolic orbits are also good tests because the close approaches will stress even good integrators.

 

[D H]

The guys at NASA, NSA, and whoever, should have devoted quit a few man-years at this. So they may have found a different method than interating t, and updating V_i and X_i. There's the Hamiltonian, if you've heard of it. It's an alternative, equivalent to Newtonian mechanics.

NSA, as far as I know, doesn't address this issue. OTOH, NASA certainly does. Hamiltonian-based techniques are called symplectic integrators; google this phrase for more info. The advantage is that these techniques inherently do a better job at conserving energy than do non-symplectic integrators (e.g., the Euler method and the standard RK4 method). The basic verlet method is a symplectic integrator, for example, and it does a much better job than does the Euler method even though both evaluate state derivatives once for each propagation interval. I believe JPL uses symplectic integration techniques for developing ephemerides.

There are many, many other orbital propagation techniques. One approach is to use Keplerian elements (e.g., $a,e,i,\Omega,\omega,M$) rather than position and velocity. In the case of two bodies, all Keplerian elements but the mean anomaly are constant and the mean anomaly is a linear function of time. Propagation is easy and exact. Making sense out of the elements (e.g., plot positions on a computer screen) takes a little work, but not too much. Now consider a stellar system which comprises a single star and several planets. The planets' orbits can still be described in terms of Keplerian elements, but these now exhibit time-varying behavior. The variations are small, so perturbation techniques can be put into play. NASA used such an approach during the Apollo era, and still uses it in some arenas today.

 

[Phrak]

... There are many, many other orbital propagation techniques. One approach is to use Keplerian elements (e.g., $a,e,i,\Omega,\omega,M$) rather than position and velocity. In the case of two bodies, all Keplerian elements but the mean anomaly are constant and the mean anomaly is a linear function of time. Propagation is easy and exact. Making sense out of the elements (e.g., plot positions on a computer screen) takes a little work, but not too much. Now consider a stellar system which comprises a single star and several planets. The planets' orbits can still be described in terms of Keplerian elements, but these now exhibit time-varying behavior. The variations are small, so perturbation techniques can be put into play. NASA used such an approach during the Apollo era, and still uses it in some arenas today.

I hadn't considered perturbations from elliptical orbits, but then again, I haven't considered orbital simulations in a long while.

An interesting possibility to avoid close-approach divergence for any pair of bodies, might be to switch to elliptical orbits + perturbations (for the pair only) based upon some nearness condition, then switch back. It may even be possible to survive a direct hit where the mass centers pass through one another.

 

[Phrak]

In the above I failed to mention the possibility of close encounders of hyperbolic and parabolic trajectories where these might also be degenerate (head on collisions).

 

[lzkelley]

I'm also trying to make an n-body orbital simulator, its going quite well - I've recently upgraded from euler to fourth order runga-kutta.
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!

 

[csprof2000]

Regardless of what everyone is telling you, your method will work fine. Just be prepared to use very small time steps if you are using Euler's method.

Besides inaccuracies that can occur due to the approximate nature of the solution, the near-collision "explosion" can be solved by simply defining the force piecewise. Calculate the distance between each pair of bodies (this must be done anyway) and see if it's less than some minimum value (this can be as simple or as complicated as you like). Then decide whether or not you want to "do" the force there, and act accordingly. Otherwise, you could use a similar method to use conservation of momentum to have them collide (elastically or inelastically) like billiard balls. Or have them lump together and run until all masses have lumped together.

There's a lot of potential for fun and good programming in what you're talking about.Good graphics can be done in Java, C, C++, C#, Visual Basic, etc. For 3-D, you may even look into DirectX or OpenGL.

The programming won't be hard, but there are some things you can do to make it easy on yourself. Look for libraries with functions / objects you can use. I recommend taking the opportunity to explore OO approaches, if you're not already familiar with these. You can make objects for vectors and massy particles and this will greatly ease the notation. Have fun with it.

 

[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