#include <G4INCLRootFinder.hh>

Static Public Member Functions
static G4bool	solve (RootFunctor const *const f, const G4double x0)
	Numerically solve a one-dimensional equation.

static std::pair< G4double, G4double > const &	getSolution ()
	Get the solution of the last call to solve().

Protected Member Functions
	RootFinder ()

	~RootFinder ()

Detailed Description

Definition at line 66 of file G4INCLRootFinder.hh.

Constructor & Destructor Documentation

◆ RootFinder()

G4INCL::RootFinder::RootFinder ( )

inlineprotected

Definition at line 110 of file G4INCLRootFinder.hh.

110{};

◆ ~RootFinder()

G4INCL::RootFinder::~RootFinder ( )

inlineprotected

Definition at line 111 of file G4INCLRootFinder.hh.

111{};

Member Function Documentation

◆ getSolution()

static std::pair< G4double, G4double > const & G4INCL::RootFinder::getSolution ( )

inlinestatic

Get the solution of the last call to solve().

Returns: the solution, as an (x,y) pair

Definition at line 85 of file G4INCLRootFinder.hh.

85{ return RootFinder::solution; }

Referenced by G4INCL::InteractionAvatar::enforceEnergyConservation().

◆ solve()

G4bool G4INCL::RootFinder::solve	(	RootFunctor const *const	f,
		const G4double	x0
	)

static

Numerically solve a one-dimensional equation.

Numerically solves the equation f(x)==0. This implementation uses the false-position method.

If a root is found, it can be retrieved using the getSolution() method,

Parameters

f	pointer to a RootFunctor
x0	initial value of the function argument

Returns: true if a root was found

Definition at line 60 of file G4INCLRootFinder.cc.

                                                                         {
    // If we already have the solution, do nothing
    const G4double y0 = (*f)(x0);
    if( std::abs(y0) < toleranceY ) {
      solution = std::make_pair(x0,y0);
      return true;
    }
 
    // Bracket the root and set the initial values
    std::pair<G4double,G4double> bracket = bracketRoot(f,x0);
    G4double x1 = bracket.first;
    G4double x2 = bracket.second;
    // If x1>x2, it means that we could not bracket the root. Return false.
    if(x1>x2) {
      // Maybe zero is a good solution?
      G4double y_at_zero = (*f)(0.);
      if(std::abs(y_at_zero)<=toleranceY) {
        f->cleanUp(true);
        solution = std::make_pair(0.,y_at_zero);
        return true;
      } else {
        WARN("Root-finding algorithm could not bracket the root." << std::endl);
        f->cleanUp(false);
        return false;
      }
    }
 
    G4double y1 = (*f)(x1);
    G4double y2 = (*f)(x2);
    G4double x = x1;
    G4double y = y1;
 
    /* ********************************
     * Start of the false-position loop
     * ********************************/
 
    // Keep track of the last updated interval end (-1=left, 1=right)
    G4int lastUpdated = 0;
 
    for(G4int iterations=0; std::abs(y) > toleranceY; iterations++) {
 
      if(iterations > maxIterations) {
        WARN("Root-finding algorithm did not converge." << std::endl);
        f->cleanUp(false);
        return false;
      }
 
      // Estimate the root position by linear interpolation
      x = (y1*x2-y2*x1)/(y1-y2);
 
      // Update the value of the function
      y = (*f)(x);
 
      // Update the bracketing interval
      if(Math::sign(y) == Math::sign(y1)) {
        x1=x;
        y1=y;
        if(lastUpdated==-1) y2 *= 0.5;
        lastUpdated = -1;
      } else {
        x2=x;
        y2=y;
        if(lastUpdated==1) y1 *= 0.5;
        lastUpdated = 1;
      }
    }
 
    /* ******************************
     * End of the false-position loop
     * ******************************/
 
    solution = std::make_pair(x,y);
    f->cleanUp(true);
    return true;
  }