How to simulate planets orbit curce around the sun ?

[kjensen]

How to simulate planets orbit curve around the sun ?

Hi

I'm about to write a java program where I wish to simulate planets orbit curve around the sun, so that I can set the speed of a planet and mass of planet and sun and simulate the curves around the sun (just basic simulations). I only know high-school level physics but I know university level mathematics. I want to build my programs on Newtons laws. I know that there are plenty of others programs that can do such simulations, but I wish to know the theory and math behind it so that I can build my program from ground up. I know Newton laws of gravity, but I do not understand how to apply them to a simulation, and when I fool around with the formulas I manage to get different results depending on how I interpreted applying the formulas to a simulation. I hope that someone here will be so kind to explain me the theory of doing such a simulation or tell me where I can find books, documents, or websites about it. I have tried searching on the net, but have just got more confused. Thanks.

 

[dacruick]

Hi there,

Keplers laws of planetary motion will explain to you the physics behind the orbit, which I think is probably all you need to get a Johandle on the situation. excuse the cheesy pun.

http://en.wikipedia.org/wiki/Kepler's_laws_of_planetary_motion

One thing I notice is that you say you can set the speed of the planet and the mass of the planet, but since planets travel in elliptical orbits, the speed of the planet will be a function of the distance from the sun. Setting the speed of the planet will not suffice for a stable model. There is a value called eccentricity which depicts how elliptical an orbit is. an eccentricity of 0 is a circular orbit, but you will come across all of this if you read through the link I posted.

 

[DaveC426913]

In terms of simulating:
Decide on your accuracy (some algorthms are harder than others but more accurate).
Decide if you want it 2D or 3D (recommended: 2D)
Decide if the orbiting bodies have significant mass compared to the parent (if no then your job is easer. If yes, you can simulate star clusters).

What you do is make nested iterative loops.

Plot all points in orbits, putting their xy coords in an array
Calc the force vector pulling on them from the parent (prop to mass, inverse prop to sqaure of dist)
From the force, calc the change in velocity.
From that, calc the new position.
Do this for all point, then update the display and repeat.

 

[kjensen]

One thing I notice is that you say you can set the speed of the planet and the mass of the planet, but since planets travel in elliptical orbits, the speed of the planet will be a function of the distance from the sun. Setting the speed of the planet will not suffice for a stable model.

Thanks a lot for your reply. But I thought that the shape of the curve is depending on the total mechanical energy of the system E_mek which is equal to E_kin + E_pot (E_mek = E_kin + E_pot = 1/2*G*(M*m)/r

Therefore I thought that the curve form is dependent on the speed and/or the mass of planet (and/or the mass of the sun)

My plan was to be apple to simulate all sort of curves (comets and so on) including hyperbolic and parabolic curves with the speed and mass as starting conditions. But maybe that is possible with Keplers laws, even though I found it hard to see how to make a simulation out of it!

 

[Grizzled]

In terms of simulating:
Plot all points in orbits, putting their xy coords in an array
Calc the force vector pulling on them from the parent (prop to mass, inverse prop to sqaure of dist)
From the force, calc the change in velocity.
From that, calc the new position.
Do this for all point, then update the display and repeat.

... and if you stop at that, you will soon discover that your simulated planets have wandered off, or crashed into the Sun, or did something else equally bizzare. The worst possibility (as regards consequences of the error) is that results may look credible enough to pass cursory examination but be nevertheless wrong.

The reason for the above is accumulated systemic error. Increments (distances between the points on your graph) are by necessity finite. That means that the force, acceleration etc you calculated at the start of next step is wrong for the rest of that step. Even if you are crafty and use a central spot or several of them - that's still not enough. It will only reduce the error but not eliminate it.

To make sure your planets behave, you need to ensure that (at every step of your calculations!) they comply with conservation laws of energy and momentum. Force-adjust the results at every step to make sure your model at least doesn't break the natural laws. Not too much anyway :-)

 

[kjensen]

In terms of simulating:
Decide on your accuracy (some algorthms are harder than others but more accurate).

Thanks a lot for your reply. But does that mean that some algorithms might give different shapes? I would to make it as precise as possible (only 2D for now).

What you do is make nested iterative loops.
Plot all points in orbits, putting their xy coords in an array
Calc the force vector pulling on them from the parent (prop to mass, inverse prop to sqaure of dist)
From the force, calc the change in velocity.
From that, calc the new position.
Do this for all point, then update the display and repeat.

So that should give the same result no matter the step size? If so then I'm not sure why that is, but then I have enough to make the program.

 

[kjensen]

... and if you stop at that, you will soon discover that your simulated planets have wandered off, or crashed into the Sun, or did something else equally bizzare. The worst possibility (as regards consequences of the error) is that results may look credible enough to pass cursory examination but be nevertheless wrong.

The reason for the above is accumulated systemic error. Increments (distances between the points on your graph) are by necessity finite. That means that the force, acceleration etc you calculated at the start of next step is wrong for the rest of that step. Even if you are crafty and use a central spot or several of them - that's still not enough. It will only reduce the error but not eliminate it.

To make sure your planets behave, you need to ensure that (at every step of your calculations!) they comply with conservation laws of energy and momentum. Force-adjust the results at every step to make sure your model at least doesn't break the natural laws. Not too much anyway :-)

Thanks a lot for your reply. That was also what I have been thinking according to that the step size will influence the result. Do You know anywhere (documents, books or something) where it is described how to make the best possible model? (I'm a computer programmer and would like the very best model possible. My plan is to start in 2D (with a single planet) and then move to 3D and then up with a complete solar system simulator)

 

[DaveC426913]

... and if you stop at that, you will soon discover that your simulated planets have wandered off, or crashed into the Sun, or did something else equally bizzare.

I've built these several times and they are stable over long periods.

However, what I'm building are not models of solar systems, they're simply generic orbital paths of up to 20 or so stars in a cluster. There's no accuracy to any real system, but they are stable.
[/QUOTE]

 

[DaveC426913]

I would to make it as precise as possible (only 2D for now).

Oh. In your OP, you said "basic simulations", which I interpreted as "rough".

 

[kjensen]

Oh. In your OP, you said "basic simulations", which I interpreted as "rough".

Actually I meant with a single planet orbiting around the sun, and when I get that working then with several planets in 2D and then hopefully ending up with a 3D solar system simulator (a solar system space travel program). I just haven't been able to find books and/or documents about how to make such simulations (also I would like to make a simulation where it is possible to add new objects).

 

[kjensen]

I know that there exist ephemerides but I would like to make it as a simulation as I am very interested in astronomy and would like to understand the mechanisms in the solar system.

 

[kjensen]

To make sure your planets behave, you need to ensure that (at every step of your calculations!) they comply with conservation laws of energy and momentum. Force-adjust the results at every step to make sure your model at least doesn't break the natural laws. Not too much anyway :-)

I already knew from experiments that the step size would affect the curve form. Can you (or someone) say more about how to force-adjust at every step (formulas, code fragments or so)?

Also what do you mean with not breaking the natural laws too much anyway? Is it not possible to make a perfect simulation, and if not then why not? And if it is possible to make a perfect simulation where can I find informations about it? I have tried to search for documents and at Amazon for books about it, but I haven't been able to find any documents or books about high quality computer simulations of planets orbiting the sun!

Also my questions is not based on passing any exams, but is the starting point in making a perfect planet simulation program (I will start in 2D with one planet), and therefore I would like to know everything about the physics related to it. I have a degree in computer science (and also I have taken a few courses in physics), so I know how to make graphics etc., but I need to understand the physics and mathematics of planetary computer simulations a lot better.

 

[DrStupid]

I want to build my programs on Newtons laws.

As DaveC426913 already wrote you have to decide which numerical method you want to use. This may strongly influence the structure of your program. The best choice would be a symplectic multi-step integrator for second order IDEs with variable step width. But there are a lot of simple methods that give sufficient results. I prefer a Runge-Kutta-Nyström method.

Than you should think about the representation of the data. There are several different possibilities such as putting the data into different arrays or encapsulating them into objects. I prefer objects for each body and an array for all bodies. That makes it very easy to remove or add bodies even during a running simulation.

I would not recommend to limit the program to 2D first because it might be very difficult to extend it to 3D later. As you have to use vector algebra in both cases there is actually no reason for such a limitation.

As you want to extend your simulation to a space travel program you have to implement a full simulation of Newtons gravitational law with individual mass for different bodies including the possibility of external forces (e.g. from rocket motors). But you may use a hybrid simulator with Kepler for Sun and planets and Newton for everything else.

Is it not possible to make a perfect simulation, and if not then why not?

It is not possible because there is no global solution for n-body problems with n>2. Thats why such systems can only be approximated numerical and such methods have a limited accuracy.

 

[TurtleMeister]

Hi kjensen. I think http://www.cs.princeton.edu/courses/archive/spring11/cos126/assignments/nbody.html" may be what you're looking for. I used it to write a java simulation which I later ported over to VB with added bells and whistles. I learned a lot and had a lot of fun with it.

 

[kjensen]

As DaveC426913 already wrote you have to decide which numerical method you want to use. This may strongly influence the structure of your program. The best choice would be a symplectic multi-step integrator for second order IDEs with variable step width. But there are a lot of simple methods that give sufficient results. I prefer a Runge-Kutta-Nyström method.

Thanks a lot for your reply. If I google "symplectic multi-step integrator" then I get a lot of interesting stuff (documents, books and so on). It might be what I am looking for! Do you know any books or documents where it is related to computer programming (any programming language). Thanks.

 

[kjensen]

Hi kjensen. I think http://www.cs.princeton.edu/courses/archive/spring11/cos126/assignments/nbody.html" may be what you're looking for. I used it to write a java simulation which I later ported over to VB with added bells and whistles. I learned a lot and had a lot of fun with it.

Thanks a lot for your reply. The link looks really great and even though I will build my program from ground up with a model as precise as possible then it might be worth figuring out the exercise (but it looks rather similar to another programming experiment I have been doing). But I think that the curve form will depend on the step size doesn't it?

 

[tony873004]

If you want to learn to walk before you learn to run, I'd start with Euler's method. It's the method DaveC described in post 3. It will give a very intiutive understanding of what you're doing. Leapfrog and RK4 are fine, but more difficult to code. At this point high accuracy is not an issure. You're creating a single planet moving around a stationary star, just to get a feel as to how to do it. So who cares if your simulation isn't accurate enough to predict solar eclipses decades into the future?

Actually Euler's method can be accurate enough to predict solar elipses decades into the future if you take a small enough time step. (delta t). The higher order integrators allow you to achieve good accuracy at larger timesteps, meaning it will take your computer less time to simulate your system decades into the future.

Later, you can upgrade to a higher order integrator. Since an Euler step is performed in the RK4 method, you'll be building upon what you've already done.

 

[kjensen]

If you want to learn to walk before you learn to run ...

Thanks a lot for your reply. Of course you are right, and I will design my program so that I easily can change the 'calculation engines' along the way. But now when I have started a thread about it then I would like to know what the best possible method is, so that I don't have to start another thread later on.

 

[DrStupid]

Do you know any books or documents where it is related to computer programming (any programming language).

Unfortunately not. Maybe you will find corresponding sources in the references of your google-results. But you should also be able to write your program without literature. The program code will mainly result from the selected integrator and the representation the physical parameters. For example http://www.drstupid.de/Newton.html" is based on this Runge-Kutta-Nytröm integrator:

$\begin{array}{l} a_1 = a\left( {x_n } \right) \cdot \Delta t \\ a_2 = a\left( {x_n + {\textstyle{{\Delta t} \over 2}}v_n + {\textstyle{{\Delta t} \over 8}}a_1 ,v_n + {\textstyle{1 \over 2}}a_1 } \right) \cdot \Delta t \\ a_3 = a\left( {x_n + {\textstyle{{\Delta t} \over 2}}v_n + {\textstyle{{\Delta t} \over 8}}a_1 ,v_n + {\textstyle{1 \over 2}}a_2 } \right) \cdot \Delta t \\ a_4 = a\left( {x_n + \Delta t \cdot v_n + {\textstyle{{\Delta t} \over 2}}a_3 ,v_n + a_3 } \right) \cdot \Delta t \\ x_{n + 1} = x_n + v_n \cdot \Delta t + {\textstyle{{\Delta t} \over 6}}\left( {a_1 + a_2 + a_3 } \right) \\ v_{n + 1} = v_n + {\textstyle{1 \over 6}}\left( {a_1 + 2a_2 + 2a_3 + a_4 } \right) \\ \end{array}$

[http://theory.gsi.de/~vanhees/faq/gravitation/node62.html" ]

This algorithm already tells you something about the corresponding program:
As it is a 4-step algorithm intermediate data must be handled. For every step 4 accelerations (for different positions) must be calculated for every body and kept in memory. I solved this problem by an array of 4 accelerations for each body:

function Body(Mass,sx,sy,sz,vx,vy,vz,Name) {
  this.m = Mass ; this.s = new Vector3D(sx,sy,sz) ; this.v = new Vector3D(vx,vy,vz) ;
  this.F = new Vector3D(0,0,0) ; this.a = new Array ; for (var k=0;k<4;k++) this.a.push(new Vector3D(0,0,0)) ;
  [...]
}

Within each time step the accelerations a[0] - a[3] are calculated using the first four equations of the algorithms. Than the velocities s and v are updated according to the last two equations using the accelerations calculated before. Each step of the algorithm must be calculated for all bodies at once. In doing so the force F is accumulated for each pair of bodies (in order to calculate each force only once).

With other algorithms and other data management you will get different code. If you select a numerical method and your favorite data management you may already have an idea about the structure of your program.

 

[DrStupid]

If you want to learn to walk before you learn to run, I'd start with Euler's method. It's the method DaveC described in post 3. It will give a very intiutive understanding of what you're doing. Leapfrog and RK4 are fine, but more difficult to code.

You are right. But as a compromise I would suggest http://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet". It is almost as simple to understand and to implement as Euler but it conserves the energy.

 

[kjensen]

Thank you DrStupid. Your code looks really great. I will give that approach a try!

 

[D H]

Thanks a lot for your reply. If I google "symplectic multi-step integrator" then I get a lot of interesting stuff (documents, books and so on). It might be what I am looking for! Do you know any books or documents where it is related to computer programming (any programming language). Thanks.

Symplectic multi-step methods are a big step for someone just starting out.

I suggest you start with something simple such as Euler's method or symplectic Euler, but design your system so that the most of the details of the integration technique are hidden from the uses of those techniques. This way you can start with something simple and step up to more advanced techniques as you learn.

Besides, accurate ephemeris models typically do not use symplectic methods. JPL uses a variable step-size, variable-order, Adams method to compute their Development Ephemerides (DE series; the latest of which is the DE421); the Russian Institute of Applied Astronomy of Saint-Petersburg uses a similar technique for their Ephemerides of the Planets and the Moon (EPM), as does the French Observatoire de Paris for their INPOP ephemerides. NORAD uses Stoermer-Cowell (aka Gauss-Jackson) to integrate the states of objects that are orbiting the Earth.There are lots of things to learn about besides integration techniques.

You'll need to learn about how to represent time. Eventually you'll learn about why we have leap seconds. As a starter, I suggest you use TAI (International Atomic Time) as your time standard. If you later want to address relativistic effects you will need to move to a relativistic time scale such as JPL's T_eph or the IERS's Barycentric Dynamic Time. If you want to relate the calculated ephemerides with what you can see from the ground you'll also have to learn about things like Universal Time.

You'll need to learn about how to represent positions in space. There are many choices of reference frames. I suggest the IERS's International Celestial Reference System. If you want to relate the calculated ephemerides with what you can see from the ground you'll also have to learn about things like Terrestrial Reference Frames, and a whole lot about the Earth's rotation.

As a starter, I suggest you use TAI as the dynamic time for your system with coordinates in ICRF, and ignore all of the subtleties. You can add the subtleties as you learn more (and as you need them).

Some key websites:
US Naval Observatory Astronomical Information Center
http://www.usno.navy.mil/USNO/astronomical-applications/astronomical-information-center

Standards of Fundamental Astronomy:
http://www.iausofa.org/

International Earth Rotation and Reference Systems Service:
http://www.iers.orgAs far as integration techniques go, I suggest that you start with something very simple. The two simplest techniques are basic Euler and symplectic Euler. Basic Euler treats position and velocity as a six-vector and integrates according to
$$\vec u(t+\Delta t) = \vec u(t) + \Delta t\frac{d}{dt}(\vec u(t))$$
This is equivalent to propagating position and velocity according to
$$\begin{aligned} \vec x(t+\Delta t) &= \vec x(t) + \Delta t\,\vec v(t) \\ \vec v(t+\Delta t) &= \vec v(t) + \Delta t\,\vec a(t) \end{aligned}$$
Symplectic Euler merely reverses the order of the calculation:
$$\begin{aligned} \vec v(t+\Delta t) &= \vec v(t) + \Delta t\,\vec a(t) \\ \vec x(t+\Delta t) &= \vec x(t) + \Delta t\,\vec v(t+\Delta t) \end{aligned}$$

Basic Euler is incredibly inaccurate and unstable. Symplectic Euler is fairly inaccurate but is stable. The only reason to learn about basic Euler is that it is the basis for a lot of other integration techniques. The Adams family of integrators (which includes Stoermer-Cowell) and the Runge-Kutta family of integrators use basic Euler as the underlying mechanism for integration.

A couple of techniques that are a bit better than the Euler techniques are leapfrog, verlet integration, and Beeman's algorithm. A lot of galactic simulations, simulations in chemistry, and robotic simulations use these techniques or similar ones. Note very well: They don't use these techniques because they are good. They use these techniques because anything more advanced would kill the computer while dropping down to Euler integration would kill accuracy.

Beyond these, you'll find two key families, the Adams and Runge-Kutta families of integrators, that predominate in this field. (There are plenty more techniques than these, but these are the dominant ones). Neither Adams nor RK-type integrators are symplectic. The losses are small, particularly if you go to higher order integrators. One problem with symplectic integrators is that it is very, very hard to get beyond second order. It is fairly easy to come up with Adams and Runge-Kutta integrators that are well beyond second order.

Perturbation techniques open the door to a completely different style of integration. The Apollo program used an Encke-Nystrom integrator to accurately propagate the state of the Apollo spacecraft . An orbiting body will follow something that is fairly close to an elliptical path. An Encke-Nystrom integrator takes advantage of this fact, propagating the path as a small deviation from a simple elliptical path. The integrator has to re-anchor the ellipse to on occasion to keep the integration accurate.

Integrating something other than cartesian position and velocity opens the door to yet another style of integration. Variation of parameters techniques integrate Lagrange's planetary equations (or some variant of those planetary equations). In these techniques, it is orbital elements that are integrated rather than position and velocity. The planetary equations provide the partial derivatives of these elements with respect to the perturbing forces. A somewhat related technique is to use the Karhunen-Loève transform. An elliptical path falls out directly from this transform.

 

[tony873004]

Thanks a lot for your reply. Of course you are right, and I will design my program so that I easily can change the 'calculation engines' along the way. But now when I have started a thread about it then I would like to know what the best possible method is...

The best possible method might be one that you invent :)
I've done exactly what you've proposed: design code so that I can easily change 'calculation engines'.

The code for these integrators is available on my website ( www.gravitysimulator.com ). If you have a Windows computer, you can download it and compare the different routines to see their strengths and weaknesses.

 

[kjensen]

Thank you very much D H. Your answer seems very deep and profound, and it looks as a very good 'recipe' for what I need to make the program. I think that I better make the program in C/C++ as it otherwise will get to slow when I extend it to higher precision! Thanks.

 

[kjensen]

The code for these integrators is available on my website ( www.gravitysimulator.com ). If you have a Windows computer, you can download it and compare the different routines to see their strengths and weaknesses.

Your program looks great and also it runs under Linux with Wine even though i couldn't find the Integrator menu that you have shown on the picture. Also I couldn't find any source code on your website. I personally prefer open source systems and also my program when it is finished one day will be released as open source.

 

[TurtleMeister]

But I think that the curve form will depend on the step size doesn't it?

If you mean accuracy then yes, it will depend on Δt. Smaller Δt will give greater accuracy but it will require more computation for the same animation speed or run time.

 

[tony873004]

Sorry, I forgot to mention that the Integrator menu exists on the latest Beta version which is available here :
http://www.orbitsimulator.com/cgi-bin/yabb/YaBB.pl?num=1176774875/47#47

RK4 source code is here:
http://www.orbitsimulator.com/cgi-bin/yabb/YaBB.pl?num=1176774875/30

Euler-Cromer is here:
http://www.orbitsimulator.com/cgi-bin/yabb/YaBB.pl?num=1160874415/6#6

I can't find Verlet or Leapfrog posted, so I'll find them and post them soon on my website's forum.

Most of Gravity Simulator's math code is posted on the forum. Anybody who does any research with it is going to want to know what it's doing. Anybody who wants to write their own 'calculation engine' can do so, compile it into a .dll file and drop it into Gravity Simulator. So in that way its kind of open source. The vast majority of the unpublished code is just simple math (how to compute a diameter from a radius, etc), and instructions for how to make Graphic User Interfaces (GUIs).

Edit: I just posted Verlet and Leapfrog here: http://www.orbitsimulator.com/cgi-bin/yabb/YaBB.pl?num=1306976185/0

 

[D H]

If you mean accuracy then yes, it will depend on Δt. Smaller Δt will give greater accuracy but it will require more computation for the same animation speed or run time.

This is only true if you are computing with infinite precision, but computing with infinite precision would require an infinite amount of time to make even the smallest of steps.

With finite precision arithmetic (e.g., the single and double precision numbers in Fortran, or the floats and doubles in C/C++), this is not true. One problem is that 1.0 + 1e-16 is exactly 1.0 with double precision numbers. Integrate with way too small a Δt value and x(t)+v(t)Δt is just x(t). Position and velocity simply freeze. Things start moving for Δt values that are above this minimal value, but the motion is just plain lousy. Increase Δt even more and the quality/accuracy of the integration improves -- to an extent.

Now let's look at the opposite extreme of very large Δt values. Here the integration technique itself is at fault. Consider simple Euler, x(t+Δt)=x(t)+v(t)Δt. With a Δt equal to many times the orbital period, that first step can take the object right out of the solar system. Decrease Δt somewhat and the path will start to follow a curve as opposed to a straight line, but the quality of the integration is still terrible. Decrease Δt even more and the quality/accuracy of the integration improves -- to an extent. Eventually you'll get to the point where numerical accuracy starts to degrade accuracy.

The end result is a bathtub curve. Every integration technique follows such a bathtub curve. There is some optimal Δt that gives the best accuracy. This optimal Δt depends on the problem at hand, on the numerical precision of the underlying number representation scheme, and on the integration technique. Adaptive integration techniques try to find this optimal Δt value.

 

[kjensen]

With finite precision arithmetic (e.g., the single and double precision numbers in Fortran, or the floats and doubles in C/C++), this is not true. One problem is that 1.0 + 1e-16 is exactly 1.0 with double precision numbers. Integrate with way too small a Δt value and x(t)+v(t)Δt is just x(t). Position and velocity simply freeze. Things start moving for Δt values that are above this minimal value, but the motion is just plain lousy. Increase Δt even more and the quality/accuracy of the integration improves -- to an extent.

Again thanks a lot for your deep insights. I agree to what you saying, but actually it is possible in for example C/C++ to extent to arbitrary precision as you can see here:

http://en.wikipedia.org/wiki/Arbitrary-precision_arithmetic

My experience is that if you use the newest processors and compile for example Linux to your newest processor (or processor cluster), then you can still achieve good performance in very high precision. Also if you write the calculation engine in assembler using parallel processing techniques then it is possible on new processors to get very good performance.

 

[kjensen]

@tony873004

It is a really cool project you are doing, and your calculation engines looks very interesting. Did you know that with wxWidgets then it is possible rather easily to write a program with a 'killer' GUI that can compile and run very fast on (almost) all operating systems!

You can check it out here:

http://www.wxwidgets.org/

I personally think that wxWidgets is one of the coolest computer projects I have ever seen!

 

[kjensen]

This is only true if you are computing with infinite precision, but computing with infinite precision would require an infinite amount of time to make even the smallest of steps.

I'm thinking about if it is possible to get (almost) infinite precision by adding some sort of factor to the result from the finite calculations?

 

[D H]

Again thanks a lot for your deep insights. I agree to what you saying, but actually it is possible in for example C/C++ to extent to arbitrary precision as you can see here:

http://en.wikipedia.org/wiki/Arbitrary-precision_arithmetic

Sure, you can extend to some greater degree of precision than nominal, but you cannot extend to infinite precision. That would require infinite memory and an infinitely fast processor. For any representation of the reals on a digital computer there will always exist a number δ>0 such that 1+ε=1 for all |ε|<δ.

My experience is that if you use the newest processors and compile for example Linux to your newest processor (or processor cluster), then you can still achieve good performance in very high precision.

My experience is otherwise. Extended precision arithmetic is very expensive computationally, typically 10 times slower (very optimistic) to 1000 times slower (or worse!) than using native floating point. This is an expense that the typical application of orbit propagation programs cannot endure. Besides, there is little need for extended precision arithmetic here. Sophisticated integration techniques can achieve a relative precision of 10^-14 for a long span of time using native doubles; simpler but still good techniques can achieve a relative precision of about 10^-12 (but only for a shorter span of time). At 60 AU, these levels of precision correspond to 9 centimeters and 9 meters respectively.

Very few applications in physics need extended precision arithmetic for the simple reason that most physical measurements aren't good to anywhere close to 16 decimal digits of accuracy.

Also if you write the calculation engine in assembler using parallel processing techniques then it is possible on new processors to get very good performance.

Assembly? Not usually, especially not for sophisticated numerical integration techniques. A good compiler will typically do a better job. Parallel computing? Numerical solution of the N-body gravitational problem as applied to the solar system is a bit tough to parallelize. Parallel algorithms work quite nicely for modeling a bunch of stars where behaviors rather than accuracy is what is important. Those galactic simulations typically use simple techniques such as leapfrog and involve a huge number of interacting bodies. A solar system simulation involves a small number of interacting bodies and accuracy takes on greater importance. Variable step-size, variable-order Adams methods are a bit tough to parallelize. These factors make it much easier to write a highly parallel solver for a galactic simulation than for the solar system.

The people who do go the parallel computing route inevitable do so using native floating point arithmetic. Using extended precision arithmetic would defeat the purpose. A somewhat recent development is to perform those parallel computations on a computer's graphics processor -- using native floating point arithmetic of course.

 

[kjensen]

My experience is otherwise. Extended precision arithmetic is very expensive computationally, typically 10 times slower (very optimistic) to 1000 times slower (or worse!) than using native floating point.

But according to the benchmarks then this GNU library can perform arbitrary precision pretty fast:

http://gmplib.org/

Also I think that a real good assembly programmer can beat any compiler by far or at least make code that performs just as fast, especially if using multiple cores in parallel clusters (but I admit that it will be very tedious to program such algorithms).

Also to get extreme speed then VHDL hardware programmed algorithms can be used (extremely fast performance if asynchronous fpga's is used)!

But I admit that I easily get carried away when it comes to computers, and for my first one planet orbiting the sun in 2D experiment it might not be necessary with multi cores CPU's connected in clusters supported by asynchronous fpga's calculation engines. But sometimes it is nice to dream :-) Maybe later if I extent my program to simulating the whole universe it might become handy. Anyway I appreciate your deep and profound replies and obviously you know a lot more about physics simulations than I do.

 

[DrStupid]

Parallel computing? Numerical solution of the N-body gravitational problem as applied to the solar system is a bit tough to parallelize.

Most of the time is used for the calculation of accelerations and this is not that difficult to parallelize. For an n-body system there are n·(n-1)/2 connections. With k processors just assign n·(n-1)/(2·k) to every processor. Simultaneous memory access can be avoided by mirroring position and mass data k times. It might be also useful to parallelize parts of the numerical integrator but everything else does not need to be optimized.

 

[D H]

But according to the benchmarks then this GNU library can perform arbitrary precision pretty fast:

http://gmplib.org/

Pretty fast compared to other bignum packages, but still very slow compared to using native floating point. The native floating point calculations are done on a separate processor that is specially architected to perform floating point calculations.

Also I think that a real good assembly programmer can beat any compiler by far, especially if using multiple cores in parallel clusters (but I admit that it will be very tedious to program such algorithms).

Assembly programming is a a niche field today, mostly relegated to area of firmware. To quote John Moore (of Moore's law fame), "He who hasn't hacked assembly language as a youth has no heart. He who does as an adult has no brain."

Assembly is pretty much dead in the field of scientific programming, and has been for a couple of decades. It made a very brief comeback when people first started using GPUs for scientific computing, but even that pretty much ended when compiler writers added those targets to their toolsets. Assembly code isn't portable from one machine to another. Finely-tuned assembly isn't even portable across machines with the same underlying machine language; those fine-tunings that make the function run so fast on machine X makes it run slow on machine Y. Lack of portability is one reason people stay away from assembly nowadays. An even bigger problem: Who will maintain the code? Assembly is such a niche skill. The job market for it is rather small. Finding someone who can do math, physics and general purpose computing well is hard enough. Finding some who can do all that and can also write assembly well: It's easier to find one needle randomly buried in one of a very large number of haystacks.