CppNoddy  0.85
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Groups Pages
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.
void iterate (DenseVector< _Type > &x)
 The Newton iteration method.
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.
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.
void arclength_solve (DenseVector< _Type > &x)
 Arc-length solve the system.
- 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.
double & ds ()
 Return a handle to the arclength step.
double & arcstep_multiplier ()
 Used to set the multiplication constant used when increasing or decreasing the arclength step.
bool & rescale_theta ()
 Handle to the RESCALE_THETA flag.
double & theta ()
 Set the arclength theta parameter.
double & desired_arc_proportion ()
 Handle to the desired proportion of the parameter to be used in the arc length solver.

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.
double arclength_residual (const DenseVector< _Type > &x) const
 The extra constraint that is to be used to replace the unknown 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.
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.
- Protected Attributes inherited from CppNoddy::ArcLength_base< _Type >
_Type * p_PARAM
 pointer to the parameter in arclength solves
DenseVector< _Type > LAST_X
 state variable at the last computed solution
DenseVector< _Type > X_DERIV_S
 derivative of the state variable w.r.t. arc length
_Type LAST_PARAM
 parameter value at the last computed solution
_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
bool INITIALISED
 for the arc-length solver - to show it has been initialised
double THETA
 the arclength theta

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.

: ArcLength_base< _Type >(),
TOL( tolerance ),
MAX( max_steps ),
p_RESIDUAL( ptr_to_residual_object )
{
p_RESIDUAL -> delta() = derivative_step;
DELTA = derivative_step;
MONITOR_DET = false;
LAST_DET_SIGN = 1;
}

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::Utility::dot(), CppNoddy::LinearSystem_base::get_det_sign(), CppNoddy::DenseVector< _Type >::inf_norm(), CppNoddy::LinearSystem_base::set_monitor_det(), CppNoddy::DenseLinearSystem< _Type >::solve(), and CppNoddy::Example::z().

Referenced by main().

{
#ifdef PARANOID
if ( !this -> INITIALISED )
{
std::string problem;
problem = "The Newton.arclength_solve method has been called, but \n";
problem += " you haven't called the arc_init method first. This means \n";
problem += " that a starting solution & derivatives thereof are unknown. \n";
problem += " Please initialise things appropriately! ";
throw ExceptionRuntime( problem );
}
#endif
#ifdef DEBUG
std::cout.precision( 6 );
std::cout << "[ DEBUG ] : Entering arclength_solve of the Newton class with\n";
std::cout << "[ DEBUG ] : a parameter of " << *( this -> p_PARAM ) << "\n";
#endif
// backup the state/system in case we fail
DenseVector<_Type> backup_state( x );
_Type backup_parameter( *( this -> p_PARAM ) );
int det_sign( 1 );
// init some local vars
bool step_succeeded( false );
unsigned itn( 0 );
// make a guess at the next solution
*( this -> p_PARAM ) = this -> LAST_PARAM + this -> PARAM_DERIV_S * this -> DS;
x = this -> LAST_X + this -> X_DERIV_S * this -> DS;
// the order of the system
unsigned N = this -> p_RESIDUAL -> get_order();
DenseMatrix<_Type> J( N, N, 0.0 );
DenseVector<_Type> R1( N, 0.0 );
DenseVector<_Type> R2( N, 0.0 );
DenseVector<_Type> dR_dp( N, 0.0 );
//
do
{
// update the residual object to the current guess
this -> p_RESIDUAL -> update( x );
// get the Jacobian of the system
J = this -> p_RESIDUAL -> jacobian();
R1 = this -> p_RESIDUAL -> residual();
// get the residual of the EXTRA arclength residual
double E1 = this -> arclength_residual( x );
// compute derivatives w.r.t the parameter
*( this -> p_PARAM ) += this -> DELTA;
this -> p_RESIDUAL -> residual_fn( x, R2 );
double E2 = this -> arclength_residual( x );
*( this -> p_PARAM ) -= this -> DELTA;
dR_dp = ( R2 - R1 ) / this -> DELTA;
_Type dE_dp = ( E2 - E1 ) / this -> DELTA;
// bordering algorithm
DenseVector<_Type> y( -R1 );
DenseVector<_Type> z( -dR_dp );
#ifdef LAPACK
DenseLinearSystem<_Type> sys1( &J, &y, "lapack" );
#else
DenseLinearSystem<_Type> sys1( &J, &y, "native" );
#endif
sys1.set_monitor_det( MONITOR_DET );
try
{
sys1.solve();
det_sign = sys1.get_det_sign();
}
catch ( ExceptionExternal )
{
break;
}
// the solve will have overwritten the LHS, so we need to replace it
J = this -> p_RESIDUAL -> jacobian();
#ifdef LAPACK
DenseLinearSystem<_Type> sys2( &J, &z, "lapack" );
#else
DenseLinearSystem<_Type> sys2( &J, &z, "native" );
#endif
try
{
sys2.solve();
}
catch ( ExceptionExternal )
{
break;
}
DenseVector<_Type> JacE( this -> Jac_arclength_residual( x ) );
_Type delta_p = - ( E1 + Utility::dot( JacE, y ) ) /
( dE_dp + Utility::dot( JacE, z ) );
DenseVector<_Type> delta_x = y + z * delta_p;
double max_correction = std::max( delta_x.inf_norm(), std::abs( delta_p ) );
if ( max_correction < this -> TOL )
{
step_succeeded = true;
break;
}
// add the corrections to the state variables
x += delta_x;
*( this -> p_PARAM ) += delta_p;
++itn;
if ( itn > MAX )
{
step_succeeded = false;
break;
}
}
while ( true );
// is this a successful step?
if ( !step_succeeded )
{
#ifdef DEBUG
std::cout << "[ DEBUG ] : REJECTING STEP \n";
#endif
// if not a successful step then restore things
x = backup_state;
*( this -> p_PARAM ) = backup_parameter;
// restore the residual object to its initial state
this -> p_RESIDUAL -> update( x );
// reduce our step length
this -> DS /= this -> ARCSTEP_MULTIPLIER;
}
else
{
// update the variables needed for arc-length continuation
this -> update( x );
if ( LAST_DET_SIGN * det_sign < 0 )
{
LAST_DET_SIGN = det_sign;
std::string problem;
problem = "[ INFO ] : Determinant monitor has changed signs in the Newton class.\n";
problem += "[ INFO ] : Bifurcation detected.\n";
throw ExceptionBifurcation( problem );
}
else
{
LAST_DET_SIGN = det_sign;
}
#ifdef DEBUG
std::cout << "[ DEBUG ] : Number of iterations = " << itn << "\n";
std::cout << "[ DEBUG ] : Parameter p = " << *( this -> p_PARAM )
<< "; arclength DS = " << this -> DS << "\n";
#endif
if ( itn >= 7 )
{
// converging too slowly, so decrease DS
this -> DS /= this -> ARCSTEP_MULTIPLIER;
#ifdef DEBUG
std::cout << "[ DEBUG ] : I decreased DS to " << this -> DS << "\n";
#endif
}
if ( itn <= 2 )
{
if ( std::abs( this -> DS * this -> ARCSTEP_MULTIPLIER ) < this -> MAX_DS )
{
// converging too quickly, so increase DS
this -> DS *= this -> ARCSTEP_MULTIPLIER;
#ifdef DEBUG
std::cout << "[ DEBUG ] : I increased DS to " << this -> DS << "\n";
#endif
}
}
}
}
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(), and CppNoddy::DenseLinearSystem< _Type >::solve().

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

{
// length of the vector
unsigned N = x.size();
// Jacobian
DenseMatrix<_Type> J( N, N, 0.0 );
// NVectors
DenseVector<_Type> oldFn( N, 0.0 ), newFn( N, 0.0 );
// linear system object - native because no LAPACK complex at the moment
DenseLinearSystem<_Type> system( &J, &oldFn, "native" );
system.set_monitor_det( MONITOR_DET );
// iteration counter
unsigned itn = 0;
do
{
// increment the counter
++itn;
// evaluate the residuals and jacobian at the current state
p_RESIDUAL -> update( x );
// get the residuals
oldFn = p_RESIDUAL -> residual();
#ifdef DEBUG
std::cout << " DEBUG: starting with |Residuals| = " << oldFn.inf_norm() << "\n";
#endif
if ( ( std::abs( oldFn.inf_norm() ) < TOL ) || ( itn == MAX ) )
{
break;
}
// retrieve the current jacobian
J = p_RESIDUAL -> jacobian();
// linear solver
// \todo LU interface to LAPACK for complex matrices
system.solve();
#ifdef DEBUG
std::cout << " DEBUG: Iteration number = " << itn << "\n";
std::cout << " DEBUG: |Newton correction| = " << oldFn.inf_norm() << "\n";
#endif
// must *subtract* delta
x -= oldFn;
}
while ( true );
// More the 'MAX' iterations currently triggers a failure.
if ( itn == MAX )
{
std::string problem;
problem = " The Newton.iterate method took too many iterations. \n";
problem += " At the moment, this is set as a failure. \n";
throw ExceptionItn( problem, itn, oldFn.inf_norm() );
}
LAST_DET_SIGN = system.get_det_sign();
if ( MONITOR_DET )
{
if ( system.get_det_sign() * LAST_DET_SIGN < 0 )
{
std::string problem;
problem = "[ INFO ] : Determinant monitor has changed signs in ODE_BVP.\n";
problem += "[ INFO ] : Bifurcation detected.\n";
throw ExceptionBifurcation( problem );
}
}
}
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().

{
MONITOR_DET = flag;
}
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 >::iterate().

{
iterate( x );
}

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

© 2012

R.E. Hewitt