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

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