G4INCL::RootFinder Namespace Reference

Data Structures
class	Solution

Functions
Solution	solve (RootFunctor const *const f, const G4double x0)
	Numerically solve a one-dimensional equation. More...

Function Documentation

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

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

a Solution object describing the root, if it was found

Definition at line 117 of file G4INCLRootFinder.cc.

References G4INCL::RootFunctor::cleanUp(), INCL_DEBUG, G4INCL::Math::sign(), and test::x.

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

                                                                   {
      // If we already have the solution, do nothing
      const G4double y0 = (*f)(x0);
      if( std::abs(y0) < toleranceY ) {
        return Solution(x0,y0);
      }

      // 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);
          return Solution(0.,y_at_zero);
        } else {
          INCL_DEBUG("Root-finding algorithm could not bracket the root." << std::endl);
          f->cleanUp(false);
          return Solution();
        }
      }

      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) {
          INCL_DEBUG("Root-finding algorithm did not converge." << std::endl);
          f->cleanUp(false);
          return Solution();
        }

        // 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
       * ******************************/

      f->cleanUp(true);
      return Solution(x,y);
    }