All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines
Template class to search roots of a 1D function. More...
#include <RootSolver.h>
Classes | |
| class | AbsoluteDownIterator |
| class | AbsoluteUpIterator |
| class | RelativeDownIterator |
| class | RelativeUpIterator |
Public Member Functions | |
| void | halving () |
| void | inversePolarity () |
| RealType | lower () const |
| void | neville () |
| void | newtonRaphson () |
| RealType | precision () const |
| RootSolverTemplate (FunctionClass *values) | |
| bool | searchDown (RealType from, RealType to) |
| bool | searchDown (RealType from, RealType min, RealType max) |
| bool | searchDownSlope (RealType from, RealType min, RealType max) |
| RealType | searchStep () const |
| bool | searchUp (RealType from, RealType to) |
| bool | searchUp (RealType from, RealType min, RealType max) |
| bool | searchUpSlope (RealType from, RealType min, RealType max) |
| void | setAbsoluteStep (RealType step) |
| bool | setInterval (RealType min, RealType max) |
| void | setPolarity (RealType x0) |
| void | setPrecision (RealType prec) |
| void | setRelativeStep (RealType step) |
| RealType | upper () const |
| FunctionClass * | values () |
Template class to search roots of a 1D function.
The roots are searched in order provided that the search step is sufficiently small. The search step can be adaptative as it can defined by a ratio. For searching a root arround zero, the search step must be absolute. After bracketing a root, it can be refined by various algorithms: halving, Neville, or Newton-Raphson.
The function values are calculated by the function class provided to the constructor.
Before any call to a searching function, the search step and the polarity of the starting point must be set. The polarity is calculated by setPolarity(), inversePolarity(). The step is accessed with searchStep(), setAbsoluteStep(), and setRelativeStep().
After a successfull search, the current bracket() can be retrieved and passed to RootRefine to refine the root location down to an arbitrary precision.
| QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::RootSolverTemplate | ( | FunctionClass * | values | ) |
Constructor
The search step is set by default to one tenth of the tail value of the bracket (relative step). The default polarity is 1. The precision is the relative precision that will be achieved by halving(), neville(), or newtonRaphson(). The default value is 1e-7.
{
_values=values;
_dxType=Relative;
_dx=0.1;
_polarity=1.0;
_precision=1e-7;
}
| void QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::halving | ( | ) |
Refines the current bracket interval down to precision set by setPrecision(). The method is the most robust one. The interval is iteratively divided in two parts.
References QGpCoreTools::endl(), ROOTSOLVER_MAX_ITERATIONS, and QGpCoreTools::App::stream().
{
// First halving
RealType x3=0.5 * (_x1 + _x2);
RealType y3=_values->y(x3);
// Relative Convergence criteria
RealType tmpAbs1;
tmpAbs1=_x1 + _x2;
if(tmpAbs1 < 0) tmpAbs1=-tmpAbs1;
RealType maxDiff=tmpAbs1 * _precision;
RealType minDiff=-maxDiff;
for(register int i=1;i < ROOTSOLVER_MAX_ITERATIONS;i++ ) {
if(y3==0.0) {
// just move the X3 by one tenth of actuel range
x3 -= (_x2 - _x1) * 0.1;
// recalculate the Y3, it will always be different from 0
y3=_values->y(x3);
}
if(isFound(_y1, y3)) {
_x2=x3;
_y2=y3;
} else {
_x1=x3;
_y1=y3;
}
// Check for convergence with a relative criteria, if so exit
RealType diff=_x1 - _x2;
if(minDiff <= diff && diff <= maxDiff) {
return;
}
x3=0.5 * (_x1 + _x2);
y3=_values->y(x3);
}
App::stream() << "Halving maximum iteration reached" << endl;
}
| void QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::inversePolarity | ( | ) | [inline] |
Simply inverses the sign of the polarity, used to calculate the next root after a first successful solve().
{_polarity=-_polarity;}
| RealType QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::lower | ( | ) | const [inline] |
Returns the lower value of the current bracket interval.
{return _x1<_x2 ? _x1 : _x2;}
| void QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::neville | ( | ) |
Refines the current bracket interval down to precision set by setPrecision(). Interval halving is used where other schemes would be inefficient. Once a suitable region is found the Neville's algorithm (polynomial fit of x=x(y) at y=0) is used to find the root. The procedure alternates between the interval halving and Neville techniques using whichever is most efficient.
The scheme of this function comes from FORTRAN 77 code written by D. R. Russell, R. B. Herrmann (Department of Earth and Atmospheric Sciences Saint Louis University 221 North Grand Boulevard St. Louis, Missouri 63103 U. S. A. COPYRIGHT 1986, 1991). Recent optimizations to take advantage of the full power of the polynomial fit are due to Marc Wathelet. The main feature is a voluntary shift of the polynomial interpolation to maintain the order of magnitude of the function value at the limits of the current interval in the same range.
See Numerical Recipes in C, chapter 3.1 for more details on Neville
References QGpCoreTools::endl(), ROOTSOLVER_MAX_ITERATIONS, and QGpCoreTools::App::stream().
{
//#define neville_DEBUG
//#define neville_STAT
#if defined(neville_DEBUG) || defined(neville_STAT)
printf("RootSolverTemplate::neville() @ precision=%lg\n", _precision);
#endif
// Variables for polynomial fit of Neville
int m=0;
RealType p[ 20 ], nevY[ 20 ];
// First halving
RealType x3=0.5 * (_x1 + _x2);
RealType y3=_values->y(x3);
// Relative Convergence criteria
RealType maxDiff=fabs(x3)*_precision;
for (register int i=1; i<ROOTSOLVER_MAX_ITERATIONS; i++) {
// define new bounds according to signs of y1 and y3
// If y1 and y3 have different signs then the root is between
// y1 and y3, thus forget (x2,y2) else forget (x1,y1).
// this case is not implemented in Herrmann's code
// In our case, it is very important to closely bracket the root
// X1 and X2 must ALWAYS stay on the same sides of the root during all
// Neville's iterations.
#ifdef neville_DEBUG
//printf("1/_x1=%.20lf _y1=%.20lg 1/_x2=%.20lf _y2=%.20lg\n",1/_x1,_y1,1/_x2,_y2);
printf("_x1=%.20lf _y1=%.20lg _x2=%.20lf _y2=%.20lg x3=%.20lf y3=%.20lg diff=%.20lg\n", _x1, _y1, _x2, _y2, x3, y3, fabs(_x1-_x2));
#endif
if(y3==0.0) {
// Move the x3 towards the furthest limit
if(fabs(_x1-x3)>fabs(_x2-x3)) {
x3=0.5*(_x1+x3);
} else {
x3=0.5*(_x2+x3);
}
// recalculate the Y3, it will always be different from 0
y3=_values->y(x3);
#ifdef neville_DEBUG
printf("y3 is null: x3=%.20lf y3=%.20lg\n", x3, y3);
#endif
}
if(isFound(_y1, y3)) {
_x2=x3;
_y2=y3;
} else {
_x1=x3;
_y1=y3;
}
// Check for convergence with a relative criteria, if so exit
if(fabs(_x1-_x2)<=maxDiff) {
#if defined(neville_DEBUG) || defined(neville_STAT)
printf( "Number of iterations=%i \n", i);
printf( "Neville poly m=%i \n", m);
printf( "Convergence achieved\n" );
#endif
return;
}
// Check the slopes between x1,x3 and x3,x2
// If they are different, the function is strongly non linear,
// thus use only halving, not Neville
if(isFound(_y1 - y3, y3 - _y2)) {
#if defined(neville_DEBUG) || defined(neville_STAT)
printf( "Neville poly m=%i \n", m);
printf( "Inverse slopes: " );
#endif
x3=0.5* (_x1 + _x2);
y3=_values->y(x3);
m=0;
} else {
// Until 2009-02-03, this test was performed, found to be useless
// If y1 and y2 differ by more than a factor of 100
// use only halving to avoid poor behavior of polynomial fit
#if 0
RealType ay1=_y1;
if(ay1 < 0) ay1=-ay1;
RealType ay2=_y2;
if(ay2 < 0) ay2=-ay2;
if(ay1 > 100 * ay2 || ay2 > 100 * ay1) {
#if defined(neville_DEBUG) || defined(neville_STAT)
printf( "Neville poly m=%i \n", m);
printf( "Factor of more than 100: " );
#endif
x3=0.5* (x1 + x2);
y3=_values->y(x3);
m=0;
} else {
#endif
#ifdef neville_DEBUG
printf( "Neville poly m=%i \n", m);
#endif
// if(m>10) cout << "**** WARNING: m=" << m << endl;
if(m>0) {
p[ m + 1 ]=x3;
nevY[ m + 1 ]=y3;
} else {
p[ 0 ]=_x1;
nevY[ 0 ]=_y1;
p[ 1 ]=_x2;
nevY[ 1 ]=_y2;
}
// perform Neville iteration. Instead of generating y(x)
// we interchange the x and y of formula to solve for x(y)
// when y=0.
// A new point is added to the pyramidal contruction.
// The new point is always added at the base, adding a new
// slice at each iteration, each time larger
// As pyramide grows, the degree of the polynome increase
// thus the number of point taken into account increase also.
for(int j=m;j >= 0;j-- ) {
RealType denom=nevY[ m + 1 ] - nevY[ j ];
RealType tmpAbs1=nevY[ m + 1 ];
if(tmpAbs1 < 0) tmpAbs1=-tmpAbs1;
RealType tmpAbs2=denom;
if(tmpAbs2 < 0) tmpAbs2=-tmpAbs2;
if(tmpAbs2 < 1.0e-10 * tmpAbs1) {
#if defined(neville_DEBUG) || defined(neville_STAT)
printf( "Neville poly m=%i \n", m);
printf( "Denom too small (%lf), halving\n", denom);
#endif
p[ 0 ]=0.5 * (_x1 + _x2);
m=0;
break;
} else p[ j ]=(nevY[ m + 1 ] * p[ j ] - nevY[ j ] * p[ j + 1 ] )/denom;
}
x3=p[ 0 ];
// Just move the X3 by one tenth of actual range towards the limit with the highest y
RealType tmpAbs1=_y1<0 ? -_y1 : _y1;
RealType tmpAbs2=_y2<0 ? -_y2 : _y2;
if(tmpAbs1 > tmpAbs2)
x3 -= 0.1 * (x3 - _x1);
else
x3 -= 0.1 * (x3 - _x2);
//printf("Function estimation\n");
m++;
if(m > 18) m=1;
// Make sure new estimate is inside the previous values
// if not perform interval halving
if(_x1 <= _x2) {
if(x3 < _x1 || x3 > _x2) {
#if defined(neville_DEBUG) || defined(neville_STAT)
printf( "Neville poly m=%i \n", m);
printf( "New estimate out of range, halving\n" );
#endif
x3=0.5* (_x1 + _x2);
m=1;
}
} else {
if(x3 < _x2 || x3 > _x1) {
#if defined(neville_DEBUG) || defined(neville_STAT)
printf( "Neville poly m=%i \n", m);
printf( "New estimate out of range, halving\n" );
#endif
x3=0.5* (_x1 + _x2);
m=1;
}
}
y3=_values->y(x3);
}
}
App::stream() << "Neville maximum iteration reached" << endl;
}
| void QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::newtonRaphson | ( | ) |
Refines the current bracket interval down to precision set by setPrecision(). The method is based upon the derivative. The derivative is computed by the function class under derivative().
References QGpCoreTools::endl(), ROOTSOLVER_MAX_ITERATIONS, and QGpCoreTools::App::stream().
{
// First halving
RealType x3=0.5 * (_x1 + _x2);
RealType y3=_values->y(x3);
// Relative Convergence criteria
RealType tmpAbs1;
tmpAbs1=_x1 + _x2;
if(tmpAbs1 < 0) tmpAbs1=-tmpAbs1;
RealType maxDiff=tmpAbs1 * _precision;
RealType minDiff=-maxDiff;
for(register int i=1;i < ROOTSOLVER_MAX_ITERATIONS;i++ ) {
if(y3==0.0) {
// just move the X3 by one tenth of actuel range
x3 -= (_x2 - _x1) * 0.1;
// recalculate the Y3, it will always be different from 0
y3=_values->y(x3);
}
if(isFound(_y1, y3)) {
_x2=x3;
_y2=y3;
} else {
_x1=x3;
_y1=y3;
}
// Check for convergence with a relative criteria, if so exit
RealType diff=_x1 - _x2;
if(minDiff <= diff && diff <= maxDiff) {
return;
}
// Compute the derivative estimation from the lowest y
// Just move the X3 by one tenth of actual range towards the limit with the highest y
RealType tmpAbs1=_y1<0 ? -_y1 : _y1;
RealType tmpAbs2=_y2<0 ? -_y2 : _y2;
if(tmpAbs1 > tmpAbs2) {
x3=_x2 - _y2/_values->derivative(_x2);
x3 -= 0.1 * (x3 - _x1);
} else {
x3=_x1 - _y1/_values->derivative(_x1);
x3 -= 0.1 * (x3 - _x2);
}
// Make sure new estimate is inside the previous values
// if not perform interval halving
if(_x1 <= _x2) {
if(x3 < _x1 || x3 > _x2) {
x3=0.5* (_x1 + _x2);
}
} else {
if(x3 < _x2 || x3 > _x1) {
x3=0.5* (_x1 + _x2);
}
}
y3=_values->y(x3);
}
App::stream() << "Newton-Raphson maximum iteration reached" << endl;
}
| double QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::precision | ( | ) | const [inline] |
Return the current relative precision.
{return _precision;}
| bool QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::searchDown | ( | RealType | from, |
| RealType | to | ||
| ) | [inline] |
Searches for roots starting from from. If no root is found, it stops at to. from must be greater than to.
Returns true if a root is effectively found between from and to.
{
if(from>to) return false;
if(_dxType==Absolute) {
AbsoluteDownIterator it(_dx);
return search(from, to, it);
} else {
RelativeDownIterator it(_dx);
return search(from, to, it);
}
}
| bool QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::searchDown | ( | RealType | from, |
| RealType | min, | ||
| RealType | max | ||
| ) | [inline] |
{
if(_dxType==Absolute) {
AbsoluteDownIterator it(_dx);
return search(from, min, max, it);
} else {
RelativeDownIterator it(_dx);
return search(from, min, max, it);
}
}
| bool QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::searchDownSlope | ( | RealType | from, |
| RealType | min, | ||
| RealType | max | ||
| ) | [inline] |
Overload function for oscillating functions. A root is searched by decreasing the initial value from.
Returns true if a root is effectively found between min and max.
{
if(_dxType==Absolute) {
AbsoluteDownIterator it(_dx);
return searchSlope(from, min, max, it);
} else {
RelativeDownIterator it(_dx);
return searchSlope(from, min, max, it);
}
}
| RealType QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::searchStep | ( | ) | const [inline] |
Returns the current step (ratio or absolute according to current type).
{return _dx;}
| bool QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::searchUp | ( | RealType | from, |
| RealType | to | ||
| ) | [inline] |
Searches for roots starting from from. If no root is found, it stops at to. from must be less than to.
Returns true if a root is effectively found between from and to.
{
if(from>to) return false;
if(_dxType==Absolute) {
AbsoluteUpIterator it(_dx);
return search(from, to, it);
} else {
RelativeUpIterator it(_dx);
return search(from, to, it);
}
}
| bool QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::searchUp | ( | RealType | from, |
| RealType | min, | ||
| RealType | max | ||
| ) | [inline] |
{
if(_dxType==Absolute) {
AbsoluteUpIterator it(_dx);
return search(from, min, max, it);
} else {
RelativeUpIterator it(_dx);
return search(from, min, max, it);
}
}
| bool QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::searchUpSlope | ( | RealType | from, |
| RealType | min, | ||
| RealType | max | ||
| ) | [inline] |
Overload function for oscillating functions. A root is searched by increasing the initial value from.
If the initial sampling step is large, a minimum or a maximum may be missed, looking regularly at the slope of the curve helps: if the sign of the slope changes, the minimum or maximum is refined and the presence of hidden root is checked.
Returns true if a root is effectively found between min and max.
{
if(_dxType==Absolute) {
AbsoluteUpIterator it(_dx);
return searchSlope(from, min, max, it);
} else {
RelativeUpIterator it(_dx);
return searchSlope(from, min, max, it);
}
}
| void QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::setAbsoluteStep | ( | RealType | step | ) | [inline] |
Set the current step.
{_dx=fabs(step); _dxType=Absolute;}
| bool QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::setInterval | ( | RealType | min, |
| RealType | max | ||
| ) |
If the root is already localized, there is no need to search for a root before refining it. Set directly the interval limits and returns true if there is effectively a root in the range.
{
_x1=min;
_y1=_values->y(_x1);
_x2=max;
_y2=_values->y(_x2);
return isFound(_y1, _y2);
}
| void QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::setPolarity | ( | RealType | x0 | ) |
Calculate the initial polarity. x0 must be chosen so that the value of the function at x0 is representative of the sign of the function at positive infinity.
{
_polarity=_values->y(x0);
}
| void QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::setPrecision | ( | RealType | prec | ) | [inline] |
Set the current relative precision used as a completion criterium for halving(), neville(), and newtonRaphson().
{_precision=prec;}
| void QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::setRelativeStep | ( | RealType | step | ) | [inline] |
Set the current step ratio. The absolute search step is defined by the product of the step and the current lower bound of the bracket interval.
{_dx=fabs(step); _dxType=Relative;}
| RealType QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::upper | ( | ) | const [inline] |
Returns the upper value of the current bracket interval.
{return _x1<_x2 ? _x2 : _x1;}
| FunctionClass* QGpCoreTools::RootSolverTemplate< RealType, FunctionClass >::values | ( | ) | [inline] |
{return _values;}
