Logo ROOT   6.08/07
Reference Guide
RichardsonDerivator.cxx
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Authors: David Gonzalez Maline 01/2008
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2006 , LCG ROOT MathLib Team *
7  * *
8  * *
9  **********************************************************************/
10 
12 #include "Math/IFunctionfwd.h"
13 #include <cmath>
14 #include <limits>
15 
16 #ifndef ROOT_Math_Error
17 #include "Math/Error.h"
18 #endif
19 
20 namespace ROOT {
21 namespace Math {
22 
24  fFunctionCopied(false),
25  fStepSize(h),
26  fLastError(0),
27  fFunction(0)
28 {
29  // Default Constructor.
30 }
32  fFunctionCopied(copyFunc),
33  fStepSize(h),
34  fLastError(0),
35  fFunction(0)
36 {
37  // Constructor from a function and step size
38  if (copyFunc) fFunction = f.Clone();
39  else fFunction = &f;
40 }
41 
42 
44 {
45  // destructor
46  if ( fFunction != 0 && fFunctionCopied )
47  delete fFunction;
48 }
49 
51 {
52  // copy constructor
53  // copy constructor (deep copy or not depending on fFunctionCopied)
54  fStepSize = rhs.fStepSize;
55  fLastError = rhs.fLastError;
57  SetFunction(*rhs.fFunction);
58  }
59 
61 {
62  // Assignment operator
63  if (&rhs == this) return *this;
65  fStepSize = rhs.fStepSize;
66  fLastError = rhs.fLastError;
67  SetFunction(*rhs.fFunction);
68  return *this;
69 }
70 
72 {
73  // set function
74  if (fFunctionCopied) {
75  if (fFunction) delete fFunction;
76  fFunction = f.Clone();
77  }
78  else fFunction = &f;
79 }
80 
81 double RichardsonDerivator::Derivative1 (const IGenFunction & function, double x, double h)
82 {
83  const double keps = std::numeric_limits<double>::epsilon();
84 
85 
86  double xx;
87  xx = x+h; double f1 = (function)(xx);
88 
89  xx = x-h; double f2 = (function)(xx);
90 
91  xx = x+h/2; double g1 = (function)(xx);
92  xx = x-h/2; double g2 = (function)(xx);
93 
94  //compute the central differences
95  double h2 = 1/(2.*h);
96  double d0 = f1 - f2;
97  double d2 = g1 - g2;
98  double deriv = h2*(8*d2 - d0)/3.;
99  // // compute the error ( to be improved ) this is just a simple truncation error
100  // fLastError = kC1*h2*0.5*(f1+f2); //compute the error
101 
102  // compute the error ( from GSL deriv implementation)
103 
104  double e0 = (std::abs( f1) + std::abs(f2)) * keps;
105  double e2 = 2* (std::abs( g1) + std::abs(g2)) * keps + e0;
106  double delta = std::max( std::abs( h2*d0), std::abs( deriv) ) * std::abs( x)/h * keps;
107 
108  // estimate the truncation error from d2-d0 which is O(h^2)
109  double err_trunc = std::abs( deriv - h2*d0 );
110  // rounding error due to cancellation
111  double err_round = std::abs( e2/h) + delta;
112 
113  fLastError = err_trunc + err_round;
114  return deriv;
115 }
116 
117 double RichardsonDerivator::DerivativeForward (const IGenFunction & function, double x, double h)
118 {
119  const double keps = std::numeric_limits<double>::epsilon();
120 
121 
122  double xx;
123  xx = x+h/4.0; double f1 = (function)(xx);
124  xx = x+h/2.0; double f2 = (function)(xx);
125 
126  xx = x+(3.0/4.0)*h; double f3 = (function)(xx);
127  xx = x+h; double f4 = (function)(xx);
128 
129  //compute the forward differences
130 
131  double r2 = 2.0*(f4 - f2);
132  double r4 = (22.0 / 3.0) * (f4 - f3) - (62.0 / 3.0) * (f3 - f2) +
133  (52.0 / 3.0) * (f2 - f1);
134 
135  // Estimate the rounding error for r4
136 
137  double e4 = 2 * 20.67 * (fabs (f4) + fabs (f3) + fabs (f2) + fabs (f1)) * keps;
138 
139  // The next term is due to finite precision in x+h = O (eps * x)
140 
141  double dy = std::max (fabs (r2 / h), fabs (r4 / h)) * fabs (x / h) * keps;
142 
143  // The truncation error in the r4 approximation itself is O(h^3).
144  // However, for safety, we estimate the error from r4-r2, which is
145  // O(h). By scaling h we will minimise this estimated error, not
146  // the actual truncation error in r4.
147 
148  double result = r4 / h;
149  double abserr_trunc = fabs ((r4 - r2) / h); // Estimated truncation error O(h)
150  double abserr_round = fabs (e4 / h) + dy;
151 
152  fLastError = abserr_trunc + abserr_round;
153 
154  return result;
155 }
156 
157 
158  double RichardsonDerivator::Derivative2 (const IGenFunction & function, double x, double h)
159 {
160  const double kC1 = 4*std::numeric_limits<double>::epsilon();
161 
162  double xx;
163  xx = x+h; double f1 = (function)(xx);
164  xx = x; double f2 = (function)(xx);
165  xx = x-h; double f3 = (function)(xx);
166 
167  xx = x+h/2; double g1 = (function)(xx);
168  xx = x-h/2; double g3 = (function)(xx);
169 
170  //compute the central differences
171  double hh = 1/(h*h);
172  double d0 = f3 - 2*f2 + f1;
173  double d2 = 4*g3 - 8*f2 +4*g1;
174  fLastError = kC1*hh*f2; //compute the error
175  double deriv = hh*(4*d2 - d0)/3.;
176  return deriv;
177 }
178 
179 double RichardsonDerivator::Derivative3 (const IGenFunction & function, double x, double h)
180 {
181  const double kC1 = 4*std::numeric_limits<double>::epsilon();
182 
183  double xx;
184  xx = x+2*h; double f1 = (function)(xx);
185  xx = x+h; double f2 = (function)(xx);
186  xx = x-h; double f3 = (function)(xx);
187  xx = x-2*h; double f4 = (function)(xx);
188  xx = x; double fx = (function)(xx);
189  xx = x+h/2; double g2 = (function)(xx);
190  xx = x-h/2; double g3 = (function)(xx);
191 
192  //compute the central differences
193  double hhh = 1/(h*h*h);
194  double d0 = 0.5*f1 - f2 +f3 - 0.5*f4;
195  double d2 = 4*f2 - 8*g2 +8*g3 - 4*f3;
196  fLastError = kC1*hhh*fx; //compute the error
197  double deriv = hhh*(4*d2 - d0)/3.;
198  return deriv;
199 }
200 
201 } // end namespace Math
202 
203 } // end namespace ROOT
double Derivative3(double x)
Returns the third derivative of the function at point x, computed by Richardson&#39;s extrapolation metho...
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:133
This namespace contains pre-defined functions to be used in conjuction with TExecutor::Map and TExecu...
Definition: StringConv.hxx:21
double DerivativeForward(double x)
Computation of the first derivative using a forward formula.
TH1 * h
Definition: legend2.C:5
double Derivative2(double x)
Returns the second derivative of the function at point x, computed by Richardson&#39;s extrapolation meth...
Double_t x[n]
Definition: legend1.C:17
void SetFunction(const IGenFunction &f)
Set function for derivative calculation (copy the function if option has been enabled in the construc...
double Derivative1(double x)
Returns the first derivative of the function at point x, computed by Richardson&#39;s extrapolation metho...
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
~RichardsonDerivator()
Destructor: Removes function if needed.
REAL epsilon
Definition: triangle.c:617
RichardsonDerivator(double h=0.001)
Default Constructor.
double f(double x)
Namespace for new Math classes and functions.
RichardsonDerivator & operator=(const RichardsonDerivator &rhs)
Assignment operator.
double f2(const double *x)
TF1 * f1
Definition: legend1.C:11
virtual IBaseFunctionOneDim * Clone() const =0
Clone a function.
double result[121]
User class for calculating the derivatives of a function.
unsigned int r2[N_CITIES]
Definition: simanTSP.cxx:322