429 lines
15 KiB
C++
429 lines
15 KiB
C++
/**
|
|
* @file RootFind.h Header file for implicit nonlinear solver of a one
|
|
* dimensional function (see \ref numerics and class \link
|
|
* Cantera::RootFind RootFind\endlink).
|
|
*/
|
|
|
|
/*
|
|
* Copyright 2004 Sandia Corporation. Under the terms of Contract
|
|
* DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
|
|
* retains certain rights in this software.
|
|
* See file License.txt for licensing information.
|
|
*/
|
|
|
|
#ifndef CT_ROOTFIND_H
|
|
#define CT_ROOTFIND_H
|
|
/**
|
|
* @defgroup solverGroup Solvers for Equation Systems
|
|
*/
|
|
|
|
#include "ResidEval.h"
|
|
|
|
namespace Cantera
|
|
{
|
|
|
|
//@{
|
|
/// @name Constant which determines the return integer from the routine
|
|
|
|
//! This means that the root solver was a success
|
|
#define ROOTFIND_SUCCESS 0
|
|
//! This return value means that the root finder resolved a solution in the x coordinate
|
|
//! However, convergence in F was not achieved.
|
|
/*!
|
|
* A common situation for this to happen is that f(x) is discontinuous about f(x) = f_0,
|
|
* where we seek the x where the function is equal to f_0. f(x) spans the
|
|
* f_0 while not being equal to f_0 anywhere.
|
|
*/
|
|
#define ROOTFIND_SUCCESS_XCONVERGENCEONLY 1
|
|
//! This means that the root solver failed to achieve convergence
|
|
#define ROOTFIND_FAILEDCONVERGENCE -1
|
|
//! This means that the input to the root solver was defective
|
|
#define ROOTFIND_BADINPUT -2
|
|
//! This means that the rootfinder believes the solution is lower than xmin
|
|
#define ROOTFIND_SOLNLOWERTHANXMIN -3
|
|
//! This means that the rootfinder believes the solution is higher than xmax
|
|
#define ROOTFIND_SOLNHIGHERTHANXMAX -4
|
|
//@}
|
|
|
|
//! Root finder for 1D problems
|
|
/*!
|
|
* The root finder solves a single nonlinear equation described below.
|
|
*
|
|
* \f[
|
|
* f(x) = f_0
|
|
* \f]
|
|
*
|
|
* \f$ f(x) \f$ is assumed to be single valued as a function of x.\f$ f(x) \f$ is not assumed to be continuous nor is
|
|
* its derivative assumed to be well formed.
|
|
*
|
|
* Root finders are significantly different in the sense that do not have to rely
|
|
* solely on Newton's method to find the answer to the problem. Instead they use a method to bound
|
|
* the solution between high and low values and then use a method to refine that bound. The eventual
|
|
* solution to the problem is presented as x_best and as a bound, delta_X, on the solution
|
|
* component. Because of this, they are far more stable for functions and Jacobians that have discontinuities
|
|
* or noise associated with them.
|
|
*
|
|
* The algorithm is a convolution of a local Secant method with an approach of finding a straddle in x.
|
|
* The Jacobian is never required.
|
|
*
|
|
* There is a general breakdown of the algorithm into stages. The first stage seeks to find a straddle of the
|
|
* function. The second stage seeks to reduce the bounds in x and f in order to satisfy the specification of the
|
|
* stopping criteria. In the last stage the algorithm seeks to find the base value of x that satisfies the
|
|
* original equation given what it current knows about the function.
|
|
*
|
|
* Globalization strategy
|
|
*
|
|
* Specifying the General Changes in x
|
|
*
|
|
* Supplying Hints with General Function Behavior Flags
|
|
*
|
|
* Stopping Criteria
|
|
*
|
|
* Specification of the Stopping Criteria
|
|
*
|
|
* Additional constraints
|
|
*
|
|
* Bounds Criteria For the Routine
|
|
*
|
|
* Example
|
|
*
|
|
* @code
|
|
* // Define a residual. The definition of a residual involves a lot more work than is shown here.
|
|
* ResidEval * ec;
|
|
* // Instantiate the root finder with the residual to be solved, ec.
|
|
* RootFind rf(&ec);
|
|
* // Set the relative and absolute tolerancess for f and x.
|
|
* rf.setTol(1.0E-5, 1.0E-10, 1.0E-5, 1.0E-11);
|
|
* // Give a hint about the function's dependence on x. This is needed, for example, if the function has
|
|
* // flat regions.
|
|
* rf.setFuncIsGenerallyIncreasing(true);
|
|
* rf.setDeltaX(0.01);
|
|
* // Supply an initial guess for the solution
|
|
* double xbest = phiM;
|
|
* double oldP = printLvl_;
|
|
* // Set the print level for the solver. Zero produces no output. Two produces a summary table of each iteration.
|
|
* rf.setPrintLvl(2);
|
|
* // Define a minimum and maximum for the independent variable.
|
|
* double phimin = 1.3;
|
|
* double phimax = 2.2;
|
|
* // Define a maximum iteration number
|
|
* int itmax = 100;
|
|
* // Define the f_0 value, and on return will contain the actual value of f(x) obtained
|
|
* double currentObtained;
|
|
* // Call the solver
|
|
* status = rf.solve(phimin, phimax, 100, currentObtained, &xbest);
|
|
* if (status == 0) {
|
|
* if (printLvl_ > 1) {
|
|
* printf("Electrode::integrateConstantCurrent(): Volts (%g amps) = %g\n", currentObtained, xbest);
|
|
* }
|
|
* } else {
|
|
* if (printLvl_) {
|
|
* printf("Electrode::integrateConstantCurrent(): bad status = %d Volts (%g amps) = %g\n",
|
|
* status, currentObtained, xbest);
|
|
* }
|
|
* }
|
|
* @endcode
|
|
*
|
|
* @todo Noise
|
|
* @todo General Search to be done when all else fails
|
|
*/
|
|
class RootFind
|
|
{
|
|
public:
|
|
//! Constructor for the object
|
|
/*!
|
|
* @param resid Pointer to the residual function to be used to calculate f(x)
|
|
*/
|
|
RootFind(ResidEval* resid);
|
|
|
|
//! Copy constructor
|
|
RootFind(const RootFind& r);
|
|
|
|
~RootFind() {}
|
|
|
|
//! Assignment operator
|
|
RootFind& operator=(const RootFind& right);
|
|
|
|
private:
|
|
//! Calculate a deltaX from an input value of x
|
|
/*!
|
|
* This routine ensure that the deltaX will be greater or equal to DeltaXNorm_
|
|
* or 1.0E-14 x
|
|
*
|
|
* @param x1 input value of x
|
|
*/
|
|
doublereal delXNonzero(doublereal x1) const;
|
|
|
|
//! Calculate a deltaX from an input value of x
|
|
/*!
|
|
* This routine ensure that the deltaX will be greater or equal to DeltaXNorm_
|
|
* or 1.0E-14 x or deltaXConverged_.
|
|
*
|
|
* @param x1 input value of x
|
|
*/
|
|
doublereal delXMeaningful(doublereal x1) const;
|
|
|
|
//! Calculate a controlled, nonzero delta between two numbers
|
|
/*!
|
|
* The delta is designed to be greater than or equal to delXMeaningful(x) defined above
|
|
* with the same sign as the original delta. Therefore if you subtract it from either
|
|
* of the two original numbers, you get a different number.
|
|
*
|
|
* @param x2 first number
|
|
* @param x1 second number
|
|
*/
|
|
doublereal deltaXControlled(doublereal x2, doublereal x1) const;
|
|
|
|
//! Function to decide whether two real numbers are the same or not
|
|
/*!
|
|
* A comparison is made between the two numbers to decide whether they
|
|
* are close to one another. This is defined as being within factor * delXMeaningful() of each other.
|
|
*
|
|
* The basic premise here is that if the two numbers are too close, the noise
|
|
* will prevent an accurate calculation of the function and its slope.
|
|
*
|
|
* @param x1 First number
|
|
* @param x2 second number
|
|
* @param factor Multiplicative factor to multiple deltaX with
|
|
* @return Returns a boolean indicating whether the two numbers are the same or not.
|
|
*/
|
|
bool theSame(doublereal x2, doublereal x1, doublereal factor = 1.0) const;
|
|
|
|
public:
|
|
//! Using a line search method, find the root of a 1D function
|
|
/*!
|
|
* This routine solves the following equation.
|
|
*
|
|
* \f[
|
|
* R(x) = f(x) - f_o = 0
|
|
* \f]
|
|
*
|
|
* @param xmin Minimum value of x to be used.
|
|
* @param xmax Maximum value of x to be used
|
|
* @param itmax maximum number of iterations. Usually, it can be less than 50.
|
|
* @param funcTargetValue
|
|
* Value of \f$ f_o \f$ in the equation.
|
|
* On return, it contains the value of the function actually obtained.
|
|
* @param xbest Returns the x that satisfies the function
|
|
* On input, xbest should contain the best estimate of the solution.
|
|
* An attempt to find the solution near xbest is made.
|
|
* @return:
|
|
* 0 = ROOTFIND_SUCCESS Found function
|
|
* -1 = ROOTFIND_FAILEDCONVERGENCE Failed to find the answer
|
|
* -2 = ROOTFIND_BADINPUT Bad input was detected
|
|
*/
|
|
int solve(doublereal xmin, doublereal xmax, int itmax, doublereal& funcTargetValue, doublereal* xbest);
|
|
|
|
//! Return the function value
|
|
/*!
|
|
* This routine evaluates the following equation.
|
|
*
|
|
* \f[
|
|
* R(x) = f(x) - f_o = 0
|
|
* \f]
|
|
*
|
|
* @param x Value of the independent variable
|
|
*
|
|
* @return The routine returns the value of \f$ R(x) \f$
|
|
*/
|
|
doublereal func(doublereal x);
|
|
|
|
//! Set the tolerance parameters for the rootfinder
|
|
/*!
|
|
* These tolerance parameters are used on the function value and the independent value
|
|
* to determine convergence
|
|
*
|
|
* @param rtolf Relative tolerance. The default is 10^-5
|
|
* @param atolf absolute tolerance. The default is 10^-11
|
|
* @param rtolx Relative tolerance. The default is 10^-5
|
|
* Default parameter is 0.0, in which case rtolx is set equal to rtolf
|
|
* @param atolx absolute tolerance. The default is 10^-11
|
|
* Default parameter is 0.0, in which case atolx is set equal to atolf
|
|
*/
|
|
void setTol(doublereal rtolf, doublereal atolf, doublereal rtolx = 0.0, doublereal atolx = 0.0);
|
|
|
|
//! Set the print level from the rootfinder
|
|
/*!
|
|
* 0 -> absolutely nothing is printed for a single time step.
|
|
* 1 -> One line summary per solve_nonlinear call
|
|
* 2 -> short description, points of interest: Table of nonlinear solve - one line per iteration
|
|
* 3 -> Table is included -> More printing per nonlinear iteration (default) that occurs during the table
|
|
* 4 -> Summaries of the nonlinear solve iteration as they are occurring -> table no longer printed
|
|
* 5 -> Algorithm information on the nonlinear iterates are printed out
|
|
* 6 -> Additional info on the nonlinear iterates are printed out
|
|
* 7 -> Additional info on the linear solve is printed out.
|
|
* 8 -> Info on a per iterate of the linear solve is printed out.
|
|
*
|
|
* @param printLvl integer value
|
|
*/
|
|
void setPrintLvl(int printLvl);
|
|
|
|
//! Set the function behavior flag
|
|
/*!
|
|
* If this is true, the function is generally an increasing function of x.
|
|
* In particular, if the algorithm is seeking a higher value of f, it will look
|
|
* in the positive x direction.
|
|
*
|
|
* This type of function is needed because this algorithm must deal with regions of f(x) where
|
|
* f is not changing with x.
|
|
*
|
|
* @param value boolean value
|
|
*/
|
|
void setFuncIsGenerallyIncreasing(bool value);
|
|
|
|
//! Set the function behavior flag
|
|
/*!
|
|
* If this is true, the function is generally a decreasing function of x.
|
|
* In particular, if the algorithm is seeking a higher value of f, it will look
|
|
* in the negative x direction.
|
|
*
|
|
* This type of function is needed because this algorithm must deal with regions of f(x) where
|
|
* f is not changing with x.
|
|
*
|
|
* @param value boolean value
|
|
*/
|
|
void setFuncIsGenerallyDecreasing(bool value);
|
|
|
|
//! Set the minimum value of deltaX
|
|
/*!
|
|
* This sets the value of deltaXNorm_
|
|
*
|
|
* @param deltaXNorm
|
|
*/
|
|
void setDeltaX(doublereal deltaXNorm);
|
|
|
|
//! Set the maximum value of deltaX
|
|
/*!
|
|
* This sets the value of deltaXMax_
|
|
*
|
|
* @param deltaX
|
|
*/
|
|
void setDeltaXMax(doublereal deltaX);
|
|
|
|
//! Print the iteration history table
|
|
void printTable();
|
|
|
|
public:
|
|
//! Pointer to the residual function evaluator
|
|
ResidEval* m_residFunc;
|
|
|
|
//! Target value for the function. We seek the value of f that is equal to this value
|
|
doublereal m_funcTargetValue;
|
|
|
|
//! Absolute tolerance for the value of f
|
|
doublereal m_atolf;
|
|
|
|
//! Absolute tolerance for the value of x
|
|
doublereal m_atolx;
|
|
|
|
//! Relative tolerance for the value of f and x
|
|
doublereal m_rtolf;
|
|
|
|
//! Relative tolerance for the value of x
|
|
doublereal m_rtolx;
|
|
|
|
//! Maximum number of step sizes
|
|
doublereal m_maxstep;
|
|
|
|
protected:
|
|
//! Print level
|
|
/*!
|
|
* 0 No printing of any kind
|
|
* 1 Single print line indicating success or failure of the routine.
|
|
* 2 Summary table printed at the end of the routine, with a convergence history
|
|
* 3 Printouts during the iteration are added. Summary table is printed out at the end.
|
|
* if writeLogAllowed_ is turned on, a file is written out with the convergence history.
|
|
*/
|
|
int printLvl;
|
|
|
|
public:
|
|
//! Boolean to turn on the possibility of writing a log file.
|
|
bool writeLogAllowed_;
|
|
|
|
protected:
|
|
//! Delta X norm. This is the nominal value of deltaX that will be used by the program
|
|
doublereal DeltaXnorm_;
|
|
|
|
//! Boolean indicating whether DeltaXnorm_ has been specified by the user or not
|
|
int specifiedDeltaXnorm_;
|
|
|
|
//! Delta X Max. This is the maximum value of deltaX that will be used by the program
|
|
/*!
|
|
* Sometimes a large change in x causes problems.
|
|
*/
|
|
doublereal DeltaXMax_;
|
|
|
|
//! Boolean indicating whether DeltaXMax_ has been specified by the user or not
|
|
int specifiedDeltaXMax_;
|
|
|
|
//! Boolean indicating whether the function is an increasing with x
|
|
bool FuncIsGenerallyIncreasing_;
|
|
|
|
//! Boolean indicating whether the function is decreasing with x
|
|
bool FuncIsGenerallyDecreasing_;
|
|
|
|
//! Value of delta X that is needed for convergence
|
|
/*!
|
|
* X will be considered as converged if we are within deltaXConverged_ of the solution
|
|
* The default is zero.
|
|
*/
|
|
doublereal deltaXConverged_;
|
|
|
|
//! Internal variable tracking largest x tried.
|
|
doublereal x_maxTried_;
|
|
|
|
//! Internal variable tracking f(x) of largest x tried.
|
|
doublereal fx_maxTried_;
|
|
|
|
//! Internal variable tracking smallest x tried.
|
|
doublereal x_minTried_;
|
|
|
|
//! Internal variable tracking f(x) of smallest x tried.
|
|
doublereal fx_minTried_;
|
|
|
|
//! Structure containing the iteration history
|
|
struct rfTable {
|
|
//@{
|
|
int its;
|
|
int TP_its;
|
|
double slope;
|
|
double xval;
|
|
double fval;
|
|
int foundPos;
|
|
int foundNeg;
|
|
double deltaXConverged;
|
|
double deltaFConverged;
|
|
double delX;
|
|
|
|
std::string reasoning;
|
|
|
|
void clear() {
|
|
its = 0;
|
|
TP_its = 0;
|
|
slope = -1.0E300;
|
|
xval = -1.0E300;
|
|
fval = -1.0E300;
|
|
reasoning = "";
|
|
};
|
|
|
|
rfTable() :
|
|
its(-2),
|
|
TP_its(0),
|
|
slope(-1.0E300),
|
|
xval(-1.0E300),
|
|
fval(-1.0E300),
|
|
foundPos(0),
|
|
foundNeg(0),
|
|
deltaXConverged(-1.0E300),
|
|
deltaFConverged(-1.0E300),
|
|
delX(-1.0E300),
|
|
reasoning("") {
|
|
};
|
|
//@}
|
|
};
|
|
|
|
//! Vector of iteration histories
|
|
std::vector<struct rfTable> rfHistory_;
|
|
};
|
|
}
|
|
#endif
|