Computational Physics - "Numerical Integration": RK4, Verlet.. "n-body problem"

[abaraba]

hello everyone,


this is about an algorithm,
time stepping function aka "iterative solver" - it is the one and the same algorithm no mater if you're simulating planets, moons, galaxies.. atoms, molecules, ping-pong balls or maybe even 'strange attractors' and chaos itself..




PLEASE,
can someone confirm or prove this wrong, im happy with either - THANK YOU!



summary:
- the 1st basic problem is that Wikipedia and the rest of the internet claim that there is some ERROR in "Euler method", but there is not. the error exist only with nonuniform rate of change in acceleration - this caused practical implementation of algorithm to perform "integration" even if there was no need for it, which in turn caused further multiple mistakes and errors in implementation

- the 2nd basic problem has to do with limiting of number of steps per main loop iteration, but this is not really important for non real-time applications, tho its one of the factors in producing the error under certain situation coupled with 1st problem


BOTH of these "bugs" are present in every single software package i came across, that include all of the commercial and Open Source physics libraries - Havok, PhysX, Bullet, ODE, Ogre, Newton.. but it also seem to include scientific simulators; "Software for molecular mechanics modeling": AGM, AMBER, Ascalaph, BOSS, CHARMM, GROMACS, LAMMPS, MCCCS, MOE, MOIL, MOLDY, NAB, QMOL, RasMol, TINKER, VMD + NAMD, YASARA, Zodiac... (from Wikipedia)


and who knows where else..
i hope, for example, at NASA they actually dont have this buggy algorithm, eh?



ok, let me now copy/paste the message that i hope explains all the subjects that i would like to discuss here, and of course i can go into more details depending on what turns out to be most interesting to public here...



+++COPY/PASTE:
//=== OVERTURE =======================================

Numerical Integration: RK4, Verlet.. WIKIPEDIA is WRONG!

KEYWORDS:
Constraint Algorithm, Iterative Solver, Interpolation, Extrapolation, Verlet, Runge Kutta, Taylor, Sampling, Multistep; Monte Carlo, SHAKE, RATTLE, SETTLE algorithms


do the search and you may find that:
- Wikipedia and every other info on the Internet regarding these matters is simply FALSE in the context of computer simulation, but no one seem to listen to me and im not allowed to make changes.. as 'unverified source'

So, there you go.. all these studies, books, papers, articles.. WRONG?!

//==================================================

//--- FUGUE..

to whom it may concern...
+++ i believe i found the heart of the whole problem:


"numerical integration",
as i suggested before, is completely misplaced in its practical implementation



--- the cause ---

before applying mathematics to physics, and before implementing the solution in computer simulation application there was a failure to recognize TWO very different types of "simulation", i'll call them 'FIELD FORCES' and 'CLASSIC' simulation



A.) 'FIELD FORCES simulation', i.e. NON-UNIFORM "external" ACCELERATION
- deals with oscillations, gravity fields, planet or atom orbits..

B.) 'CLASSIC simulation', i.e CONSTANT "external" ACCELERATION
- everything that is not in the first category.. majority of simulations that deal with objects under constant gravity (constant altitude)



--- the problem ---

1.) with constant acceleration - basic kinematic equations will actually give EXACT solution, i.e. the next POSITION is exact solution as function of time, previous position and acceleration/velocity
(what was measured as error was actually "bug" in implementation)

ergo,
we do not need any APPROXIMATION... no INTERPOLATION is necessary with constant acceleration - we do not need any fancy mathematics, we already have EXACT solution, and so addition of any other computation will only hinder precision and speed


2.) ..plus, if you have this popular implementation that uses MAX steps, as everyone does, then i hope everyone can see now where all the problems came from



--- the solution ---

huh.. this thing surely turned much deeper than i expected, and im still tumbling down the rabbit hole.. but i know this - there are many problems here and solutions that i already proposed do work.. tho, the general and complete solution has to do with defining the problems of the "FIELD FORCE simulation" and basically just define stuff more clearly - which is what im trying to do, and i would appreciate some help







--- the general and complete solution ---

"FIELD FORCE simulation" - NON-UNIFORM ACCELERATION
- so you doing simulation of planetary or maybe atomic orbits..


then, i think its important to realize this sort of simulation is to CLASSIC simulation what ray tracing is to vector graphics


basically, all of the objects are under COLLISION all the time, they all interacting with each and every other object 'ALL the time', in NON-UNIFORM acceleration topology


now we have two problems:

1.) physics libraries and all the spatial partitioning and speed-ups are of no use, the whole "constant collision" is a huge performance problem in its own and solution is to buy faster computer or build cluster of CPUs


2.) the solution to "constant collision" algorithm has actually to do with "n-body problem", "Kepler problem", maybe even "Coulomb's law" and stuff like that - basically, a problem is to find a position as function of time when given initial positions, masses and velocities


one of the ways to go about it is by ITERATIONS, i.e every time step we calculate new AVERAGE ACCELERATION, then calculate velocities and positions as usual




and so,
we finally arrive to the point where "numerical interpolation" makes sense - its nothing else but one of the ways too approximate that *AVERAGE ACCELERATION*, thats the whole point - to compensate for the rate of change in acceleration and try to estimate it given the elapsed time period (time step)


obviously, with enough small time step and maybe some "simple average" the error caused by this would probably be minimal compared to other numerical instabilities already present in computer hardware... and so, even in this case RK4, Verlet.. might be less than desirable solution


*ONLY* "EXTERNAL" NON-UNIFORM ACCELERATION "needs" some sort of APPROXIMATION, there is no need to interpolate position or velocity



of course,
exact solution is much better than some approximation... (to be concluded)



am i right or wrong?

//------------------------------------------------------------
+++END COPY/PASTE





thank you all

[cjameshuff]

summary:
- the 1st basic problem is that Wikipedia and the rest of the internet claim that there is some ERROR in "Euler method", but there is not. the error exist only with nonuniform rate of change in acceleration - this caused practical implementation of algorithm to perform "integration" even if there was no need for it, which in turn caused further multiple mistakes and errors in implementation

You contradict yourself. Yes, there is an error in the result of Euler integration when acceleration is not uniform...so the various sources out there are correct.

- the 2nd basic problem has to do with limiting of number of steps per main loop iteration, but this is not really important for non real-time applications, tho its one of the factors in producing the error under certain situation coupled with 1st problem

I'm not sure what you refer to here. However, computer arithmetic is done to finite precision, and it is not very useful to continue reducing step size when the limiting factor is the numeric precision.

BOTH of these "bugs" are present in every single software package i came across, that include all of the commercial and Open Source physics libraries - Havok, PhysX, Bullet, ODE, Ogre, Newton.. but it also seem to include scientific simulators; "Software for molecular mechanics modeling": AGM, AMBER, Ascalaph, BOSS, CHARMM, GROMACS, LAMMPS, MCCCS, MOE, MOIL, MOLDY, NAB, QMOL, RasMol, TINKER, VMD + NAMD, YASARA, Zodiac... (from Wikipedia)

and who knows where else..
i hope, for example, at NASA they actually dont have this buggy algorithm, eh?

Carefully and precisely describe exactly what you think the "bugs" are.

A.) 'FIELD FORCES simulation', i.e. NON-UNIFORM "external" ACCELERATION
- deals with oscillations, gravity fields, planet or atom orbits..

B.) 'CLASSIC simulation', i.e CONSTANT "external" ACCELERATION
- everything that is not in the first category.. majority of simulations that deal with objects under constant gravity (constant altitude)

The usual division is between numeric simulations and analytical ones. Many of what you call "field forces simulations" can be done analytically.

1.) with constant acceleration - basic kinematic equations will actually give EXACT solution, i.e. the next POSITION is exact solution as function of time, previous position and acceleration/velocity
(what was measured as error was actually "bug" in implementation)

ergo,
we do not need any APPROXIMATION... no INTERPOLATION is necessary with constant acceleration - we do not need any fancy mathematics, we already have EXACT solution, and so addition of any other computation will only hinder precision and speed

This does not follow. If acceleration is constant, d = v*t + 0.5*a*t^2 is indeed ideal, giving exact results to the precision of the arithmetic...but acceleration is usually not constant, or even varying in a way conducive to analytical solution.

then, i think its important to realize this sort of simulation is to CLASSIC simulation what ray tracing is to vector graphics

I think you mean scanline rendering. And a better analogy would be a raster image to a vector image.

1.) physics libraries and all the spatial partitioning and speed-ups are of no use, the whole "constant collision" is a huge performance problem in its own and solution is to buy faster computer or build cluster of CPUs

This is not true. At a distance, the effects of the individual objects in a small volume of space are often insignificant, and so they can be grouped together and treated as a single object as far as interactions with distant objects are concerned. The maximum error produced by this simplification can be quantified. It is a valuable technique for handling complex systems without requiring excessive amounts of processing power.

*ONLY* "EXTERNAL" NON-UNIFORM ACCELERATION "needs" some sort of APPROXIMATION, there is no need to interpolate position or velocity

Pretty much. What's your point? Analytical solutions are almost always preferred over numeric ones, when they are possible. There are a few cases where exact analytical solutions are just too unwieldy for practical use, or the precision isn't worth the trouble of working the equations out, but as a general rule, they are preferred. Those physics packages will generally not do analytical solutions, the assumption is that if you need to use one, you are doing something that is complex enough to need numeric solving. If "d = v*t + 0.5*a*t^2" is sufficient for your calculations, why would you be using a complicated physics library?

So...what exactly is the problem you see? Or the question you have?

[abaraba]

cjameshuff,


are you proving me wrong?
..or just saying "stuff" with no reason, purpose or meaning?


if you go on pretending that you know - you will never learn




>>"You contradict yourself."
>>"...so the various sources out there are correct."

- no, i do not, but you do..
please make yourself familiar with "various sources"





>>"I'm not sure what you refer to here."

- of course you're not,
and yet you find appropriate to be giving advice?





>>"Carefully and precisely describe exactly what you think the "bugs" are."


- i do not think, assume or guess - i KNOW


the 2nd bug is mostly related to real-time simulation, it has to do with smooth animation and FPS, but since i had the most terrible time trying to inform public of this, i think that i should keep this "bug-fix" for myself ..you dont even know that you have bugs, so its all the same to you





>>"The usual division is between numeric simulations and analytical ones. Many of what you call "field forces simulations" can be done analytically."


- you should be asking questions and not pretending to be giving advice, its silly





>>"This does not follow. If acceleration is constant, d = v*t + 0.5*a*t^2 is indeed ideal, giving exact results to the precision of the arithmetic...but acceleration is usually not constant, or even varying in a way conducive to analytical solution."


- of course it follows,
thats why i first made that classification and so i was clearly talking about:
>"1.) with constant acceleration - ...with constant acceleration"


so it is you who does not follow,

you say:
>>"but acceleration is usually not constant" - even if my first three words clearly stated what case im talking about, then i repeated it once more, just in case, but that wasn't enough apparently - whats your problem?





>>"I think you mean scanline rendering. And a better analogy would be a raster image to a vector image."


- stop insulting yourself,
you clearly know nothing about any of these then





>>"This is not true. At a distance, the effects of the individual objects.."

- you can simulate your objects as you like,
i said i was describing GENERAL solution, which also means unoptimized, or yet in other words, it also means 'as precise as possible'

..you can optimize as you like, but dont complain and go write studies about some "phantom errors", mysterious "energy leaks" and other "phenomenon" ..its just bugs

good luck!




>>"Pretty much. What's your point?"

- my points is that everyone has wrong implementation, that Wikipedia and the rest of the internet have wrong explanations which most likely was actual cause that everyone have same wrong implementation


to understand that you need to be looking at source code of all these software packages and libraries i mentioned or whichever software you're using, then you might see bugs im talking about ..and then you might search internet and you will find out where the bugs came from




>>"So...what exactly is the problem you see? Or the question you have? "

- are you blind?
--------------------------------------------------------------------------
PLEASE,
can someone confirm or prove this wrong, im happy with either - THANK YOU!
--------------------------------------------------------------------------


..for that you need to use ARGUMENTS and you need to know what are you talking about



actually when i look at it now,
when you're not plain wrong or confused, you mostly DO agree with me.. thanks

[cjameshuff]

If this is the tone you've been using elsewhere, it's no wonder people don't pay any attention to you.

are you proving me wrong?

You haven't really given enough of an argument to prove wrong. There were some factual errors, I covered most of those in my previous post.

>>"I'm not sure what you refer to here."

- of course you're not,
and yet you find appropriate to be giving advice?

I don't know what you think is a problem, because you haven't stated it clearly. I do appear to be equipped to give you advice. Have you ever even taken a numeric analysis class? Or written an orbit simulator?

>>"Carefully and precisely describe exactly what you think the "bugs" are."


- i do not think, assume or guess - i KNOW


the 2nd bug is mostly related to real-time simulation, it has to do with smooth animation and FPS, but since i had the most terrible time trying to inform public of this, i think that i should keep this "bug-fix" for myself ..you dont even know that you have bugs, so its all the same to you

This is not a careful and precise description of what you think the problem is. Have you considered that if you want people to be informed about it...you might have to give them some information?

>>"The usual division is between numeric simulations and analytical ones. Many of what you call "field forces simulations" can be done analytically."

- you should be asking questions and not pretending to be giving advice, its silly

I did ask a question. And here's another: can you point out anything wrong with what I just said?

- of course it follows,
thats why i first made that classification and so i was clearly talking about:
>"1.) with constant acceleration - ...with constant acceleration"


so it is you who does not follow,

you say:
>>"but acceleration is usually not constant" - even if my first three words clearly stated what case im talking about, then i repeated it once more, just in case, but that wasn't enough apparently - whats your problem?

Okay...so you're talking about problems involving constant acceleration. Why claim a problem with orbit simulations, then? Why talk about numeric integration packages, when you're talking about something with a closed form analytical solution?

>>"I think you mean scanline rendering. And a better analogy would be a raster image to a vector image."

- stop insulting yourself,
you clearly know nothing about any of these then

I've written raytracers, scanline renderers, and programs using raster and vector graphics. I've contributed a fair amount of code to the POV-Ray raytracer.
Raytracing is to vector images as apples are to oranges. One's a way of representing 2D images with polygons and curves, the other is a way of generating a 2D image of a 3D scene by performing intersection tests on rays going into the scene. Scanline rendering instead projects polygons onto a 2D plane, and raster images represent an image with an array of pixels.

>>"This is not true. At a distance, the effects of the individual objects.."

- you can simulate your objects as you like,
i said i was describing GENERAL solution, which also means unoptimized, or yet in other words, it also means 'as precise as possible'

Then why'd you bring it up?

>>"Pretty much. What's your point?"

- my points is that everyone has wrong implementation, that Wikipedia and the rest of the internet have wrong explanations which most likely was actual cause that everyone have same wrong implementation

In what way is it wrong?

>>"So...what exactly is the problem you see? Or the question you have? "

- are you blind?

That would have made it rather difficult to test the renderers and simulations I've written...
Above, you got rather upset about me "giving advice" rather than asking questions. In case you didn't recognize it, that was a question. Please answer it.

..for that you need to use ARGUMENTS and you need to know what are you talking about

If you want people to pay attention to what you're saying, you might try that yourself...

actually when i look at it now,
when you're not plain wrong or confused, you mostly DO agree with me.. thanks

I'm certainly confused, your rant has not exactly been the most coherent of arguments. As for the other...where am I wrong?

[Tinaa]

abaraba take time to read http://www.bautforum.com/about-baut/...ing-board.html

Pay especial attention to 2. Civility and Decorum

Politeness is the top rule here. Of course, we expect to have spirited debates! That’s fine, as long as the people involved extend one another basic respect. Disagreements are inevitable, but even in those situations you must still be nice.

Attack the ideas, not the person(s) presenting them. If you've got concerns with what someone is saying, feel free dismantle their arguments, but do not resort to ad hominem or personal attacks. Be mindful and respectful of others' feelings. If you feel that someone has crossed the line and insulted you, please contact one of the moderators, preferably via the reporting mechanism described here, or by PM or email. Don't write scathing posts in the forum to try and humiliate people publicly.

If these guidelines are not followed, the administrators/moderators will take swift and appropriate action, so please behave accordingly.

[abaraba]

cjameshuff,

- please copy/paste your time stepping function so i can show you on concrete and practical example..

would that be enough of a proof?






//--------------------------------------------------------

Tinaa,


>>"Politeness is the top rule here. Of course, we expect to have spirited debates!"

- why do yo saying this to me?

is it not obvious that it was 'cjameshuff' who was pointless, trying to put me down with no reason and with no arguments, skewing my sentences and lecturing about stuff he clearly showed he has no knowledge of



>>"Attack the ideas, not the person(s) presenting them."

- clearly it was me who was attacked here first, with no reason and NO ARGUMENTS, i simply do not back down and no one will **** over me, i hope that is clear



>>"If these guidelines are not followed, the administrators/moderators will take swift and appropriate action, so please behave accordingly. "


- i fixed many problems that concern almost everyone... there is no reason, what so ever, for anyone to try and put me down


i dont need silly lectures and i can do whatever i please.. you, on the other hand, have these choices:

- you can say "thank you"
- you can ask anything
- you can agree
- you can disagree and we can discuss it until one is proven wrong

but what ever you say, you MUST QUOTE what are you referring to and use ARGUMENTS to support your opinion, NOT GUESSES, NOT ASSUMPTIONS, so we can make it clear for everyone and continue the discussion, everything else is rude, unclear or unnecessary




//----------------------------------------------------------
-"All truth goes through three stages.
First it is ridiculed;
then it is violently opposed;
finally it is accepted as self-evident." (Schopenhauer)

[ToSeek]

i dont need silly lectures and i can do whatever i please.. you, on the other hand, have these choices:

- you can say "thank you"
- you can ask anything
- you can agree
- you can disagree and we can discuss it until one is proven wrong

There is one other possibility here that you need to take into consideration: we can suspend or terminate your posting privileges here. And this is exactly what's going to happen unless you start treating the rules and other members of this forum with respect and consideration.

I see nothing objectionable in cjameshuff's posts so far and a great deal objectionable in yours. People are trying to help you out and understand what you're talking about here - you should take that into account rather than treating any criticism as hostility.

ToSeek
BAUT Forum Moderator

[mugaliens]

...but what ever you say, you MUST QUOTE what are you referring to and use ARGUMENTS to support your opinion, NOT GUESSES, NOT ASSUMPTIONS, so we can make it clear for everyone and continue the discussion, everything else is rude, unclear or unnecessary.

Have you considered pursuing your explorations via a graduate level course in philosophy, where such mandatory rules apply?

Here on these message forums, collectively known as Bad Astronomy/Universe Today, we take more of an eclectic approach, and quite often, guesses, assumptions, and other "soft" approaches (some might call them "touchy-feely") are better suited for the many and varied individuals who call this online community "home."

You rules sound reasonable, however, for such a discussion as you propose. Have you considered starting your own message forum website, wherein you would be free to impose your own rules as you see fit?

I do not say this in jest, but as a mod and admin on other boards, a few of which I've started from scratch. Sometimes, which certain approaches are not feasible on one board, it's time to start another, with a different set of rules under which they would be feasible.

On the other hand, you can accomplish much the same here by simply choosing to adhere to those rules yourself while mentioning that you're doing, without requiring that others do so. Then, you're free to respond to whatever portion of anyone's comment you wish, and if you find someone's response not to your liking, you're also free to simply ignore it.

I hope you find this bit of constructive commentary useful.

Good luck!

[abaraba]

i apologize to everyone..

i suppose its all matter of opinion, i thought this concerns everyone

Computer Simulation,
used in Chemistry, Biology, Medicine.. to develop new drugs, to study diseases, in genetic engineering.. its used to help us predict earthquake, tsunamis and other natural disasters.. in industry to build automobiles and planes.. its used to help us understand the world around us and model interactions of atomic and subatomic particles.. its even used in development and testing of nuclear weapons.. to model planets, stars, galaxies, used by NASA for planning trajectories, orbits, exit and entry points ...it used everywhere


i thought it was important enough, so i felt responsibility to share this with the rest of the world.. in return, i got only stones in my face and fingers in my eyes ..interesting

[Veeger]

This particular discussion is very interesting to me, but it definitely got off on the wrong foot, and I think mainly because of the opening post. Perhaps because of the posting style, or whatever, I found it difficult to ascertain what your point was, abaraba. For example, you stated 1) there was no inherent error in the Euler method for uniform acceleration 2) there is no need to limit step size on non realtime computations; 3) both of these bugs are present on all software packages.

(If I incorrectly summarized your premises, forgive me)

What bugs?

I think cjameshuff was right on track when he stated that (more or less - I'm paraphrasing) that the biggest limitation of digital computation models is the round-off error inherent in the computation system. After many interations, the error accumulates and soon the error exceeds the threshold of what is acceptable.

I presently have several implementations of integrators (in source code format) and have been looking into them more deeply and while each approach seems to be according to the whim and skill of the programmer, at the end of the day, if the simulation yields acceptable results then it works. (disclaimer - I am not a mathematician, but I have written a few programs in my day ).

[cjameshuff]

What bugs?

I'm still wondering that myself. It'd be nice if he answered some of the questions and gave us some idea what he was talking about, rather than howling about unjust persecution in response to every reply...

I think cjameshuff was right on track when he stated that (more or less - I'm paraphrasing) that the biggest limitation of digital computation models is the round-off error inherent in the computation system. After many interations, the error accumulates and soon the error exceeds the threshold of what is acceptable.

To be more precise, I was pointing out that for a given solver, problem, and numeric precision, there is a point where reducing the size of the time-steps can not improve results because the numeric precision of the computations becomes the limiting factor. One of his "bugs" is apparently a practice of limiting the number of steps you can specify...I've not come across such a limit myself, and would assume that NASA, etc. are doing the analysis to determine what step size meets their requirements.

The other "bug", as near as I can tell, is that numeric solvers will solve problems numerically even if it is possible to do them analytically. I would put forth that if you don't want a numeric solution, you shouldn't use a numeric solver. In molecular dynamics, orbital trajectories, etc, you generally don't have any other options, though.

abaraba : I have written simulations using RK4, velocity Verlet, and Beeman's method, and a few variations of Euler's method. However, none of them are really conducive to discussion of the algorithms involved...too much unrelated code intertwined. There are many minimal examples of these algorithms out there...why not pick one, and demonstrate the "bugs" and your fix?

[frankuitaalst]

In the case of integration of the n-body gravitational problem there also exists code which uses algorithms based on Taylor-series . This approach is quite different of the methods used by Euler , Leap Frog , RK4...
The only limitation in this approach seems to be the numerical limitation
The method works quite well , fast and accurate , especially for low numbers of bodies .
This method approaches the analytical solution of the n-body problem , as a Taylor series does .

[Veeger]

What initially struck my interest as a programmer, was the fact I was looking at Aldo Vitagliano's early implementation of a Stormer 8th order integration yesterday and then boom, this thread appears. I have never successfully written a program which uses first order forward Euler that remains stable. Perhaps the algorithm is flawed, perhaps the initial conditions, I don't know, but I have just started exploring this interesting topic. I, for one, am ready to learn something.

[abaraba]

everyone,

feel free to copy/paste your algorithms, your time stepping functions..

i'll show all the problems and then fix them on concrete and practical examples, then you can report your findings and compare with previous results, you also might compare the speed of algorithms, not only precision/correctness


..hope that should be convincing enough

[Van Rijn]

everyone,

feel free to copy/paste your algorithms, your time stepping functions..

i'll show all the problems and then fix them on concrete and practical examples, then you can report your findings and compare with previous results, you also might compare the speed of algorithms, not only precision/correctness


..hope that should be convincing enough

Why don't you show us an example you think is relevant (but not your own code, not something you came up with yourself) and show us the problems you think you've found in it?

[Veeger]

Just to see where this leads, the snippet of code below was from one my first attempts at testing the feasibility of a simple algorithm. The code is written in VB6 which I usually use for prototyping before converting to C/C++.

I am posting the parts which calculate accelerations and which would normally integrate them. The loops shown below are intended to be inside an outer loop which is basically the step loop.

For i = 0 To nObjects - 1
    pi(0) = Obj(i).px               ' position of the i'th body
    pi(1) = Obj(i).py
    pi(2) = Obj(i).pz
    tax = 0                         ' start with no accelerations
    tay = 0
    taz = 0
    For j = 0 To nObjects - 1
        DoEvents                    ' yield to other processes while looping
        ' i cannot equal j
        If i <> j Then
            pj(0) = Obj(j).px       ' position of the j'th body
            pj(1) = Obj(j).py
            pj(2) = Obj(j).pz
            ' Find the acceleration caused by
            ' the gravity of j on i
            dx = pj(0) - pi(0)      ' derive a vector from j to i
            dy = pj(1) - pi(1)
            dz = pj(2) - pi(2)
            dist = Sqr(dx * dx + dy * dy + dz * dz)
            tax = tax + tstep * (G * Obj(j).m / (dist ^ 3)) * dx
            tay = tay + tstep * (G * Obj(j).m / (dist ^ 3)) * dy
            taz = taz + tstep * (G * Obj(j).m / (dist ^ 3)) * dz
        
        End If
    Next


    Obj(i).ax = Obj(i).ax + tax         ' accumulate the accelerations
    Obj(i).ay = Obj(i).ay + tay
    Obj(i).az = Obj(i).az + taz
    
    Obj(i).px = Obj(i).px + Obj(i).ax   ' update the object positions
    Obj(i).py = Obj(i).py + Obj(i).ay
    Obj(i).pz = Obj(i).pz + Obj(i).az
Next

When putting in the masses, starting positions, and accelerations of the Mercury through Neptune and the earth's moon, this code worked and produced decent orbits. The earth orbit, however, soon decayed, and the moon broke free into a solar orbit; obviously flawed. I tried several variants of this code but none allowed me to use a time slice of any value other than one.

[cjameshuff]

In the case of integration of the n-body gravitational problem there also exists code which uses algorithms based on Taylor-series . This approach is quite different of the methods used by Euler , Leap Frog , RK4...
The only limitation in this approach seems to be the numerical limitation
The method works quite well , fast and accurate , especially for low numbers of bodies .
This method approaches the analytical solution of the n-body problem , as a Taylor series does .

Even Euler's method is closely related to Taylor series. For example:

vnew = vold + a*t

That's your simple, first-order Euler's method, and also the first two terms of a Taylor series. It'll give exact results, aside from roundoff error, as long as a is constant. If it is not constant, the error is bounded by the third term of the Taylor series, a'*t^2/2. The position computation is in this situation, and using Euler's method with the initial or final velocity will give very poor results, but if you instead average the velocities at the beginning and end of the step:

pnew = pold + (vold + vnew)/2*t
     = pold + (vold + vold + a*t)/2*t
     = pold + vold*t + a*t^2/2

Which should look familiar...both as the exact solution for motion under constant acceleration, and as the first three terms of a Taylor series. This can be applied further if more derivatives are available, but that can become rather unwieldy. If you infer such derivatives from past accelerations, velocities, or positions, you end up with something that can be shown to be analogous to one of the higher order methods. To my knowledge, *all* such techniques can be analyzed in terms of Taylor series (after all, they are all solving the same problem), and some of them, like Runge-Kutta, are easily extensible to higher orders...RK4 is simply the 4th order version.

[abaraba]

hey cjameshuff, thank you!

you actually did it,
you confirmed the basic premise from which most of the problems spawned

let me take back everything i said ..can i do that?

thanks again..

i'll soon post some practical examples of source code i have here

but let me tell you there are TWO main bugs and there are about 4-5 other "bugs" or common mistakes that go closely related to this - so, although very simple in nature, this may appear complicated just because it interconnects somewhat unrelated stuff...

[edit]
ok, so the code below is ACTUAL code from Open Source Bullet physics library... as much as i know so far, i can tell you the similar/same implementation exist in Intel Havok, NVIDIA PhysX, ODE, Newton...

i've also seen it in molecular dynamics software that is used by universities and large scientific institutions
actually, as far as i know - EVERYONE has some similar implementation that has 1, 2, 3.. or more of these bugs

in this particular case i will concentrate on 2nd bug - "MAX SUBSTEPS BUG"

this is original:

//-----------------------------------------------------------------------
int	btDiscreteDynamicsWorld::stepSimulation( btScalar timeStep,int maxSubSteps, btScalar fixedTimeStep)
{
	startProfiling(timeStep);

	BT_PROFILE("stepSimulation");

	int numSimulationSubSteps = 0;

	if (maxSubSteps)
	{
		//fixed timestep with interpolation
		m_localTime += timeStep;
		if (m_localTime >= fixedTimeStep)
		{
			numSimulationSubSteps = int( m_localTime / fixedTimeStep);
			m_localTime -= numSimulationSubSteps * fixedTimeStep;
		}
	} else
	{
		//variable timestep
		fixedTimeStep = timeStep;
		m_localTime = timeStep;
		if (btFuzzyZero(timeStep))
		{
			numSimulationSubSteps = 0;
			maxSubSteps = 0;
		} else
		{
			numSimulationSubSteps = 1;
			maxSubSteps = 1;
		}
	}

	//process some debugging flags
	if (getDebugDrawer())
	{
		gDisableDeactivation = (getDebugDrawer()->getDebugMode() & btIDebugDraw::DBG_NoDeactivation) != 0;
	}
	if (numSimulationSubSteps)
	{

		saveKinematicState(fixedTimeStep);

		applyGravity();

		//clamp the number of substeps, to prevent simulation grinding spiralling down to a halt
		int clampedSimulationSteps = (numSimulationSubSteps > maxSubSteps)? maxSubSteps : numSimulationSubSteps;

		for (int i=0;i<clampedSimulationSteps;i++)
		{
			internalSingleStepSimulation(fixedTimeStep);
			synchronizeMotionStates();
		}

	} 

	synchronizeMotionStates();

	clearForces();

#ifndef BT_NO_PROFILE
	CProfileManager::Increment_Frame_Counter();
#endif //BT_NO_PROFILE
	
	return numSimulationSubSteps;
}
//-----------------------------------------------------------------------

this completely replaces that thing from above, and fixes 1st & 2nd bug (in most cases):

//-----------------------------------------------------------------------
int btDiscreteDynamicsWorld::stepSimulation(
btScalar timeStep,int maxSubSteps, btScalar fixedTimeStep)
{
    m_localTime+= (timeStep < fixedTimeStep*2)? timeStep : fixedTimeStep*2;
    saveKinematicState(fixedTimeStep); applyGravity(); maxSubSteps= 0;
    while( m_localTime >= fixedTimeStep )
    {
       internalSingleStepSimulation(fixedTimeStep);
       m_localTime-= fixedTimeStep; maxSubSteps++;
    }
    synchronizeMotionStates(); clearForces();
  return maxSubSteps;
}
//-----------------------------------------------------------------------

this is actual, practical example implementation and you can try it in a matter of minutes: http://www.bulletphysics.com/Bullet/phpBB3/


+++symbolic explanation:

--- PROBLEM, "NUMERICAL INTEGRATION" ---

//-------------------------- *** ORIGINAL algo - capping the NUMBER
cnt= 0;
time= 10;

maxSteps= 5;

while(time > 0 && cnt < maxSteps)
{
step(); cnt++;
time--;
}
//--------------------------
Iterations, cnt= 5
"Time left", time= 5 <<<--- TIME REMAINDER

combination of bug1 & bug2 (integration + maxSubSteps)
"integration" wrongly trying to deal with this reminder of time caused by MAX sybSteps- in all of the implementations of this algorithm i've seen.. all physics libraries, ALL ..and some molecular dynamics software.. SAME!

this mostly has to do with smooth animation and FPS under loaded CPU with real-time applications, but if implemented in scientific software, scientist could, in effect, observe some "energy drift" not knowing there was a bug inside the software

this bug would heavily depend on the definition of "MAX ITERATIONS", "MAX NUM of SUBSTEPS" or similarly named variable and the size of time-step (fixed time-step, deltaTime)

[pzkpfw]

Could you edit your posts using [code]...code...[/code] tags to make them readable?

For one thing, it stops HTML-style space compression, so you can use indenting.

[frankuitaalst]

[QUOTE=cjameshuff;1346589]Even Euler's method is closely related to Taylor series. For example:

vnew = vold + a*t

That's your simple, first-order Euler's method, and also the first two terms of a Taylor series. It'll give exact results, aside from roundoff error, as long as a is constant. If it is not constant, the error is bounded by the third term of the Taylor series, a'*t^2/2. The position computation is in this situation, and using Euler's method with the initial or final velocity will give very poor results, but if you instead average the velocities at the beginning and end of the step:

pnew = pold + (vold + vnew)/2*t
     = pold + (vold + vold + a*t)/2*t
     = pold + vold*t + a*t^2/2

QUOTE]
Yes indeed you're right . This is some kind of Taylor expansion .
But the method a referred to writes down the velocities and positions of each body as a time power series up to ...fi. power 22 .
So ri(t) is represented by ri1+ri1*t+ri2*t^2.....rin*t^n
( same for the velocities ) .
Together with the accelerations formulas one obtains a set of algebraic equations which can be solved for ri0, ri1....rin .

The great advantage of this approach is that the motion can be described as a function of time and can be interpolated , and ...the desired accurancy per timestep can be given .
Accurancies of fi. 10m or even less can be obtained over an interval of time of 120 hours ( which is in this case the timestep) .

If interested I can provide some references about this method .
The method , described by Parker et al . is used fi. by Carles Simo

[abaraba]

let me argue these two:

"NUMERICAL INTEGRATION": Interpolation, Extrapolation, Verlet, Runge Kutta, Taylor, Modified Euler...


A.) "NUMERICAL INTEGRATION" is completely misplaced in its practical implementation, time stepping algorithm

B.) "NUMERICAL INTEGRATION" is completely NOT NECESSARY in this algorithm, and as such, it is actually a cause of inaccuracy and slowdown



"Euler method" is all you need.. but maybe i can argue 3rd point as well

C.) Not even "Euler method" is required, all you need is basic kinematic equations and you solve the problem in easy three steps like you did in high school:


STEP 1: calculate INSTANTANEOUS acceleration (this is the ONLY topic for discussion, i think)


then calculate position at the end of the step as usual...

STEP 2: v=a*dt
STEP 3: s=v*dt

--->REPEAT


i dont see why would that be called "Euler method" and why would there be any kind of "numerical integration"

sure, since it is iterative loop it is some sort of "integration" in its nature, but calling it fancy names while it is a matter of simple kinematic equations is misleading, yes?



//--------------------------------
...and once we agree,
we have a chance of making a really nice, general and universal algorithm... having in mind that in case where acceleration was uniform during the time-step the INSTANTANEOUS acceleration is nothing more but simple AVERAGE

so, the REAL QUESTION is - what is the "AVERAGE" of INVERSE SQUARE change of rate in acceleration (jolt, surge)



--- THE TRICK ---
because you not quite sure what is the acceleration at the end of the step - you are not quite sure about the "average" ..and even worse you are not quite sure even about acceleration at the begging of the time-step, since you calculated it with this same uncertainty the last time as well

..kind of like Schrödinger's cat and thing about determining the position of an electron in atom orbit, you can not know velocity, mass and position in the same time.. interestingly, in analogy we could achieve COMPLETE SPATIAL ACCURACY, by giving up the certainty of TIME, but sure without knowing "what time it is" during the simulation, we cant quite iterate to the next step and we would not know what time it is "around" other objects ..ah damn quantum mechanics!

all i wanted is to make the bloody Moon spin around the Earth ..sheesh! can human not do anything on this planet without bumping into some PARADOX?!




let me repeat this now,
- this is COMPLETE and EXACT solution with constant gravity (constant altitude), simply because we KNOW the exact acceleration at ANY time and so we can calculate that average acceleration easy: (a1+a2)/2



so again,
the REAL QUESTION.. The Ultimate Question of Life, Universe and Everything:

*** what is the "AVERAGE" of INVERSE SQUARE change of rate in acceleration (jolt, surge) ? ***


42?


answer that question,
and, in return, i'll simulate 'LIFE' inside the computer.. virus or bacteria perhaps, on atomic or molecular level


how does that sound?

[cjameshuff]

A.) "NUMERICAL INTEGRATION" is completely misplaced in its practical implementation, time stepping algorithm

B.) "NUMERICAL INTEGRATION" is completely NOT NECESSARY in this algorithm, and as such, it is actually a cause of inaccuracy and slowdown

The time stepping process, accumulating a quantity over time, is numerical integration.

"Euler method" is all you need.. but maybe i can argue 3rd point as well

Due to the error involved when the function being integrated is not constant...no, it isn't. About the only thing it's useful for is a teaching tool, and rough approximations where processing power and code size are very limited.

i dont see why would that be called "Euler method" and why would there be any kind of "numerical integration"

It's called the "Euler method" because Leonhard Euler developed much of the math related to it. It's a method for performing numerical integration.

sure, since it is iterative loop it is some sort of "integration" in its nature, but calling it fancy names while it is a matter of simple kinematic equations is misleading, yes?

No. It is integration. That's not a technicality, or an obscure and trivial analogy...it's what it fundamentally is. Where do you think those kinematic equations came from?

so, the REAL QUESTION is what is the "AVERAGE" of SQUARE INVERSE change of rate in acceleration (jolt, surge)

That depends on your trajectory through the field. It's more than tricky...it requires numerical algorithms to compute for arbitrary numbers of objects.

Also, there are often other effects...Lorentz forces that depend on velocity, drag which depends on both velocity and the surrounding medium (in the case of planets, on winds and density profiles), photon and solar wind pressure which depend on solar weather, variations of planets from perfect spheres, the measured performance of rocket engines...

*** what is the "AVERAGE" of SQUARE INVERSE change of rate in acceleration (jolt, surge) ***

Use a numerical approximation to compute it. It's generally not (a1 + a2)/2, because acceleration is usually not changing at a constant rate.

[abaraba]

i dont understand what is the point of your message


- you do not prove wrong any of my statements

- you do basically agree with all i said

- you just rearranged the words and made fancy labels


[edit:]
>"sure, since it is iterative loop it is some sort of "integration" in its nature, but calling it fancy names while it is a matter of simple kinematic equations is misleading, yes?"

your answer to this question is therefore - YES
(but are you doing it on purpose to confuse everyone or you just confused yourself..)




please, have in mind that i can actually CONFIRM in PRACTICE that im telling the truth, its very pointless to argue about some empty meanings

i expect now some of moderators will give you a lecture.. >>"Politeness is the top rule here."

i think if you plan on blindly attacking me, you should consider that you failed to provide your algorithm... anyway, why so angry?

[abaraba]

ooops sorry,
i have to address this, since its very, very.. terribly, awfully wrong.. wrong, wrong, wrong, wrong ...FALSE!


ME:
>"Euler method" is all you need..

YOU:
>>"Due to the error involved when the function being integrated is not constant...no, it isn't."


while(TRUE != FALSE) print( " WRONG, WRONG, WRONG !!! " );

once you realize what that "function" you're talking about actually is, you will agree with me, even more


YOU:
>>"About the only thing it's useful for is a teaching tool, and rough approximations where processing power and code size are very limited."

while(TRUE != FALSE) print( " WRONG, WRONG, WRONG !!! " );


my algorithm will give EXACT and PERFECT solution with ALL,
but ONE SPECIAL CASE - when acceleration was nonuniform for the given time interval, just like division has special case - "division by zero"

BUT,
even with this special case, with nonuniform acceleration, my algorithm will be "MORE CORRECT" and FASTER than whatever algorithm you have there based on these wrong assumptions - wanna bet? give me your algorithm!


(the keyword here is actually: "constant COLLISION DETECTION" )

[Veeger]

Are we talking about two-body or n-body solutions? I am not understanding how acceleration can be constant in an n-body system when the effective acceleration is the resultant of all of the other acceleration vectors also acting in the system.

[abaraba]

>>"Are we talking about two-body or n-body solutions? "

- am talking about Life, Universe and Everything, so about n-body, which includes two-body


but im also talking about all other everyday simulations, all the games, soccer, billiard, ping-pong, racing, flying, falling, shooting cannon, all those simulations are under CONSTANT GRAVITY, constant altitude

so, EULER WILL DO PERFECTLY

but to answer THE QUESTION, we then talk about n-body problem, we talk about nonuniform acceleration, or more precisely we talk about INVERSE SQUARE since its the same thing with planets and atoms...


>>"I am not understanding how acceleration can be constant in an n-body system when the effective acceleration is the resultant of all of the other acceleration vectors also acting in the system."

- its not, but it is where it is, on Earth surface for example which will satisfy the most, and again Euler is COMPLETE and PERFECT solution for "everyday" simulations on Earth surface, or some other "constant" altitude - constant "external' acceleration


As for for 'n-body problem' and nonuniform acceleration - so we can simulate planets and atoms - we need to address this:


The Ultimate Question of Life, Universe and Everything:
- what is the "AVERAGE" of INVERSE SQUARE change of rate in acceleration?

answer that question,
and, in return, i'll simulate 'LIFE' inside the computer.. virus or bacteria perhaps, on atomic or molecular level



...or maybe even on planetary level.. have you considered the possibility that our galaxy might be just a molecule inside some huge interstellar virus?


but don't forget that our equations are flawed as they are,
masses are not point masses and size does matter, not really in relative way, such that atoms could easily be the size of planets, that's possible, size does not matter in that sense - but it matters as a property of "distance to mass" or some sort of "density parameter" - this actually has quite to do with practical case, because its something you can start with in order to make decision on the size of time-step (granularity) ..everything is about how to "control the error" and recognize how errors form, how they accumulate with time and how decimal precision of computer hardware plays its role in all this...

[Veeger]

Sorry, I am not interested in philosophy or science fiction. I am interested in n-body simulations which is what I thought I understood from the title of the opening post.

[abaraba]

did i not answer your question before?

>>"Are we talking about two-body or n-body solutions? "
- ..so about n-body, which includes two-body

>>"I am interested in n-body simulations.."

- yes, this is one of the better ways to go about solving the n-body problem



>>"Sorry, I am not interested in philosophy or science fiction."

- im sorry, can't help you there,
its really unfortunate since it is all interconnected..

[cjameshuff]

my algorithm will give EXACT and PERFECT solution with all but ONE SPECIAL CASE - when acceleration was nonuniform for the given time interval, just like division has special case - "division by zero"

That is not a special case, it is the general rule. Euler's method only gives exact results when the derivative is constant...that is a special case.

even with this special case, with nonuniform acceleration, my algorithm will be MORE CORRECT and FASTER than whatever algorithm you have there based on these wrong assumptions - wanna bet? give me your algorithm!

Velocity Verlet? Or Beeman's algorithm? They're simpler to implement than RK4. I'm not claiming to have any super-special secret algorithm that gives better results. I'm claiming everybody uses the more complex algorithms because they need to in order to get good results.

Have you actually tried this at all? It's extremely clear exactly why Euler's method is inaccurate. And in fact, an orbit simulation using Euler's method will very quickly diverge from anything resembling reality, even with a huge number of steps. Higher order methods can give stable results for many orbits with far larger steps. The other methods do require more processing per step, but can achieve equivalent accuracy with far fewer steps, and so are much faster overall, not to mention avoiding roundoff error.

[Veeger]

but don't forget that our equations are flawed as they are,
masses are not point masses and size does matter, not really in relative way, such that atoms could easily be the size of planets, that's possible, size does not matter in that sense - but it matters as a property of "distance to mass" or some sort of "density parameter" - this actually has quite to do with practical case, because its something you can start with in order to make decision on the size of time-step (granularity) ..everything is about how to "control the error" and recognize how errors form, how they accumulate with time and how decimal precision of computer hardware plays its role in all this...

If I'm understanding you well (and I doubt I am), considering the size of object and the sum of its relative gravity as anything other than a point-source, adds a whole other layer of complexity to the equations. Even if the goal is to only consider size in computationally trying to determine and control error accumulation. Consider when two masses are very close together and their collective gravities are acting on a third mass.