Code-Eli  0.3.6
newton_raphson_system_method.hpp
Go to the documentation of this file.
1 /*********************************************************************************
2 * Copyright (c) 2013 David D. Marshall <ddmarsha@calpoly.edu>
3 *
4 * All rights reserved. This program and the accompanying materials
5 * are made available under the terms of the Eclipse Public License v1.0
6 * which accompanies this distribution, and is available at
7 * http://www.eclipse.org/legal/epl-v10.html
8 *
9 * Contributors:
10 * David D. Marshall - initial code and implementation
11 ********************************************************************************/
12 
13 #ifndef eli_mutil_nls_newton_raphson_system_method_hpp
14 #define eli_mutil_nls_newton_raphson_system_method_hpp
15 
16 #include "eli/code_eli.hpp"
17 
19 
20 namespace eli
21 {
22  namespace mutil
23  {
24  namespace nls
25  {
26  template<typename data__, size_t N__, size_t NSOL__=1>
27  class newton_raphson_system_method : public iterative_system_root_base<data__, N__, NSOL__>
28  {
29  public:
30  static const int hit_constraint = 101;
31 
32  public:
34  : iterative_system_root_base<data__, N__, NSOL__>()
35  {
36  x0.setConstant(static_cast<data__>(0));
37  }
38 
40  : iterative_system_root_base<data__, N__, NSOL__>(nrm), x0(nrm.x0)
41  {
42  }
43 
45  {
46  }
47 
49  {
50  x0=xg;
51  }
52 
54  {
55  return x0;
56  }
57 
58  template<typename f__, typename g__>
60  {
61  typename iterative_system_root_base<data__, N__, NSOL__>::solution_matrix dx, x(x0), fx(fun(x0)), eval1, eval2;
63  data__ abs_tol_norm, rel_tol_norm;
65 
66  // calculate the function evaluated at the initial location
67  eval1=fx-f0;
68  abs_tol_norm=this->calculate_norm(eval1);
69  eval2=(fx-f0).array()/f0.array();
70  eval2.setConstant(1);
71  rel_tol_norm=this->calculate_norm(eval2);
72  if (this->test_converged(0, rel_tol_norm, abs_tol_norm))
73  {
74  root=x;
75  return this->converged;
76  }
77 
78  bool all_zero(false);
79  count=0;
80  while (!this->test_converged(count, rel_tol_norm, abs_tol_norm) && !all_zero)
81  {
82  // Don't have any easy (efficient) way of determining if matrix in invertible
83  // if (fpx==0)
84  // return iterative_root_base<data__>::no_root_found;
85 
86  dx=-fpx.lu().solve(eval1);
87  dx=calculate_delta_factor(x, dx);
88  x+=dx;
89  fx=fun(x);
90  fpx=fprime(x);
91  eval1=fx-f0;
92  abs_tol_norm=this->calculate_norm(eval1);
93  bool nonzero(false);
94  all_zero=true;
95  for (size_t i=0; i<N__; ++i)
96  {
97  // check if stuck and cannot move x anymore
98  if (std::abs(dx(i))>std::numeric_limits<data__>::epsilon())
99  {
100  all_zero=false;
101  }
102 
103  if (std::abs(f0(i))<=std::numeric_limits<data__>::epsilon())
104  eval2(i)=std::numeric_limits<data__>::epsilon();
105  else
106  {
107  nonzero=true;
108  eval2(i)=eval1(i)/f0(i);
109  }
110  }
111  if (nonzero)
112  rel_tol_norm=this->calculate_norm(eval2);
113  else
114  rel_tol_norm=static_cast<data__>(0);
115 
116  ++count;
117  }
118 
119  root=x;
120  if (this->max_iteration_reached(count))
121  return this->max_iteration; // could not converge
122  if (all_zero)
123  return this->hit_constraint; // constraints limited convergence
124 
125  return this->converged;
126  }
127 
128  private:
132  {
133  return dx;
134  }
135 
136  private:
138  };
139  }
140  }
141 }
142 #endif
Definition: math.hpp:20
bool max_iteration_reached(const iteration_type &it) const
Definition: iterative_root_base.hpp:261
const iterative_system_root_base< data__, N__, NSOL__ >::solution_matrix & get_initial_guess() const
Definition: newton_raphson_system_method.hpp:53
newton_raphson_system_method()
Definition: newton_raphson_system_method.hpp:33
Eigen::Matrix< data__, N__, NSOL__ > solution_matrix
Definition: iterative_system_root_base.hpp:28
Eigen::Matrix< data__, N__, N__ > jacobian_matrix
Definition: iterative_system_root_base.hpp:29
Definition: iterative_system_root_base.hpp:25
bool test_converged(const iteration_type &it, const tolerance_type &relv, const tolerance_type &absv) const
Definition: iterative_root_base.hpp:256
newton_raphson_system_method(const newton_raphson_system_method< data__, N__, NSOL__ > &nrm)
Definition: newton_raphson_system_method.hpp:39
void set_initial_guess(const typename iterative_system_root_base< data__, N__, NSOL__ >::solution_matrix &xg)
Definition: newton_raphson_system_method.hpp:48
data__ calculate_norm(const solution_matrix &mat) const
Definition: iterative_system_root_base.hpp:66
virtual iterative_system_root_base< data__, N__, NSOL__ >::solution_matrix calculate_delta_factor(const typename iterative_system_root_base< data__, N__, NSOL__ >::solution_matrix &, const typename iterative_system_root_base< data__, N__, NSOL__ >::solution_matrix &dx) const
Definition: newton_raphson_system_method.hpp:130
~newton_raphson_system_method()
Definition: newton_raphson_system_method.hpp:44
static const int hit_constraint
Definition: newton_raphson_system_method.hpp:30
iterative_system_root_base< data__, N__, NSOL__ >::solution_matrix x0
Definition: newton_raphson_system_method.hpp:137
static const int max_iteration
Definition: iterative_root_base.hpp:155
max_iteration_type::data_type iteration_type
Definition: iterative_root_base.hpp:161
Definition: newton_raphson_system_method.hpp:27
static const int converged
Definition: iterative_root_base.hpp:154
int find_root(typename iterative_system_root_base< data__, N__, NSOL__ >::solution_matrix &root, const f__ &fun, const g__ &fprime, const typename iterative_system_root_base< data__, N__, NSOL__ >::solution_matrix &f0) const
Definition: newton_raphson_system_method.hpp:59