Add python and C implementation of dop853

[henrysky]

This PR adds python and C implementation of DOP853 integrator to galpy and test. DOP853 is an explicit Runge-Kutta method of order 8(5,3) by Dormand & Prince with dense output.

Here is a simple testing I have done so far to see the performance and accuracy of leapfrog and dop853 in galpy

from galpy.potential import MWPotential2014
from galpy.orbit import Orbit
import numpy as np
import time
from astropy import units as u

# Orbit Integration using galpy for the Sun
op = Orbit([0., 0., 0., 0., 0., 0.], radec=True, ro=8., vo=220.)
ts = np.linspace(0, 20, 10000) * u.Gyr
start_time = time.time()
op.integrate(ts, MWPotential2014, method="method_here")
print(time.time() - start_time, "s")
E = op.E(ts)
print("Energy accurancy: ", np.std(E)/np.fabs(np.median(E)))

For leapfrog:

27.408s
Energy accurancy: 4.894e-09

For dop853:

8.446s
Energy accurancy: 1.361e-13

For dop853_c:

0.082s
Energy accurancy: 1.177e-12

For dopr54_c:

0.199s
Energy accurancy: 2.625e-10

For symplec4_c:

0.220s
Energy accurancy: 4.528e-12

[jobovy]

travis python 2.7 builds are failing with

from galpy.util.leung_dop853 import dop853
     File "/home/travis/build/jobovy/galpy/galpy/util/leung_dop853.py", line 467
       return np.concatenate((x[len(x) // 2:], func(x[:len(x) // 2], *args, t)))
   SyntaxError: only named arguments may follow *expression

which seems to work in python 3, but not in python 2.7.

[codecov-io]

Codecov Report

Merging #363 into master will increase coverage by <.01%.
The diff coverage is 100%.

@@            Coverage Diff            @@
##           master    #363      +/-   ##
=========================================
+ Coverage    99.8%   99.8%   +<.01%     
=========================================
  Files         144     147       +3     
  Lines       23079   23578     +499     
=========================================
+ Hits        23034   23533     +499     
  Misses         45      45

Impacted Files	Coverage Δ
galpy/orbit/orbit_c_ext/integrateLinearOrbit.c	100% <100%> (ø)	⬆️
galpy/util/leung_dop853.c	100% <100%> (ø)
galpy/orbit/linearOrbit.py	100% <100%> (ø)	⬆️
galpy/util/leung_dop853.py	100% <100%> (ø)
galpy/util/leung_dop853.h	100% <100%> (ø)
galpy/orbit/Orbit.py	100% <100%> (ø)	⬆️
galpy/orbit/integratePlanarOrbit.py	100% <100%> (ø)	⬆️
galpy/orbit/orbit_c_ext/integrateFullOrbit.c	100% <100%> (ø)	⬆️
galpy/orbit/orbit_c_ext/integratePlanarOrbit.c	100% <100%> (ø)	⬆️
galpy/orbit/FullOrbit.py	99.68% <100%> (ø)	⬆️
... and 5 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 451d8c7...94bd4f4. Read the comment docs.

[jobovy]

I went through all but the leung_dop853.py code, here are some minor changes that are required.

[jobovy]

Started testing this, looking at how the actual orbit (like o.R(ts)) compares to other integrators rather than just energy. It seems like there's something wrong with the C implementation, because the radius has large deviations compared to other integrators, which the python implementation does not have. Example:

import numpy
from galpy.potential import MWPotential2014
from galpy.orbit import Orbit
%pylab inline
o= Orbit()
ts= numpy.linspace(0.,100.,1001)
o.integrate(ts,MWPotential2014,method='dop853')
od= o()
od.integrate(ts,MWPotential2014,method='dop853_c')
os= o()
os.integrate(ts,MWPotential2014,method='symplec4_c')
plot(ts,od.R(ts,use_physical=False)-os.R(ts,use_physical=False),label='dop853_c - symplec4_c')
plot(ts,o.R(ts,use_physical=False)-os.R(ts,use_physical=False),label='dop853 - symplec4_c')
xlabel(r'$t$')
ylabel(r'$\Delta R$')
legend();

gives the following:

The appveyor build also fails because the orbit tests time out (they take an hour, used to take 25 minutes), which happened twice, so likely something to do with the new C code (travis is fine though...).

[henrysky]

Started testing this, looking at how the actual orbit (like o.R(ts)) compares to other integrators rather than just energy. It seems like there's something wrong with the C implementation, because the radius has large deviations compared to other integrators, which the python implementation does not have. Example:

I gave dense output in c implementation the wrong time which leads to wrong result. Now it behaves exactly like dop853 python implementation. I also changed rtol and atol from 1e-14 to whatever galpy gives. i did not know that is exponential but fixed now

The appveyor build also fails because the orbit tests time out (they take an hour, used to take 25 minutes), which happened twice, so likely something to do with the new C code (travis is fine though...).

SIGINT test gets stuck. It happens on my Windows too. Seems like the issue is python cannot get exit code correctly on Windows. CTRL+C works in practice tho. So I removed them from test for now until I figure it out whats the deal on Windows.

[jobovy]

Great, that test seems to pass now!

You should also add dop853_c to the 'liouville' test in test_orbit.

There are a few lines in leung_dop853.c not being touched by the tests now:

• There are two interruption blocks? Because there is only a single interrupt in the tests only one of these gets hit. Maybe you can remove one? We don't care that the interruption happens immediately, so perhaps it can just be in the outer loop?
• A special case, probably difficult to hit with a test...

[jobovy]

Started looking at this again and this all looks quite good to me now (I also did an explicit test of the Python DOP853 implementation with a differential equation with an analytical solution, which passed). But I did one more tests that gives somewhat odd results (although may be nothing...):

import numpy
from galpy.potential import MWPotential2014
from galpy.orbit import Orbit
%pylab inline
o= Orbit([1.3,-0.5,1.,0.])
ts= numpy.linspace(0.,300.,1001)
o.integrate(ts,[MWPotential2014],method='dop853')
od= o()
od.integrate(ts,[MWPotential2014],method='dop853_c')
od.plotE()
o.plotE(overplot=True)

gives

Apologies for the lack of labels, but the green line is the Python implementation of DOP853 and the blue line the C implementation: the Python version has a much more stable energy error. Do you understand why this is?

Finally, could you add an entry in HISTORY.txt about this addition and add yourself in the Direct contributions to the code base: section of AUTHORS.txt about this as well?

[henrysky]

I found a bug that the time step is being reset to hmax and the steps are rejected and scaled down until acceptable error. dop853_c is now a bit faster due to time step not being reset to hmax in every loop. But this does not seems to change the result. Still not sure why C implementation seems worse than python implementation

[jobovy]

If anything, with this latest change things seem a little worse still (but again, these are still very small errors):

[henrysky]

I try to trace every variable and found python and C version differ even in very early in initial step estimation.

And I check with the initial condition array given by galpy, the initial condition array in c-implementation is [1.3 0. 0.5 1.] and in python is [1.3 0.5 0. 0.76923077]. Is it supposed to happen??

[jobovy]

Yes, that's expected. In C the equation of motion is integrated in rectangular coordinates, so the initial conditions are x,y,vx,vy (I think you mean -0.5); in Python the equation of motion is integrated in cylindrical coordinates and the initial conditions are R,vR,phi,vT/R.

[jobovy]

But this is a good point and I think the difference in the equations of motion might explain the difference. Because the potential is axisymmetric and angular momentum is conserved, integrating the system in cylindrical coordinates is more stable. So there is likely nothing wrong with the DOP853 integrator!

[henrysky]

I thought the C version also integrate in different coordinate between symplectic and non-symplectic like the python one you told me? Because I see the calls to C integrators are different between RK/DOP54 and symplectic integrators. I do find lowering the error tolerance of dop853_c to something like 1e-14 can bring it much closer to the python implementation and I do think the python and c implementation is working the same way. I cant find any difference anymore.

[henrysky]

The cause is clearly coordinates now. I use my standalone python implementation of dop853 in solar system. The integration using rectangular and polar coordinates show the exact same trend.

import pylab as plt
import numpy as np
from leung_dop853 import dop853
from scipy.integrate import odeint


G = 4*np.pi**2 # define G for convenience.

# Set initial position and speed of satellite.
# Compute and display some values for a circular orbit to help decide on initial conditions.
Ro = 1.0    # Radius of the orbit 
M  = 1.0    # mass of the central mass

v = np.sqrt(G*M/Ro) # v assuming circular orbit

# Set initial conditions and define needed array.
x0 = 1.
vy0 = 5.
t0 = 0.
tf = 10.
tau = 1e-4

X0 = [x0, 0, 0, vy0]       # set initial state of the system
t = np.arange(t0,tf,tau)   # create time array starting at t0, ending at tf with a spacing tau

# Define function to return f(X,t)
def f_func(state, time):
    f = np.zeros(4)
    f[0] = state[2]
    f[1] = state[3]
    r = np.sqrt(state[0]**2 + state[1]**2)
    f[2] = -(G*M*state[0])/r**3
    f[3] = -(G*M*state[1])/r**3
    return f

X_dop853 = dop853(func=f_func, x=X0, t=t, rtol=1e-12)
x_dop853 = X_dop853[:,0]  # Giving a ':' as the index specifies all of the elements for that index so
y_dop853 = X_dop853[:,1]  # x, y, vx, and vy are arrays that specify their 
vx_dop853 = X_dop853[:,2] # values at any given time index
vy_dop853 = X_dop853[:,3]

X0 = [x0, 0, 0, vy0]       # set initial state of the system, same as retengular but in [r, phi, vr, vphi]

# Define function to return f(X,t) in polar
def f_func_polar(state, time):
    f = np.zeros(4)
    f[0] = state[2]
    f[1] = state[3]
    f[2] = -(G*M)/state[0]**2 + state[0]*state[3]**2
    f[3] = (-2 * state[2]*state[3])/state[0]
    return f

X_dop853_polar = dop853(func=f_func_polar, x=X0, t=t, rtol=1e-12)

# Plot the results
plt.figure(figsize=(10, 10))
plt.clf()

# calculate KE + PE
eps_dop853 = 0.5*(vx_dop853**2 + vy_dop853**2) - G*M/np.sqrt(x_dop853**2 + y_dop853**2)
eps_dop853_polar = 0.5*(X_dop853_polar[:,2]**2 + X_dop853_polar[:,0]**2*X_dop853_polar[:,3]**2) - G*M/X_dop853_polar[:,0]

# Plot the energy
plt.subplot(2,1,2)
plt.title("Energy Plot", fontsize=15)
plt.plot(t, eps_dop853, label='dop853 in rectangular')
plt.plot(t, eps_dop853_polar, label='dop853 in polar')
# plt.axis([min(t),max(t),eps_plot_min,0])
# plt.ylim((-31.4784, -31.47845))
plt.xlabel('Time', fontsize=15)
plt.ylabel('Normalized Energy', fontsize=15)
plt.legend(loc="best", fontsize=15)
plt.tight_layout()
plt.savefig("standalone_integration.png")

[jobovy]

Okay, great, thanks for checking this in detail! I think all is good then and will merge this.

Yes, the integration in C is different for symplectic and non-symplectic integrators, because they use very different algorithms, so they can't be the same. But I think they are always in rectangular coordinates, even though the calls are different.