CppNoddy  0.90
Public Member Functions | List of all members
CppNoddy::Newton< _Type > Class Template Reference

A vector NEWTON iteration class. More...

#include <Newton.h>

Inheritance diagram for CppNoddy::Newton< _Type >:
CppNoddy::ArcLength_base< _Type > CppNoddy::Uncopyable

Public Member Functions

 Newton (Residual< _Type > *ptr_to_residual_object, unsigned max_steps=20, double tolerance=1.e-8, double derivative_step=1.e-8)
 Constructor. More...
 
void iterate (DenseVector< _Type > &x)
 The Newton iteration method. More...
 
void set_monitor_det (bool flag)
 If set then the system will monitor the sign of determinant of the Jacobian matrix and cause an ExceptionBifurcation when it changes sign. More...
 
void solve (DenseVector< _Type > &x)
 Solve the system for an initial guess by Newton iteration, this method is inherited from the ArcLength_base class and this simply points it to the iteration method. More...
 
void arclength_solve (DenseVector< _Type > &x)
 Arc-length solve the system. More...
 
- Public Member Functions inherited from CppNoddy::ArcLength_base< _Type >
 ArcLength_base ()
 
virtual ~ArcLength_base ()
 
void init_arc (DenseVector< _Type > x, _Type *p, const double &length, const double &max_length)
 Initialise the class ready for arc-length continuation. More...
 
double & ds ()
 Return a handle to the arclength step. More...
 
double & arcstep_multiplier ()
 Used to set the multiplication constant used when increasing or decreasing the arclength step. More...
 
bool & rescale_theta ()
 Handle to the RESCALE_THETA flag. More...
 
double & theta ()
 Set the arclength theta parameter. More...
 
double & desired_arc_proportion ()
 Handle to the desired proportion of the parameter to be used in the arc length solver. More...
 

Additional Inherited Members

- Protected Member Functions inherited from CppNoddy::ArcLength_base< _Type >
void update (const DenseVector< _Type > &x)
 A method called by arclength_solve and init_arc which stores the current converged state and parameter and hence computes the derivatives w.r.t the arc-length. More...
 
double arclength_residual (const DenseVector< _Type > &x) const
 The extra constraint that is to be used to replace the unknown arc length. More...
 
DenseVector< _Type > Jac_arclength_residual (DenseVector< _Type > &x) const
 The derivative of the arclength_residual function with respect to each of the state variables. More...
 
void update_theta (const DenseVector< _Type > &x)
 Automatically update the Keller THETA such that the proportion of the arclength obtained from the parameter is the desired value. More...
 
- Protected Attributes inherited from CppNoddy::ArcLength_base< _Type >
_Type * p_PARAM
 pointer to the parameter in arclength solves More...
 
DenseVector< _Type > LAST_X
 state variable at the last computed solution More...
 
DenseVector< _Type > X_DERIV_S
 derivative of the state variable w.r.t. arc length More...
 
_Type LAST_PARAM
 parameter value at the last computed solution More...
 
_Type PARAM_DERIV_S
 derivative of the parameter w.r.t arc length More...
 
double DS
 size of the arc length step More...
 
double MAX_DS
 maximum arc length step to be taken More...
 
double ARCSTEP_MULTIPLIER
 step change multiplier More...
 
bool INITIALISED
 for the arc-length solver - to show it has been initialised More...
 
double THETA
 the arclength theta More...
 

Detailed Description

template<typename _Type>
class CppNoddy::Newton< _Type >

A vector NEWTON iteration class.

This allows for Newton iteration to be performed for a vector function of a vector unknown. Use templates to allow double or complex.

Definition at line 26 of file Newton.h.

Constructor & Destructor Documentation

template<typename _Type >
CppNoddy::Newton< _Type >::Newton ( Residual< _Type > *  ptr_to_residual_object,
unsigned  max_steps = 20,
double  tolerance = 1.e-8,
double  derivative_step = 1.e-8 
)
explicit

Constructor.

Parameters
ptr_to_residual_objectA pointer to an inherited Residual object
max_stepsThe maximum number of iteration steps.
toleranceA tolerance used as a convergence criterion.
derivative_stepA step length used to compute derivatives.

Definition at line 20 of file Newton.cpp.

24  : ArcLength_base< _Type >(),
25  TOL( tolerance ),
26  MAX_STEPS( max_steps ),
27  p_RESIDUAL( ptr_to_residual_object )
28  {
29  p_RESIDUAL -> delta() = derivative_step;
30  DELTA = derivative_step;
31  MONITOR_DET = false;
32  LAST_DET_SIGN = 1;
33  }

Member Function Documentation

template<typename _Type >
void CppNoddy::Newton< _Type >::arclength_solve ( DenseVector< _Type > &  x)

Arc-length solve the system.

Before this can be called the arc_init method should have been called in order to ensure we know a solution and have derivatives w.r.t. the arc-length parameter.

Definition at line 107 of file Newton.cpp.

References CppNoddy::ArcLength_base< _Type >::arclength_residual(), CppNoddy::ArcLength_base< _Type >::ARCSTEP_MULTIPLIER, CppNoddy::Utility::dot(), CppNoddy::ArcLength_base< _Type >::DS, CppNoddy::LinearSystem_base::get_det_sign(), CppNoddy::DenseVector< _Type >::inf_norm(), CppNoddy::ArcLength_base< _Type >::INITIALISED, CppNoddy::ArcLength_base< _Type >::Jac_arclength_residual(), CppNoddy::ArcLength_base< _Type >::LAST_PARAM, CppNoddy::ArcLength_base< _Type >::LAST_X, CppNoddy::ArcLength_base< _Type >::MAX_DS, CppNoddy::ArcLength_base< _Type >::p_PARAM, CppNoddy::ArcLength_base< _Type >::PARAM_DERIV_S, CppNoddy::LinearSystem_base::set_monitor_det(), CppNoddy::DenseLinearSystem< _Type >::solve(), CppNoddy::ArcLength_base< _Type >::update(), CppNoddy::ArcLength_base< _Type >::X_DERIV_S, and CppNoddy::Example::z().

Referenced by main(), and CppNoddy::Newton< _Type >::solve().

108  {
109 #ifdef PARANOID
110  if ( !this -> INITIALISED )
111  {
112  std::string problem;
113  problem = "The Newton.arclength_solve method has been called, but \n";
114  problem += " you haven't called the arc_init method first. This means \n";
115  problem += " that a starting solution & derivatives thereof are unknown. \n";
116  problem += " Please initialise things appropriately! ";
117  throw ExceptionRuntime( problem );
118  }
119 #endif
120 #ifdef DEBUG
121  std::cout.precision( 6 );
122  std::cout << "[ DEBUG ] : Entering arclength_solve of the Newton class with\n";
123  std::cout << "[ DEBUG ] : a parameter of " << *( this -> p_PARAM ) << "\n";
124 #endif
125  // backup the state/system in case we fail
126  DenseVector<_Type> backup_state( x );
127  _Type backup_parameter( *( this -> p_PARAM ) );
128 
129  int det_sign( 1 );
130  // init some local vars
131  bool step_succeeded( false );
132  unsigned itn( 0 );
133  // make a guess at the next solution
134  *( this -> p_PARAM ) = this -> LAST_PARAM + this -> PARAM_DERIV_S * this -> DS;
135  x = this -> LAST_X + this -> X_DERIV_S * this -> DS;
136 
137  // the order of the system
138  unsigned N = this -> p_RESIDUAL -> get_order();
139  DenseMatrix<_Type> J( N, N, 0.0 );
140  DenseVector<_Type> R1( N, 0.0 );
141  DenseVector<_Type> R2( N, 0.0 );
142  DenseVector<_Type> dR_dp( N, 0.0 );
143  //
144  do
145  {
146  // update the residual object to the current guess
147  this -> p_RESIDUAL -> update( x );
148  // get the Jacobian of the system
149  J = this -> p_RESIDUAL -> jacobian();
150  R1 = this -> p_RESIDUAL -> residual();
151  // get the residual of the EXTRA arclength residual
152  double E1 = this -> arclength_residual( x );
153  // compute derivatives w.r.t the parameter
154  *( this -> p_PARAM ) += this -> DELTA;
155  this -> p_RESIDUAL -> residual_fn( x, R2 );
156  double E2 = this -> arclength_residual( x );
157  *( this -> p_PARAM ) -= this -> DELTA;
158  dR_dp = ( R2 - R1 ) / this -> DELTA;
159  _Type dE_dp = ( E2 - E1 ) / this -> DELTA;
160  // bordering algorithm
161  DenseVector<_Type> y( -R1 );
162  DenseVector<_Type> z( -dR_dp );
163 #ifdef LAPACK
164  DenseLinearSystem<_Type> sys1( &J, &y, "lapack" );
165 #else
166  DenseLinearSystem<_Type> sys1( &J, &y, "native" );
167 #endif
168  sys1.set_monitor_det( MONITOR_DET );
169  try
170  {
171  sys1.solve();
172  det_sign = sys1.get_det_sign();
173  }
174  catch ( ExceptionExternal )
175  {
176  break;
177  }
178  // the solve will have overwritten the LHS, so we need to replace it
179  J = this -> p_RESIDUAL -> jacobian();
180 #ifdef LAPACK
181  DenseLinearSystem<_Type> sys2( &J, &z, "lapack" );
182 #else
183  DenseLinearSystem<_Type> sys2( &J, &z, "native" );
184 #endif
185  try
186  {
187  sys2.solve();
188  }
189  catch ( ExceptionExternal )
190  {
191  break;
192  }
193  DenseVector<_Type> JacE( this -> Jac_arclength_residual( x ) );
194  _Type delta_p = - ( E1 + Utility::dot( JacE, y ) ) /
195  ( dE_dp + Utility::dot( JacE, z ) );
196  DenseVector<_Type> delta_x = y + z * delta_p;
197  double max_correction = std::max( delta_x.inf_norm(), std::abs( delta_p ) );
198  if ( max_correction < this -> TOL )
199  {
200  step_succeeded = true;
201  break;
202  }
203  // add the corrections to the state variables
204  x += delta_x;
205  *( this -> p_PARAM ) += delta_p;
206  ++itn;
207  if ( itn > MAX_STEPS )
208  {
209  step_succeeded = false;
210  break;
211  }
212  }
213  while ( true );
214 
215  // is this a successful step?
216  if ( !step_succeeded )
217  {
218 #ifdef DEBUG
219  std::cout << "[ DEBUG ] : REJECTING STEP \n";
220 #endif
221  // if not a successful step then restore things
222  x = backup_state;
223  *( this -> p_PARAM ) = backup_parameter;
224  // restore the residual object to its initial state
225  this -> p_RESIDUAL -> update( x );
226  // reduce our step length
227  this -> DS /= this -> ARCSTEP_MULTIPLIER;
228  }
229  else
230  {
231  // update the variables needed for arc-length continuation
232  this -> update( x );
233  if ( LAST_DET_SIGN * det_sign < 0 )
234  {
235  LAST_DET_SIGN = det_sign;
236  std::string problem;
237  problem = "[ INFO ] : Determinant monitor has changed signs in the Newton class.\n";
238  problem += "[ INFO ] : Bifurcation detected.\n";
239  throw ExceptionBifurcation( problem );
240  }
241  else
242  {
243  LAST_DET_SIGN = det_sign;
244  }
245 #ifdef DEBUG
246  std::cout << "[ DEBUG ] : Number of iterations = " << itn << "\n";
247  std::cout << "[ DEBUG ] : Parameter p = " << *( this -> p_PARAM )
248  << "; arclength DS = " << this -> DS << "\n";
249 #endif
250  if ( itn >= 7 )
251  {
252  // converging too slowly, so decrease DS
253  this -> DS /= this -> ARCSTEP_MULTIPLIER;
254 #ifdef DEBUG
255  std::cout << "[ DEBUG ] : I decreased DS to " << this -> DS << "\n";
256 #endif
257  }
258  if ( itn <= 2 )
259  {
260  if ( std::abs( this -> DS * this -> ARCSTEP_MULTIPLIER ) < this -> MAX_DS )
261  {
262  // converging too quickly, so increase DS
263  this -> DS *= this -> ARCSTEP_MULTIPLIER;
264 #ifdef DEBUG
265  std::cout << "[ DEBUG ] : I increased DS to " << this -> DS << "\n";
266 #endif
267  }
268  }
269  }
270  }
double z(const double &x)
Topography shape.
_Type dot(const DenseVector< _Type > &X, const DenseVector< _Type > &Y)
Templated dot product.
Definition: Utility.h:381
DenseVector< _Type > LAST_X
state variable at the last computed solution
double arclength_residual(const DenseVector< _Type > &x) const
The extra constraint that is to be used to replace the unknown arc length.
_Type PARAM_DERIV_S
derivative of the parameter w.r.t arc length
double DS
size of the arc length step
double MAX_DS
maximum arc length step to be taken
double ARCSTEP_MULTIPLIER
step change multiplier
_Type LAST_PARAM
parameter value at the last computed solution
_Type * p_PARAM
pointer to the parameter in arclength solves
void update(const DenseVector< _Type > &x)
A method called by arclength_solve and init_arc which stores the current converged state and paramete...
bool INITIALISED
for the arc-length solver - to show it has been initialised
DenseVector< _Type > X_DERIV_S
derivative of the state variable w.r.t. arc length
DenseVector< _Type > Jac_arclength_residual(DenseVector< _Type > &x) const
The derivative of the arclength_residual function with respect to each of the state variables...
template<typename _Type >
void CppNoddy::Newton< _Type >::iterate ( DenseVector< _Type > &  x)

The Newton iteration method.

Parameters
xAn initial guess vector and returns the solution via this too.

Definition at line 37 of file Newton.cpp.

References CppNoddy::LinearSystem_base::get_det_sign(), CppNoddy::LinearSystem_base::set_monitor_det(), CppNoddy::DenseVector< _Type >::size(), CppNoddy::DenseLinearSystem< _Type >::solve(), and CppNoddy::ArcLength_base< _Type >::update().

Referenced by main(), and CppNoddy::Newton< _Type >::solve().

38  {
39  // length of the vector
40  unsigned N = x.size();
41  // Jacobian
42  DenseMatrix<_Type> J( N, N, 0.0 );
43  // NVectors
44  DenseVector<_Type> oldFn( N, 0.0 ), newFn( N, 0.0 );
45  // linear system object - native because no LAPACK complex at the moment
46  DenseLinearSystem<_Type> system( &J, &oldFn, "native" );
47  system.set_monitor_det( MONITOR_DET );
48  // iteration counter
49  unsigned itn = 0;
50  do
51  {
52  // increment the counter
53  ++itn;
54 
55  // evaluate the residuals and jacobian at the current state
56  p_RESIDUAL -> update( x );
57  // get the residuals
58  oldFn = p_RESIDUAL -> residual();
59 #ifdef DEBUG
60  std::cout << " DEBUG: starting with |Residuals| = " << oldFn.inf_norm() << "\n";
61 #endif
62 
63  if ( ( std::abs( oldFn.inf_norm() ) < TOL ) || ( itn == MAX_STEPS ) )
64  {
65  break;
66  }
67  // retrieve the current jacobian
68  J = p_RESIDUAL -> jacobian();
69 
70  // linear solver
71  // \todo LU interface to LAPACK for complex matrices
72  system.solve();
73 
74 #ifdef DEBUG
75  std::cout << " DEBUG: Iteration number = " << itn << "\n";
76  std::cout << " DEBUG: |Newton correction| = " << oldFn.inf_norm() << "\n";
77 #endif
78 
79  // must *subtract* delta
80  x -= oldFn;
81  }
82  while ( true );
83 
84  // More the 'MAX_STEPS' iterations currently triggers a failure.
85  if ( itn == MAX_STEPS )
86  {
87  std::string problem;
88  problem = " The Newton.iterate method took too many iterations. \n";
89  problem += " At the moment, this is set as a failure. \n";
90  throw ExceptionItn( problem, itn, oldFn.inf_norm() );
91  }
92  LAST_DET_SIGN = system.get_det_sign();
93  if ( MONITOR_DET )
94  {
95  if ( system.get_det_sign() * LAST_DET_SIGN < 0 )
96  {
97  std::string problem;
98  problem = "[ INFO ] : Determinant monitor has changed signs in ODE_BVP.\n";
99  problem += "[ INFO ] : Bifurcation detected.\n";
100  throw ExceptionBifurcation( problem );
101  }
102  }
103  }
void update(const DenseVector< _Type > &x)
A method called by arclength_solve and init_arc which stores the current converged state and paramete...
template<typename _Type>
void CppNoddy::Newton< _Type >::set_monitor_det ( bool  flag)
inline

If set then the system will monitor the sign of determinant of the Jacobian matrix and cause an ExceptionBifurcation when it changes sign.

Parameters
flagThe value to be set.

Definition at line 49 of file Newton.h.

Referenced by main().

50  {
51  MONITOR_DET = flag;
52  }
template<typename _Type>
void CppNoddy::Newton< _Type >::solve ( DenseVector< _Type > &  x)
inlinevirtual

Solve the system for an initial guess by Newton iteration, this method is inherited from the ArcLength_base class and this simply points it to the iteration method.

Parameters
xAn initial guess vector, the solution overwrites it.

Implements CppNoddy::ArcLength_base< _Type >.

Definition at line 58 of file Newton.h.

References CppNoddy::Newton< _Type >::arclength_solve(), and CppNoddy::Newton< _Type >::iterate().

59  {
60  iterate( x );
61  }
void iterate(DenseVector< _Type > &x)
The Newton iteration method.
Definition: Newton.cpp:37

The documentation for this class was generated from the following files:

© 2012

R.E. Hewitt