From d78bdaf6e878010b1a4773241718725eabdadecc Mon Sep 17 00:00:00 2001
From: Norihiro Watanabe <norihiro.watanabe@ufz.de>
Date: Fri, 7 Mar 2014 14:39:11 +0100
Subject: [PATCH] add NewtonRaphson
---
MathLib/Nonlinear/NewtonRaphson-impl.h | 141 +++++++++++++++++++++++++
MathLib/Nonlinear/NewtonRaphson.h | 123 +++++++++++++++++++++
2 files changed, 264 insertions(+)
create mode 100644 MathLib/Nonlinear/NewtonRaphson-impl.h
create mode 100644 MathLib/Nonlinear/NewtonRaphson.h
diff --git a/MathLib/Nonlinear/NewtonRaphson-impl.h b/MathLib/Nonlinear/NewtonRaphson-impl.h
new file mode 100644
index 00000000000..4a466b3c657
--- /dev/null
+++ b/MathLib/Nonlinear/NewtonRaphson-impl.h
@@ -0,0 +1,141 @@
+/**
+ * \author Norihiro Watanabe
+ * \date 2012-06-25
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <limits>
+
+#include "logog/include/logog.hpp"
+
+namespace MathLib
+{
+
+namespace Nonlinear
+{
+
+NewtonRaphson::NewtonRaphson()
+: _normType(VecNormType::INFINITY_N),
+ _r_abs_tol(std::numeric_limits<double>::max()), _r_rel_tol(1e-6), _dx_rel_tol(.0),
+ _max_itr(25), _printErrors(true), _n_iterations(0), _r_abs_error(.0), _r_rel_error(.0), _dx_rel_error(.0)
+{
+}
+
+template<class F_RESIDUAL, class F_DX, class T_VALUE>
+bool NewtonRaphson::solve(F_RESIDUAL &f_residual, F_DX &f_dx, const T_VALUE &x0, T_VALUE &x_new)
+{
+ T_VALUE r(x0), dx(x0);
+ r = .0;
+ dx = .0;
+
+ const bool checkAbsResidualError = (_r_abs_tol<std::numeric_limits<double>::max());
+ const bool checkRelResidualError = (_r_rel_tol<std::numeric_limits<double>::max());
+ const bool checkRelDxError = (_dx_rel_tol>.0);
+ const bool needXNorm = (checkRelResidualError || checkRelDxError);
+
+ INFO("------------------------------------------------------------------");
+ INFO("*** NEWTON-RAPHSON nonlinear solver");
+ INFO("-> iteration started");
+
+ double x_norm = -1.;
+ double dx_norm = -1.;
+ double r_norm = -1.;
+ bool converged = false;
+
+ // evaluate initial residual
+ x_new = x0;
+ f_residual(x_new, r);
+ // check convergence
+ r_norm = norm(r, _normType);
+ dx_norm = norm(dx, _normType);
+ if (needXNorm)
+ x_norm = norm(x_new, _normType);
+ converged = ((r_norm < _r_abs_tol && r_norm < _r_rel_tol*x_norm)
+ || (checkRelDxError && dx_norm < _dx_rel_tol*x_norm));
+ if (_printErrors)
+ INFO("-> %d: ||r||=%1.3e, ||dx||=%1.3e, ||x||=%1.3e, ||dx||/||x||=%1.3e", 0, r_norm, dx_norm, x_norm, x_norm==0 ? dx_norm : dx_norm/x_norm);
+
+ std::size_t itr_cnt = 0;
+ if (!converged) {
+ for (itr_cnt=1; itr_cnt<_max_itr; itr_cnt++) {
+ // solve dx=-J^-1*r
+ f_dx(x_new, r, dx);
+ x_new += dx;
+ // evaluate residual
+ f_residual(x_new, r);
+#ifdef NDEBUG
+ printout(std::cout, itr_cnt, x_new, r, dx);
+#endif
+ // check convergence
+ r_norm = norm(r, _normType);
+ dx_norm = norm(dx, _normType);
+ if (needXNorm)
+ x_norm = norm(x_new, _normType);
+ converged = ((r_norm < _r_abs_tol && r_norm < _r_rel_tol*x_norm)
+ || (checkRelDxError && dx_norm < _dx_rel_tol*x_norm));
+ if (_printErrors)
+ INFO("-> %d: ||r||=%1.3e, ||dx||=%1.3e, ||x||=%1.3e, ||dx||/||x||=%1.3e", itr_cnt, r_norm, dx_norm, x_norm, x_norm==0 ? dx_norm : dx_norm/x_norm);
+
+ if (converged)
+ break;
+ }
+ }
+
+ INFO("-> iteration finished");
+ if (_max_itr==1) {
+ INFO("status : iteration not required");
+ } else {
+ INFO("status : %s", (converged ? "CONVERGED" : "***DIVERGED***"));
+ }
+ INFO("iteration : %d/%d", itr_cnt, _max_itr);
+ if (checkAbsResidualError)
+ INFO("abs. res. : %1.3e (tolerance=%1.3e)", r_norm, _r_abs_tol);
+ if (checkRelResidualError)
+ INFO("rel. res. : %1.3e (tolerance=%1.3e)", x_norm==0?r_norm:r_norm/x_norm, _r_rel_tol);
+ if (checkRelDxError)
+ INFO("dx : %1.3e (tolerance=%1.3e)", x_norm==0?dx_norm:dx_norm/x_norm, _dx_rel_tol);
+ INFO("norm type : %s", convertVecNormTypeToString(_normType).c_str());
+ INFO("------------------------------------------------------------------");
+
+ this->_n_iterations = itr_cnt;
+ this->_r_abs_error = r_norm;
+ this->_r_rel_error = (x_norm==0 ? r_norm : r_norm / x_norm);
+ this->_dx_rel_error = (x_norm==0 ? dx_norm : dx_norm / x_norm);
+
+ return converged;
+}
+
+
+#ifdef NDEBUG
+template<class T_VALUE>
+inline void NewtonRaphson::printout(std::ostream& os, std::size_t i, T_VALUE& x_new, T_VALUE& r, T_VALUE& dx)
+{
+ os << "-> " << i <<": x=(";
+ for (std::size_t i=0; i<x_new.size(); i++)
+ os << x_new[i] << " ";
+ os << "), r=(";
+ for (std::size_t i=0; i<dx.size(); i++)
+ os << r[i] << " ";
+ os << "), dx=(";
+ for (std::size_t i=0; i<dx.size(); i++)
+ os << dx[i] << " ";
+ os << ")\n";
+}
+
+// in case of double
+template<>
+inline void NewtonRaphson::printout(std::ostream& os, std::size_t i, double& x_new, double& r, double& dx)
+{
+ os << "-> " << i <<": x=" << x_new << ", r=" << r << ", dx=" << dx << "\n";
+}
+#endif
+
+}
+
+} //end
diff --git a/MathLib/Nonlinear/NewtonRaphson.h b/MathLib/Nonlinear/NewtonRaphson.h
new file mode 100644
index 00000000000..feb44f3c4c8
--- /dev/null
+++ b/MathLib/Nonlinear/NewtonRaphson.h
@@ -0,0 +1,123 @@
+/**
+ * \author Norihiro Watanabe
+ * \date 2012-06-25
+ *
+ * \copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef MATHLIB_NONLINEAR_NEWTONRAPHSON_H_
+#define MATHLIB_NONLINEAR_NEWTONRAPHSON_H_
+
+#include "LinAlg/VectorNorms.h"
+
+namespace MathLib
+{
+
+namespace Nonlinear
+{
+
+/**
+ * \brief Newton-Raphson method
+ */
+class NewtonRaphson
+{
+public:
+ /// Default constructor with norm type INFINITY_N, relative tolerance 1e-6
+ /// for residual, and the maximum number of iterations 25
+ NewtonRaphson();
+
+ /// set a vector norm type
+ void setNormType(VecNormType normType) {_normType = normType;}
+
+ /// set the maximum number of iterations
+ void setMaxIterations(std::size_t max_itr) {_max_itr = max_itr;}
+
+ /// set the absolute residual tolerance used by the stopping criteria
+ void setAbsResidualTolerance(double abs_tol) {_r_abs_tol = abs_tol;}
+
+ /// set the relative residual tolerance used by the stopping criteria
+ void setRelResidualTolerance(double rel_tol) {_r_rel_tol = rel_tol;}
+
+ /// set the relative solution increment tolerance used by the stopping criteria
+ void setRelDxTolerance(double rel_tol) {_dx_rel_tol = rel_tol;}
+
+ /// print errors during iterations
+ void printErrors(bool flag) {_printErrors = flag;}
+
+ /**
+ * solve a nonlinear problem
+ *
+ * \tparam F_RESIDUAL Function object returning a residual of an equation.
+ * The object should have an operator ()(\f$x_k\f$, \f$r_k}\f$).
+ * \f$x_k\f$ is current solution. \f$r_k\f$ is calculated residual and an output
+ * of this function.
+ * \tparam F_DX Function object returning a solution increment.
+ * The object should have an operator ()(\f$x_k\f$, \f$r_k, \f$\Delta x_{k+1}}\f$).
+ * \f$\Delta x_{k+1}}\f$ is a calculated solution increment and an output of
+ * this function.
+ * \tparam T_VALUE Data type of \f$x_k\f$
+ * Both scalar and vector types are available as far as the following conditions
+ * are satisfied
+ * - T_VALUE has a copy constructor (for non-primitive data types)
+ * - MathLib::norm(T_VALUE) exists
+ * \param f_residual Residual function object \f$r(x)\f$
+ * \param f_dx Solution increment function object \f$dx(x)=J^{-1}r\f$
+ * \param x0 Initial guess
+ * \param x_new Solution
+ * \return true if converged
+ */
+ template<class F_RESIDUAL, class F_DX, class T_VALUE>
+ bool solve(F_RESIDUAL &f_residual, F_DX &f_dx, const T_VALUE &x0, T_VALUE &x_new);
+
+ /// return the number of iterations
+ std::size_t getNIterations() const {return _n_iterations; }
+
+ /// return absolute error in the last iteration
+ double getAbsResidualError() const {return _r_abs_error; }
+
+ /// return relative error in the last iteration
+ double getRelResidualError() const {return _r_rel_error; }
+
+ /// return relative error in the last iteration
+ double getRelDxError() const {return _dx_rel_error; }
+
+private:
+#ifdef NDEBUG
+ /// print out for debugging
+ template<class T_VALUE>
+ inline void printout(std::ostream& os, std::size_t i, T_VALUE& x_new, T_VALUE& r, T_VALUE& dx);
+#endif
+
+ /// vector norm type used to evaluate errors
+ VecNormType _normType;
+ /// absolute tolerance for residual
+ double _r_abs_tol;
+ /// relative tolerance for residual
+ double _r_rel_tol;
+ /// relative tolerance for dx
+ double _dx_rel_tol;
+ /// the maximum allowed iteration number
+ std::size_t _max_itr;
+ /// print iteration errors or not
+ bool _printErrors;
+ /// the number of iterations in the last calculation
+ std::size_t _n_iterations;
+ /// absolute residual error in the last calculation
+ double _r_abs_error;
+ /// relative residual error in the last calculation
+ double _r_rel_error;
+ /// relative dx error in the last calculation
+ double _dx_rel_error;
+};
+
+} // Nonlinear
+} // MathLib
+
+#include "NewtonRaphson-impl.h"
+
+#endif // MATHLIB_NONLINEAR_NEWTONRAPHSON_H_
--
GitLab