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
19namespace ROOT {
20namespace 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;
57 }
58
60{
61 // Assignment operator
62 if (&rhs == this) return *this;
64 fStepSize = rhs.fStepSize;
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
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
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
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
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
#define f(i)
Definition: RSha256.hxx:104
#define h(i)
Definition: RSha256.hxx:106
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:135
User class for calculating the derivatives of a function.
double Derivative2(double x)
Returns the second derivative of the function at point x, computed by Richardson's extrapolation meth...
~RichardsonDerivator()
Destructor: Removes function if needed.
RichardsonDerivator(double h=0.001)
Default Constructor.
double DerivativeForward(double x)
Computation of the first derivative using a forward formula.
RichardsonDerivator & operator=(const RichardsonDerivator &rhs)
Assignment operator.
double Derivative3(double x)
Returns the third derivative of the function at point x, computed by Richardson's extrapolation metho...
double Derivative1(double x)
Returns the first derivative of the function at point x, computed by Richardson's extrapolation metho...
void SetFunction(const IGenFunction &f)
Set function for derivative calculation (copy the function if option has been enabled in the construc...
Double_t x[n]
Definition: legend1.C:17
TF1 * f1
Definition: legend1.C:11
Namespace for new Math classes and functions.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:151
VSD Structures.
Definition: StringConv.hxx:21
REAL epsilon
Definition: triangle.c:617