Logo ROOT  
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 
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 
93  xx = x-h/2; double g2 = (function)(xx);
94 
95  //compute the central differences
96  double h2 = 1/(2.*h);
97  double d0 = f1 - f2;
98  double d2 = g1 - g2;
99  double deriv = h2*(8*d2 - d0)/3.;
100  // // compute the error ( to be improved ) this is just a simple truncation error
101  // fLastError = kC1*h2*0.5*(f1+f2); //compute the error
102 
103  // compute the error ( from GSL deriv implementation)
104 
105  double e0 = (std::abs( f1) + std::abs(f2)) * keps;
106  double e2 = 2* (std::abs( g1) + std::abs(g2)) * keps + e0;
107  double delta = std::max( std::abs( h2*d0), std::abs( deriv) ) * std::abs( x)/h * keps;
108 
109  // estimate the truncation error from d2-d0 which is O(h^2)
110  double err_trunc = std::abs( deriv - h2*d0 );
111  // rounding error due to cancellation
112  double err_round = std::abs( e2/h) + delta;
113 
114  fLastError = err_trunc + err_round;
115  return deriv;
116 }
117 
118 double RichardsonDerivator::DerivativeForward (const IGenFunction & function, double x, double h)
119 {
120  const double keps = std::numeric_limits<double>::epsilon();
121 
122 
123  double xx;
124  xx = x+h/4.0; double f1 = (function)(xx);
125  xx = x+h/2.0; double f2 = (function)(xx);
126 
127  xx = x+(3.0/4.0)*h; double f3 = (function)(xx);
128  xx = x+h; double f4 = (function)(xx);
129 
130  //compute the forward differences
131 
132  double r2 = 2.0*(f4 - f2);
133  double r4 = (22.0 / 3.0) * (f4 - f3) - (62.0 / 3.0) * (f3 - f2) +
134  (52.0 / 3.0) * (f2 - f1);
135 
136  // Estimate the rounding error for r4
137 
138  double e4 = 2 * 20.67 * (fabs (f4) + fabs (f3) + fabs (f2) + fabs (f1)) * keps;
139 
140  // The next term is due to finite precision in x+h = O (eps * x)
141 
142  double dy = std::max (fabs (r2 / h), fabs (r4 / h)) * fabs (x / h) * keps;
143 
144  // The truncation error in the r4 approximation itself is O(h^3).
145  // However, for safety, we estimate the error from r4-r2, which is
146  // O(h). By scaling h we will minimise this estimated error, not
147  // the actual truncation error in r4.
148 
149  double result = r4 / h;
150  double abserr_trunc = fabs ((r4 - r2) / h); // Estimated truncation error O(h)
151  double abserr_round = fabs (e4 / h) + dy;
152 
153  fLastError = abserr_trunc + abserr_round;
154 
155  return result;
156 }
157 
158 
159  double RichardsonDerivator::Derivative2 (const IGenFunction & function, double x, double h)
160 {
161  const double kC1 = 4*std::numeric_limits<double>::epsilon();
162 
163  double xx;
164  xx = x+h; double f1 = (function)(xx);
165  xx = x; double f2 = (function)(xx);
166  xx = x-h; double f3 = (function)(xx);
167 
168  xx = x+h/2; double g1 = (function)(xx);
169  xx = x-h/2; double g3 = (function)(xx);
170 
171  //compute the central differences
172  double hh = 1/(h*h);
173  double d0 = f3 - 2*f2 + f1;
174  double d2 = 4*g3 - 8*f2 +4*g1;
175  fLastError = kC1*hh*f2; //compute the error
176  double deriv = hh*(4*d2 - d0)/3.;
177  return deriv;
178 }
179 
180 double RichardsonDerivator::Derivative3 (const IGenFunction & function, double x, double h)
181 {
182  const double kC1 = 4*std::numeric_limits<double>::epsilon();
183 
184  double xx;
185  xx = x+2*h; double f1 = (function)(xx);
186  xx = x+h; double f2 = (function)(xx);
187  xx = x-h; double f3 = (function)(xx);
188  xx = x-2*h; double f4 = (function)(xx);
189  xx = x; double fx = (function)(xx);
190  xx = x+h/2; double g2 = (function)(xx);
191  xx = x-h/2; double g3 = (function)(xx);
192 
193  //compute the central differences
194  double hhh = 1/(h*h*h);
195  double d0 = 0.5*f1 - f2 +f3 - 0.5*f4;
196  double d2 = 4*f2 - 8*g2 +8*g3 - 4*f3;
197  fLastError = kC1*hhh*fx; //compute the error
198  double deriv = hhh*(4*d2 - d0)/3.;
199  return deriv;
200 }
201 
202 } // end namespace Math
203 
204 } // end namespace ROOT
f
#define f(i)
Definition: RSha256.hxx:122
ROOT::Math::RichardsonDerivator::Derivative2
double Derivative2(double x)
Returns the second derivative of the function at point x, computed by Richardson's extrapolation meth...
Definition: RichardsonDerivator.h:172
IFunctionfwd.h
ROOT::Math::RichardsonDerivator::DerivativeForward
double DerivativeForward(double x)
Computation of the first derivative using a forward formula.
Definition: RichardsonDerivator.h:125
RichardsonDerivator.h
x
Double_t x[n]
Definition: legend1.C:17
ROOT::Math::RichardsonDerivator
User class for calculating the derivatives of a function.
Definition: RichardsonDerivator.h:55
ROOT::Math::RichardsonDerivator::SetFunction
void SetFunction(const IGenFunction &f)
Set function for derivative calculation (copy the function if option has been enabled in the construc...
Definition: RichardsonDerivator.cxx:81
ROOT::Math::RichardsonDerivator::Derivative1
double Derivative1(double x)
Returns the first derivative of the function at point x, computed by Richardson's extrapolation metho...
Definition: RichardsonDerivator.h:116
ROOT::Math::RichardsonDerivator::fStepSize
double fStepSize
Definition: RichardsonDerivator.h:235
ROOT::Math::RichardsonDerivator::fFunctionCopied
bool fFunctionCopied
Definition: RichardsonDerivator.h:234
ROOT::Math::fabs
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:131
ROOT::Math::RichardsonDerivator::fFunction
const IGenFunction * fFunction
Definition: RichardsonDerivator.h:237
ROOT::Math::RichardsonDerivator::~RichardsonDerivator
~RichardsonDerivator()
Destructor: Removes function if needed.
Definition: RichardsonDerivator.cxx:53
h
#define h(i)
Definition: RSha256.hxx:124
Error.h
epsilon
REAL epsilon
Definition: triangle.c:617
ROOT::Math::RichardsonDerivator::RichardsonDerivator
RichardsonDerivator(double h=0.001)
Default Constructor.
Definition: RichardsonDerivator.cxx:33
ROOT::Math::IBaseFunctionOneDim
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:135
f1
TF1 * f1
Definition: legend1.C:11
ROOT::Math::RichardsonDerivator::fLastError
double fLastError
Definition: RichardsonDerivator.h:236
ROOT::Math::RichardsonDerivator::operator=
RichardsonDerivator & operator=(const RichardsonDerivator &rhs)
Assignment operator.
Definition: RichardsonDerivator.cxx:70
ROOT::Math::RichardsonDerivator::Derivative3
double Derivative3(double x)
Returns the third derivative of the function at point x, computed by Richardson's extrapolation metho...
Definition: RichardsonDerivator.h:211
ROOT
VSD Structures.
Definition: StringConv.hxx:21
Math