Problem understanding NBody simulation physics

[Bekos]

Hello everyone! :)

This is my first post in this forum! :) I am a computer programmer but the last weeks I am trying to combine my field with astronomy and astrophysics. I decided to create my own nBody simulator but I have problems understanding the calculation of the velocity and speed during the integration. My skills in physics are not that good so please bear with me. Let's start with a simple example. Assume that we have a set of bodies with some initial velocity. I thought that the position of each body is calculated like this:

for each integration frame (or whatever you can call it) do the following:
X = Xo + Uo*dt + 0.5 * A * dt * dt;
U = Uo + A * dt;
Uo = U;

X is the position at the end of the interval. Xo is the position at the start of the interval.
Uo is the velocity at the start of the interval. dt is the time interval. A is the acceleration at the start of the interval. What I do above is to calculate the position at the end of the interval using the position, velocity and acceleration at the start of the interval. Then I calculate the velocity at the end of the interval in order to use it for the next physics-frame/integration. Isn't this right?

After checking some NBody simulation examples from nVidia and Microsoft SDK I realized that they use the following to calculate the position of each particle:

U = Uo + A*dt; // This will Calculate the velocity at the end of the interval. Right?
X = Xo + U*dt;
Uo = U;

From my point of view this is wrong because here the position at the end of the interval is calculated using the velocity at the end of the interval and not the velocity at the start of the interval. Apparently, using the second method, the nbody system looks "nice" in nVidia's and Microsoft's demos. But if I replace these 3 lines of code with my solution, the nbody won't look "nice". Using the second solution I get a "galaxy-looking" result. With a lot of particles at the center of the "galaxy". With my solution, the "galaxy" still have more particles on its center than on the edges but the difference is not very noticeable like it is when using the 2nd solution.

What is the difference between these two solutions? And why my solution is wrong?
Thanks a lot for your time.

Cheers,
Bekos

 

[Jano L.]

Hello Bekos,

welcome to physicsforums:)

In fact if you want to simulate the motion of a many particle system, you have to start from the exact equations of motion and then devise an algorithm which would solve it numerically.

Your first set of equations is called Euler forward method, and unfortunately is not very accurate. The trajectories of the particles will usually diverge from the correct ones very quickly. Typical behaviour is that the energy of the system grows without limit and the particles will run away. Maybe that's why you have more particles on the rim if you use this method. What sort of force do you use to keep the particles together?

The second method is called leapfrog or Stoermer-Verlet method, if I remember correctly. It is much better for solving equations of motions in mechanics.
As you pointed out, the interesting question is, why

x_n+1 = x_n + v_n+1 *t

is better than

x_n+1 = x_n + v_n *t ?

Here is how it can understand it intuitively. Draw a graph of the function x(t) = t^2 (t is time and x is the distance from the origin). The plot is a record of the actual motion of the particle.

Now draw two points on the curve, A and B, in times t_n, t_n+1, that are 4cm apart. Then attach two arrows of length 1cm to both A and B so that they are tangential to the curve. The arrows represent velocities at times t_n, t_n+1 (derivatives of the position with respect to time).

Now, if you used the Euler formula, you would in fact use the first arrow to get from A to B: draw a dotted line extending the arrow A. You see, because the graph is curved in general, you won't get at the point B in time t_n+1 but you will come into another point, let's designate it B', which for our curve is under B.

But if you transported paralelly the second arrow from B to A, extended it as before into a dotted line, the line will come much closer to the point B, let's designate it B''.

If you choose small enough time step, for most cases, B'' is closer to B than B' and therefore will give you much better estimate of the next position x_n+1.

If you are serious about these calculations, there are much better algorithms, but they are also much more complicated... For example, try to search N-Body integrator, symplectic integrators or Runge-Kutta methods on the Internet.

 

[D H]

What is the difference between these two solutions? And why my solution is wrong?

Welcome to the forums, and welcome to the wonderful world of numerical integration.

The second method is variously called semi-implicit Euler's method, Euler-Cromer method, symplectic Euler method, plus some others. That last name, symplectic Euler, indicates why it is used on occasion. Symplectic integration techniques conserve energy (or come close to it). Your method does not.

There are boatload upon boatload of numerical integration techniques. How to evaluate them? A naive look at the equations doesn't cut it. Some techniques (e.g., symplectic Euler) just look wrong. You need to look to the results. How accurate are the results over the short term and over the long term? Is the technique stable? How sensitive is it to numerical errors? One can also look at techniques in terms of ease of programming, memory utilization, and CPU utilization. Regarding the latter, CPU utilization: Some techniques such as symplectic Euler require very small time steps to achieve even a modicum of accuracy. Better techniques can use comparatively huge time steps and still yield accuracies that could never be obtained with symplectic Euler without going to arbitrary precision arithmetic. (Note that using arbitrary precision arithmetic instead of hardware floating point would increase CPU utilization by orders of magnitude.)

A slight modification to your technique does yield something that is useful. Simply update the velocity using the average of the acceleration at times t and t+Δt:

[tex]\begin{aligned}
x(t+\Delta t) = x(t) + v(t)\Delta t + 0.5a(t)\Delta t^2 \\
v(t+\Delta t) = v(t) + 0.5(a(t) + a(t+\Delta t))
\end{aligned}[/tex]
This is the velocity verlet. While it is much better than symplectic Euler, it is not a very good integrator. People use low order integrators such as these because they don't know better, they don't care, or because they have no choice because of things like collisions. If you want accuracy and precision you have no choice but to use better integration techniques.

 

[Bekos]

Ohh my god guys... :)

Seriously, I am surprised with the almost instant help! I have a much better understanding now. Thank you both guys!
If I understood correctly, I could get precise results if I had a constant acceleration right?
Thank you very much for you help guys! I am going to have a look at some other integration methods now.

Cheers!

 

[pmacias]

Hello Bekos,

First of all, it's great to see someone from a different field trying to combine their skills with astronomy and astrophysics. The nBody simulation is a complex and challenging subject, so it's understandable that you are facing difficulties.

From what I understand, your main concern is about the calculation of velocity and position during the integration process. Let me try to explain the difference between the two solutions you mentioned.

In the first solution, you are using the velocity and acceleration at the start of the interval to calculate the position at the end of the interval. This is the correct approach because the position at the end of the interval should depend on the initial conditions (position, velocity, and acceleration) at the start of the interval.

In the second solution, the velocity at the end of the interval is calculated first, using the initial velocity and acceleration. Then, this new velocity is used to calculate the position at the end of the interval. This approach can also work, but it might result in some errors in the simulation. This is because the velocity at the end of the interval is not directly dependent on the initial conditions at the start of the interval.

It's difficult to say why the second solution is giving a "galaxy-looking" result without looking at your code and the specific values you are using. However, it's important to note that nBody simulations are highly sensitive to initial conditions and even small changes can result in drastically different outcomes. So, it's possible that the difference in results is due to some small variations in the initial conditions.

In conclusion, your first solution is the correct approach for calculating the position and velocity during the integration process in an nBody simulation. However, it's important to carefully consider the initial conditions and make sure they are accurate to get the desired results.

I hope this helps and good luck with your nBody simulation!