NumA/root_finding/NewtonMD.h

0001 /**
0002  * @file NewtonMD.h
0003  * @author Nabil CHOUIKA (SPhN / CEA Saclay)
0004  * @date Feb 5 2016
0005  * @version 1.0
0006  */
0007
0008 #ifndef NEWTONMD_H_
0009 #define NEWTONMD_H_
0010
0011 #include <stddef.h>
0012
0013 #include "../linear_algebra/LinAlgUtils.h"
0014 #include "../linear_algebra/matrix/MatrixD.h"
0015 #include "../linear_algebra/vector/VectorD.h"
0016
0017 namespace NumA {
0018
0019 /**
0020  * @class NewtonMD
0021  *
0022  * @brief Newton method generalized to any dimension N:
0023  * Solves the system G(X) = 0,
0024  * where G is a regular but non-linear application \f$ R^N \rightarrow R^N \f$ .
0025  *
0026  * The class can be used either through templates with the solve method
0027  * providing G and the Jacobian matrix functions pointers, or manually
0028  * through instantiating an object NewtonMD and iterating the method iterate
0029  * by giving G(X) and the Jacobian matrix at each step, to avoid using functions.
0030  *
0031  * For linear G (matrices), use the LinearSystem class.
0032  */
0033 class NewtonMD {
0034 public:
0035     NewtonMD(int N = 1, LinAlgUtils::Method linMethod = LinAlgUtils::ColQR);
0036     virtual ~NewtonMD();
0037
0038     int getN() const;
0039     void setN(int n);
0040     LinAlgUtils::Method getLinMethod() const;
0041     void setLinMethod(LinAlgUtils::Method linMethod);
0042
0043     template<class OBJ, typename FUNC, typename JACOB>
0044     VectorD solve(OBJ* object, FUNC G,
0045             const std::vector<double> &funcParameters, VectorD X0,
0046             JACOB J_G = 0, const double absTolerance = 1.e-4,
0047             const double relTolerance = 1.e-4, const size_t maxIteration = 50) {
0048
0049         bool noConvergence = true;
0050         double absDiff = 0, relDiff = 0.;
0051         VectorD X(X0.size());
0052
0053         for (unsigned int n = 0; n < maxIteration && noConvergence; n++) {
0054
0055             VectorD G_X0 = ((*object).*G)(X0, funcParameters);
0056             MatrixD J_G_X0 = ((*object).*J_G)(X0, funcParameters);
0057
0058             X = iterate(X0, G_X0, J_G_X0);
0059
0060             absDiff = (X - X0).norm();
0061             relDiff = absDiff / (X0.norm() + 1.e-16);
0062             noConvergence = absDiff > absTolerance || relDiff > relTolerance;
0063
0064             X0 = X;
0065         }
0066
0067         return X;
0068     }
0069
0070     /** Produces a new iteration by solving the linear system
0071      * J^G(X0) * (X0 - X1) = G(X0)
0072      * and returns X1
0073      *
0074      * @param X0 : VectorD of size N
0075      * @param G_X0 : VectorD of size N
0076      * @param J_G_X0 : MatrixD of size NxN
0077      * @return X1 : New solution after the new iteration.
0078      */
0079     VectorD iterate(const VectorD & X0, const VectorD & G_X0,
0080             const MatrixD & J_G_X0);
0081
0082 private:
0083     int m_N; ///< Dimension of the problem
0084     LinAlgUtils::Method m_linMethod; ///< Method used for solving linear systems
0085 };
0086
0087 } /* namespace NumA */
0088
0089 #endif /* NEWTONMD_H_ */