deriv_cern.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 deriv_cern.h
24  \brief File defining \ref o2scl::deriv_cern
25 */
26 #ifndef O2SCL_DERIV_CERN_H
27 #define O2SCL_DERIV_CERN_H
28 
29 #include <o2scl/deriv.h>
30 #include <o2scl/funct.h>
31 #include <o2scl/string_conv.h>
32 #include <o2scl/err_hnd.h>
33 
34 #ifndef DOXYGEN_NO_O2NS
35 namespace o2scl {
36 #endif
37 
38  /** \brief Numerical differentiation routine (CERNLIB)
39 
40  This uses Romberg extrapolation to compute the
41  derivative with the finite-differencing formula
42  \f[
43  f^{\prime}(x) = [f(x+h)-f(x-h)]/(2h)
44  \f]
45 
46  If \ref deriv_base::verbose is greater than zero, then each iteration
47  prints out the extrapolation table, and if \ref deriv_base::verbose
48  is greater than 1, then a keypress is required at the end of
49  each iteration.
50 
51  For sufficiently difficult functions, the derivative computation
52  can fail, and will call the error handler and return zero with
53  zero error.
54 
55  Based on the CERNLIB routine DERIV, which was
56  based on \ref Rutishauser63 and is documented at
57  http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/d401/top.html
58 
59  An example demonstrating the usage of this class is
60  given in <tt>examples/ex_deriv.cpp</tt> and the \ref ex_deriv_sect .
61 
62  If \ref deriv_base::verbose is greater than zero, at each iteration
63  this class prints something similar to
64  \verbatim
65  deriv_cern, iteration: 1
66  (hh, a[k], derivative) list:
67  -4.193459e-05 4.387643e-14 8.775286e-01
68  -2.995402e-05 4.387792e-14 8.775585e-01
69  -1.048405e-05 4.387845e-14 8.775690e-01
70  -7.488654e-06 4.387882e-14 8.775765e-01
71  -2.621038e-06 4.387895e-14 8.775791e-01
72  -1.872173e-06 4.387905e-14 8.775810e-01
73  -6.552611e-07 4.387908e-14 8.775817e-01
74  -4.680438e-07 4.387910e-14 8.775821e-01
75  -1.638153e-07 4.387911e-14 8.775823e-01
76  \endverbatim
77  If \ref deriv_base::verbose is greater than 1, a keypress is required
78  after each iteration.
79 
80  \note Second and third derivatives are computed by naive nested
81  applications of the formula for the first derivative.
82  No uncertainty for these derivatives is provided.
83 
84  \future All of the coefficients appear to be fractions which
85  could be replaced with exact representation?
86  \future Record the number of function calls?
87 
88  \comment
89  - Maybe we should consider moving the table size to a template
90  parameter? (1/29/07 - Probably not, as we'd have to re-derive
91  the coefficients for sizes other than 10)
92  \endcomment
93  */
94  template<class func_t=funct11> class deriv_cern :
95  public deriv_base<func_t> {
96 
97  public:
98 
99  /// A scaling factor (default 1.0)
100  double delta;
101 
102  /// Extrapolation tolerance (default is \f$ 5 \times 10^{14} \f$)
103  double eps;
104 
105  deriv_cern() {
106  dx[0]=0.0256;
107  dx[1]=0.0192;
108  dx[2]=0.0128;
109  dx[3]=0.0096;
110  dx[4]=0.0064;
111  dx[5]=0.0048;
112  dx[6]=0.0032;
113  dx[7]=0.0024;
114  dx[8]=0.0016;
115  dx[9]=0.0012;
116 
117  w[1][1]=1.3333333333333333;
118  w[3][1]=1.0666666666666667;
119  w[5][1]=1.0158730158730159;
120  w[7][1]=1.0039215686274510;
121 
122  w[2][1]=3.3333333333333333e-1;
123  w[4][1]=6.6666666666666667e-2;
124  w[6][1]=1.5873015873015873e-2;
125  w[8][1]=3.9215686274509804e-3;
126 
127  w[0][2]=2.2857142857142857;
128  w[2][2]=1.1636363636363636;
129  w[4][2]=1.0364372469635628;
130  w[6][2]=1.0088669950738916;
131  w[8][2]=1.0022021042329337;
132 
133  w[1][2]=1.2857142857142857;
134  w[3][2]=1.6363636363636364e-1;
135  w[5][2]=3.6437246963562753e-2;
136  w[7][2]=8.8669950738916256e-3;
137  w[9][2]=2.2021042329336922e-3;
138 
139  w[0][3]=1.8000000000000000;
140  w[2][3]=1.1250000000000000;
141  w[4][3]=1.0285714285714286;
142  w[6][3]=1.0069930069930070;
143  w[8][3]=1.0017391304347826;
144 
145  w[1][3]=8.0000000000000000e-1;
146  w[3][3]=1.2500000000000000e-1;
147  w[5][3]=2.8571428571428571e-2;
148  w[7][3]=6.9930069930069930e-3;
149  w[9][3]=1.7391304347826087e-3;
150 
151  delta=1.0;
152  eps=5.0e-14;
153  }
154 
155  /** \brief Calculate the first derivative of \c func w.r.t. x and the
156  uncertainty
157  */
158  virtual int deriv_err(double x, func_t &func,
159  double &dfdx, double &err) {
160  return deriv_tlate<func_t>(x,func,dfdx,err);
161  }
162 
163  /// Return string denoting type ("deriv_cern")
164  virtual const char *type() { return "deriv_cern"; }
165 
166  protected:
167 
168 #ifndef DOXYGEN_INTERNAL
169 
170  /** \brief Internal template version of the derivative function
171  */
172  template<class func2_t> int deriv_tlate(double x, func2_t &func,
173  double &dfdx, double &err) {
174 
175  double t[10][10], a[10], del, hh;
176  bool lev[10]={1,0,1,0,1,0,1,0,1,0}, lmt;
177  int is, k, m;
178 
179  del=10.0*fabs(delta);
180  is=10;
181 
182  do {
183  is--;
184  del=del/10.0;
185 
186  if (is==0 || x+del*dx[9]==x) {
187  delta=del;
188  dfdx=0.0;
189  err=0.0;
190  std::string str="Calculation of derivative failed (is="+
191  itos(is)+" and del*dx[9]="+dtos(del*dx[9])+
192  ") in deriv_cern::deriv_tlate().";
193  O2SCL_ERR(str.c_str(),exc_efailed);
194  }
195 
196  for(k=0;k<=9;k++) {
197  hh=del*dx[k];
198  t[k][0]=(func(x+hh)-func(x-hh))/(hh+hh);
199  a[k]=t[k][0];
200  }
201 
202  if (a[0]>=a[9]) {
203  for(k=0;k<=9;k++) a[k]=-a[k];
204  }
205 
206  lmt=true;
207  for(k=1;k<=9;k++) {
208  hh=a[k-1]-a[k];
209  lmt=(lmt && (hh<=0.0 || fabs(hh)<=eps*fabs(a[k])));
210  }
211 
212  if (this->verbose>0) {
213  std::cout << "deriv_cern, iteration: " << 10-is << std::endl;
214  std::cout << "(hh, a[k], derivative) list: "
215  << std::endl;
216  for(k=1;k<=9;k++) {
217  std::cout << a[k-1]-a[k] << " " << eps*fabs(a[k]) << " "
218  << t[k][0] << std::endl;
219  }
220  std::cout << "Converged: " << lmt << std::endl;
221  if (this->verbose>1) {
222  char ch;
223  std::cin >> ch;
224  }
225  }
226 
227  } while (lmt==false);
228 
229  for(m=1;m<=9;m++) {
230  for(k=0;k<=9-m;k++) {
231  if (lev[m]) {
232  t[k][m]=w[m-1][1]*t[k+1][m-1]-w[m][1]*t[k][m-1];
233  } else if (lev[k]) {
234  t[k][m]=w[m-1][2]*t[k+1][m-1]-w[m][2]*t[k][m-1];
235  } else {
236  t[k][m]=w[m-1][3]*t[k+1][m-1]-w[m][3]*t[k][m-1];
237  }
238  }
239  }
240  dfdx=t[0][9];
241  if (dfdx!=0.0) err=(dfdx-t[0][8])/dfdx;
242  else err=0.0;
243  delta=del;
244 
245  return 0;
246  }
247 
248  /** \brief Calculate the first derivative of \c func w.r.t. x
249 
250  This is an internal version of deriv() which is used in
251  computing second and third derivatives
252  */
253  virtual int deriv_err_int(double x, funct11 &func, double &dfdx,
254  double &err) {
255  return deriv_tlate<>(x,func,dfdx,err);
256  }
257 
258  /// \name Storage for the fixed coefficients
259  //@{
260  double dx[10];
261  double w[10][4];
262  //@}
263 
264 #endif
265 
266  };
267 
268 #ifndef DOXYGEN_NO_O2NS
269 }
270 #endif
271 
272 #endif
273 
274 
virtual const char * type()
Return string denoting type ("deriv_cern")
Definition: deriv_cern.h:164
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
double delta
A scaling factor (default 1.0)
Definition: deriv_cern.h:100
std::function< double(double)> funct11
One-dimensional function typedef.
Definition: funct.h:44
Numerical differentiation routine (CERNLIB)
Definition: deriv_cern.h:94
generic failure
Definition: err_hnd.h:61
int deriv_tlate(double x, func2_t &func, double &dfdx, double &err)
Internal template version of the derivative function.
Definition: deriv_cern.h:172
double eps
Extrapolation tolerance (default is )
Definition: deriv_cern.h:103
Numerical differentiation base [abstract base].
Definition: deriv.h:63
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
int verbose
Output control.
Definition: deriv.h:152
virtual int deriv_err(double x, func_t &func, double &dfdx, double &err)
Calculate the first derivative of func w.r.t. x and the uncertainty.
Definition: deriv_cern.h:158
virtual int deriv_err_int(double x, funct11 &func, double &dfdx, double &err)
Calculate the first derivative of func w.r.t. x.
Definition: deriv_cern.h:253
std::string itos(int x)
Convert an integer to a string.

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