Logo ROOT   6.14/05
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 #include <algorithm>
16 
17 #include "Math/Error.h"
18 
19 namespace ROOT {
20 namespace Math {
21 
23  fFunctionCopied(false),
24  fStepSize(h),
25  fLastError(0),
26  fFunction(0)
27 {
28  // Default Constructor.
29 }
31  fFunctionCopied(copyFunc),
32  fStepSize(h),
33  fLastError(0),
34  fFunction(0)
35 {
36  // Constructor from a function and step size
37  if (copyFunc) fFunction = f.Clone();
38  else fFunction = &f;
39 }
40 
41 
43 {
44  // destructor
45  if ( fFunction != 0 && fFunctionCopied )
46  delete fFunction;
47 }
48 
50 {
51  // copy constructor
52  // copy constructor (deep copy or not depending on fFunctionCopied)
53  fStepSize = rhs.fStepSize;
54  fLastError = rhs.fLastError;
56  SetFunction(*rhs.fFunction);
57  }
58 
60 {
61  // Assignment operator
62  if (&rhs == this) return *this;
64  fStepSize = rhs.fStepSize;
65  fLastError = rhs.fLastError;
66  SetFunction(*rhs.fFunction);
67  return *this;
68 }
69 
71 {
72  // set function
73  if (fFunctionCopied) {
74  if (fFunction) delete fFunction;
75  fFunction = f.Clone();
76  }
77  else fFunction = &f;
78 }
79 
80 double RichardsonDerivator::Derivative1 (const IGenFunction & function, double x, double h)
81 {
82  const double keps = std::numeric_limits<double>::epsilon();
83 
84 
85  double xx;
86  xx = x+h; double f1 = (function)(xx);
87 
88  xx = x-h; double f2 = (function)(xx);
89 
90  xx = x+h/2; double g1 = (function)(xx);
91  xx = x-h/2; double g2 = (function)(xx);
92 
93  //compute the central differences
94  double h2 = 1/(2.*h);
95  double d0 = f1 - f2;
96  double d2 = g1 - g2;
97  double deriv = h2*(8*d2 - d0)/3.;
98  // // compute the error ( to be improved ) this is just a simple truncation error
99  // fLastError = kC1*h2*0.5*(f1+f2); //compute the error
100 
101  // compute the error ( from GSL deriv implementation)
102 
103  double e0 = (std::abs( f1) + std::abs(f2)) * keps;
104  double e2 = 2* (std::abs( g1) + std::abs(g2)) * keps + e0;
105  double delta = std::max( std::abs( h2*d0), std::abs( deriv) ) * std::abs( x)/h * keps;
106 
107  // estimate the truncation error from d2-d0 which is O(h^2)
108  double err_trunc = std::abs( deriv - h2*d0 );
109  // rounding error due to cancellation
110  double err_round = std::abs( e2/h) + delta;
111 
112  fLastError = err_trunc + err_round;
113  return deriv;
114 }
115 
116 double RichardsonDerivator::DerivativeForward (const IGenFunction & function, double x, double h)
117 {
118  const double keps = std::numeric_limits<double>::epsilon();
119 
120 
121  double xx;
122  xx = x+h/4.0; double f1 = (function)(xx);
123  xx = x+h/2.0; double f2 = (function)(xx);
124 
125  xx = x+(3.0/4.0)*h; double f3 = (function)(xx);
126  xx = x+h; double f4 = (function)(xx);
127 
128  //compute the forward differences
129 
130  double r2 = 2.0*(f4 - f2);
131  double r4 = (22.0 / 3.0) * (f4 - f3) - (62.0 / 3.0) * (f3 - f2) +
132  (52.0 / 3.0) * (f2 - f1);
133 
134  // Estimate the rounding error for r4
135 
136  double e4 = 2 * 20.67 * (fabs (f4) + fabs (f3) + fabs (f2) + fabs (f1)) * keps;
137 
138  // The next term is due to finite precision in x+h = O (eps * x)
139 
140  double dy = std::max (fabs (r2 / h), fabs (r4 / h)) * fabs (x / h) * keps;
141 
142  // The truncation error in the r4 approximation itself is O(h^3).
143  // However, for safety, we estimate the error from r4-r2, which is
144  // O(h). By scaling h we will minimise this estimated error, not
145  // the actual truncation error in r4.
146 
147  double result = r4 / h;
148  double abserr_trunc = fabs ((r4 - r2) / h); // Estimated truncation error O(h)
149  double abserr_round = fabs (e4 / h) + dy;
150 
151  fLastError = abserr_trunc + abserr_round;
152 
153  return result;
154 }
155 
156 
157  double RichardsonDerivator::Derivative2 (const IGenFunction & function, double x, double h)
158 {
159  const double kC1 = 4*std::numeric_limits<double>::epsilon();
160 
161  double xx;
162  xx = x+h; double f1 = (function)(xx);
163  xx = x; double f2 = (function)(xx);
164  xx = x-h; double f3 = (function)(xx);
165 
166  xx = x+h/2; double g1 = (function)(xx);
167  xx = x-h/2; double g3 = (function)(xx);
168 
169  //compute the central differences
170  double hh = 1/(h*h);
171  double d0 = f3 - 2*f2 + f1;
172  double d2 = 4*g3 - 8*f2 +4*g1;
173  fLastError = kC1*hh*f2; //compute the error
174  double deriv = hh*(4*d2 - d0)/3.;
175  return deriv;
176 }
177 
178 double RichardsonDerivator::Derivative3 (const IGenFunction & function, double x, double h)
179 {
180  const double kC1 = 4*std::numeric_limits<double>::epsilon();
181 
182  double xx;
183  xx = x+2*h; double f1 = (function)(xx);
184  xx = x+h; double f2 = (function)(xx);
185  xx = x-h; double f3 = (function)(xx);
186  xx = x-2*h; double f4 = (function)(xx);
187  xx = x; double fx = (function)(xx);
188  xx = x+h/2; double g2 = (function)(xx);
189  xx = x-h/2; double g3 = (function)(xx);
190 
191  //compute the central differences
192  double hhh = 1/(h*h*h);
193  double d0 = 0.5*f1 - f2 +f3 - 0.5*f4;
194  double d2 = 4*f2 - 8*g2 +8*g3 - 4*f3;
195  fLastError = kC1*hhh*fx; //compute the error
196  double deriv = hhh*(4*d2 - d0)/3.;
197  return deriv;
198 }
199 
200 } // end namespace Math
201 
202 } // 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:135
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
double DerivativeForward(double x)
Computation of the first derivative using a forward formula.
#define f(i)
Definition: RSha256.hxx:104
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
#define h(i)
Definition: RSha256.hxx:106
RichardsonDerivator(double h=0.001)
Default Constructor.
Namespace for new Math classes and functions.
RichardsonDerivator & operator=(const RichardsonDerivator &rhs)
Assignment operator.
TF1 * f1
Definition: legend1.C:11
virtual IBaseFunctionOneDim * Clone() const =0
Clone a function.
User class for calculating the derivatives of a function.