root.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file root.h
24  \brief File for one-dimensional solver base class \ref o2scl::root
25 */
26 #ifndef O2SCL_ROOT_H
27 #define O2SCL_ROOT_H
28 
29 #include <iostream>
30 #include <cmath>
31 #include <o2scl/err_hnd.h>
32 #include <o2scl/funct.h>
33 #include <o2scl/misc.h>
34 
35 #ifndef DOXYGEN_NO_O2NS
36 namespace o2scl {
37 #endif
38 
39  /** \brief One-dimensional solver [abstract base]
40 
41  See the \ref onedsolve_subsect section of the User's guide for
42  general information about \o2 solvers.
43 
44  \future Maybe consider allowing the user to specify
45  the stream to which 'verbose' information is sent.
46  */
47  template<class func_t=funct11, class dfunc_t=func_t> class root {
48 
49  public:
50 
51  root() {
52  ntrial=100;
53  tol_rel=1.0e-8;
54  verbose=0;
55  tol_abs=1.0e-12;
56  last_ntrial=0;
57  err_nonconv=true;
58  }
59 
60  virtual ~root() {}
61 
62  /** \brief The maximum value of the functions for success
63  (default \f$ 10^{-8} \f$ )
64  */
65  double tol_rel;
66 
67  /** \brief The minimum allowable stepsize
68  (default \f$ 10^{-12} \f$ )
69  */
70  double tol_abs;
71 
72  /// Output control (default 0)
73  int verbose;
74 
75  /// Maximum number of iterations (default 100)
76  int ntrial;
77 
78  /** \brief If true, call the error handler if the solver does
79  not converge (default true)
80  */
82 
83  /// The number of iterations used in the most recent solve
85 
86  /// Return the type, \c "root".
87  virtual const char *type() { return "root"; }
88 
89  /** \brief Print out iteration information.
90 
91  Depending on the value of the variable verbose, this prints
92  out the iteration information. If verbose=0, then no
93  information is printed, while if verbose>1, then after each
94  iteration, the present values of \c x and \c y are
95  output to std::cout along with the iteration number. If
96  verbose>=2 then each iteration waits for a character before
97  continuing.
98  */
99  virtual int print_iter(double x, double y, int iter, double value=0.0,
100  double limit=0.0, std::string comment="") {
101  if (verbose<=0) return success;
102 
103  char ch;
104 
105  std::cout << comment << " Iteration: " << iter << std::endl;
106  if (x<0) std::cout << x << " ";
107  else std::cout << " " << x << " ";
108  if (y<0) std::cout << y << " ";
109  else std::cout << " " << y << " ";
110  if (value<0) std::cout << value << " ";
111  else std::cout << " " << value << " ";
112  if (limit<0) std::cout << limit << std::endl;
113  else std::cout << " " << limit << std::endl;
114  if (verbose>1) {
115  std::cout << "Press a key and type enter to continue. ";
116  std::cin >> ch;
117  }
118 
119  return success;
120  }
121 
122  /** \brief Solve \c func using \c x as an initial guess
123  */
124  virtual int solve(double &x, func_t &func)=0;
125 
126  /** \brief Solve \c func in region \f$ x_1<x<x_2 \f$
127  returning \f$ x_1 \f$ .
128  */
129  virtual int solve_bkt(double &x1, double x2, func_t &func) {
130  return solve(x1,func);
131  }
132 
133  /** \brief Solve \c func using \c x as an initial
134  guess using derivatives \c df.
135  */
136  virtual int solve_de(double &x, func_t &func, dfunc_t &df) {
137  return solve(x,func);
138  }
139 
140  };
141 
142  /** \brief One-dimensional bracketing solver [abstract base]
143 
144  */
145  template<class func_t=funct11, class dfunc_t=func_t> class root_bkt :
146  public root<func_t,dfunc_t> {
147 
148  public:
149 
150  root_bkt() {
151  bracket_step=0.0;
152  bracket_min=0.0;
153  bracket_iters=20;
154  }
155 
156  virtual ~root_bkt() {}
157 
158  /** \brief The relative stepsize for automatic bracketing
159  (default value is zero)
160 
161  If this is less than or equal to zero, then
162  \f$ 10^{-4} \f$, will be used.
163  */
164  double bracket_step;
165 
166  /** \brief The minimum stepsize for automatic bracketing
167  (default zero)
168  */
169  double bracket_min;
170 
171  /// The number of iterations in attempt to bracket root (default 10)
173 
174  /// Return the type, \c "root_bkt".
175  virtual const char *type() { return "root_bkt"; }
176 
177  /** \brief Solve \c func in region \f$ x_1<x<x_2 \f$
178  returning \f$ x_1 \f$ .
179  */
180  virtual int solve_bkt(double &x1, double x2, func_t &func)=0;
181 
182  /** \brief Solve \c func using \c x as an initial guess
183 
184  This has been updated to automatically bracket functions
185  which are undefined away from the initial guess, but there
186  are still various ways in which the bracketing algorithm
187  can fail.
188 
189  \future Return early if the bracketing procedure finds a root
190  early?
191  */
192  virtual int solve(double &x, func_t &func) {
193 
194  // Internal version of bracket step
195  double bstep_int;
196 
197  // Set value of bstep_int
198  if (bracket_step<=0.0) {
199  bstep_int=fabs(x)*1.0e-4;
200  if (bstep_int<bracket_min) bstep_int=bracket_min;
201  if (bstep_int<=0.0) bstep_int=1.0e-4;
202  } else {
203  bstep_int=fabs(x)*bracket_step;
204  if (bstep_int<bracket_min) bstep_int=bracket_min;
205  if (bstep_int<=0.0) bstep_int=bracket_step;
206  }
207 
208  double x2=0.0, df, fx, fx2;
209  size_t i=0;
210  bool done=false;
211 
212  // Use function to try to bracket a root
213  while(done==false && i<bracket_iters) {
214 
215  // -----------------------------------------------------------
216  // Phase 1: Find a step which gives a finite function value and
217  // a non-zero derivative
218 
219  fx=func(x);
220  if (!std::isfinite(fx)) {
221  O2SCL_ERR2("First function value not finite in ",
222  "root_bkt::solve().",exc_ebadfunc);
223  }
224 
225  x2=x+bstep_int;
226  fx2=func(x2);
227  df=(fx2-fx)/bstep_int;
228  if (this->verbose>0) {
229  std::cout << "P1: x,x2,fx,fx2,df: " << x << " " << x2 << " "
230  << fx << " " << fx2 << " " << df << std::endl;
231  }
232 
233  // If the stepsize made the function not finite or the
234  // derivative zero, adjust accordingly
235  size_t j=0;
236  double step_phase1=bstep_int;
237  while ((!std::isfinite(fx2) || df==0.0) && j<bracket_iters*2) {
238  // Alternate between flipping the sign and making it smaller
239  if (j%2==0) step_phase1*=-1.0;
240  else step_phase1/=2.0;
241  x2=x+step_phase1;
242  fx2=func(x2);
243  df=(fx2-fx)/step_phase1;
244  if (this->verbose>1) {
245  std::cout << "P1: x,x2,fx,fx2,df: " << x << " " << x2 << " "
246  << fx << " " << fx2 << " " << df << std::endl;
247  }
248  j++;
249  }
250 
251  // If we failed to compute the derivative
252  if (df==0.0) {
253  O2SCL_CONV_RET("Failed to bracket (df==0) in root_bkt::solve().",
255  }
256  if (!std::isfinite(fx2)) {
257  O2SCL_CONV2_RET("Failed to bracket (f2 not finite, phase 1) in ",
258  "root_bkt::solve().",
260  }
261 
262  // Output the current iteration information
263  if (this->verbose>0) {
264  std::cout << "root_bkt::solve(): Iteration " << i << " of "
265  << bracket_iters << "." << std::endl;
266  std::cout << "\tx1: " << x << " f(x1): " << fx << std::endl;
267  std::cout << "\tx2: " << x2 << " f(x2): " << fx2 << std::endl;
268  }
269 
270  // -----------------------------------------------------------
271  // Phase 2: Create a new value which may provide a bracket
272 
273  double step_phase2=2.0;
274  x2=x-step_phase2*fx/df;
275  fx2=func(x2);
276  if (this->verbose>0) {
277  std::cout << "P2: x,x2,fx,fx2: " << x << " " << x2 << " "
278  << fx << " " << fx2 << std::endl;
279  }
280 
281  // Adjust if the function value is not finite
282  size_t k=0;
283  while (!std::isfinite(fx2) && k<bracket_iters) {
284  step_phase2/=2.0;
285  x2=x-step_phase2*fx/df;
286  fx2=func(x2);
287  k++;
288  if (this->verbose>0) {
289  std::cout << "P2: x,x2,fx,fx2: " << x << " " << x2 << " "
290  << fx << " " << fx2 << std::endl;
291  }
292  }
293 
294  // Throw if we failed
295  if (!std::isfinite(fx2)) {
296  O2SCL_CONV2_RET("Failed to bracket (f2 not finite, phase 2) in ",
297  "root_bkt::solve().",
299  }
300 
301  // Continue verbose output
302  if (this->verbose>0) {
303  std::cout << "\tx2: " << x2 << " f(x2): " << fx2 << std::endl;
304  if (this->verbose>1) {
305  std::cout << "Press a key and type enter to continue. ";
306  char ch;
307  std::cin >> ch;
308  }
309  }
310 
311  // -----------------------------------------------------------
312  // See if we're done
313 
314  if (fx*fx2<0.0) {
315  done=true;
316  } else {
317  // The step didn't cause a sign flip. Update 'x' to move
318  // closer to the purported root.
319  x=(x2+x)/2.0;
320  }
321 
322  // Continue for the next bracketing iteration
323  i++;
324  }
325 
326  // Throw if we failed
327  if (done==false) {
328  O2SCL_CONV_RET("Failed to bracket (iters>max) in root_bkt::solve().",
330  }
331 
332  // Go to the bracketing solver
333  if (this->verbose>0) {
334  std::cout << "root_bkt::solve(): Going to solve_bkt()."
335  << std::endl;
336  }
337  return solve_bkt(x,x2,func);
338  }
339 
340  /** \brief Solve \c func using \c x as an initial
341  guess using derivatives \c df.
342  */
343  virtual int solve_de(double &x, func_t &func, dfunc_t &df) {
344 
345  double x2=0.0, dx, fx, fx2;
346  int i=0;
347  bool done=false;
348 
349  // Use derivative information to try to bracket a root
350  while(done==false && i<10) {
351 
352  fx=func(x);
353  dx=df(x);
354 
355  x2=x-2.0*fx/dx;
356 
357  fx2=func(x2);
358 
359  if (fx*fx2<0.0) {
360  done=true;
361  } else {
362  x=(x2+x)/2.0;
363  }
364  i++;
365  }
366 
367  if (done==false) {
368  O2SCL_CONV_RET("Failed to bracket function in root_bkt::solve_de().",
370  }
371 
372  return solve_bkt(x,x2,func);
373  }
374 
375  };
376 
377  /** \brief One-dimensional with solver with derivatives [abstract base]
378 
379  \note At the moment, the functions solve() and solve_bkt()
380  are not implemented for derivative solvers.
381 
382  \future Implement the functions solve() and solve_bkt()
383  for derivative solvers.
384  */
385  template<class func_t=funct11, class dfunc_t=func_t>
386  class root_de : public root<func_t> {
387 
388  public:
389 
390  root_de() {
391  }
392 
393  virtual ~root_de() {}
394 
395  /// Return the type, \c "root_de".
396  virtual const char *type() { return "root_de"; }
397 
398  /** \brief Solve \c func in region \f$ x_1<x<x_2 \f$
399  returning \f$ x_1 \f$ .
400  */
401  virtual int solve_bkt(double &x1, double x2, func_t &func) {
402  O2SCL_ERR("Function solve_bkt() not implemented.",exc_eunimpl);
403  return 0;
404  }
405 
406  /** \brief Solve \c func using \c x as an initial guess
407  */
408  virtual int solve(double &x, func_t &func) {
409  O2SCL_ERR("Function solve() not implemented.",exc_eunimpl);
410  return 0;
411  }
412 
413  /** \brief Solve \c func using \c x as an initial
414  guess using derivatives \c df.
415  */
416  virtual int solve_de(double &x, func_t &func, dfunc_t &df)=0;
417 
418  };
419 
420 #ifndef DOXYGEN_NO_O2NS
421 }
422 #endif
423 
424 #endif
425 
virtual int solve(double &x, func_t &func)
Solve func using x as an initial guess.
Definition: root.h:408
int ntrial
Maximum number of iterations (default 100)
Definition: root.h:76
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
Definition: err_hnd.h:292
virtual const char * type()
Return the type, "root_de".
Definition: root.h:396
double tol_abs
The minimum allowable stepsize (default )
Definition: root.h:70
virtual int solve(double &x, func_t &func)=0
Solve func using x as an initial guess.
exceeded max number of iterations
Definition: err_hnd.h:73
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
Definition: err_hnd.h:298
double bracket_min
The minimum stepsize for automatic bracketing (default zero)
Definition: root.h:169
One-dimensional with solver with derivatives [abstract base].
Definition: root.h:386
One-dimensional bracketing solver [abstract base].
Definition: root.h:145
virtual int solve_bkt(double &x1, double x2, func_t &func)
Solve func in region returning .
Definition: root.h:129
requested feature not (yet) implemented
Definition: err_hnd.h:99
double tol_rel
The maximum value of the functions for success (default )
Definition: root.h:65
bool err_nonconv
If true, call the error handler if the solver does not converge (default true)
Definition: root.h:81
virtual int solve_de(double &x, func_t &func, dfunc_t &df)
Solve func using x as an initial guess using derivatives df.
Definition: root.h:343
One-dimensional solver [abstract base].
Definition: root.h:47
virtual int print_iter(double x, double y, int iter, double value=0.0, double limit=0.0, std::string comment="")
Print out iteration information.
Definition: root.h:99
virtual const char * type()
Return the type, "root".
Definition: root.h:87
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
size_t bracket_iters
The number of iterations in attempt to bracket root (default 10)
Definition: root.h:172
virtual int solve(double &x, func_t &func)
Solve func using x as an initial guess.
Definition: root.h:192
virtual int solve_de(double &x, func_t &func, dfunc_t &df)
Solve func using x as an initial guess using derivatives df.
Definition: root.h:136
virtual int solve_bkt(double &x1, double x2, func_t &func)
Solve func in region returning .
Definition: root.h:401
problem with user-supplied function
Definition: err_hnd.h:69
int last_ntrial
The number of iterations used in the most recent solve.
Definition: root.h:84
double bracket_step
The relative stepsize for automatic bracketing (default value is zero)
Definition: root.h:164
virtual const char * type()
Return the type, "root_bkt".
Definition: root.h:175
static const double x2[5]
Definition: inte_qng_gsl.h:66
static const double x1[5]
Definition: inte_qng_gsl.h:48
int verbose
Output control (default 0)
Definition: root.h:73
Success.
Definition: err_hnd.h:47

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).