deriv_gsl.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 /* deriv/deriv.c
24  *
25  * Copyright (C) 2004, 2007 Brian Gough
26  *
27  * This program is free software; you can redistribute it and/or modify
28  * it under the terms of the GNU General Public License as published by
29  * the Free Software Foundation; either version 3 of the License, or (at
30  * your option) any later version.
31  *
32  * This program is distributed in the hope that it will be useful, but
33  * WITHOUT ANY WARRANTY; without even the implied warranty of
34  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
35  * General Public License for more details.
36  *
37  * You should have received a copy of the GNU General Public License
38  * along with this program; if not, write to the Free Software
39  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
40  * 02110-1301, USA.
41  */
42 
43 #ifndef O2SCL_DERIV_GSL_H
44 #define O2SCL_DERIV_GSL_H
45 
46 /** \file deriv_gsl.h
47  \brief File defining \ref o2scl::deriv_gsl
48 */
49 
50 #include <iostream>
51 #include <cmath>
52 #include <limits>
53 
54 #include <gsl/gsl_deriv.h>
55 #include <gsl/gsl_errno.h>
56 
57 #include <o2scl/deriv.h>
58 #include <o2scl/funct.h>
59 #include <o2scl/err_hnd.h>
60 
61 #ifndef DOXYGEN_NO_O2NS
62 namespace o2scl {
63 #endif
64 
65  /** \brief Numerical differentiation (GSL)
66 
67  This class computes the numerical derivative of a function. The
68  stepsize \ref h should be specified before use. If similar
69  functions are being differentiated in succession, the user may
70  be able to increase the speed of later derivatives by setting
71  the new stepsize equal to the optimized stepsize from the
72  previous differentiation, by setting \ref h to \ref h_opt.
73 
74  The results will be incorrect for sufficiently difficult
75  functions or if the step size is not properly chosen.
76 
77  Some successive derivative computations can be made more
78  efficient by using the optimized stepsize in \ref
79  deriv_gsl::h_opt , which is set by the most recent last
80  derivative computation.
81 
82  If the function returns a non-finite value, or if \ref func_max
83  is greater than zero and the absolute value of the function is
84  larger than \ref func_max, then this class attempts to decrease
85  the step size by a factor of 10 in order to compute the
86  derivative. The class gives up after 20 reductions of the
87  step size.
88 
89  If \ref h is negative or zero, the initial step size is chosen
90  to be \f$ 10^{-4} |x| \f$ or if \f$x=0\f$, then the initial step
91  size is chosen to be \f$ 10^{-4} \f$ .
92 
93  Setting \ref deriv_base::verbose to a number greater than zero
94  results in output for each call to \ref central_deriv() which
95  looks like:
96  \verbatim
97  deriv_gsl:
98  step: 1.000000e-04
99  abscissas: 4.999500e-01 4.999000e-01 5.000500e-01 5.001000e-01
100  ordinates: 4.793377e-01 4.793816e-01 4.794694e-01 4.795132e-01
101  res: 8.775825e-01 trc: 1.462163e-09 rnd: 7.361543e-12
102  \endverbatim
103  where the last line contains the result (<tt>res</tt>), the
104  truncation error (<tt>trc</tt>) and the rounding error
105  (<tt>rnd</tt>). If \ref deriv_base::verbose is greater than 1, a
106  keypress is required after each iteration.
107 
108  If the function always returns a finite value, then computing
109  first derivatives requires either 1 or 2 calls to \ref
110  central_deriv() and thus either 4 or 8 function calls.
111 
112  \note Second and third derivatives are computed by naive nested
113  applications of the formula for the first derivative. No
114  uncertainty for these derivatives is provided.
115 
116  An example demonstrating the usage of this class is given in
117  <tt>examples/ex_deriv.cpp</tt> and the \ref ex_deriv_sect .
118 
119  \future Include the forward and backward GSL derivatives.
120  These would be useful for EOS classes which run in to
121  trouble for negative densities.
122  */
123  template<class func_t=funct11> class deriv_gsl :
124  public deriv_base<func_t> {
125 
126  public:
127 
128  deriv_gsl() {
129  h=0.0;
130  h_opt=0.0;
131  func_max=-1.0;
132  }
133 
134  virtual ~deriv_gsl() {}
135 
136  /** \brief Initial stepsize
137 
138  This should be specified before a call to deriv() or
139  deriv_err(). If it is less than or equal to zero, then \f$ x
140  10^{-4} \f$ will used, or if \c x is zero, then \f$ 10^{-4} \f$
141  will be used.
142  */
143  double h;
144 
145  /** \brief Maximum absolute value of function, or
146  a negative value for no maximum (default -1)
147  */
148  double func_max;
149 
150  /** \brief The last value of the optimized stepsize
151 
152  This is initialized to zero in the constructor and set by
153  deriv_err() to the most recent value of the optimized stepsize.
154  */
155  double h_opt;
156 
157  /** \brief Calculate the first derivative of \c func w.r.t. x and
158  uncertainty
159  */
160  virtual int deriv_err(double x, func_t &func, double &dfdx, double &err) {
161  return deriv_tlate<func_t>(x,func,dfdx,err);
162  }
163 
164  /// Return string denoting type ("deriv_gsl")
165  virtual const char *type() { return "deriv_gsl"; }
166 
167 #ifndef DOXYGEN_INTERNAL
168 
169  protected:
170 
171  /** \brief Internal template version of the derivative function
172  */
173  template<class func2_t> int deriv_tlate(double x, func2_t &func,
174  double &dfdx, double &err) {
175  double hh;
176  if (h<=0.0) {
177  if (x==0.0) hh=1.0e-4;
178  else hh=1.0e-4*fabs(x);
179  } else {
180  hh=h;
181  }
182 
183  double r_0, round, trunc, error;
184 
185  size_t it_count=0;
186  bool fail=true;
187  while (fail && it_count<20) {
188 
189  fail=false;
190 
191  int cret=central_deriv(x,hh,r_0,round,trunc,func);
192  if (cret!=0) fail=true;
193 
194  error=round+trunc;
195 
196  if (fail==false && round < trunc && (round > 0 && trunc > 0)) {
197  double r_opt, round_opt, trunc_opt, error_opt;
198 
199  /* Compute an optimised stepsize to minimize the total error,
200  using the scaling of the truncation error (O(h^2)) and
201  rounding error (O(1/h)). */
202 
203  h_opt=hh*pow(round/(2.0*trunc),1.0/3.0);
204  cret=central_deriv(x,h_opt,r_opt,round_opt,trunc_opt,func);
205  if (cret!=0) fail=true;
206  error_opt=round_opt+trunc_opt;
207 
208  /* Check that the new error is smaller, and that the new derivative
209  is consistent with the error bounds of the original estimate. */
210 
211  if (fail==false && error_opt < error &&
212  fabs (r_opt-r_0) < 4.0*error) {
213  r_0=r_opt;
214  error=error_opt;
215  }
216  }
217 
218  it_count++;
219  if (fail==true) {
220  hh/=10.0;
221  if (this->verbose>0) {
222  std::cout << "Function deriv_gsl::deriv_tlate out of range. "
223  << "Decreasing step." << std::endl;
224  }
225  }
226  }
227 
228  if (fail==true || it_count>=20) {
229  if (this->err_nonconv) {
230  O2SCL_ERR2("Failed to find finite derivative in ",
231  "deriv_gsl::deriv_tlate<>.",o2scl::exc_efailed);
232  }
233  return o2scl::exc_efailed;
234  }
235 
236  dfdx=r_0;
237  err=error;
238 
239  return 0;
240  }
241 
242  /** \brief Internal version of calc_err() for second
243  and third derivatives
244  */
245  virtual int deriv_err_int
246  (double x, funct11 &func, double &dfdx, double &err) {
247  return deriv_tlate<>(x,func,dfdx,err);
248  }
249 
250  /** \brief Compute derivative using 5-point rule
251 
252  Compute the derivative using the 5-point rule (x-h, x-h/2, x,
253  x+h/2, x+h) and the error using the difference between the
254  5-point and the 3-point rule (x-h,x,x+h). Note that the
255  central point is not used for either.
256 
257  This must be a class template because it is used by
258  both deriv_err() and deriv_err_int().
259  */
260  template<class func2_t>
261  int central_deriv(double x, double hh, double &result,
262  double &abserr_round, double &abserr_trunc,
263  func2_t &func) {
264 
265  double fm1, fp1, fmh, fph;
266 
267  double eps=std::numeric_limits<double>::epsilon();
268 
269  fm1=func(x-hh);
270  fp1=func(x+hh);
271 
272  fmh=func(x-hh/2);
273  fph=func(x+hh/2);
274 
275  if (this->verbose>0) {
276  std::cout << "deriv_gsl: " << std::endl;
277  std::cout << "step: " << hh << std::endl;
278  std::cout << "abscissas: " << x-hh/2 << " " << x-hh << " "
279  << x+hh/2 << " " << x+hh << std::endl;
280  std::cout << "ordinates: " << fm1 << " " << fmh << " " << fph << " "
281  << fp1 << std::endl;
282  }
283 
284  if (!std::isfinite(fm1) ||
285  !std::isfinite(fp1) ||
286  !std::isfinite(fmh) ||
287  !std::isfinite(fph) ||
288  (func_max>0.0 && (fabs(fm1)>func_max ||
289  fabs(fp1)>func_max ||
290  fabs(fmh)>func_max ||
291  fabs(fph)>func_max))) {
292  return 1;
293  }
294 
295  double r3=0.5*(fp1-fm1);
296  double r5=(4.0/3.0)*(fph-fmh)-(1.0/3.0)*r3;
297 
298  double e3=(fabs(fp1)+fabs(fm1))*eps;
299  double e5=2.0*(fabs(fph)+fabs(fmh))*eps+e3;
300 
301  /* The next term is due to finite precision in x+h=O (eps*x) */
302 
303  double dy=GSL_MAX(fabs(r3/hh),fabs(r5/hh))*fabs(x/hh)*eps;
304 
305  /* The truncation error in the r5 approximation itself is O(h^4).
306  However, for safety, we estimate the error from r5-r3, which is
307  O(h^2). By scaling h we will minimise this estimated error, not
308  the actual truncation error in r5.
309  */
310 
311  result=r5/hh;
312  /* Estimated truncation error O(h^2) */
313  abserr_trunc=fabs((r5-r3)/hh);
314  /* Rounding error (cancellations) */
315  abserr_round=fabs(e5/hh)+dy;
316 
317  if (this->verbose>0) {
318  std::cout << "res: " << result << " trc: " << abserr_trunc
319  << " rnd: " << abserr_round << std::endl;
320  if (this->verbose>1) {
321  char ch;
322  std::cin >> ch;
323  }
324  }
325 
326  return 0;
327  }
328 
329 #endif
330 
331  };
332 
333 #ifndef DOXYGEN_NO_O2NS
334 }
335 #endif
336 
337 #endif
338 
339 
340 
double h_opt
The last value of the optimized stepsize.
Definition: deriv_gsl.h:155
int central_deriv(double x, double hh, double &result, double &abserr_round, double &abserr_trunc, func2_t &func)
Compute derivative using 5-point rule.
Definition: deriv_gsl.h:261
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
std::function< double(double)> funct11
One-dimensional function typedef.
Definition: funct.h:44
virtual int deriv_err(double x, func_t &func, double &dfdx, double &err)
Calculate the first derivative of func w.r.t. x and uncertainty.
Definition: deriv_gsl.h:160
generic failure
Definition: err_hnd.h:61
double h
Initial stepsize.
Definition: deriv_gsl.h:143
Numerical differentiation base [abstract base].
Definition: deriv.h:63
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
bool err_nonconv
If true, call the error handler if the routine does not "converge".
Definition: deriv.h:97
double func_max
Maximum absolute value of function, or a negative value for no maximum (default -1) ...
Definition: deriv_gsl.h:148
Numerical differentiation (GSL)
Definition: deriv_gsl.h:123
virtual const char * type()
Return string denoting type ("deriv_gsl")
Definition: deriv_gsl.h:165
int verbose
Output control.
Definition: deriv.h:152
int deriv_tlate(double x, func2_t &func, double &dfdx, double &err)
Internal template version of the derivative function.
Definition: deriv_gsl.h:173
virtual int deriv_err_int(double x, funct11 &func, double &dfdx, double &err)
Internal version of calc_err() for second and third derivatives.
Definition: deriv_gsl.h:246

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