A vector NEWTON iteration class. More...

#include <Newton.h>

Inheritance diagram for CppNoddy::Newton< _Type >:

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...

virtual void	solve (DenseVector< _Type > &x)=0
	Compute a solution for that state & parameter variables that are an arc-length 'ds' from the current state. 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 25 of file Newton.h.

Constructor & Destructor Documentation

◆ Newton()

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_object	A pointer to an inherited Residual object
max_steps	The maximum number of iteration steps.
tolerance	A tolerance used as a convergence criterion.
derivative_step	A step length used to compute derivatives.

Definition at line 19 of file Newton.cpp.

    : ArcLength_base< _Type >(),
      TOL(tolerance),
      MAX_STEPS(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

◆ arclength_solve()

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 98 of file Newton.cpp.

                                                           {
#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 ( const ExceptionExternal &error ) {
        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( const ExceptionExternal &error ) {
        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_STEPS) {
        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
        }
      }
    }
  }

References CppNoddy::Utility::dot(), CppNoddy::DenseLinearSystem< _Type >::get_det_sign(), CppNoddy::DenseVector< _Type >::inf_norm(), CppNoddy::DenseLinearSystem< _Type >::set_monitor_det(), and CppNoddy::DenseLinearSystem< _Type >::solve().

Referenced by main().

◆ iterate()

template<typename _Type >

void CppNoddy::Newton< _Type >::iterate ( DenseVector< _Type > & x )

The Newton iteration method.

Parameters

x	An initial guess vector and returns the solution via this too.

Definition at line 35 of file Newton.cpp.

                                                   {
    // 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_STEPS)) {
        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_STEPS' iterations currently triggers a failure.
    if(itn == MAX_STEPS) {
      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);
      }
    }
  }

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

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

◆ set_monitor_det()

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

flag	The value to be set.

Definition at line 47 of file Newton.h.

                                    {
      MONITOR_DET = flag;
    }

Referenced by main().

◆ solve()

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

x	An initial guess vector, the solution overwrites it.

Implements CppNoddy::ArcLength_base< _Type >.

Definition at line 55 of file Newton.h.

                                      {
      iterate(x);
    }

References CppNoddy::Newton< _Type >::iterate().

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

include/Newton.h
src/Newton.cpp

Public Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ Newton()

Member Function Documentation

◆ arclength_solve()

◆ iterate()

◆ set_monitor_det()

◆ solve()