hannorein/rebound

trying to use variational equations to get Jacobian but get an error...

JonathanHammer2 opened this issue · 6 comments

JonathanHammer2 commented

Hello again! Thanks for your answer to my last question. What I'm trying to do here is to get the Jacobian of the Earth's eccentricity with respect to the cartesian positions and velocities of the other planets. So H = [de_earth/dx_mercury, de_earth/dy_mecury, de_earth/dz_mecury, de_earth/dvx_mecury... de_earth/dvz_mecury... de_earth/dx_neptune... etc.

I tried to set this up as follows, but got this error: ValueError: Variational particles can only be initializes using the derivatives with respect to one of the following: m, a, e, inc, omega, Omega, f, k, h, lambda, ix, iy.

So I'm wondering if there is a straightforward way to do what I want?

my code:
for j in range(1, 8):
if j != 2:
var_dx = sim.add_variation()
var_dx.vary(j, "x")
var_dy = sim.add_variation()
var_dy.vary(j, "y")
var_dz = sim.add_variation()
var_dz.vary(j, "z")
var_dxd = sim.add_variation()
var_dxd.vary(j, "vx")
var_dyd = sim.add_variation()
var_dyd.vary(j, "vy")
var_dzd = sim.add_variation()
var_dzd.vary(j, "vz")

hannorein commented

The syntax is:

var_dx = sim.add_variation()
var_dx.particles[j].x = 1
JonathanHammer2 commented
hannorein commented

she was concerned that the variational equations implementation in Rebound only work with celestial variables, not with Cartesian variables.

That's not correct. The "vary" function is a convenience function that allows you to vary orbital elements. The function exists because it's non-trivial to map these variations to cartesian coordinates. Varying cartesian coordinates on the other hand is trivial. See below.

I’m wondering what the 1 signifies in initializing var_dx.particles[j].x = 1?

This is varying the x coordinate of the j-th particle. This is what you want for the Jacobian you mentioned. The "1" is just a normalization. You can use another other non-zero number. That would give you the same result if you divide it out again at the end. This is just coming from the fact that variational equations are invariant under multiplications with a constant. But really, there is no good reason to use a number other than 1 here.

then I should be able to get the derivate of earth eccentricity

Your setup looks good. But you want to keep a reference to all the var_dx, var_dy, ... variables. Right now you are overwriting them in the loop. Use an array and set var_dx[j]. After the integration you can get the matrix elements:

$$ \begin{bmatrix} \frac{dx_i}{dx_j} & \frac{dy_i}{dx_j} & \ldots \\ \frac{dx_i}{dy_j} & \frac{dy_i}{dy_j} & \\ \ldots & & \end{bmatrix} $$

using the syntax

var_dx[j].particles[i].x        var_dx[j].particles[i].y     ...
var_dy[j].particles[i].x        var_dy[j].particles[i].y
...

From that you can calculate

$$ \begin{bmatrix} \frac{de_i}{dx_j} & \frac{de_i}{dx_j} & \ldots \\ \frac{de_i}{dy_j} & \frac{de_i}{dy_j} & \\ \ldots & & \end{bmatrix} $$

using the chain rule. Just write down $e$ as a function of $x, y, z, \ldots$. The take the derivative of $e$ with respect to all the cartesian coordinates, multiply each derivate with the relevant matrix element, and add them up.

I hope that helps.

JonathanHammer2 commented
hannorein commented

I'll close this issue for know, but feel free to re-open if you have more questions.

JonathanHammer2 commented