centered_difference.h
Go to the documentation of this file.00001 // -*- C++ -*- 00002 00008 #if !defined(__numerical_centered_difference_h__) 00009 #define __numerical_centered_difference_h__ 00010 00011 #include "../defs.h" 00012 00013 #include "../../ads/array/FixedArray.h" 00014 00015 #include <functional> 00016 #include <limits> 00017 00018 #include <cmath> 00019 00020 BEGIN_NAMESPACE_NUMERICAL 00021 00022 //----------------------------------------------------------------------------- 00030 // @{ 00031 00033 template <class Functor> 00034 inline 00035 void 00036 derivative_centered_difference 00037 ( const Functor& f, 00038 const typename Functor::argument_type x, 00039 typename Functor::result_type& deriv, 00040 const typename Functor::argument_type 00041 delta = std::pow( std::numeric_limits<typename Functor::argument_type>:: 00042 epsilon(), 1.0 / 3.0 ) ) 00043 { 00044 const typename Functor::argument_type p = x + delta; 00045 const typename Functor::argument_type n = x - delta; 00046 deriv = ( f( p ) - f( n ) ) / ( p - n ); 00047 } 00048 00050 template <class Functor> 00051 inline 00052 typename Functor::result_type 00053 derivative_centered_difference 00054 ( const Functor& f, 00055 const typename Functor::argument_type x, 00056 const typename Functor::argument_type 00057 delta = std::pow( std::numeric_limits<typename Functor::argument_type>:: 00058 epsilon(), 1.0 / 3.0 ) ) 00059 { 00060 typename Functor::result_type deriv; 00061 derivative_centered_difference( f, x, deriv, delta ); 00062 return deriv; 00063 } 00064 00066 template <class Functor> 00067 class DerivativeCenteredDifference 00068 : public std::unary_function<typename Functor::argument_type, 00069 typename Functor::result_type> 00070 { 00071 // 00072 // Types 00073 // 00074 00075 private: 00076 00077 typedef std::unary_function<typename Functor::argument_type, 00078 typename Functor::result_type> base_type; 00079 00080 public: 00081 00083 typedef Functor function_type; 00085 typedef typename base_type::argument_type argument_type; 00087 typedef typename base_type::result_type result_type; 00088 00089 private: 00090 00091 // 00092 // Member data. 00093 // 00094 00095 const Functor& _f; 00096 00097 argument_type _delta; 00098 00099 // 00100 // Not implemented. 00101 // 00102 00103 // Default constructor not implemented. 00104 DerivativeCenteredDifference(); 00105 00106 // Assignment operator not implemented. 00107 DerivativeCenteredDifference& 00108 operator=( const DerivativeCenteredDifference& x ); 00109 00110 public: 00111 00113 DerivativeCenteredDifference 00114 ( const function_type& f, 00115 const argument_type delta = 00116 std::pow( std::numeric_limits<argument_type>::epsilon(), 1.0 / 3.0 ) ) : 00117 _f( f ), 00118 _delta( delta ) 00119 {} 00120 00122 DerivativeCenteredDifference( const DerivativeCenteredDifference& x ) : 00123 _f( x._f ), 00124 _delta( x._delta ) 00125 {} 00126 00128 void 00129 operator()( const argument_type x, result_type& deriv ) const 00130 { 00131 derivative_centered_difference( _f, x, deriv, _delta ); 00132 } 00133 00135 result_type 00136 operator()( const argument_type x ) const 00137 { 00138 return derivative_centered_difference( _f, x, _delta ); 00139 } 00140 00142 argument_type 00143 delta() const 00144 { 00145 return _delta; 00146 } 00147 00149 void 00150 set_delta( const argument_type delta ) 00151 { 00152 _delta = delta; 00153 } 00154 }; 00155 00156 // @} 00157 00158 00159 //----------------------------------------------------------------------------- 00165 // @{ 00166 00168 template <int N, class Functor, typename T> 00169 inline 00170 void 00171 gradient_centered_difference 00172 ( const Functor& f, 00173 typename Functor::argument_type x, 00174 ads::FixedArray<N,typename Functor::result_type>& gradient, 00175 const T delta = std::pow( std::numeric_limits<T>::epsilon(), 1.0 / 3.0 ) ) 00176 { 00177 T temp, xp, xn; 00178 typename Functor::result_type fp, fn; 00179 for ( int n = 0; n != N; ++n ) { 00180 temp = x[n]; 00181 // Positive offset. 00182 xp = x[n] += delta; 00183 fp = f( x ); 00184 // Negative offset. 00185 x[n] = temp; 00186 xn = x[n] -= delta; 00187 fn = f( x ); 00188 // Reset 00189 x[n] = temp; 00190 // Centered difference. 00191 gradient[n] = ( fp - fn ) / ( xp - xn ); 00192 } 00193 } 00194 00196 00200 template <int N, class Functor, typename T> 00201 inline 00202 ads::FixedArray<N,typename Functor::result_type> 00203 gradient_centered_difference 00204 ( const Functor& f, 00205 typename Functor::argument_type x, 00206 const T delta = std::pow( std::numeric_limits<T>::epsilon(), 1.0 / 3.0 ) ) 00207 { 00208 ads::FixedArray<N,typename Functor::result_type> gradient; 00209 gradient_centered_difference( f, x, gradient, delta ); 00210 return gradient; 00211 } 00212 00214 template <int N, 00215 class Functor, 00216 typename T = typename Functor::argument_type::value_type> 00217 class GradientCenteredDifference : 00218 public std::unary_function 00219 < typename Functor::argument_type, 00220 ads::FixedArray<N,typename Functor::result_type> > 00221 { 00222 private: 00223 typedef std::unary_function 00224 < typename Functor::argument_type, 00225 ads::FixedArray<N,typename Functor::result_type> > 00226 base_type; 00227 typedef T number_type; 00228 00229 private: 00230 00231 // 00232 // Member data. 00233 // 00234 00235 const Functor& _f; 00236 00237 number_type _delta; 00238 00239 // 00240 // Not implemented. 00241 // 00242 00243 // Default constructor not implemented. 00244 GradientCenteredDifference(); 00245 00246 // Assignment operator not implemented. 00247 GradientCenteredDifference& 00248 operator=( const GradientCenteredDifference& x ); 00249 00250 public: 00251 00253 typedef Functor function_type; 00255 typedef typename base_type::argument_type argument_type; 00257 typedef typename base_type::result_type result_type; 00258 00260 GradientCenteredDifference 00261 ( const function_type& f, 00262 const number_type delta = 00263 std::pow( std::numeric_limits<T>::epsilon(), 1.0 / 3.0 ) ) : 00264 _f( f ), 00265 _delta( delta ) 00266 {} 00267 00269 GradientCenteredDifference( const GradientCenteredDifference& x ) : 00270 _f( x._f ), 00271 _delta( x._delta ) 00272 {} 00273 00275 void 00276 operator()( const argument_type x, result_type& deriv ) const 00277 { 00278 gradient_centered_difference( _f, x, deriv, _delta ); 00279 } 00280 00282 result_type 00283 operator()( const argument_type x ) const 00284 { 00285 return gradient_centered_difference<N>( _f, x, _delta ); 00286 } 00287 00289 number_type 00290 delta() const 00291 { 00292 return _delta; 00293 } 00294 00296 void 00297 set_delta( const number_type delta ) 00298 { 00299 _delta = delta; 00300 } 00301 }; 00302 00303 // @} 00304 00305 END_NAMESPACE_NUMERICAL 00306 00307 #endif
Generated on Thu Jun 30 02:14:59 2016 for Numerical Algorithms Package by
1.6.3