summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorTavian Barnes <tavianator@gmail.com>2010-10-06 23:56:11 -0400
committerTavian Barnes <tavianator@gmail.com>2010-10-06 23:56:11 -0400
commit2b438504e0fc68ea8224e88675247e555ee6a6e6 (patch)
tree68f832949da9d0eb6201f0fab1a87ecd7679bf4e
parentae956bffbf796005e1f59a9c430cef69e21e13f6 (diff)
downloadvz-2b438504e0fc68ea8224e88675247e555ee6a6e6.tar.xz
Handle edge case in AdaptiveIntegrator::step().
-rw-r--r--src/vZ/Adaptive.hpp16
1 files changed, 10 insertions, 6 deletions
diff --git a/src/vZ/Adaptive.hpp b/src/vZ/Adaptive.hpp
index 256bda6..05fc439 100644
--- a/src/vZ/Adaptive.hpp
+++ b/src/vZ/Adaptive.hpp
@@ -76,30 +76,34 @@ namespace vZ
GenericAdaptiveIntegrator<Y>::step()
{
static const Scalar S = Scalar(19)/Scalar(20); // Arbitrary saftey factor
- Scalar newH;
+ Scalar newH = this->h();
Y y;
// Attempt the integration step in a loop
- bool rejected = true;
- while (rejected) {
+ while (true) {
KVector k = calculateK(m_a);
y = calculateY(k, m_b);
Y yStar = calculateY(k, m_bStar);
// Get an error estimate
+
using std::abs;
using std::pow;
Scalar delta = abs(y - yStar);
- Scalar scale = m_atol + std::max(abs(y), abs(this->y()))*m_rtol;
- newH = S*this->h()*pow(scale/delta, Scalar(1)/m_order);
+ if (delta == Scalar(0)) {
+ break;
+ }
+
+ Scalar scale = m_atol + std::max(abs(y), abs(this->y()))*m_rtol;
+ newH = S*this->h()*pow(scale/delta, Scalar(1)/m_order);
if (delta > scale) {
// Reject the step
this->h(newH);
++m_rejections;
} else {
- rejected = false;
+ break;
}
}