Skip to content
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

ENH: Implement DOP853 ODE integrator #10173

Merged
merged 20 commits into from Aug 15, 2019
Merged

ENH: Implement DOP853 ODE integrator #10173

merged 20 commits into from Aug 15, 2019

Conversation

@nmayorov
Copy link
Contributor

nmayorov commented May 12, 2019

I tried to figure out DOP853 algorithm in order to finish #9290

Eventually I decided to make a RungeKutta subclass implementation trying to make it more straightforward.

The main features of DOP853 compared to RK45 and RK23:

  1. There are two error estimators. Both estimators are implemented in the same way --- by a linear combination of computed RK stages. Then the final error is computed by a certain combination. This requires very little modification of the original code, two sets of coefficients needs to be defined.
  2. Three extra stages are required to compute the interpolant. This again is relatively straightforward and done inside _dense_output_impl.
  3. The suggested interpolant is not expressed in the polynomial basis. Of course it can be expressed as such, but to use the original coefficients we need to use the suggested evaluation scheme. This is implemented inside Dop853DenseOutput.
  4. Enhanced step size control is suggested with parameter beta. I haven't implemented this yet and put it as a low priority.

Comparison of RK23, RK45 and DOP853 on a celestial mechanics problem:

import numpy as np
from scipy.integrate import solve_ivp


mu = 0.0122771
nu = 1 - mu

t0 = 0
tf = 17.0652165601579625588917206249
y0 = [0.994, 0, 0, -2.00158510637908252240537862224]

rtol=1e-10
atol=1e-14


def fun(x, y):
    D1 = ((y[0] + mu)**2 + y[1]**2) ** 1.5
    D2 = ((y[0] - nu)**2 + y[1]**2) ** 1.5

    return np.array([
        y[2],
        y[3],
        y[0] + 2 * y[3] - nu * (y[0] + mu) / D1 - mu * (y[0] - nu) / D2,
        y[1] - 2 * y[2] - nu * y[1] / D1 - mu * y[1] / D2
    ])


np.set_printoptions(10)
for method in ['RK23', 'RK45', 'DOP853']:
    res = solve_ivp(fun, [t0, tf], y0, method=method, rtol=rtol, atol=atol)
    print(method)
    print("Solution at tf: {}".format(res.y[:, -1]))
    print("Number of function evaluations: {}".format(res.nfev))
    print("--------------------------------------------------")

The result:

RK23
Solution at tf: [ 0.9929909741 -0.0024338004 -0.4547767983 -2.0295309619]
Number of function evaluations: 99674
--------------------------------------------------
RK45
Solution at tf: [ 0.9929909646 -0.0024338146 -0.4547799205 -2.0295307335]
Number of function evaluations: 7622
--------------------------------------------------
DOP853
Solution at tf: [ 0.9929909675 -0.0024338098 -0.4547789079 -2.0295308438]
Number of function evaluations: 4046
--------------------------------------------------

Dense output illustration with rtol=1e-8 (same atol):
traj

Overall it seems to work correctly and perform more efficient with low tolerances.

However currently it does ~20% more function evaluations compared to the Fortran code. It must have something to do with error estimation and step size control / step rejection. At least I want to determine the reason for it. But I certainly don't plan to make it behave exactly as the Fortran code --- hard to achieve, creates strange code and doesn't give any real benefits.


@Juanlu001 @nikita-astronaut

I hope you guys are still interested in making it part of scipy and don't mind that I took a different approach to the implementation. Your code was certainly useful for me.

@Juanlu001

This comment has been minimized.

Copy link
Contributor

Juanlu001 commented May 13, 2019

Hi @nmayorov, thanks a lot for taking this over! I'm so happy that it's on good hands :) And glad my original attempt was helpful!

Regarding the number of function evaluations:

At least I want to determine the reason for it. But I certainly don't plan to make it behave exactly as the Fortran code --- hard to achieve, creates strange code and doesn't give any real benefits.

I was thinking of trying a side-by-side debugging, but perhaps that's not the best strategy if the code is not going to be exactly the same. Still, I will take a look at it this week.

On the other hand, we will try this method in poliastro and report here how it went.

@Juanlu001

This comment has been minimized.

Copy link
Contributor

Juanlu001 commented May 13, 2019

Pinging @helgee as well in case he wants to take a look.

@jorgepiloto

This comment has been minimized.

Copy link

jorgepiloto commented May 16, 2019

@nmayorov This is amazing! Thank you for this contribution! As @Juanlu001 said we will be testing this in poliastro 🚀

I will take a closer look to your implementation for DOP853 and study it in depth since we are also working on Verner78 method. Some time ago I opened an issue #9698 because I got stacked with the implementation of this method from Runge-Kutta class. This is the official page for Verner's Methods, from which we took the coefficients.

Again, thank you for this pull request! 😄

@nmayorov

This comment has been minimized.

Copy link
Contributor Author

nmayorov commented May 16, 2019

Hey guys! Good news --- I have figured out what the difference between this code and Fortran code is and will make an update soon.

It comes down to 2 main things:

  1. Error estimation is done somewhat different to what I implemented. The denominator (err5**2 + err3**2 * 0.01) ** 0.5 is computed not component-wise but as a scalar using norms of the corresponding errors (with scale involved).
  2. After a step was rejected the original code does not allow for the step size to be increased after a single successful step.

After I implemented these changes I got the nfev higher by exactly the number of rejected steps. Why is that? Because I always compute f_new = fun(t_new, y_new) (for generality), but DOP853 requires this only when the step was accepted (i.e. the error estimation does not require f_new).

My suggestion is not to change this place in order not to complicate the code.

@nmayorov

This comment has been minimized.

Copy link
Contributor Author

nmayorov commented May 16, 2019

The updated result on the same problem:

RK23
Solution at tf: [ 0.9929909729 -0.0024338022 -0.4547771854 -2.0295309264]
Number of function evaluations: 100658
--------------------------------------------------
RK45
Solution at tf: [ 0.9929909649 -0.0024338142 -0.4547798365 -2.0295307436]
Number of function evaluations: 7094
--------------------------------------------------
DOP853
Solution at tf: [ 0.9929909685 -0.0024338082 -0.4547785701 -2.0295308818]
Number of function evaluations: 3506
--------------------------------------------------

The remaining difference of the implementation:

  1. Minimum and maximum values for the step size change factors. I don't think these are terribly important, probably we can keep the current values.
  2. Initial step size selection. It seems to make a bigger difference, I will study what is implemented in DOP853 Fortran code and perhaps transfer this implementation to Python. I would prefer to have a single implementation of the step size selection (at least for RK methods).
@nmayorov

This comment has been minimized.

Copy link
Contributor Author

nmayorov commented May 17, 2019

Re-evaluated the situation quite a bit. I did compare the results to the MATLAB version posted on the same page with Fortran code.

And considering MATLAB version --- the codes work similar, the nfev is higher by the number of rejected steps which I explained above. However, the Fortran code (called through scipy) still works a bit different. But this code is painful to study in microscopic details.

As for the initial step selection --- original code uses norms without division by the number of components, which I believe was some arbitrary decision by the authors (perhaps small "mistake" on their part, at least it seems to contradict the description from the book). But I will think about it some more.


My current perspective is to keep the current version of the code as a legit implementation of DOP853. It looks conceptually correct, has right properties and works very similar to the MATLAB version and quite similar to the original Fortran version.

Please report the result of your evaluation within poliastro.

@nikita-astronaut

This comment has been minimized.

Copy link

nikita-astronaut commented May 17, 2019

Hi @nmayorov ,

thanks a lot for pushing it through!

@Juanlu001

This comment has been minimized.

Copy link
Contributor

Juanlu001 commented May 18, 2019

All our tests passed, and I did some quick tests and I'd say that this implementation is noticeably faster than ours, although I didn't measure it. Looking forward to see it merged, too bad we'll have to wait for SciPy 1.4 :)

@nmayorov nmayorov force-pushed the nmayorov:dop853 branch 3 times, most recently from 95e626b to 6eb294d Jun 30, 2019
@nmayorov nmayorov force-pushed the nmayorov:dop853 branch 2 times, most recently from 6fe7358 to 5d982ef Jul 28, 2019
@nmayorov nmayorov force-pushed the nmayorov:dop853 branch from 5d982ef to b7e16d3 Aug 13, 2019
@nmayorov

This comment has been minimized.

Copy link
Contributor Author

nmayorov commented Aug 13, 2019

@rgommers

Not sure how to proceed with it. Could you suggest something please? I think we can merge it.

@rgommers rgommers added this to the 1.4.0 milestone Aug 13, 2019
@rgommers

This comment has been minimized.

Copy link
Member

rgommers commented Aug 13, 2019

It looks like @Juanlu001 and @nikita-astronaut are happy with this, so I think we should get this merged indeed. Let's see if CI is completely happy. I'm not much of an expert, but happy to have a look over this PR this week and merge it (if anyone wants to beat me to it, please feel free!).

@nmayorov

This comment has been minimized.

Copy link
Contributor Author

nmayorov commented Aug 14, 2019

@rgommers sounds great, thanks.

Copy link
Member

rgommers left a comment

Looks great, merging. Thanks @nmayorov and everyone else who pitched in!

@rgommers rgommers merged commit b1bcb12 into scipy:master Aug 15, 2019
9 checks passed
9 checks passed
ci/circleci: build_docs Your tests passed on CircleCI!
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
scipy.scipy Build #20190814.15 succeeded
Details
scipy.scipy (Linux_Python_36_32bit_full) Linux_Python_36_32bit_full succeeded
Details
scipy.scipy (Windows Python35-64bit-full) Windows Python35-64bit-full succeeded
Details
scipy.scipy (Windows Python36-32bit-full) Windows Python36-32bit-full succeeded
Details
scipy.scipy (Windows Python36-64bit-full) Windows Python36-64bit-full succeeded
Details
scipy.scipy (Windows Python37-64bit-full) Windows Python37-64bit-full succeeded
Details
@rgommers

This comment has been minimized.

Copy link
Member

rgommers commented Aug 15, 2019

@nmayorov if you'd like to add a release note in https://github.com/scipy/scipy/wiki/Release-note-entries-for-SciPy-1.4.0, that would be useful

@Juanlu001

This comment has been minimized.

Copy link
Contributor

Juanlu001 commented Aug 15, 2019

Thanks @nmayorov and everyone! We're so happy to see this in SciPy ♥️

@nmayorov

This comment has been minimized.

Copy link
Contributor Author

nmayorov commented Aug 15, 2019

@rgommers updated release notes about DOP853.

@Juanlu001 happy to be useful!

@nmayorov nmayorov deleted the nmayorov:dop853 branch Oct 18, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Linked issues

Successfully merging this pull request may close these issues.

None yet

7 participants
You can’t perform that action at this time.