LeapfrogIntegrator will reverse the sign of velocity which may lead to incorrect result

[azz147]

Hi, thank you for the great package. I realized the below code in the LeapfrogIntegrator may change the sign of the velocity when _dt is negative, which may change the result of the force function because it could depend on the velocity. In my case I add a dynamical friction term in my force function and LeapfrogIntegrator will give me incorrect results.

gala/gala/integrate/pyintegrators/leapfrog.py

Lines 146 to 150 in 782a8b1

	if _dt < 0.:
	v0 *= -1.
	dt = np.abs(_dt)
	else:
	dt = _dt

[adrn]

Hi @azz147 - can you say more about what you would like to do? In my use cases where I have implemented dynamical friction, if I am integrating backwards in time, I do want to reverse the sign of the friction force. Are you sure you don't want this?

[azz147]

Hi @adrn , if v_i and dt are reversed then I need to keep the sign of dynamical friction (calculated by the original v_i). My friction force in F will produce a force at the opposite direction of the velocity, directly pass a reversed velocity argument will produce a force to slow down reversed velocity, but when integrating backwards in time I expect it will accelerate the reversed velocity.

[adrn]

Ah, yes, but note that in line L148 in the code you linked above the timestep is turned positive dt = np.abs(dt)

[azz147]

I think the current code tries to reverse the direction of velocity and timestep to do integrating backwards. Based on the below substitution (not the Leapfrog formula but they should be the same):

So for the force function F, it should receive the original velocity instead of the reversed one

[adrn]

I'm at a conference now so can't look closely until later this week, but I think the integrator is behaving as I would expect. For example, here's a toy: A Kepler potential with a simple/fake drag force that acts opposite the velocity vector. Integrating forwards in time, the orbit decays. Taking the final position and integrating backwards (i.e. just setting dt to a negative value), the orbit grows back to the starting position. https://gist.github.com/62138226e1daf734db1aa9cf7c12781d

[azz147]

Thanks for the demo, that is very helpful. I realized the code I referred above actually comes from an old version (I use search in GitHub to find the code of that function and somehow I get into an old version code), and 51278ee already fixed this in v1.5. My installed version stays at v1.4.1 due to python 3.7, but I guess I could try to upgrade python to get v1.5.