search

headmyshoulder / odeint-v2

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

Altering error tolerance and step allotment Bulirsch-Stoer dense_output #202

Closed ceoconnor1996 closed 7 years ago

ceoconnor1996 commented 7 years ago

Dear Developers:

I have been developing a relatively straightforward simulation using the integrator bulirsch_stoer_dense_out. For moderate error tolerances (down to 1e-10 or so), the program runs without incident; but the results are somewhat less accurate than I would prefer (I monitor a conserved quantity to check this). However, when I attempt to enforce a more stringent error tolerance (1e-12, say), the integrator eventually throws an overflow error.

My intuition initially suggested that the integrator was hitting its maximal allowed number of steps (which I assume to be max_count, on line 302 of bulirsch_stoer_dense_out.hpp). When I increased this quantity, though, nothing appeared to change, even for augmentations by a factor of 1e+3 or more.

Am I approaching this the wrong way? Or might this error be the result of something subtler?

Thank you for your time! This package already has been immensely helpful.

headmyshoulder commented 7 years ago

Your approach sounds correct. It is difficult to debug this kind of problem from here. What error is thrown exactly? Is it an exception from odeint or a real overflow error from the standard library. Maybe you could add some logging to the integrator to see more

ceoconnor1996 commented 7 years ago

In hindsight, "overflow error" is the wrong term; my apologies. The error message is a simple "vector subscript out of range," citing something in the vector library.

My debugger takes me to line 525 in the dense output header, where indeed some elements of a matrix (m_diffs) are called. In addition, there are questionable values for the variables fac and sign: approximately 6.5e-107 and -9.3e+61, respectively. Those seem like nonsense numbers that come up when something hasn't been initialized correctly, but it's difficult for me to see how that could have happened.

mariomulansky commented 7 years ago

this sounds like either your rhs implementation is not correct, or you are not using the correct resizing scheme. You should validate the rhs function. Also the fact that you have bad convergence hints towards a bug in the rhs implementation.

mariomulansky commented 7 years ago

Please feel free to reopen if this is still an issue.