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