From 2b438504e0fc68ea8224e88675247e555ee6a6e6 Mon Sep 17 00:00:00 2001 From: Tavian Barnes Date: Wed, 6 Oct 2010 23:56:11 -0400 Subject: Handle edge case in AdaptiveIntegrator::step(). --- src/vZ/Adaptive.hpp | 16 ++++++++++------ 1 file 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::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; } } -- cgit v1.2.3