diff --git a/include/GIntegral.hpp b/include/GIntegral.hpp index 8fe1fbc..6d11e08 100644 --- a/include/GIntegral.hpp +++ b/include/GIntegral.hpp @@ -4,10 +4,18 @@ * copyright (C) 2010-2011 by Jurgen Knodlseder * * ----------------------------------------------------------------------- * * * - * This program is free software; you can redistribute it and/or modify * - * it under the terms of the GNU General Public License as published by * - * the Free Software Foundation; either version 2 of the License, or * - * (at your option) any later version. * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * * * ***************************************************************************/ /** @@ -64,6 +72,7 @@ public: void integrand(GIntegrand* integrand) { m_integrand=integrand; } const GIntegrand* integrand(void) const { return m_integrand; } double romb(double a, double b, int k = 5); + double trapzd(double a, double b, int n = 1, double result = 0.0); std::string print(void) const; protected: @@ -72,7 +81,6 @@ protected: void copy_members(const GIntegral& integral); void free_members(void); double polint(double* xa, double* ya, int n, double x, double *dy); - double trapzd(double a, double b, int n, double result); // Protected data area GIntegrand* m_integrand; //!< Pointer to integrand diff --git a/src/numerics/GIntegral.cpp b/src/numerics/GIntegral.cpp index 0314907..963a444 100644 --- a/src/numerics/GIntegral.cpp +++ b/src/numerics/GIntegral.cpp @@ -196,63 +196,157 @@ double GIntegral::romb(double a, double b, int k) double result = 0.0; double ss = 0.0; double dss = 0.0; + + // Continue only if integration range is valid + if (b > a) { + + // Allocate temporal storage + double* s = new double[m_max_iter+2]; + double* h = new double[m_max_iter+2]; + + // Initialise step size + h[1] = 1.0; + s[0] = 0.0; + + // Iterative loop + for (m_iter = 1; m_iter <= m_max_iter; ++m_iter) { + + // Integration using Trapezoid rule + s[m_iter] = trapzd(a, b, m_iter, s[m_iter-1]); + + // Compile option: Check for NaN/Inf + #if defined(G_NAN_CHECK) + if (isnotanumber(s[m_iter]) || isinfinite(s[m_iter])) { + std::cout << "*** ERROR: GIntegral::romb"; + std::cout << "(a=" << a << ", b=" << b << ", k=" << k << "):"; + std::cout << " NaN/Inf encountered"; + std::cout << " (s[" << m_iter << "]=" << s[m_iter] << ")"; + std::cout << std::endl; + } + #endif + + // Starting from iteration k on, use polynomial interpolation + if (m_iter >= k) { + ss = polint(&h[m_iter-k], &s[m_iter-k], k, 0.0, &dss); + if (std::abs(dss) <= m_eps * std::abs(ss)) { + converged = true; + result = ss; + break; + } + } - // Allocate temporal storage - double* s = new double[m_max_iter+2]; - double* h = new double[m_max_iter+2]; - - // Initialise step size - h[1] = 1.0; - s[0] = 0.0; + // Reduce step size + h[m_iter+1]= 0.25 * h[m_iter]; - // Loop - for (m_iter = 1; m_iter <= m_max_iter; ++m_iter) { + } // endfor: iterative loop - // Integration using Trapezoid rule - s[m_iter] = trapzd(a, b, m_iter, s[m_iter-1]); + // Free temporal storage + delete [] s; + delete [] h; - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (std::isnan(s[m_iter]) || std::isinf(s[m_iter])) { - std::cout << "*** ERROR: GIntegral::romb"; - std::cout << "(a=" << a << ", b=" << b << ", k=" << k << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (s[" << m_iter << "]=" << s[m_iter] << ")"; - std::cout << std::endl; - } - #endif - - // Starting from iteration k on, use polynomial interpolation - if (m_iter >= k) { - ss = polint(&h[m_iter-k], &s[m_iter-k], k, 0.0, &dss); - if (std::abs(dss) <= m_eps * std::abs(ss)) { - converged = true; - result = ss; - break; + // Dump warning + if (!m_silent) { + if (!converged) { + std::cout << "*** WARNING: GIntegral::romb: "; + std::cout << "Integration did not converge "; + std::cout << "(iter=" << m_iter; + std::cout << ", result=" << ss; + std::cout << ", d=" << std::abs(dss); + std::cout << " > " << m_eps * std::abs(ss) << ")"; + std::cout << std::endl; + throw; } } + + } // endif: integration range was valid + + // Return result + return result; +} - // Reduce step size - h[m_iter+1]= 0.25 * h[m_iter]; +/***********************************************************************//** + * @brief Perform Trapezoidal integration + * + * @param[in] a Left integration boundary. + * @param[in] b Right integration boundary. + * @param[in] n Number of steps. + * @param[in] result Result from a previous trapezoidal integration step. + * + * The original Numerical Recipes function had result declared as a static + * variable, yet this led to some untrackable integration problems. For this + * reason, previous results are now passed using an argument. + * Result initialisation is done if n=1. + ***************************************************************************/ +double GIntegral::trapzd(double a, double b, int n, double result) +{ + // Handle case of identical boundaries + if (a == b) { + result = 0.0; } + + // ... otherwise use trapeziodal rule + else { + + // Case A: Only a single step is requested + if (n == 1) { + + // Evaluate integrand at boundaries + double y_a = m_integrand->eval(a); + double y_b = m_integrand->eval(b); + + // Compute result + result = 0.5*(b-a)*(m_integrand->eval(a) + m_integrand->eval(b)); + + } // endif: only a single step was requested - // Free temporal storage - delete [] s; - delete [] h; - - // Dump warning - if (!m_silent) { - if (!converged) { - std::cout << "*** WARNING: GIntegral::romb: "; - std::cout << "Integration did not converge "; - std::cout << "(iter=" << m_iter; - std::cout << ", result=" << ss; - std::cout << ", d=" << std::abs(dss); - std::cout << " > " << m_eps * std::abs(ss) << ")"; - std::cout << std::endl; + // Case B: More than a single step is requested + else { + // Compute step level 2^(n-1) + int it = 1; + for (int j = 1; j < n-1; ++j) + it <<= 1; + + // Verify that step level is valid + if (it == 0) { + std::cout << "*** ERROR: GIntegral::trapzd("; + std::cout << "a=" << a << ", b=" << b << ", n=" << n; + std::cout << ", result=" << result << "): it=" << it << ". "; + std::cout << "Looks like an overflow?"; + std::cout << std::endl; + } + + // Set step size + double tnm = double(it); + double del = (b-a)/tnm; + + // Verify that step is >0 + if (del == 0) { + std::cout << "*** ERROR: GIntegral::trapzd("; + std::cout << "a=" << a << ", b=" << b << ", n=" << n; + std::cout << ", result=" << result << "): del=" << del << ". "; + std::cout << "Step is too small to make sense."; + std::cout << std::endl; + } + + // Sum up values + double x = a + 0.5*del; + double sum = 0.0; + for (int j = 0; j < it; ++j, x+=del) { + + // Evaluate integrand + double y = m_integrand->eval(x); + + // Add integrand + sum += y; + + } // endfor: looped over steps + + // Set result + result = 0.5*(result + (b-a)*sum/tnm); } - } + + } // endelse: trapeziodal rule was applied // Return result return result; @@ -411,94 +505,6 @@ double GIntegral::polint(double* xa, double* ya, int n, double x, double *dy) } -/***********************************************************************//** - * @brief Perform Trapezoidal integration - * - * @param[in] a Left integration boundary. - * @param[in] b Right integration boundary. - * @param[in] n Number of steps. - * @param[in] result Result from a previous trapezoidal integration step. - * - * The original Numerical Recipes function had result declared as a static - * variable, yet this led to some untrackable integration problems. For this - * reason, previous results are now passed using an argument. - * Result initialisation is done if n=1. - ***************************************************************************/ -double GIntegral::trapzd(double a, double b, int n, double result) -{ - // Handle case of identical boundaries - if (a == b) { - result = 0.0; - } - - // ... otherwise use trapeziodal rule - else { - - // Case A: Only a single step is requested - if (n == 1) { - - // Evaluate integrand at boundaries - double y_a = m_integrand->eval(a); - double y_b = m_integrand->eval(b); - - // Compute result - result = 0.5*(b-a)*(m_integrand->eval(a) + m_integrand->eval(b)); - - } // endif: only a single step was requested - - // Case B: More than a single step is requested - else { - // Compute step level 2^(n-1) - int it = 1; - for (int j = 1; j < n-1; ++j) - it <<= 1; - - // Verify that step level is valid - if (it == 0) { - std::cout << "*** ERROR: GIntegral::trapzd("; - std::cout << "a=" << a << ", b=" << b << ", n=" << n; - std::cout << ", result=" << result << "): it=" << it << ". "; - std::cout << "Looks like an overflow?"; - std::cout << std::endl; - } - - // Set step size - double tnm = double(it); - double del = (b-a)/tnm; - - // Verify that step is >0 - if (it == 0) { - std::cout << "*** ERROR: GIntegral::trapzd("; - std::cout << "a=" << a << ", b=" << b << ", n=" << n; - std::cout << ", result=" << result << "): del=" << del << ". "; - std::cout << "Step is too small to make sense."; - std::cout << std::endl; - } - - // Sum up values - double x = a + 0.5*del; - double sum = 0.0; - for (int j = 0; j < it; ++j, x+=del) { - - // Evaluate integrand - double y = m_integrand->eval(x); - - // Add integrand - sum += y; - - } // endfor: looped over steps - - // Set result - result = 0.5*(result + (b-a)*sum/tnm); - } - - } // endelse: trapeziodal rule was applied - - // Return result - return result; -} - - /*========================================================================== = = = Friends =