CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: FunctionNumDeriv.cc,v 1.4 2003/10/10 17:40:39 garren Exp $ 00003 // --------------------------------------------------------------------------- 00004 00005 #include "CLHEP/GenericFunctions/FunctionNumDeriv.hh" 00006 #include <assert.h> 00007 #include <cmath> // for pow() 00008 00009 namespace Genfun { 00010 FUNCTION_OBJECT_IMP(FunctionNumDeriv) 00011 00012 FunctionNumDeriv::FunctionNumDeriv(const AbsFunction *arg1, unsigned int index): 00013 _arg1(arg1->clone()), 00014 _wrtIndex(index) 00015 { 00016 } 00017 00018 FunctionNumDeriv::FunctionNumDeriv(const FunctionNumDeriv & right): 00019 _arg1(right._arg1->clone()), 00020 _wrtIndex(right._wrtIndex) 00021 { 00022 } 00023 00024 00025 FunctionNumDeriv::~FunctionNumDeriv() 00026 { 00027 delete _arg1; 00028 } 00029 00030 unsigned int FunctionNumDeriv::dimensionality() const { 00031 return _arg1->dimensionality(); 00032 } 00033 00034 #define ROBUST_DERIVATIVES 00035 #ifdef ROBUST_DERIVATIVES 00036 00037 double FunctionNumDeriv::f_x (double x) const { 00038 return (*_arg1)(x); 00039 } 00040 00041 00042 double FunctionNumDeriv::f_Arg (double x) const { 00043 _xArg [_wrtIndex] = x; 00044 return (*_arg1)(_xArg); 00045 } 00046 00047 00048 double FunctionNumDeriv::operator ()(double x) const 00049 { 00050 assert (_wrtIndex==0); 00051 return numericalDerivative ( &FunctionNumDeriv::f_x, x ); 00052 } 00053 00054 double FunctionNumDeriv::operator ()(const Argument & x) const 00055 { 00056 assert (_wrtIndex<x.dimension()); 00057 _xArg = x; 00058 double xx = x[_wrtIndex]; 00059 return numericalDerivative ( &FunctionNumDeriv::f_Arg, xx ); 00060 } 00061 00062 00063 double FunctionNumDeriv::numericalDerivative 00064 ( double (FunctionNumDeriv::*f)(double)const, double x ) const { 00065 00066 const double h0 = 5 * pow(2.0, -17); 00067 00068 const double maxErrorA = .0012; // These are the largest errors in steps 00069 const double maxErrorB = .0000026; // A, B consistent with 8-digit accuracy. 00070 00071 const double maxErrorC = .0003; // Largest acceptable validation discrepancy. 00072 00073 // This value of gives 8-digit accuracy for 1250 > curvature scale < 1/1250. 00074 00075 const int nItersMax = 6; 00076 int nIters; 00077 double bestError = 1.0E30; 00078 double bestAns = 0; 00079 00080 const double valFactor = pow(2.0, -16); 00081 00082 const double w = 5.0/8; 00083 const double wi2 = 64.0/25.0; 00084 const double wi4 = wi2*wi2; 00085 00086 double size = fabs((this->*f)(x)); 00087 if (size==0) size = pow(2.0, -53); 00088 00089 const double adjustmentFactor[nItersMax] = { 00090 1.0, 00091 pow(2.0, -17), 00092 pow(2.0, +17), 00093 pow(2.0, -34), 00094 pow(2.0, +34), 00095 pow(2.0, -51) }; 00096 00097 for ( nIters = 0; nIters < nItersMax; ++nIters ) { 00098 00099 double h = h0 * adjustmentFactor[nIters]; 00100 00101 // Step A: Three estimates based on h and two smaller values: 00102 00103 double A1 = ((this->*f)(x+h) - (this->*f)(x-h))/(2.0*h); 00104 // size = max(fabs(A1), size); 00105 if (fabs(A1) > size) size = fabs(A1); 00106 00107 double hh = w*h; 00108 double A2 = ((this->*f)(x+hh) - (this->*f)(x-hh))/(2.0*hh); 00109 // size = max(fabs(A2), size); 00110 if (fabs(A2) > size) size = fabs(A2); 00111 00112 hh *= w; 00113 double A3 = ((this->*f)(x+hh) - (this->*f)(x-hh))/(2.0*hh); 00114 // size = max(fabs(A3), size); 00115 if (fabs(A3) > size) size = fabs(A3); 00116 00117 if ( (fabs(A1-A2)/size > maxErrorA) || (fabs(A1-A3)/size > maxErrorA) ) { 00118 continue; 00119 } 00120 00121 // Step B: Two second-order estimates based on h h*w, from A estimates 00122 00123 double B1 = ( A2 * wi2 - A1 ) / ( wi2 - 1 ); 00124 double B2 = ( A3 * wi2 - A2 ) / ( wi2 - 1 ); 00125 if ( fabs(B1-B2)/size > maxErrorB ) { 00126 continue; 00127 } 00128 00129 // Step C: Third-order estimate, from B estimates: 00130 00131 double ans = ( B2 * wi4 - B1 ) / ( wi4 - 1 ); 00132 double err = fabs ( ans - B1 ); 00133 if ( err < bestError ) { 00134 bestError = err; 00135 bestAns = ans; 00136 } 00137 00138 // Validation estimate based on much smaller h value: 00139 00140 hh = h * valFactor; 00141 double val = ((this->*f)(x+hh) - (this->*f)(x-hh))/(2.0*hh); 00142 if ( fabs(val-ans)/size > maxErrorC ) { 00143 continue; 00144 } 00145 00146 // Having passed both apparent accuracy and validation, we are finished: 00147 break; 00148 } 00149 00150 return bestAns; 00151 00152 } 00153 #endif // ROBUST_DERIVATIVES 00154 00155 00156 00157 #ifdef SIMPLER_DERIVATIVES 00158 double FunctionNumDeriv::operator ()(double x) const 00159 { 00160 assert (_wrtIndex==0); 00161 const double h=1.0E-6; 00162 return ((*_arg1)(x+h) - (*_arg1)(x-h))/(2.0*h); 00163 } 00164 00165 double FunctionNumDeriv::operator ()(const Argument & x) const 00166 { 00167 assert (_wrtIndex<x.dimension()); 00168 const double h=1.0E-6; 00169 Argument x1=x, x0=x; 00170 x1[_wrtIndex] +=h; 00171 x0[_wrtIndex] -=h; 00172 return ((*_arg1)(x1) - (*_arg1)(x0))/(2.0*h); 00173 } 00174 #endif // SIMPLER_DERIVATIVES 00175 00176 } // namespace Genfun