Opened 3 years ago
Closed 3 years ago
#28707 closed enhancement (fixed)
More control on the numerical ODE solver for integrated curves and geodesics
Reported by:  Eric Gourgoulhon  Owned by:  

Priority:  major  Milestone:  sage9.0 
Component:  geometry  Keywords:  manifolds, integrated_curve, geodesic 
Cc:  ghFlorentinJ, Karim VanAelst  Merged in:  
Authors:  Eric Gourgoulhon  Reviewers:  Florentin Jaffredo 
Report Upstream:  N/A  Work issues:  
Branch:  aca08a8 (Commits, GitHub, GitLab)  Commit:  aca08a898e22586d2d9fcb4f29e830cc788f165f 
Dependencies:  Stopgaps: 
Description (last modified by )
This ticket allows to pass extra parameters to the method IntegratedCurve.solve()
in order to control the numerical ODE solver. In particular, this allows to fix an issue with the odeint
solver by increasing its precision via some tolerance parameters (see comment:1 for details).
Besides, this ticket makes the odeint
solver (introduced in #25936) the default one, instead of rk4_maxima
, since the former is much faster than
the latter, as already noticed in #25936. The method name 'ode_int'
has been changed to 'odeint'
, in order to match with desolve_odeint()
and the SciPy name (for backward compatibility, method='ode_int'
is still accepted though).
Moreover, this ticket performs some improvement of the documentation of integrated curves.
Attachments (2)
Change History (13)
Changed 3 years ago by
Attachment:  test_geod_bad.png added 

Changed 3 years ago by
Attachment:  test_geod.png added 

comment:2 Changed 3 years ago by
Branch:  → public/manifolds/geod_integrator_options 

Commit:  → 511e26d6b637c8fb85c0112c8a5185d24c105672 
Description:  modified (diff) 
New commits:
7b44ef4  Add control parameters to IntegratedCurve.solve()

e31f6cb  Expose methods of scipy.integrate.ode in IntegratedCurve.solve()

0cc2077  Change default method in IntegratedCurve.solve() from rk4_maxima to odeint

bd1392b  Improve documentation of integrated curves

511e26d  More improvements in the documentation of integrated curves

comment:3 Changed 3 years ago by
Cc:  ghFlorentinJ Karim VanAelst added 

Status:  new → needs_review 
comment:4 Changed 3 years ago by
Commit:  511e26d6b637c8fb85c0112c8a5185d24c105672 → 2e918f44c02dcf185126ae3865c6769f446ceddc 

Branch pushed to git repo; I updated commit sha1. New commits:
2e918f4  Correct some doctests after the change of default method in IntegratedCurve.solve()

comment:5 Changed 3 years ago by
Commit:  2e918f44c02dcf185126ae3865c6769f446ceddc → aca08a898e22586d2d9fcb4f29e830cc788f165f 

Branch pushed to git repo; I updated commit sha1. New commits:
aca08a8  Merge branch public/manifolds/geod_integrator_options into Sage 9.0.beta5

comment:6 Changed 3 years ago by
The above commit fixes a merge conflict with Sage 9.0.beta5, due to #28669.
comment:7 followup: 9 Changed 3 years ago by
Reviewers:  → ghFlorentinJ 

Status:  needs_review → positive_review 
Seems perfect to me. Many thanks Éric.
(The failed plugins are only there because of the merge with #28669, no issue here)
comment:8 Changed 3 years ago by
Reviewers:  ghFlorentinJ → Florentin Jaffredo 

comment:9 Changed 3 years ago by
Replying to ghFlorentinJ:
Seems perfect to me. Many thanks Éric.
Thank you for the review!
(The failed plugins are only there because of the merge with #28669, no issue here)
I think it's rather an issue with this patchbot: there are many errors and all reports say that they happen "on 0 nonempty lines". The "0" is highly suspicious here...
comment:10 Changed 3 years ago by
The diff for the failed plugin only contain 2 lines:
 9.0.beta5 +++ 9.0.beta5 + #28707
Isn't that the reason it failed ?
comment:11 Changed 3 years ago by
Branch:  public/manifolds/geod_integrator_options → aca08a898e22586d2d9fcb4f29e830cc788f165f 

Resolution:  → fixed 
Status:  positive_review → closed 
In Sage 8.9, when integrating a null geodesic in Schwarzschild spacetime with an impact parameter close to critical (all details are in this notebook), the solution obtained with
solve(method='ode_int')
(green in the figure below) differs significantly from those obtained by means of other solvers, for instance therkf45
solver (red in the figure below):The geodesic arises from the upper right ((x,y) approx. (10, 5)), circles twice around the black hole (the grey area) and then departs away to the lower left (
rkf45
"exact" solution) or to the lower right (ode_int
bad solution). This issue is due to the default values of 1.49012e8 for the parametersrtol
andatol
of scipy.integrate.odeint being not sufficient in the present case. Now, with the Sage 8.9 implementation ofIntegratedCurve.solve()
, there is no way to change the values ofrtol
andatol
. Thanks to the new argument**control_param
introduced here this becomes possible. Morevoer, the default values ofrtol
andatol
are changed to 1.e10. Accordingly, the curves obtained from theodeint
andrkf45
become almost undistinguishable:See the notebook for a detailed study of the agreement between
odeint
and the other solvers when using the code of this ticket branch.