43 #ifndef O2SCL_DERIV_GSL_H 44 #define O2SCL_DERIV_GSL_H 54 #include <gsl/gsl_deriv.h> 55 #include <gsl/gsl_errno.h> 57 #include <o2scl/deriv.h> 58 #include <o2scl/funct.h> 59 #include <o2scl/err_hnd.h> 61 #ifndef DOXYGEN_NO_O2NS 123 template<
class func_t=funct,
class fp_t=
double>
class deriv_gsl :
160 virtual int deriv_err(fp_t x, func_t &func, fp_t &dfdx, fp_t &err) {
161 return deriv_tlate<func_t>(x,func,dfdx,err);
165 virtual const char *
type() {
return "deriv_gsl"; }
167 #ifndef DOXYGEN_INTERNAL 174 fp_t &dfdx, fp_t &err) {
177 if (x==0.0) hh=1.0e-4;
178 else hh=1.0e-4*std::abs(x);
183 fp_t r_0, round, trunc, error;
187 while (fail && it_count<20) {
192 if (cret!=0) fail=
true;
196 if (fail==
false && round < trunc && (round > 0 && trunc > 0)) {
197 fp_t r_opt, round_opt, trunc_opt, error_opt;
203 h_opt=hh*std::pow(round/(2.0*trunc),1.0/3.0);
205 if (cret!=0) fail=
true;
206 error_opt=round_opt+trunc_opt;
211 if (fail==
false && error_opt < error &&
212 std::abs(r_opt-r_0) < 4.0*error) {
222 std::cout <<
"Function deriv_gsl::deriv_tlate out of range. " 223 <<
"Decreasing step." << std::endl;
228 if (fail==
true || it_count>=20) {
230 O2SCL_ERR2(
"Failed to find finite derivative in ",
246 (fp_t x,
funct &func, fp_t &dfdx, fp_t &err) {
247 return deriv_tlate<>(x,func,dfdx,err);
260 template<
class func2_t>
262 fp_t &abserr_round, fp_t &abserr_trunc,
265 fp_t fm1, fp1, fmh, fph;
267 fp_t eps=std::numeric_limits<fp_t>::epsilon();
276 std::cout <<
"deriv_gsl: " << std::endl;
277 std::cout <<
"step: " << hh << std::endl;
278 std::cout <<
"abscissas: " << x-hh/2.0 <<
" " << x-hh <<
" " 279 << x+hh/2.0 <<
" " << x+hh << std::endl;
280 std::cout <<
"ordinates: " << fm1 <<
" " << fmh <<
" " << fph <<
" " 284 if (!std::isfinite(fm1) ||
285 !std::isfinite(fp1) ||
286 !std::isfinite(fmh) ||
287 !std::isfinite(fph) ||
288 (func_max>0.0 && (std::abs(fm1)>func_max ||
289 std::abs(fp1)>func_max ||
290 std::abs(fmh)>func_max ||
291 std::abs(fph)>func_max))) {
295 fp_t r3=0.5*(fp1-fm1);
296 fp_t r5=(4.0/3.0)*(fph-fmh)-(1.0/3.0)*r3;
298 fp_t e3=(std::abs(fp1)+std::abs(fm1))*eps;
299 fp_t e5=2.0*(std::abs(fph)+std::abs(fmh))*eps+e3;
303 fp_t dy=std::max(std::abs(r3/hh),std::abs(r5/hh))*std::abs(x/hh)*eps;
313 abserr_trunc=std::abs((r5-r3)/hh);
315 abserr_round=std::abs(e5/hh)+dy;
318 std::cout <<
"res: " << result <<
" trc: " << abserr_trunc
319 <<
" rnd: " << abserr_round << std::endl;
333 #ifndef DOXYGEN_NO_O2NS The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
fp_t h_opt
The last value of the optimized stepsize.
virtual const char * type()
Return string denoting type ("deriv_gsl")
bool err_nonconv
If true, call the error handler if the routine does not "converge".
int deriv_tlate(fp_t x, func2_t &func, fp_t &dfdx, fp_t &err)
Internal template version of the derivative function.
virtual int deriv_err_int(fp_t x, funct &func, fp_t &dfdx, fp_t &err)
Internal version of calc_err() for second and third derivatives.
Numerical differentiation base [abstract base].
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
int central_deriv(fp_t x, fp_t hh, fp_t &result, fp_t &abserr_round, fp_t &abserr_trunc, func2_t &func)
Compute derivative using 5-point rule.
virtual int deriv_err(fp_t x, func_t &func, fp_t &dfdx, fp_t &err)
Calculate the first derivative of func w.r.t. x and uncertainty.
fp_t func_max
Maximum absolute value of function, or a negative value for no maximum (default -1) ...
Numerical differentiation (GSL)
int verbose
Output control.
std::function< double(double)> funct
One-dimensional function typedef.