root_stef.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 /* roots/steffenson.c
24  *
25  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Reid Priedhorsky,
26  * Brian Gough
27  *
28  * This program is free software; you can redistribute it and/or modify
29  * it under the terms of the GNU General Public License as published by
30  * the Free Software Foundation; either version 3 of the License, or (at
31  * your option) any later version.
32  *
33  * This program is distributed in the hope that it will be useful, but
34  * WITHOUT ANY WARRANTY; without even the implied warranty of
35  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
36  * General Public License for more details.
37  *
38  * You should have received a copy of the GNU General Public License
39  * along with this program; if not, write to the Free Software
40  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
41  * 02110-1301, USA.
42  */
43 
44 #ifndef O2SCL_ROOT_STEF_H
45 #define O2SCL_ROOT_STEF_H
46 
47 /** \file root_stef.h
48  \brief File defining \ref o2scl::root_stef
49 */
50 
51 #include <string>
52 
53 #include <o2scl/misc.h>
54 #include <o2scl/root.h>
55 
56 #ifndef DOXYGEN_NO_O2NS
57 namespace o2scl {
58 #endif
59 
60  /** \brief Steffenson equation solver (GSL)
61 
62  This is Newton's method with an Aitken "delta-squared"
63  acceleration of the iterates. This can improve the convergence
64  on multiple roots where the ordinary Newton algorithm is slow.
65 
66  Defining the next iteration with
67  \f[
68  x_{i+1} = x_i - f(x_i) / f^{\prime}(x_i)
69  \f]
70  the accelerated value is
71  \f[
72  x_{\mathrm{acc},i} = x_i - (x_{i+1}-x_i)^2 / (x_{i+2} -
73  2 x_{i+1} + x_i)
74  \f]
75  We can only use the accelerated estimate after three iterations,
76  and use the unaccelerated value until then.
77 
78  This class finds a root of a function a derivative. If the
79  derivative is not analytically specified, it is most likely
80  preferable to use of the alternatives, \ref o2scl::root_brent_gsl,
81  \ref o2scl::root_bkt_cern, or \ref o2scl::root_cern. The function
82  solve_de() performs the solution automatically, and a
83  lower-level GSL-like interface with set() and iterate() is also
84  provided.
85 
86  By default, this solver compares the present value of the root
87  (\f$ \mathrm{root} \f$) to the previous value (\f$ \mathrm{x}
88  \f$), and returns success if \f$ | \mathrm{root} - \mathrm{x} |
89  < \mathrm{tol} \f$, where \f$ \mathrm{tol} = \mathrm{tol\_abs} +
90  \mathrm{tol\_rel2}~\mathrm{root} \f$ .
91 
92  If \ref test_residual is set to true, then the solver
93  additionally requires that the absolute value of the function is
94  less than \ref root::tol_rel.
95 
96  The original variable \c x_2 has been removed as it was unused
97  in the original GSL code.
98 
99  See the \ref onedsolve_subsect section of the User's guide for
100  general information about \o2 solvers.
101 
102  \future There's some extra copying here which can probably
103  be removed.
104  \future Compare directly to GSL.
105  \future This can probably be modified to shorten the step
106  if the function goes out of bounds as in exc_mroot_hybrids.
107 
108  */
109  template<class func_t=funct11, class dfunc_t=func_t> class root_stef :
110  public root_de<func_t,dfunc_t> {
111 
112  protected:
113 
114  /// Function value
115  double f;
116 
117  /// Derivative value
118  double df;
119 
120  /// Previous value of root
121  double x_1;
122 
123  /// Root
124  double x;
125 
126  /// Number of iterations
127  int count;
128 
129  /// The function to solve
130  func_t *fp;
131 
132  /// The derivative
133  dfunc_t *dfp;
134 
135  public:
136 
137  root_stef() {
138  test_residual=false;
139  tol_rel2=1.0e-12;
140  }
141 
142  /// Return the type, \c "root_stef".
143  virtual const char *type() { return "root_stef"; }
144 
145  /// The present solution estimate
146  double root;
147 
148  /** \brief The relative tolerance for subsequent solutions
149  (default \f$ 10^{-12} \f$)
150  */
151  double tol_rel2;
152 
153  /** \brief Perform an iteration
154 
155  After a successful iteration, \ref root contains the
156  most recent value of the root.
157  */
158  int iterate() {
159 
160  double x_new, f_new, df_new;
161 
162  double x_1t=x_1;
163  double xt=x;
164 
165  if (df == 0.0) {
166  O2SCL_CONV_RET("Derivative is zero in root_stef::iterate().",
168  }
169 
170  x_new=xt-(f/df);
171 
172  // It is important that the derivative be evaluated first here,
173  // because we want the last function evaluation to be an
174  // evaluation for the returned root
175  df_new=(*dfp)(x_new);
176  f_new=(*fp)(x_new);
177 
178  x_1=xt;
179  x=x_new;
180 
181  f=f_new;
182  df=df_new;
183 
184  if (!std::isfinite(f_new)) {
185  std::string str="Function not finite (returned "+dtos(f_new)+
186  ") in root_stef::iterate().";
187  O2SCL_ERR(str.c_str(),o2scl::exc_ebadfunc);
188  }
189 
190  if (count<3) {
191 
192  root=x_new;
193  count++;
194 
195  } else {
196 
197  double u=(xt-x_1t);
198  double v=(x_new-2*xt+x_1t);
199 
200  if (v == 0) {
201  // Avoid division by zero
202  root=x_new;
203  } else {
204  // Accelerated value
205  root=x_1t-u*u/v;
206  }
207  }
208 
209  if (!std::isfinite(df_new)) {
210  std::string str="Derivative not finite (returned "+dtos(df_new)+
211  ") in root_stef::iterate().";
212  O2SCL_ERR(str.c_str(),o2scl::exc_ebadfunc);
213  }
214 
215  return o2scl::success;
216  }
217 
218  /** \brief Solve \c func using \c x as an initial
219  guess using derivatives \c df.
220  */
221  virtual int solve_de(double &xx, func_t &fun, dfunc_t &dfun) {
222 
223  int status1=success, status2=gsl_continue, iter=0;
224 
225  set(fun,dfun,xx);
226 
227  while (status1==success && status2==gsl_continue &&
228  iter<this->ntrial) {
229  iter++;
230 
231  status1=iterate();
232 
233  // Compare present value to previous value
234  status2=gsl_root_test_delta(root,xx,this->tol_abs,tol_rel2);
235 
236  if (test_residual && status2==success) {
237  double y;
238  y=fun(root);
239  if (fabs(y)>=this->tol_rel) status2=gsl_continue;
240  }
241 
242  if (this->verbose>0) {
243  double fval;
244  fval=fun(root);
245  this->print_iter(root,fval,iter,fabs(root-xx),this->tol_abs*root,
246  "root_stef");
247  }
248  xx=root;
249  }
250 
251  this->last_ntrial=iter;
252 
253  if (status1!=success || status2!=success) {
254  int ret=o2scl::err_hnd->get_errno();
255  return ret;
256  }
257  if (iter>=this->ntrial) {
258  std::string str=((std::string)"Function solve_de() exceeded the ")+
259  "maximum number of iterations, "+itos(this->ntrial)+".";
260  O2SCL_CONV_RET(str.c_str(),exc_emaxiter,this->err_nonconv);
261  }
262 
263  return o2scl::success;
264  }
265 
266  /// True if we should test the residual also (default false)
268 
269  /** \brief Set the information for the solver
270 
271  Set the function, the derivative, the initial guess and
272  the parameters.
273  */
274  void set(func_t &fun, dfunc_t &dfun, double guess) {
275 
276  fp=&fun;
277  dfp=&dfun;
278  root=guess;
279 
280  df=dfun(root);
281  f=fun(root);
282 
283  x=root;
284  x_1=0.0;
285  count=1;
286 
287  return;
288  }
289 
290  };
291 
292 #ifndef DOXYGEN_NO_O2NS
293 }
294 #endif
295 
296 #endif
297 
int count
Number of iterations.
Definition: root_stef.h:127
dfunc_t * dfp
The derivative.
Definition: root_stef.h:133
double root
The present solution estimate.
Definition: root_stef.h:146
int ntrial
Maximum number of iterations (default 100)
Definition: root.h:76
tried to divide by zero
Definition: err_hnd.h:75
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
err_hnd_type * err_hnd
The global error handler pointer.
double f
Function value.
Definition: root_stef.h:115
double tol_abs
The minimum allowable stepsize (default )
Definition: root.h:70
exceeded max number of iterations
Definition: err_hnd.h:73
iteration has not converged
Definition: err_hnd.h:51
One-dimensional with solver with derivatives [abstract base].
Definition: root.h:386
int iterate()
Perform an iteration.
Definition: root_stef.h:158
virtual const char * type()
Return the type, "root_stef".
Definition: root_stef.h:143
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
bool test_residual
True if we should test the residual also (default false)
Definition: root_stef.h:267
virtual int solve_de(double &xx, func_t &fun, dfunc_t &dfun)
Solve func using x as an initial guess using derivatives df.
Definition: root_stef.h:221
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
double x_1
Previous value of root.
Definition: root_stef.h:121
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
Steffenson equation solver (GSL)
Definition: root_stef.h:109
func_t * fp
The function to solve.
Definition: root_stef.h:130
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
double df
Derivative value.
Definition: root_stef.h:118
double x
Root.
Definition: root_stef.h:124
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 tol_rel2
The relative tolerance for subsequent solutions (default )
Definition: root_stef.h:151
virtual int get_errno() const =0
Return the last error number.
std::string itos(int x)
Convert an integer to a string.
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).