Runge-Kutta for gravitational N-Body simulation - prediction of acceleration

[rph]

Runge-Kutta for gravitational N-Body simulation - "prediction" of acceleration

Hey!
I'm Raph, as you can see this is my first post, so, thus I'll start saying that it's nice to be here ;).

I've got the following question:
I'm simulating a system of N gravitationally interacting objects, and I implemented a 4th order Runge-Kutta method for the numerical integration. For example, the first two of the four slopes needed to determine the velocity are something like:
$$k_{v,1} = a (t , r)$$
$$k_{v,2} = a (t+0.5\Delta t , r + 0.5 \Delta t k_{v,1} )$$
where "a" denotes the acceleration, "\Delta t" is the step size and "r" the position vector.
Of course I intend "a" to be the total acceleration on the object I'm considering, calculated by summing up all the accelerations due to the other objcets.

Finding
$$a ( r + 0.5 \Delta t k_{v,1} )$$
would be really easy, as I just would have to move the object with position "r".
But to find
$$a (t+0.5\Delta t , r + 0.5 \Delta t k_{v,1} )$$
I have to "predict" the positions of the other objects at time "t+0.5\Delta t".

What is the most common way to do so? I could simply do an Euler-Step of 0.5 \Delta t, but then the position of the object in consideration would differ from
$$r + 0.5 \Delta t k_{v,1}$$
as requested by RK4. What do I do wrong? What didn't I understand right?

(Some more pieces of information: For the restricted three-body-problem that "prediction" was trivial, as the acceleration "a" was an easy function of "t" as the positions of the other objects could be predicted independently from the third object's position.

What is confusing me most is that "he": http://spiff.rit.edu/richmond/nbody/OrbitRungeKutta4.pdf doesn't even write the time-dependence of the acceleration in his equations for the slope-calculation... )

Would be great to find some help!
Raph

 

[AlephZero]

I think you are getting confused trying to solve a second-order equation directly with R-K and trying to figure out what the "prediction" fornulas should be.

You need to convert the second-order equation to two first-order equations, and then think of the displacement and velocity as two independent variables. In your PDF reference, that is equations (7) and (8).

That means you are solving the equation "velocity = dx/dt" numerically, not treating it as an exact mathematical formula. The numerical values of the velocity and displacement will not be perfectly consistent with each other, but don't worry about that, it is just a conseqence of the fact that the numerical solution will never be the "exact" solution of the differential equation.

 

[rph]

Thank you for your reply, AlephZero!

I hope I got you right but of course, I converted the 2nd order diff. eq. into two 1st order diff eq. :-)

The part I don't understand is the following:
The acceleration a is function of the time -- as the objects of the system are in motion.
To compute
$$a(t+0.5 \Delta t, r + 0.5 \Delta t k_1)$$
I thus have to do more than just shifting the object of consideration to position
$$0.5 \Delta t k_1$$
but I also have to predict the other objects' positions at time
$$t+0.5 \Delta t$$

How is this usually done? I think treating the problem as done in the PDF I cited would be right for static systems of masses, but not for dynamic systems.

Or I am totally confused and on the wrong track? :D

 

[AlephZero]

You do each step of the R-K algorithm for all the particles. You don't try to do all 4 steps for the first particle, then all 4 steps for the second one, etc.

So, you start by calculating (and storing) the value of k1 for each particle.
Then you predict (and store) the position of all the particles as $y_i + (h/2) k_{1i}$
Then you calculate all the values of k2 assuming those particle positions, and so on.

The notation above is a bit scrappy, because you actually calculate 2n values of k1 for the n particles, n for the displacements and n for the velocities. But I hope you get the idea.

 

[rph]

Excellent, thank you very much! That's actually what I implemented :D
-- but it felt kind of wrong, as in the textbooks I studied the four slopes are always given at once, and I thought there would have to be another "trick" to predicht the other object's state.
Well, thank you again!

 

[AlephZero]

Well, you can only draw pictures of the slopes for a problem with 1 degree of freedom, and your question doesn't arise in that case.

I agree this stuff can be confusing, especially if you think too hard about what is the "right" way to do it. Really, RK is just a nice computational way of creating several terms of a Taylor series, without the need to differenate the differential equation several times to get explicit formulas for the higher derivatives.

It's a bit misleading to think of the values used in the intermediate steps as "approximate solutions" of the DE. They are just a recipe for finding the end-poiint of the step, consistent with Taylor series expansions of all the functions involved in the DE. Unlike some other methods, RK doesn't have any concept of the succesive k values "coverging" to the "best" solution.

 

[D H]

Unlike some other methods, RK doesn't have any concept of the succesive k values "coverging" to the "best" solution.

Yes and no. There are extensions to the basic RK technique such as the widely implemented adaptive Runge-Kutta-Fehlberg 4/5 that do have this concept. Basically, with a clever arrange of the Butcher tableau one gets both an RK5 and and RK4 solution that can be used to adjust the time step to what is supposedly the optimal value.
One problem with Runge Kutta techniques is also one of their strengths: They are general purpose. In the problem at hand, the technique doesn't "know" that this six dimensional first order differential equation is in fact a three dimensional second order DE that necessarily conserves energy, linear momentum, and angular momentum. RK techniques don't take into account the geometry of the problem. They can't; RK integrators are a general purpose solution to the problem of numerically solving first order differential equations of all sorts.

Techniques that do account for the geometry (geometric integrators) can fare better than will general purpose techniques. They are also inherently special purpose. (Note however than JPL does quite nicely with its planetary ephemerides using what is essentially Adams Bashforth Moulton, another family of general purpose integrators.)

One last point on this: I've found that using fourth order Adams Bashforth as a predictor and a third order Adams Moulton as a single corrector (no iteration) has about the same accuracy as RK4 for the N-body problem, but at only half the computational burden of RK4.

 

[rph]

Thank you again for your input. I'll try out Runga-Kutta-Fehlberg.
Can you recommend me a textbook I could look up in our library covering numerical methods applied to physical problems in general, or, if no such book exists you enjoyed reading, a textbook about methods for numerical integration?
Thx, Rph

 

[D H]

Adaptive step size Runge Kutta Fehlberg is a bear to implement. Matlab's ode45, for example, has some problems getting the step size right -- and that's code written by people who know what they're doing. There is little point to using one of the Runge Kutta Fehlberg methods without making the step size adaptive.

As far as books go, some numerical methods texts have chapters on numerical solutions of ODEs. I've got two really old ones, and can't really recommend either.There's also the Numerical Recipes family of books by Press et al. This is an okay starting point for lots of numerical programming techniques. Keep in mind that it is a starting point. The algorithms are frequently overly elementary, oftentimes out of date, sometimes flawed, and the code is essentially Fortran 77, despite the language on the cover. Despite my criticism, I still turn to my old and battered Numerical Recipes on occasion.

The chapter in Numerical Recipes on ODE initial value problems covers a couple of the Runge Kutta techniques, including standard RK4 and a full section on the RK Fehlberg techniques. The book (in my opinion) puts way too much emphasis on Bulirsch–Stoer. That technique is nowhere near as accurate, fast, or stable as the authors like to think it is. The book downplays predictor corrector methods, which are still in widespread use.

Missing from the book: Low order techniques such as Verlet and its variants (this is a big sin, IMHO), high order techniques (although those are admittedly specialized), variation of parameters, transform techniques, and no mention at all of symplectic or geometric integration (this is another sin). Do a very large N-body simulation and you have no choice but to use a low order technique. You'll still want that simulation to obey some basic laws of physics such as conservation of energy. That means you'll want something like symplectic Euler or velocity Verlet.One last point: Since geometric integration has been around for about ten years, books are starting to come out on the topic. (Use google books.) However, unless you know what manifolds are, the math will get over your head right quick.