search

headmyshoulder / odeint-v2

odeint - solving ordinary differential equations in c++ v2
http://headmyshoulder.github.com/odeint-v2/
Other
341 stars 101 forks source link

Fixed incorrect constant in Dormand–Prince Runge–Kutta Butcher tableau. #235

Closed thecount2a closed 5 years ago

thecount2a commented 5 years ago

I already submitted this to boostorg/odeint but was told I should submit it here instead.

This constant somehow has the incorrect sign when compared to the actual Dormand-Prince Runge-Kutta butcher tableau. This causes the ODE solver to be a bit misguided as it tries to guess steps and estimate error and in the end, while it still works, it does not work as efficiently as it ought to. The entire butcher tableau was created in an internally balanced way but this coefficient being incorrect means that that the boost implementation is much less effective at solving correctly, when compared to another solver with correct coefficients, for the same solver tolerances.

Here are a few references showing that the 1/40 coefficient ought to be positive (look at the lower-right corner of the table):

https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method

Here is a paper that clearly has this coefficient as positive on page 1 equation 9:

http://depa.fquim.unam.mx/amyd/archivero/DormandPrince_19856.pdf

Here is a link to the original Dormand Prince paper published in 1980 which I was able to find freely available on sciencedirect.com (Page 5 of the PDF, table 2, again in the lower-right corner of the table):

https://www.sciencedirect.com/science/article/pii/0771050X80900133?via%3Dihub

I hope these references are sufficient to show that the coefficient is incorrect as it currently stands. We have used this ODE solver and found the solving to work much better with the coefficient corrected. This library has been very useful for us, thank you for all of the development work!

mariomulansky commented 5 years ago

Thanks a lot! I will propagate this back into boost.

mariomulansky commented 5 years ago

Tests were failing, i will have a closer look, but for now I need to revert.

mariomulansky commented 5 years ago

Ok I think I understand what's happening. The butcher tableau lists the values d_1...d_7, however in the code we need dc_n = c_n - d_n. Now for n=7 we have c_7=0, hence dc_7 = -d_7 = -1/40. So I believe the current implementation is correct.