2
$\begingroup$

This question is about techniques to calculate orbits.

The Github post On the dynamical stability of Principia's modified Jool system (cited in this answer) says:

Principia computes the trajectory of the celestial bodies by integrating the equations of motion1; as a result, if the system is unstable, it may break down in-game. This is in fact the case of the stock system: while the specifics depend on how KSP's Keplerian orbital elements are translated into a Cartesian initial state, with Principia's interpretation as hierarchical Jacobi elements, the Jool system breaks down within 19 days, with a close encounter between Vall and Laythe.

Per Wikipedia:

In vector calculus, the Jacobian matrix of a vector-valued function in several variables is the matrix of all its first-order partial derivatives. When this matrix is square, that is, when the function takes the same number of variables as input as the number of vector components of its output, its determinant is referred to as the Jacobian determinant. Both the matrix and (if applicable) the determinant are often referred to simply as the Jacobian in literature.

When I do numerical integration for problems that I am currently facing I used to use the old scipy.integrate.odeint but am switching to scipy.integrate.solve_ivp and the latter says:

Parameters: fun callable: Right-hand side of the system. The calling signature is fun(t, y). Here t is a scalar, and there are two options for the ndarray y: It can either have shape (n,); then fun must return array_like with shape (n,). Alternatively, it can have shape (n, k); then fun must return an array_like with shape (n, k), i.e., each column corresponds to a single column in y. The choice between the two options is determined by vectorized argument (see below). The vectorized implementation allows a faster approximation of the Jacobian by finite differences (required for stiff solvers).

update: I think that this might be a clue; a hierarchical system might be useful when estimating interactions between objects that vary inversely with the square of distance; you might want to only evaluate interactions when orbiting bodies pass close to each other.

I've never used a Jacobian myself but I'd like to understand this more.

In the document the stiffness of the problem during a close approach is a central issue.

The matrix shown in the Wikipedia definition looks flat, not hierarchical:

$$\mathbf{J_{ij}} = \frac{\partial f_i}{\partial x_j}$$

or as the article suggests a bunch of transposed gradients.

Question: What exactly does an n-body numerical simulator's "interpretation as hierarchical Jacobi elements" mean? How would some kind of hierarchy be used in a Jacobian used by an n-body problem solver; this one or anything similar?

Would this be used in the case of a swarm or constellation of small Earth satellites who's gravitational interaction is not considered, or would the bodies need to interact with each other in some way (gravitationally or responsively) before said hierarchy has any utility?

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ I don't know from Jacobians, but hierarchy in n-body simulations with large n is used to reduce the complexity from the naive O(n^2). For example, the subsystem of Jool and its moons could be simulated in isolation, then approximated as a single point-mass body within rest of the Kerbol system; Kerbol and all its planets and their moons could be approximated as a single body to the local star cluster level of simulation, and so on. $\endgroup$
    – Russell Borogove
    Commented Jan 22, 2021 at 3:09
  • $\begingroup$ @RussellBorogove Yes this might be the kind of thing it refers to. I was re-reading the scipy.integrate.solve_ivp documentation, saw the word "approximation" and then added an update section at the same time. $\endgroup$
    – uhoh
    Commented Jan 22, 2021 at 3:13
  • 2
    $\begingroup$ I don't know n-body simulations, nor KSP, but I do know Jacobians decently well (wrote a paper on computing them more efficiently). I can say this: I am unaware of any Jacobian which is called "hierarchical." The finite differencing method uses the definition $f'(x) = (f(x+h) -f(x))/h$ to get the derivative of a multi-variate function. Maybe there's a computation method which allows and ordered & parallelized computation for finite differencing, but it should not change the result per se. Note that finite differencing is half the precision of your computer, cf AAS 19-716 conference paper. $\endgroup$
    – ChrisR
    Commented Jan 22, 2021 at 6:57
  • 2
    $\begingroup$ Indeed, 'hierarchy' is not a concept in the Python (or other) numerical libraries dealing with Jacobians. Calculating the Jacobian by either finite difference (as @ChrisR notes) or an explicit function (as determined by you symbolically differentiating and implementing it) is often an option, depending on the numerical approach you choose. There are tradeoffs to either approach, depending on the problem at hand. $\endgroup$
    – Jon Custer
    Commented Jan 22, 2021 at 16:31
  • 1
    $\begingroup$ @uhoh, with finite differencing it would be trivial to not compute specific parts of the Jacobian by simply not varying that part of the state vector. For automatic differentiation it should also be possible by setting the dual number part to zero for that component. Finally, for the analytical derivative, it is also possible to ignore some components albeit tricky if your library uses vector math. $\endgroup$
    – ChrisR
    Commented Jan 22, 2021 at 17:02

1 Answer 1

Reset to default
3
$\begingroup$

I reserve the right to be totally wrong on this one! That said, I'm throwing this into an answer rather than a comment because I hope it can be the jumping board for a solid community answer.

I very much doubt the bolded text refers to the Jacobian matrix (presumably containing all first-order partials of all the system variables against one another). Rather, given that they say they are converting KSP's definition of where celestial bodies are at epoch (provided in Keplerian elements) to a Cartesian form, I am expecting they are next defining the system in a hierarchy of Jacobi coordinates, in which positions are coordinated relative to the barycentre of a n-body system for small n. (In the bolded text, the authors use the term "Jacobi elements," which I have not encountered before, but elsewhere in the documentation they use the more standard "Jacobi coordinates")

I'm not sure if this hierarchy is hard-coded or dynamically assigned, but either way, I suspect it is used primarily to reduce floating-point error in position calculations between close objects. Unlike what some comments have suggested, this does not mean that Principia approximates n-body gravitation by only allowing nearby masses to interact--again (I suspect), this is solely a method for coordinating bodies while mitigating floating-point error in the relative positions.

For a definite answer, I suggest hitting up either one of the Leroys, either through Github or on this discord in which Principia development discussion occurs (I am not on there any more, but egg has always been kind to humor my questions about their numerical methods). Please share with the SE what you find!

Improve this answer
$\endgroup$
1
  • 2
    $\begingroup$ Oh, I think this will turn out to be the right answer, yes. I've never heard of Jacobi coordinates but after re-reading "Cartesian initial state, with Principia's interpretation as hierarchical Jacobi elements" it certainly seems like this is exactly what it's talking about. Thanks! $\endgroup$
    – uhoh
    Commented Jan 23, 2021 at 0:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.