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