CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

FunctionNumDeriv.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id: FunctionNumDeriv.cc,v 1.4 2003/10/10 17:40:39 garren Exp $
3 // ---------------------------------------------------------------------------
4 
6 #include <assert.h>
7 #include <cmath> // for pow()
8 
9 namespace Genfun {
10 FUNCTION_OBJECT_IMP(FunctionNumDeriv)
11 
12 FunctionNumDeriv::FunctionNumDeriv(const AbsFunction *arg1, unsigned int index):
13  _arg1(arg1->clone()),
14  _wrtIndex(index)
15 {
16 }
17 
19  _arg1(right._arg1->clone()),
20  _wrtIndex(right._wrtIndex)
21 {
22 }
23 
24 
26 {
27  delete _arg1;
28 }
29 
30 unsigned int FunctionNumDeriv::dimensionality() const {
31  return _arg1->dimensionality();
32 }
33 
34 #define ROBUST_DERIVATIVES
35 #ifdef ROBUST_DERIVATIVES
36 
37 double FunctionNumDeriv::f_x (double x) const {
38  return (*_arg1)(x);
39 }
40 
41 
42 double FunctionNumDeriv::f_Arg (double x) const {
43  _xArg [_wrtIndex] = x;
44  return (*_arg1)(_xArg);
45 }
46 
47 
48 double FunctionNumDeriv::operator ()(double x) const
49 {
50  assert (_wrtIndex==0);
51  return numericalDerivative ( &FunctionNumDeriv::f_x, x );
52 }
53 
54 double FunctionNumDeriv::operator ()(const Argument & x) const
55 {
56  assert (_wrtIndex<x.dimension());
57  _xArg = x;
58  double xx = x[_wrtIndex];
59  return numericalDerivative ( &FunctionNumDeriv::f_Arg, xx );
60 }
61 
62 
63 double FunctionNumDeriv::numericalDerivative
64  ( double (FunctionNumDeriv::*f)(double)const, double x ) const {
65 
66  const double h0 = 5 * pow(2.0, -17);
67 
68  const double maxErrorA = .0012; // These are the largest errors in steps
69  const double maxErrorB = .0000026; // A, B consistent with 8-digit accuracy.
70 
71  const double maxErrorC = .0003; // Largest acceptable validation discrepancy.
72 
73  // This value of gives 8-digit accuracy for 1250 > curvature scale < 1/1250.
74 
75  const int nItersMax = 6;
76  int nIters;
77  double bestError = 1.0E30;
78  double bestAns = 0;
79 
80  const double valFactor = pow(2.0, -16);
81 
82  const double w = 5.0/8;
83  const double wi2 = 64.0/25.0;
84  const double wi4 = wi2*wi2;
85 
86  double size = fabs((this->*f)(x));
87  if (size==0) size = pow(2.0, -53);
88 
89  const double adjustmentFactor[nItersMax] = {
90  1.0,
91  pow(2.0, -17),
92  pow(2.0, +17),
93  pow(2.0, -34),
94  pow(2.0, +34),
95  pow(2.0, -51) };
96 
97  for ( nIters = 0; nIters < nItersMax; ++nIters ) {
98 
99  double h = h0 * adjustmentFactor[nIters];
100 
101  // Step A: Three estimates based on h and two smaller values:
102 
103  double A1 = ((this->*f)(x+h) - (this->*f)(x-h))/(2.0*h);
104 // size = max(fabs(A1), size);
105  if (fabs(A1) > size) size = fabs(A1);
106 
107  double hh = w*h;
108  double A2 = ((this->*f)(x+hh) - (this->*f)(x-hh))/(2.0*hh);
109 // size = max(fabs(A2), size);
110  if (fabs(A2) > size) size = fabs(A2);
111 
112  hh *= w;
113  double A3 = ((this->*f)(x+hh) - (this->*f)(x-hh))/(2.0*hh);
114 // size = max(fabs(A3), size);
115  if (fabs(A3) > size) size = fabs(A3);
116 
117  if ( (fabs(A1-A2)/size > maxErrorA) || (fabs(A1-A3)/size > maxErrorA) ) {
118  continue;
119  }
120 
121  // Step B: Two second-order estimates based on h h*w, from A estimates
122 
123  double B1 = ( A2 * wi2 - A1 ) / ( wi2 - 1 );
124  double B2 = ( A3 * wi2 - A2 ) / ( wi2 - 1 );
125  if ( fabs(B1-B2)/size > maxErrorB ) {
126  continue;
127  }
128 
129  // Step C: Third-order estimate, from B estimates:
130 
131  double ans = ( B2 * wi4 - B1 ) / ( wi4 - 1 );
132  double err = fabs ( ans - B1 );
133  if ( err < bestError ) {
134  bestError = err;
135  bestAns = ans;
136  }
137 
138  // Validation estimate based on much smaller h value:
139 
140  hh = h * valFactor;
141  double val = ((this->*f)(x+hh) - (this->*f)(x-hh))/(2.0*hh);
142  if ( fabs(val-ans)/size > maxErrorC ) {
143  continue;
144  }
145 
146  // Having passed both apparent accuracy and validation, we are finished:
147  break;
148  }
149 
150  return bestAns;
151 
152 }
153 #endif // ROBUST_DERIVATIVES
154 
155 
156 
157 #ifdef SIMPLER_DERIVATIVES
158 double FunctionNumDeriv::operator ()(double x) const
159 {
160  assert (_wrtIndex==0);
161  const double h=1.0E-6;
162  return ((*_arg1)(x+h) - (*_arg1)(x-h))/(2.0*h);
163 }
164 
165 double FunctionNumDeriv::operator ()(const Argument & x) const
166 {
167  assert (_wrtIndex<x.dimension());
168  const double h=1.0E-6;
169  Argument x1=x, x0=x;
170  x1[_wrtIndex] +=h;
171  x0[_wrtIndex] -=h;
172  return ((*_arg1)(x1) - (*_arg1)(x0))/(2.0*h);
173 }
174 #endif // SIMPLER_DERIVATIVES
175 
176 } // namespace Genfun