NumA/root_finding/Newton.h

0001 #ifndef NEWTON_H
0002 #define NEWTON_H
0003
0004 /**
0005  * @file Newton.h
0006  * @author Bryan BERTHOU (SPhN / CEA Saclay)
0007  * @date 29 January 2016
0008  * @version 1.0
0009  */
0010
0011 #include <ElementaryUtils/logger/CustomException.h>
0012 #include <stddef.h>
0013 #include <cmath>
0014 #include <iostream>
0015 #include <vector>
0016
0017 namespace NumA {
0018
0019 /**
0020  * @class Newton
0021  *
0022  * @brief The code supports on the example of the python library : scipy.
0023  * Find a zero using the Newton-Raphson or secant method.
0024  *
0025  * Find a zero of the function `func` given a nearby starting point `x0`.
0026  * The Newton-Raphson method is used if the derivative `fprime` of `func`
0027  * is provided, otherwise the secant method is used.
0028  */
0029 class Newton {
0030 public:
0031     //TODO add relative tolerance ; return status message when max iteration reached and no convergence
0032     template<class OBJ, typename FUNC, typename DERIV>
0033     double solve(OBJ* object, FUNC function,
0034             std::vector<double> &funcParameters, double x0,
0035             DERIV derivative = 0, const double absTolerance = 1.48e-8,
0036             const size_t maxIteration = 50) {
0037
0038         if (absTolerance <= 0.) {
0039             throw ElemUtils::CustomException("Newton", __func__,
0040                     "Tolerance too small");
0041         }
0042         if (maxIteration < 1) {
0043             throw ElemUtils::CustomException("Newton", __func__,
0044                     "Number of maximum iteration must be greater than 0");
0045         }
0046
0047         std::vector<double> variable(1, 0.);
0048         double x = 0.;
0049         bool noConvergence = true;
0050
0051         if (derivative != 0) {
0052             // Newton-Rapheson method
0053             double result_deriv = 0., result = 0.;
0054
0055             for (size_t i = 0; (i < maxIteration) && noConvergence; i++) {
0056                 result_deriv = ((*object).*derivative)(x0, funcParameters);
0057                 if (result_deriv == 0.) {
0058                     std::cerr
0059                             << "[WARNING] (Newton::solve) Derivative was zero."
0060                             << std::endl;
0061                     return x0;
0062                 }
0063
0064                 result = ((*object).*function)(x0, funcParameters);
0065
0066                 // Newton step
0067                 x = x0 - result / result_deriv;
0068                 noConvergence = fabs(x - x0) > absTolerance;
0069
0070                 x0 = x;
0071             }
0072
0073             return x;
0074
0075         } else {
0076             // Secant method
0077             double x1 = 0.;
0078
0079             if (x0 >= 0.) {
0080                 x1 = x0 * (1 + 1e-4) + 1e-4;
0081             } else {
0082                 x1 = x0 * (1 + 1e-4) - 1e-4;
0083             }
0084
0085             double v0 = ((*object).*function)(x0, funcParameters);
0086             double v1 = ((*object).*function)(x1, funcParameters);
0087
0088             for (size_t i = 0; (i < maxIteration) && noConvergence; i++) {
0089                 if (v1 == v0) {
0090
0091                     if (x1 != x0) {
0092                         std::cerr << "[WARNING] (Newton::solve) Tolerance of "
0093                                 << (x1 - x0) << "reached" << std::endl;
0094                     }
0095
0096                     return (x1 + x0) / 2.0;
0097                 } else {
0098                     x = x1 - v1 * (x1 - x0) / (v1 - v0);
0099                 }
0100
0101                 noConvergence = fabs(x - x1) > absTolerance;
0102
0103                 x0 = x1;
0104                 v0 = v1;
0105                 x1 = x;
0106                 v1 = ((*object).*function)(x1, funcParameters);
0107             }
0108
0109             return x;
0110         }
0111     }
0112 };
0113
0114 } /* namespace NumA */
0115
0116 #endif /* NEWTON_H */