From 12e667bcc3447b810a2afc8ba288d858d5d82b12 Mon Sep 17 00:00:00 2001 From: knodlseder <> Date: Wed, 2 Feb 2011 14:02:51 +0000 Subject: [PATCH] initial version --- include/GFunction.hpp | 58 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/numerics/GDerivative.cpp | 360 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/numerics/GException_numerics.cpp | 45 +++++++++++++++++++++++++++++++++++++++++++++ src/numerics/GFunction.cpp | 161 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 624 insertions(+), 0 deletions(-) create mode 100644 include/GFunction.hpp create mode 100644 src/numerics/GDerivative.cpp create mode 100644 src/numerics/GException_numerics.cpp create mode 100644 src/numerics/GFunction.cpp diff --git a/include/GFunction.hpp b/include/GFunction.hpp new file mode 100644 index 0000000..a2ef519 --- /dev/null +++ b/include/GFunction.hpp @@ -0,0 +1,58 @@ +/*************************************************************************** + * GFunction.hpp - Single parameter function abstract base class * + * ----------------------------------------------------------------------- * + * copyright (C) 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. * + * * + ***************************************************************************/ +/** + * @file GFunction.hpp + * @brief GFunction abstract virtual base class interface definition. + * @author J. Knodlseder + */ + +#ifndef GFUNCTION_HPP +#define GFUNCTION_HPP + +/* __ Includes ___________________________________________________________ */ + + +/***********************************************************************//** + * @class GFunction + * + * @brief GFunction class interface defintion. + * + * This class implements the abstract interface for a one parameter function. + * This function is for example used for integration or numerical computation + * of derivatives. This class has no members. The only pure virtual method + * that needs to be implemented by the derived class is the eval() method + * that provides function evaluation at a given value x, e.g. y=eval(x). + ***************************************************************************/ +class GFunction { + +public: + + // Constructors and destructors + GFunction(void); + GFunction(const GFunction& func); + virtual ~GFunction(void); + + // Operators + GFunction& operator= (const GFunction& func); + + // Methods + virtual double eval(double x) = 0; + +protected: + // Protected methods + void init_members(void); + void copy_members(const GFunction& func); + void free_members(void); +}; + +#endif /* GFUNCTION_HPP */ diff --git a/src/numerics/GDerivative.cpp b/src/numerics/GDerivative.cpp new file mode 100644 index 0000000..b15b344 --- /dev/null +++ b/src/numerics/GDerivative.cpp @@ -0,0 +1,360 @@ +/*************************************************************************** + * GDerivative.cpp - Derivative class * + * ----------------------------------------------------------------------- * + * copyright (C) 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. * + * * + ***************************************************************************/ +/** + * @file GDerivative.cpp + * @brief GDerivative class implementation. + * @author J. Knodlseder + */ + +/* __ Includes ___________________________________________________________ */ +#include // Definition of NULL +#include // for DBL_MAX +#include "GDerivative.hpp" +#include "GMatrix.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_RIDDER "GDerivative::ridder(double&, double&, double&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ +#define G_VALUE_SIGNAL_BAD_TOLERANCE 1 //!< Signal bad tolerance + +/* __ Debug definitions __________________________________________________ */ +#define G_VALUE_DEBUG 0 //!< Debug value method + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GDerivative::GDerivative(void) +{ + // Initialise private members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Function constructor + * + * @param[in] func Function pointer. + * + * Construct class instance from single parameter function. Note that the + * function is not copied but just a pointer to it is stored. The function + * should thus not be destroyed before computation of derivatives is + * finished. + ***************************************************************************/ +GDerivative::GDerivative(GFunction* func) +{ + // Initialise private members + init_members(); + + // Set function + this->function(func); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] dx Derivative. + ***************************************************************************/ +GDerivative::GDerivative(const GDerivative& dx) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(dx); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GDerivative::~GDerivative(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] dx Derivative. + ***************************************************************************/ +GDerivative& GDerivative::operator= (const GDerivative& dx) +{ + // Execute only if object is not identical + if (this != &dx) { + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(dx); + + } // endif: object was not identical + + // Return + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Returns derivative + * + * @param[in] x Function value. + * + * Computes derivative using the Ridders' method. The method starts with a + * step size h=0.02*|x| (with h >= 1e-6) and iteratively decreases the + * step size until the error becomes less than 1e-6. A maximum number of + * 5 iterations is done. + ***************************************************************************/ +double GDerivative::value(const double& x) +{ + // Set constants + const int max_iter = 5; + const double min_h = 1.0e-6; + const double max_err = 1.0e-6; + + // Initialise result + double result = 0.0; + + // Set initial step size + double h = 0.02*fabs(x); + if (h < min_h) + h = min_h; + + // Initialise tolerance + double err = 1.0; + + // Loop over Ridder's method until we have an acceptable tolerance + int iter = 0; + for (; iter < max_iter; ++iter) { + + // Compute derivative using Ridder's method + result = ridder(x, h, err); + + // Debug option: Show actual results + #if G_VALUE_DEBUG + std::cout << "GDerivative::value("; + std::cout << "iter=" << iter; + std::cout << ", x=" << x; + std::cout << ", dy/dx=" << result; + std::cout << ", h=" << h; + std::cout << ", err=" << err; + std::cout << ")" << std::endl; + #endif + + // If uncertainty is below tolerance then exit now + if (err < max_err) + break; + + // ... otherwise reduce the step size + h /= 5.0; + + } // endfor: tolerance matching loop + + // Compile option: signal if we exceed the tolerance + #if G_VALUE_SIGNAL_BAD_TOLERANCE + if (err >= max_err) { + std::cout << "WARNING: GDerivative::value("; + std::cout << "iter=" << iter; + std::cout << ", x=" << x; + std::cout << ", dy/dx=" << result; + std::cout << ", h=" << h; + std::cout << ", err=" << err; + std::cout << "): error exceeds tolerance of "; + std::cout << max_err << std::endl; + } + #endif + + + // Return derivative + return result; +} + + +/***********************************************************************//** + * @brief Returns derivative by Ridders' method + * + * @param[in] x Function value. + * @param[in] h Estimated initial step size. + * @param[out] err Estimated error in the derivative. + * + * @exception GException::invalid_argument + * Step size of 0 has been specified. + * + * Returns the derivative of a function at point x by Ridders's method of + * polynomial extrapolation. The value h is input as an estimated initial + * stepsize; it needs not to be small, but rather should be an increment + * in x over which the function changes substantially. + * + * This method has been adopted from Numerical Recipes, 3rd edition, page + * 321f. + ***************************************************************************/ +double GDerivative::ridder(const double& x, const double& h, double& err) +{ + // Set constants + const int ntab = 10; // Maximum table size + const double con = 1.4; // Stepsize decreases by con at each iteration + const double con2 = con * con; + const double big = DBL_MAX; + const double safe = 2.0; // Return when error is safe worse than best so far + + // Initialise result + double dx = 0.0; + err = big; + + // Continue only if we have a valid function + if (m_func != NULL) { + + // Allocate Neville table + GMatrix a(ntab, ntab); + + // Throw exception if step size is zero + if (h <= 0) + throw GException::invalid_argument(G_RIDDER, "Step size must be positive."); + + // Initialise adaptive step size + double hh = h; + + // Perform initial function evaluation + a(0,0) = (m_func->eval(x+hh) - m_func->eval(x-hh))/(2.0*hh); + + // Perform successive computations + for (int i = 1; i < ntab; ++i) { + + // Reduce step size + hh /= con; + + // Try new smaller step size + a(0,i) = (m_func->eval(x+hh) - m_func->eval(x-hh))/(2.0*hh); + + // Compute extrapolations of various orders, requiring no new + // function evaluations + double fac = con2; + for (int j = 1; j <= i; ++j) { + + // Compute extrapolation + a(j,i) = ( a(j-1,i) * fac - a(j-1,i-1) ) / (fac - 1.0); + fac *= con2; + + // If the error is reduced then save the improved answer + double err1 = fabs(a(j,i) - a(j-1,i)); + double err2 = fabs(a(j,i) - a(j-1,i-1)); + double errt = (err1 > err2) ? err1 : err2; + if (errt <= err) { + err = errt; + dx = a(j,i); + } + + } // endfor: looped over tableau + + // If higher order is worse by a significant factor, the quit early + if (fabs( a(i,i) - a(i-1,i-1) ) >= safe*err) + break; + + } // endfor: looped over computations + + } // endif: we had a valid function + + // Return derivative + return dx; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GDerivative::init_members(void) +{ + // Initialise members + m_func = NULL; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] dx Derivative. + ***************************************************************************/ +void GDerivative::copy_members(const GDerivative& dx) +{ + // Copy attributes + m_func = dx.m_func; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GDerivative::free_members(void) +{ + // Return + return; +} + + +/*========================================================================== + = = + = Friends = + = = + ==========================================================================*/ diff --git a/src/numerics/GException_numerics.cpp b/src/numerics/GException_numerics.cpp new file mode 100644 index 0000000..0858f41 --- /dev/null +++ b/src/numerics/GException_numerics.cpp @@ -0,0 +1,45 @@ +/*************************************************************************** + * GException_numerics.cpp - Numerics exceptions * + * ----------------------------------------------------------------------- * + * copyright (C) 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. * + * * + ***************************************************************************/ +/** + * @file GException_numerics.cpp + * @brief Numerics exception class implementations + * @author J. Knodlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "GException.hpp" +#include "GTools.hpp" + + +/***********************************************************************//** + * @brief + * + * @param[in] origin Name of method that has thrown the exception. + * @param[in] message Optional error message. + ***************************************************************************/ +/* +GException::exception::exception(std::string origin, std::string message) +{ + // Set origin of exception + m_origin = origin; + + // Set exception message + m_message = "Exception. "+ message; + + // Return + return; +} +*/ \ No newline at end of file diff --git a/src/numerics/GFunction.cpp b/src/numerics/GFunction.cpp new file mode 100644 index 0000000..41d4670 --- /dev/null +++ b/src/numerics/GFunction.cpp @@ -0,0 +1,161 @@ +/*************************************************************************** + * GFunction.cpp - Single parameter function abstract base class * + * ----------------------------------------------------------------------- * + * copyright (C) 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. * + * * + ***************************************************************************/ +/** + * @file GFunction.cpp + * @brief GFunction abstract virtual base class implementation. + * @author J. Knodlseder + */ + +/* __ Includes ___________________________________________________________ */ +#include "GFunction.hpp" + +/* __ Method name definitions ____________________________________________ */ + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GFunction::GFunction(void) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] func Function. + ***************************************************************************/ +GFunction::GFunction(const GFunction& func) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(func); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GFunction::~GFunction(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] func Function. + ***************************************************************************/ +GFunction& GFunction::operator= (const GFunction& func) +{ + // Execute only if object is not identical + if (this != &func) { + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(func); + + } // endif: object was not identical + + // Return + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GFunction::init_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] func Function. + ***************************************************************************/ +void GFunction::copy_members(const GFunction& func) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GFunction::free_members(void) +{ + // Return + return; +} + + +/*========================================================================== + = = + = Friends = + = = + ==========================================================================*/ -- libgit2 0.21.4