G4INCL::RootFinder Class Reference

#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 64 of file G4INCLRootFinder.hh.

---

Constructor & Destructor Documentation

G4INCL::RootFinder::RootFinder

(

)

[inline, protected]

Definition at line 108 of file G4INCLRootFinder.hh.

{};

G4INCL::RootFinder::~RootFinder

(

)

[inline, protected]

Definition at line 109 of file G4INCLRootFinder.hh.

{};

---

Member Function Documentation

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

(

)

[inline, static]

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

Returns:

the solution, as an (x,y) pair

Definition at line 83 of file G4INCLRootFinder.hh.

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

{ return RootFinder::solution; }

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 58 of file G4INCLRootFinder.cc.

References G4INCL::RootFunctor::cleanUp(), G4INCL::Math::sign(), and WARN.

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

                                                                         {
    // 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;
  }

---

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

• G4INCLRootFinder.hh
• G4INCLRootFinder.cc

---