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
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).