Cantera  2.1.2
RootFind.h
Go to the documentation of this file.
1 /**
2  * @file RootFind.h Header file for implicit nonlinear solver of a one
3  * dimensional function (see \ref numerics and class \link
5  */
6
7 /*
8  * Copyright 2004 Sandia Corporation. Under the terms of Contract
9  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
10  * retains certain rights in this software.
11  * See file License.txt for licensing information.
12  */
13
14 #ifndef CT_ROOTFIND_H
15 #define CT_ROOTFIND_H
16 /**
17  * @defgroup solverGroup Solvers for Equation Systems
18  */
19
20 #include "ResidEval.h"
21
22 namespace Cantera
23 {
24
25 //@{
26 /// @name Constant which determines the return integer from the routine
27
28 //! This means that the root solver was a success
29 #define ROOTFIND_SUCCESS 0
30 //! This return value means that the root finder resolved a solution in the x coordinate
31 //! However, convergence in F was not achieved.
32 /*!
33  * A common situation for this to happen is that f(x) is discontinuous about f(x) = f_0,
34  * where we seek the x where the function is equal to f_0. f(x) spans the
35  * f_0 while not being equal to f_0 anywhere.
36  */
37 #define ROOTFIND_SUCCESS_XCONVERGENCEONLY 1
38 //! This means that the root solver failed to achieve convergence
39 #define ROOTFIND_FAILEDCONVERGENCE -1
40 //! This means that the input to the root solver was defective
42 //! This means that the rootfinder believes the solution is lower than xmin
43 #define ROOTFIND_SOLNLOWERTHANXMIN -3
44 //! This means that the rootfinder believes the solution is higher than xmax
45 #define ROOTFIND_SOLNHIGHERTHANXMAX -4
46 //@}
47
48 //! Root finder for 1D problems
49 /*!
50  * The root finder solves a single nonlinear equation described below.
51  *
52  * \f[
53  * f(x) = f_0
54  * \f]
55  *
56  * \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
57  * its derivative assumed to be well formed.
58  *
59  * Root finders are significantly different in the sense that do not have to rely
60  * solely on Newton's method to find the answer to the problem. Instead they use a method to bound
61  * the solution between high and low values and then use a method to refine that bound. The eventual
62  * solution to the problem is presented as x_best and as a bound, delta_X, on the solution
63  * component. Because of this, they are far more stable for functions and Jacobians that have discontinuities
64  * or noise associated with them.
65  *
66  * The algorithm is a convolution of a local Secant method with an approach of finding a straddle in x.
67  * The Jacobian is never required.
68  *
69  * There is a general breakdown of the algorithm into stages. The first stage seeks to find a straddle of the
70  * function. The second stage seeks to reduce the bounds in x and f in order to satisfy the specification of the
71  * stopping criteria. In the last stage the algorithm seeks to find the base value of x that satisfies the
72  * original equation given what it current knows about the function.
73  *
74  * Globalization strategy
75  *
76  * Specifying the General Changes in x
77  *
78  * Supplying Hints with General Function Behavior Flags
79  *
80  * Stopping Criteria
81  *
82  * Specification of the Stopping Criteria
83  *
85  *
86  * Bounds Criteria For the Routine
87  *
88  * Example
89  *
90  * @code
91  * // Define a residual. The definition of a residual involves a lot more work than is shown here.
92  * ResidEval * ec;
93  * // Instantiate the root finder with the residual to be solved, ec.
94  * RootFind rf(&ec);
95  * // Set the relative and absolute tolerancess for f and x.
96  * rf.setTol(1.0E-5, 1.0E-10, 1.0E-5, 1.0E-11);
97  * // Give a hint about the function's dependence on x. This is needed, for example, if the function has
98  * // flat regions.
99  * rf.setFuncIsGenerallyIncreasing(true);
100  * rf.setDeltaX(0.01);
101  * // Supply an initial guess for the solution
102  * double xbest = phiM;
103  * double oldP = printLvl_;
104  * // Set the print level for the solver. Zero produces no output. Two produces a summary table of each iteration.
105  * rf.setPrintLvl(2);
106  * // Define a minimum and maximum for the independent variable.
107  * double phimin = 1.3;
108  * double phimax = 2.2;
109  * // Define a maximum iteration number
110  * int itmax = 100;
111  * // Define the f_0 value, and on return will contain the actual value of f(x) obtained
112  * double currentObtained;
113  * // Call the solver
114  * status = rf.solve(phimin, phimax, 100, currentObtained, &xbest);
115  * if (status == 0) {
116  * if (printLvl_ > 1) {
117  * printf("Electrode::integrateConstantCurrent(): Volts (%g amps) = %g\n", currentObtained, xbest);
118  * }
119  * } else {
120  * if (printLvl_) {
121  * printf("Electrode::integrateConstantCurrent(): bad status = %d Volts (%g amps) = %g\n",
122  * status, currentObtained, xbest);
123  * }
124  * }
125  * @endcode
126  *
127  * @todo Noise
128  * @todo General Search to be done when all else fails
129  *
130  */
131 class RootFind
132 {
133 public:
134  //! Constructor for the object
135  /*!
136  * @param resid Pointer to the residual function to be used to calculate f(x)
137  */
138  RootFind(ResidEval* resid);
139
140  //! Copy constructor
141  RootFind(const RootFind& r);
142
143  ~RootFind();
144
145  //! Assignment operator
146  RootFind& operator=(const RootFind& right);
147
148 private:
149  //! Calculate a deltaX from an input value of x
150  /*!
151  * This routine ensure that the deltaX will be greater or equal to DeltaXNorm_
152  * or 1.0E-14 x
153  *
154  * @param x1 input value of x
155  */
156  doublereal delXNonzero(doublereal x1) const;
157
158  //! Calculate a deltaX from an input value of x
159  /*!
160  * This routine ensure that the deltaX will be greater or equal to DeltaXNorm_
161  * or 1.0E-14 x or deltaXConverged_.
162  *
163  * @param x1 input value of x
164  */
165  doublereal delXMeaningful(doublereal x1) const;
166
167  //! Calculate a controlled, nonzero delta between two numbers
168  /*!
169  * The delta is designed to be greater than or equal to delXMeaningful(x) defined above
170  * with the same sign as the original delta. Therefore if you subtract it from either
171  * of the two original numbers, you get a different number.
172  *
173  * @param x2 first number
174  * @param x1 second number
175  */
176  doublereal deltaXControlled(doublereal x2, doublereal x1) const;
177
178  //! Function to decide whether two real numbers are the same or not
179  /*!
180  * A comparison is made between the two numbers to decide whether they
181  * are close to one another. This is defined as being within factor * delXMeaningful() of each other.
182  *
183  * The basic premise here is that if the two numbers are too close, the noise
184  * will prevent an accurate calculation of the function and its slope.
185  *
186  * @param x1 First number
187  * @param x2 second number
188  * @param factor Multiplicative factor to multiple deltaX with
189  *
190  * @return Returns a boolean indicating whether the two numbers are the same or not.
191  */
192  bool theSame(doublereal x2, doublereal x1, doublereal factor = 1.0) const;
193
194 public:
195  //! Using a line search method, find the root of a 1D function
196  /*!
197  * This routine solves the following equation.
198  *
199  * \f[
200  * R(x) = f(x) - f_o = 0
201  * \f]
202  *
203  * @param xmin Minimum value of x to be used.
204  * @param xmax Maximum value of x to be used
205  * @param itmax maximum number of iterations. Usually, it can be less than 50.
206  * @param funcTargetValue
207  * Value of \f$f_o \f$ in the equation.
208  * On return, it contains the value of the function actually obtained.
209  * @param xbest Returns the x that satisfies the function
210  * On input, xbest should contain the best estimate of the solution.
211  * An attempt to find the solution near xbest is made.
212  *
213  * @return:
214  * 0 = ROOTFIND_SUCCESS Found function
215  * -1 = ROOTFIND_FAILEDCONVERGENCE Failed to find the answer
217  */
218  int solve(doublereal xmin, doublereal xmax, int itmax, doublereal& funcTargetValue, doublereal* xbest);
219
220  //! Return the function value
221  /*!
222  * This routine evaluates the following equation.
223  *
224  * \f[
225  * R(x) = f(x) - f_o = 0
226  * \f]
227  *
228  * @param x Value of the independent variable
229  *
230  * @return The routine returns the value of \f$R(x) \f$
231  */
232  doublereal func(doublereal x);
233
234  //! Set the tolerance parameters for the rootfinder
235  /*!
236  * These tolerance parameters are used on the function value and the independent value
237  * to determine convergence
238  *
239  * @param rtolf Relative tolerance. The default is 10^-5
240  * @param atolf absolute tolerance. The default is 10^-11
241  * @param rtolx Relative tolerance. The default is 10^-5
242  * Default parameter is 0.0, in which case rtolx is set equal to rtolf
243  * @param atolx absolute tolerance. The default is 10^-11
244  * Default parameter is 0.0, in which case atolx is set equal to atolf
245  */
246  void setTol(doublereal rtolf, doublereal atolf, doublereal rtolx = 0.0, doublereal atolx = 0.0);
247
248  //! Set the print level from the rootfinder
249  /*!
250  *
251  * 0 -> absolutely nothing is printed for a single time step.
252  * 1 -> One line summary per solve_nonlinear call
253  * 2 -> short description, points of interest: Table of nonlinear solve - one line per iteration
254  * 3 -> Table is included -> More printing per nonlinear iteration (default) that occurs during the table
255  * 4 -> Summaries of the nonlinear solve iteration as they are occurring -> table no longer printed
256  * 5 -> Algorithm information on the nonlinear iterates are printed out
257  * 6 -> Additional info on the nonlinear iterates are printed out
258  * 7 -> Additional info on the linear solve is printed out.
259  * 8 -> Info on a per iterate of the linear solve is printed out.
260  *
261  * @param printLvl integer value
262  */
263  void setPrintLvl(int printLvl);
264
265  //! Set the function behavior flag
266  /*!
267  * If this is true, the function is generally an increasing function of x.
268  * In particular, if the algorithm is seeking a higher value of f, it will look
269  * in the positive x direction.
270  *
271  * This type of function is needed because this algorithm must deal with regions of f(x) where
272  * f is not changing with x.
273  *
274  * @param value boolean value
275  */
276  void setFuncIsGenerallyIncreasing(bool value);
277
278  //! Set the function behavior flag
279  /*!
280  * If this is true, the function is generally a decreasing function of x.
281  * In particular, if the algorithm is seeking a higher value of f, it will look
282  * in the negative x direction.
283  *
284  * This type of function is needed because this algorithm must deal with regions of f(x) where
285  * f is not changing with x.
286  *
287  * @param value boolean value
288  */
289  void setFuncIsGenerallyDecreasing(bool value);
290
291  //! Set the minimum value of deltaX
292  /*!
293  * This sets the value of deltaXNorm_
294  *
295  * @param deltaXNorm
296  */
297  void setDeltaX(doublereal deltaXNorm);
298
299  //! Set the maximum value of deltaX
300  /*!
301  * This sets the value of deltaXMax_
302  *
303  * @param deltaX
304  */
305  void setDeltaXMax(doublereal deltaX);
306
307  //! Print the iteration history table
308  void printTable();
309
310 public:
311  //! Pointer to the residual function evaluator
313
314  //! Target value for the function. We seek the value of f that is equal to this value
315  doublereal m_funcTargetValue;
316
317  //! Absolute tolerance for the value of f
318  doublereal m_atolf;
319
320  //! Absolute tolerance for the value of x
321  doublereal m_atolx;
322
323  //! Relative tolerance for the value of f and x
324  doublereal m_rtolf;
325
326  //! Relative tolerance for the value of x
327  doublereal m_rtolx;
328
329  //! Maximum number of step sizes
330  doublereal m_maxstep;
331
332 protected:
333  //! Print level
334  /*!
335  * 0 No printing of any kind
336  * 1 Single print line indicating success or failure of the routine.
337  * 2 Summary table printed at the end of the routine, with a convergence history
338  * 3 Printouts during the iteration are added. Summary table is printed out at the end.
339  * if writeLogAllowed_ is turned on, a file is written out with the convergence history.
340  */
341  int printLvl;
342
343 public:
344  //! Boolean to turn on the possibility of writing a log file.
346
347 protected:
348  //! Delta X norm. This is the nominal value of deltaX that will be used by the program
349  doublereal DeltaXnorm_;
350
351  //! Boolean indicating whether DeltaXnorm_ has been specified by the user or not
353
354  //! Delta X Max. This is the maximum value of deltaX that will be used by the program
355  /*!
356  * Sometimes a large change in x causes problems.
357  */
358  doublereal DeltaXMax_;
359
360  //! Boolean indicating whether DeltaXMax_ has been specified by the user or not
362
363  //! Boolean indicating whether the function is an increasing with x
365
366  //! Boolean indicating whether the function is decreasing with x
368
369  //! Value of delta X that is needed for convergence
370  /*!
371  * X will be considered as converged if we are within deltaXConverged_ of the solution
372  * The default is zero.
373  */
374  doublereal deltaXConverged_;
375
376  //! Internal variable tracking largest x tried.
377  doublereal x_maxTried_;
378
379  //! Internal variable tracking f(x) of largest x tried.
380  doublereal fx_maxTried_;
381
382  //! Internal variable tracking smallest x tried.
383  doublereal x_minTried_;
384
385  //! Internal variable tracking f(x) of smallest x tried.
386  doublereal fx_minTried_;
387
388  //! Structure containing the iteration history
389  struct rfTable {
390  //@{
391  int its;
392  int TP_its;
393  double slope;
394  double xval;
395  double fval;
396  int foundPos;
397  int foundNeg;
398  double deltaXConverged;
399  double deltaFConverged;
400  double delX;
401
402  std::string reasoning;
403
404  void clear() {
405  its = 0;
406  TP_its = 0;
407  slope = -1.0E300;
408  xval = -1.0E300;
409  fval = -1.0E300;
410  reasoning = "";
411  };
412
413  rfTable() :
414  its(-2),
415  TP_its(0),
416  slope(-1.0E300),
417  xval(-1.0E300),
418  fval(-1.0E300),
419  foundPos(0),
420  foundNeg(0),
421  deltaXConverged(-1.0E300),
422  deltaFConverged(-1.0E300),
423  delX(-1.0E300),
424  reasoning("") {
425  };
426  //@}
427  };
428
429  //! Vector of iteration histories
430  std::vector<struct rfTable> rfHistory_;
431 };
432 }
433 #endif
void setDeltaX(doublereal deltaXNorm)
Set the minimum value of deltaX.
Definition: RootFind.cpp:1177
doublereal DeltaXnorm_
Delta X norm. This is the nominal value of deltaX that will be used by the program.
Definition: RootFind.h:349
doublereal fx_minTried_
Internal variable tracking f(x) of smallest x tried.
Definition: RootFind.h:386
bool FuncIsGenerallyDecreasing_
Boolean indicating whether the function is decreasing with x.
Definition: RootFind.h:367
doublereal x_maxTried_
Internal variable tracking largest x tried.
Definition: RootFind.h:377
std::vector< struct rfTable > rfHistory_
Vector of iteration histories.
Definition: RootFind.h:430
void setFuncIsGenerallyDecreasing(bool value)
Set the function behavior flag.
Definition: RootFind.cpp:1169
Structure containing the iteration history.
Definition: RootFind.h:389
int specifiedDeltaXnorm_
Boolean indicating whether DeltaXnorm_ has been specified by the user or not.
Definition: RootFind.h:352
RootFind & operator=(const RootFind &right)
Assignment operator.
Definition: RootFind.cpp:109
void printTable()
Print the iteration history table.
Definition: RootFind.cpp:1189
doublereal delXNonzero(doublereal x1) const
Calculate a deltaX from an input value of x.
Definition: RootFind.cpp:138
int printLvl
Print level.
Definition: RootFind.h:341
doublereal m_atolx
Absolute tolerance for the value of x.
Definition: RootFind.h:321
doublereal m_rtolx
Relative tolerance for the value of x.
Definition: RootFind.h:327
bool writeLogAllowed_
Boolean to turn on the possibility of writing a log file.
Definition: RootFind.h:345
doublereal m_rtolf
Relative tolerance for the value of f and x.
Definition: RootFind.h:324
doublereal DeltaXMax_
Delta X Max. This is the maximum value of deltaX that will be used by the program.
Definition: RootFind.h:358
doublereal x_minTried_
Internal variable tracking smallest x tried.
Definition: RootFind.h:383
void setTol(doublereal rtolf, doublereal atolf, doublereal rtolx=0.0, doublereal atolx=0.0)
Set the tolerance parameters for the rootfinder.
Definition: RootFind.cpp:1140
doublereal deltaXControlled(doublereal x2, doublereal x1) const
Calculate a controlled, nonzero delta between two numbers.
Definition: RootFind.cpp:157
Root finder for 1D problems.
Definition: RootFind.h:131
void setFuncIsGenerallyIncreasing(bool value)
Set the function behavior flag.
Definition: RootFind.cpp:1161
Virtual base class for DAE residual function evaluators.
Definition: ResidEval.h:33
doublereal m_funcTargetValue
Target value for the function. We seek the value of f that is equal to this value.
Definition: RootFind.h:315
bool FuncIsGenerallyIncreasing_
Boolean indicating whether the function is an increasing with x.
Definition: RootFind.h:364
doublereal fx_maxTried_
Internal variable tracking f(x) of largest x tried.
Definition: RootFind.h:380
ResidEval * m_residFunc
Pointer to the residual function evaluator.
Definition: RootFind.h:312
int solve(doublereal xmin, doublereal xmax, int itmax, doublereal &funcTargetValue, doublereal *xbest)
Using a line search method, find the root of a 1D function.
Definition: RootFind.cpp:184
void setDeltaXMax(doublereal deltaX)
Set the maximum value of deltaX.
Definition: RootFind.cpp:1183
doublereal m_maxstep
Maximum number of step sizes.
Definition: RootFind.h:330
doublereal delXMeaningful(doublereal x1) const
Calculate a deltaX from an input value of x.
Definition: RootFind.cpp:148
doublereal deltaXConverged_
Value of delta X that is needed for convergence.
Definition: RootFind.h:374
RootFind(ResidEval *resid)
Constructor for the object.
Definition: RootFind.cpp:56
void setPrintLvl(int printLvl)
Set the print level from the rootfinder.
Definition: RootFind.cpp:1156
doublereal m_atolf
Absolute tolerance for the value of f.
Definition: RootFind.h:318
int specifiedDeltaXMax_
Boolean indicating whether DeltaXMax_ has been specified by the user or not.
Definition: RootFind.h:361
bool theSame(doublereal x2, doublereal x1, doublereal factor=1.0) const
Function to decide whether two real numbers are the same or not.
Definition: RootFind.cpp:172
doublereal func(doublereal x)
Return the function value.
Definition: RootFind.cpp:1118