New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

added tests, changed rtol and numiter in kepler #362

Merged
merged 4 commits into from Apr 29, 2018

Conversation

2 participants
@nikita-astronaut
Copy link
Contributor

nikita-astronaut commented Apr 27, 2018

The kepler method suffered from not convergence on large propagation times, when xi (universal formulation parameter) is high.

It turns out, that the method converges in a physical sense, but numerically the requirement rtol < 1e-10 is not satisfied but to the floating point operations limitation.

After relaxing the requirement xi - xi_new to 1e-7, the method converges and provides physically accurate results. The other thing one had to do is to increase the maximum number of iterations done by the Newton method. The convergence check is done after each iteration, so this would not increase the propagation time for easy orbits, but will make the convergence of hard ones possible.

I added tests, checking that the orbits stated in the old issue converge (iss for 100 years and halleys for 100 years) and that their orbital elements conserve, meaning that the propagation is done right.

@Juanlu001
Copy link
Member

Juanlu001 left a comment

In general, my main concern about raising the tolerance to make the tests pass is that we are not sure whether we will find another case that requires it to raise it again. And my other main concern is that, the more we raise the tolerance, the less we know whether our results are correct.

Ideas:

  • Compare the output to similar Keplerian propagators in SPICE, GMAT or STK
  • Remove the final tolerance check, estimate the error of the solution and just return the solution and the error
  • Do a better analysis on whether f• g - g• f suffers from catastrophic cancellation or high truncation errors when xi tends to infinity

The conclusion of all of this might be "the kepler propagator is only accurate to a certain point", but then it would be nice to say something like "but it's faster/? than mean_motion". If it's only slower and worse than the other one, perhaps we just have to remove it 😉

halleys = dastcom5.orbit_from_name('1P')[0]
r_mm, v_mm = halleys.propagate(tof, method=mean_motion).rv()
r_k, v_k = halleys.propagate(tof, method=kepler).rv()
assert_quantity_allclose(r_mm, r_k, rtol=1e-6)

This comment has been minimized.

@Juanlu001

Juanlu001 Apr 27, 2018

Member

Better to hardcore the r, v values (from the days of DASTCOM5)

@@ -248,7 +246,7 @@ def _kepler(k, r0, v0, tof, numiter, rtol):
xi_new = xi + (sqrt_mu * tof - xi * xi * xi * c3_psi -
dot_r0v0 / sqrt_mu * xi * xi * c2_psi -
norm_r0 * xi * (1 - psi * c3_psi)) / norm_r
if abs(np.divide(xi_new - xi, xi_new)) < rtol or abs(xi_new - xi) < rtol:
if abs(xi_new - xi) < 1e-6:

This comment has been minimized.

@Juanlu001

Juanlu001 Apr 27, 2018

Member

Why hardcoding the tolerance here?

@nikita-astronaut

This comment has been minimized.

Copy link
Contributor Author

nikita-astronaut commented Apr 28, 2018

Hello! Thanks for the review. I would manually set 1e-7 to the _kepler function. This value is perfect: if it is smaller, say 1e-8, then we run into machine numerical problem and the algorithm won't converge. If it is less, then the method precision is not high enough. The value 1e-7 allows to reach the maximum precision in kepler allowed by the machine floating point boundaries.

@codecov

This comment has been minimized.

Copy link

codecov bot commented Apr 28, 2018

Codecov Report

Merging #362 into master will decrease coverage by 0.12%.
The diff coverage is 100%.

Impacted file tree graph

@@            Coverage Diff            @@
##           master    #362      +/-   ##
=========================================
- Coverage   82.52%   82.4%   -0.13%     
=========================================
  Files          31      31              
  Lines        1574    1574              
  Branches      123     123              
=========================================
- Hits         1299    1297       -2     
- Misses        245     246       +1     
- Partials       30      31       +1
Impacted Files Coverage Δ
src/poliastro/twobody/propagation.py 93.33% <100%> (-2.23%) ⬇️

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 fb50528...8653da6. Read the comment docs.

@@ -15,6 +15,8 @@
from poliastro.twobody.propagation import cowell, kepler, mean_motion
from poliastro.examples import iss

from poliastro.neos import dastcom5

This comment has been minimized.

@Juanlu001

Juanlu001 Apr 29, 2018

Member

This import is no longer needed?

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Apr 29, 2018

I added a minor style comment but anyway I will merge this. Thanks!

@Juanlu001 Juanlu001 merged commit ecb788f into poliastro:master Apr 29, 2018

5 checks passed

codeclimate All good!
Details
codecov/patch 100% of diff hit (target 82.52%)
Details
codecov/project Absolute coverage decreased by -0.12% but relative coverage increased by +17.47% compared to fb50528
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details

@nikita-astronaut nikita-astronaut deleted the nikita-astronaut:fixing_kepler branch Apr 29, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment