Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
20namespace ROOT {
21namespace 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;
58 }
59
61{
62 // Assignment operator
63 if (&rhs == this) return *this;
65 fStepSize = rhs.fStepSize;
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
81double 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
118double 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
180double 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
#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
virtual IBaseFunctionOneDim * Clone() const =0
Clone a function.
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)
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...