Algorithms for planet motion question

[camodude009]

Hello Dan,
I have already posted this in an old thread in the Universe Sandbox (1) discussions.
In case you do not look at those anymore, I'm making a new post here.

I am working on scientific paper about the different algorithms used for approximating the solution to the n-body problem (Euler/Verlet/RK4). I would like to could do a quick q and a on what you are using in universe sandbox 1/2 and how you deal with timestep/rounding error or any other error that accumulates. It would really help me out.

-Felix

[Greenleaf]

I will respond for Dan, since I am the one working mainly on the NBody part of the simulation.
I do it as a comment, in case others are interested in this.

Currently you can only select between explicit euler and verlet. Euler being there primarily to show how bad it is. Previously we had support for (from memory) explicit euler, semi implicit euler, rk2, rk4, RKF, adams bashforth 6th order, adams moulton 6 order, hermite 5th order, PEFRL and forrest ruth.
I made a video showcasing their individual strengths and weaknesses, which you can find at https://youtu.be/IJ2MhXUDZ6o . there may be other videos of interest at the channel.

A large NBody rewrite took most out of service, but they will be going back in shortly.

As to error, as you probably know, there is truncation error and rounding error. Truncation error from having a limited order of integration in the integrators, and rounding errors from floating point having a limited precision.
We deal with truncation error, as best we can, by letting the user pick the right integrator for the job, and setting a sensible default. If you have a fast paced chaotic simulation, computational speed is more important than a very high degree of accuracy. If you have a long running system which you want to analyze stability for, accuracy becomes more important.

Rounding error is generally not that important for this particular simulator. We are using doubles, and running on ordinary gaming computers and doing so at at generally very accelerated time speed.

That leaves time step control. How large an integration step is too large? For that, the user is given the option of setting a tolerance which is defined in units of length. For every simulation step forward in time, the integration step(s) length will be adjusted to not generate a position error larger than this tolerance.

As an example, the simulation may want to step 1 day. The integrator may then try to do this and realize that it generates a position error larger than the given tolerance. It will then calculate how long a step is in fact generating a just small enough errors, which (for the example) may be 1 hour. It then steps the 1 day in 24 sub steps of 1 hour before returning. The integration is running asynchronously from the render and input layers. Generally this generates a new physics state for every rendered frame.

The way we estimate error is by taking one step dt forward from the current state and then also take two ½dt successive steps forward from the same current state. Then those two results are compared and a good error estimate can be found using standardized methods.

It is all a bit more involved than this, but that is the essense of it.

I hope I have answered your questions. Feel free to ask more, if I left something out... which I might have. Sitting at a noisy airport on my way to Seattle right now.

[FatherPickle]

Originally posted by RLM.Camodude009:

Hello Dan,
I have already posted this in an old thread in the Universe Sandbox (1) discussions.
In case you do not look at those anymore, I'm making a new post here.

I am working on scientific paper about the different algorithms used for approximating the solution to the n-body problem (Euler/Verlet/RK4). I would like to could do a quick q and a on what you are using in universe sandbox 1/2 and how you deal with timestep/rounding error or any other error that accumulates. It would really help me out.

-Felix

Neeeerd!

[AlfredWallace]

Can you just explain me simply what is n-body ?

[lvngbth]

If you have just two things in your universe, you can use some simple maths to work out how they interact and the results will be right.

If you have three or more things in your universe, they won't be. With three, you don't just have to consider how A interacts with B, but how A interacts with C and B with C, and the more things you have, the many more interactions you need to consider. A tiny error - because of the limits of the accuracy of your calculations, say - gets rapidly enlarged.

This is 'the n-body problem'. There are a number of ways to get answers that are roughly right, some rather rougher than others...

[AlfredWallace]

Ok thanks very much for the explanation -)

[camodude009]

Hi Greenleaf, hope you see this comment. I'm the guy on steam that asked about which algorithms you use :)
Since there are a few varieties of Verlet I wanted to ask if this is the one you are using (using this is my program):

In Simstate:
public void tickVerlet(double t){
for(Planet p: thePlanets) p.updateAccelerationVerlet();
for(Planet p: thePlanets) p.updatePositionVerlet(t);
for(Planet p: thePlanets) p.updateAccelerationVerlet();
for(Planet p: thePlanets) p.updateVelocityVerlet(t);
}

In Planet:
public void updateAccelerationVerlet(){
aXOld = aX;
aYOld = aY;
aX = 0;
aY = 0;
for(Planet p : simstate.getPlanets()){
if(p != this){
double a = Sim.G*p.getM()/getDS(p);
aX += (getDX(p)/getD(p)*a);
aY += (getDY(p)/getD(p)*a);
}
}
}

public void updatePositionVerlet(double t){
x += (t*vX)+(0.5*aX*t*t);
y += (t*vY)+(0.5*aY*t*t);
}

public void updateVelocityVerlet(double t){
vX += 0.5*(aX + aXOld)*t;
vY += 0.5*(aY + aYOld)*t;
}

[Greenleaf]

Originally posted by RLM.Camodude009:

Hi Greenleaf, hope you see this comment. I'm the guy on steam that asked about which algorithms you use :)
Since there are a few varieties of Verlet I wanted to ask if this is the one you are using

We are using velocity verlet
x(t+dt)=x(t)+v(t)*dt+0.5*a(t)*dt*dt
v(t+dt)=v(t)+0.5*(a(t)+a(t+dt))*dt

So, essentially... yes ;-)

You might want to consider a higher order integrator like PEFRL, though, which I personally like, and which is not more complex to implement.
http://arxiv.org/pdf/cond-mat/0110585.pdf

[camodude009]

Thanks for the reply.
I'm trying to add 4th Order Runge-Kutta to my program and I don't quiet understand the online explanations (usually looking something like this: http://i.imgur.com/4kmpcQQ.png).
I'm at the point where I know that I need to do this for velocity and position separately but I don't know what the f(tn,xn) represents.
I would be delighted if you or someone else on this forum that understands the issue could explain this to me.
Thanks,
Felix

[Greenleaf]

Ok, I started typing out a long answer, but realized that the formatting in here doesn't really help when you try to explain something ever so slightly complicated.

Let me therefore just address your specific question, and urge you to ask more, if you are still in trouble.

tn means the time at the start of a step and tn+dt½ means the time half way through a step.
xn is the system (that you are integrating) at the start of a step as well.

By "a system", I mean the properties changing over time, the things you are trying to evolve essentially. For nbody that would generally be the positions and velocities as a whole. Always think of it as one bit system and not a number of separated properties!


When you read f(x_n) it means the derivative of the system as it is at state n. Assuming you have two bodies in a system and they are defined with some initial state: position and velocity. The f(x_n) means "the time derivative of this system as it is now". The time derivative of position is then velocity and the time derivative of velocity is acceleration. You already have the velocity, but you need to calculate the acceleration, and you do that based on the position properties of all bodies in the system at that same point in time.

Therefore the procedure is
figure out how much each property changes, given the current system properties, and multiply that by the time step length and call that k1.
Now make a new temporary system with properties equal to x_n+0.5*k1 and calculate the derivatives for that configuration. This is the f(x_n+0.5*k1) and you can get k2
Do this for k3 and k4 as well and you can do a weighted sum in the end.

The rather important thing to keep in mind is that you advance _every_ property in lock step. You do not change positions first and then use those new values to update velocities or the other way around. You define every property for a full system of bodies and act on them as a whole.

Feel free to ask for more details. Always willing to help anyone with a genuine interest in this :-)

[camodude009]

Thanks for you answer :)
I got that part. I should have phrased me question a bit differetnly.
What I meant was how do I create this temporary system with properties equal to x_n+0.5*k1
Do the planets now have 1/2 the x and y coodinates? or du I just 1/2 the speed for the k1 system?
This is my thought process atm:
when x is the initial system and t the timestep
k1 = x advanced by t
k2 = k1 advanced by 1/2*t
k3 = k2 advanced by 1/2*t
k4 = k3 advanced by t
x advanced by t = x + k1/6 + k2/3 + k3/3 + k4/6

I advance from k1 -> k2 ->... with euler and the specific timestep
in the end i add (1/6 of velocity and position of k1) + (1/3 of velocity and position of k2) + ... to the actual system to advance by t
I have a feeling I'm wrong here :/
Sorry once again for all these questions but pretty much all of these recources are hard to grasp for a highschool student without knowledge of standard notation in physics (didnt know what a dot above a letter/ dt meant until 3 weeks ago). Should have written about climate change on some island nobody knows about ^^. That wouldn't be as much fun though :)

[Greenleaf]

In US² a body has not one position, and one velocity, but multiple in an array. This is to make rk4 easy as well as to support multistep methods.

what you write about k1 is correct. x advanced by dt
k2 is however supposed to be x advanced dt using the derivatives of the system which is x+½k1. In other words, using an acceleration calculated half way in between x and x+k1

So...
you find derivatives (velocities and accelerations) for initial state x(0)
you set k1=f( x(0))*dt
Then you find the halfway state x(½)=x(0)+0.5*dt*k1
Essentially saying x(½) is where everything would be if we took a half step from x(0)

Now, to find k2, you need the derivatives at that position f( x(½) ) which is the same as f( x(0)+½k1)
so k2= dt*f( x(½)) = dt*f( x(0)+½k1)

and so on and so on.

The point is that you need to calculate temporary states and calculate derivatives for those. We do this by explicitly storing the temporary states in multiple buffers, for various reasons reasons, but you can do it in-place by just storing the k1-4 and doing a lot of addition and subtraction. In general, you do calculate states such as x(½) and send those states to the derivative calculation which returns f( x) for you to multiply with dt to get the k's.

To be very explicit. (using x(0) as the symbol for the entire system at initial state)
Assuming x(0) has a position of 1m and a velocity of 10 m/s (for one example body) and you calculate an acceleration to 5 m/s^2 and your dt is 1s.
The derivatives [ f(x(0)) ] are then 10m/s and 5m/s^2
Multiplying those with dt gives you k1 as 10m and 5m/s

That means that x(0)+½k1 will be pos=1m+½*10m and vel=10m/s+½5m/s
or pos=6m and vel=12.5m/s

That is the new temporary state for which you next have to calculate the derivatives to get k2. You calculate the new state for all bodies in your system and then those new positions and velocities are used to calculate the new derivatives : velocities and accelerations.

Then you do the same a few more times, and you are done.

I hope I didn't again answer the wrong question ;-)

[camodude009]

No :P
Ok so this is how it works:

With a timestep of t and positions x... velocities vX... and accelerations aX....
Initially I have x (postition) and vX (velocity)
I compute the acceleration aX.

I compute state k1:
the change in velocity for K1:
vXK1 = aX * t
adding to the position by old velocity + change in velocity * t
xK1 = x + (vX + vXK1) * t
aXK1 = recompute from positions

I compute k2:
the change in velocity for K1 (0.5 the derivative of K1):
vXK2 = 0.5 * t * vXK1
adding to the position by 0.5 the derivative of K1 * t
xK2 = x + (0.5 * t * vXK1)

...
is this correct?
I update the system based on the previous derivatives. When updating position I add the factor * velocity of previous system to the -> initial position? (or the one from the previous system? In the pseudo code above I used the initial position)

Thanks and I hope you can decipher this,
Felix ;)

[Greenleaf]

Ok, lets try once more.
k1 is not a state. k1 is a change over a time step.

"vXK1 = aX * t" is right
"xK1 = x + (vX + vXK1) * t" is wrong and should have been
"xK1 = vX*t"

"I update the system based on the previous derivatives"
you constantly have to find derivatives for sysyem states which are based on the initial value (of the step) plus some fraction of k1,k2,k3, so you cannot update the initial system, but need to make temporaries.

In one step you calculate the derivatives of the integrated properties.

Let us call position x and velocity v and acceleration a. Each one of these can be a 3d vector or just a 1d value. That doesn't matter. Let us say that x(0) is the initial position.
For clarity, though a little odd syntax, let us say that system(0) means [x(0),v(0)] so when you say system(someTime) you talk about the total collection of all positions and all velocities at someTime. Imagine system(0) as a structure holding one array for all x(0) and one array for all v(0).

you have x(0) and v(0), or now called system(0), and want to get k1=dt*f( system(0) )
f( system(0) ) means the derivatives of system(0) and is equal to [v(0), a(0)]

since the derivative over time of x is v and the derivative of v is a.

so calculate k1 as
k1=dt*f( system(0) )

k1 now contains change in position over dt as well as change in velocity over dt, and this is one entire block of numbers created as one.

to get k2, the change over dt using derivatives based on something slightly different from system(0), you do

k2 = dt*f( system(0)+½k1 )

system(0)+½k1 simply means that you have this temporary intermediate system inside f(..)
x(0)+½k1
v(0)+½k1

typed out even more you have this inside f(...)
x(0)+½*v(0)*dt
v(0)+½*a(0)*dt

Personally I handle this by explicitly calculating
tmp = system(0)+½k1 so I can then send tmp to my derivative calculation as f(tmp) and then multiply by dt and call the result k1.

Tmp is the temporary system with concrete values for x and v, which you need to calculate derivatives for. Again the derivative of x is v (in this case v(0)+½k1) and the derivative of v is a or in this case a( x(0)+½*v(0)*dt)

To repeat the main point...
you deal with a system here. One position and one velocity per body. Every value in here is valid at one particular point in time.
You must treat it as a whole. That means that for n bodies, you have a structure of n positions and n velocities.
When you see x(t) it means such a system as it is at time t. All of it. All positions and all velocities as a unit.
When you see f( x(t) ) is means the derivatives of an entire system as it is at time t. For n bodies, that is n position derivatives (which happen to be the same as the velocity) and n velocity derivatives, which are accelerations you calculate.
When you see k1=dt*f(x(0)), read it as
x(0) is an entire system of n positions and n velocities, stored however you like
f( x(0) ) is returning a block of n position derivatives and n velocity derivatives
dt*f(..) means that every value in the f(..) block is multiplied with dt, yielding a new number which is not a derivative, but an actual change of something over time, so you get delta position and delta velocity. As in m/s*s=m
x(0)+½k1 therefore means that every position in the system x(0) gets half its matching delta added to it and same goes for velocities.


If it still doesn't quite make sense for you, I am afraid that I have to bow out for now. I would then suggest that you find some source code. This seems usable at a glance.
http://gafferongames.com/game-physics/integration-basics/

[camodude009]

Thanks again for your answer :)
I think I figured it out now. I kept trying to do all the calculation within the planet class itself. Going to go by the gafferongames implementaion. Will post when I have it fully implemented!
Thanks for all your work,
Felix

[camodude009]

Got it working :)
Have to go now, will post pictures when I get back!