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

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