//
//
//
//
// Lab. Calcolo II - Esempio di codice
//
//
//
//
//
//
// <![CDATA[
// Example code: Derivative (Ridders algorithm) functor.
// (francesco.prelz@mi.infn.it 20060814)
#include <iostream>
#include <iomanip>
#include <functional>
#include <cmath>
template <typename DT>
class FunctionDerivative
{
public:
typedef DT value_type;
typedef DT (*function_ptr_type)(DT); // Normal - nonmember function pointer
typedef std::function<DT (DT arg)> functor_type;
// Creator
FunctionDerivative(functor_type &f);
FunctionDerivative(function_ptr_type f);
// Destructor
~FunctionDerivative();
DT operator() (DT argument);
DT derivative(DT argument);
DT set_h(DT h) { m_h = h; };
DT get_err() { return m_err; };
private:
functor_type m_f;
DT m_h;
DT m_err;
};
template <typename DT>
FunctionDerivative<DT>::FunctionDerivative(functor_type &f) : m_f(f), m_h(0.5) {}
template <typename DT>
FunctionDerivative<DT>::FunctionDerivative(function_ptr_type f) : m_f(f), m_h(0.5) {}
template <typename DT>
FunctionDerivative<DT>::~FunctionDerivative() {}
template <typename DT>
DT FunctionDerivative<DT>::operator() (DT argument)
{
return derivative(argument);
}
template <typename DT>
DT FunctionDerivative<DT>::derivative(DT x)
{
const DT con=1.4; // Stepsize is decreased by con at each iteration
const DT con2=(con*con);
const DT big=1.0e50;
const int ntab=10; // Sets maximum size of tableau
const DT safe=2.0; // Return when error is SAFE worse than the best so far
// returns the derivative of a function func at a point x by Ridders'
// method of polynomial extrapolation. The value h is input as an estimated
// initial stepsize; it need not be small, but rather should be an
// increment in x over which func changes substantially.
// An estimate of the error in the derivative is returned as err */
int i,j;
double errt, fac, hh, a[ntab+1][ntab+1], ans;
double test1,test2;
if (m_h == 0)
{
std::cerr << "Increment parameter for Ridders' algorithm must be nonzero" << std::endl;
return 0;
}
hh = m_h;
// Note: this routine assumes a to be a one-based matrix
a[1][1]=(m_f(x+hh)-m_f(x-hh))/(2.0*hh);
m_err = big;
for (i=2;i<=ntab;i++)
{
// Successive columns in the Neville tableau will go to smaller
// stepsizes and higher orders of extrapolation.
hh /= con;
// Try new, smaller stepsize
a[1][i]=(m_f(x+hh)-m_f(x-hh))/(2.0*hh);
fac=con2;
for(j=2;j<=i;j++)
{
// Compute extrapolations of various orders, requiring no
// new function evaluations
a[j][i] = (a[j-1][i]*fac-a[j-1][i-1])/(fac-1.0);
fac = con2*fac;
test1 = std::fabs(a[j][i]-a[j-1][i]);
test2 = std::fabs(a[j][i]-a[j-1][i-1]);
if (test1>test2) errt = test1;
else errt = test2;
// The error strategy is to compare each new extrapolation to one
// order lower, both at the present stepsize ans the previous one.
if (errt <= m_err)
{
// If error is decreased, save the improved answer
m_err=errt;
ans=a[j][i];
}
}
// If higher order is worse by a significant factor SAFE, then quit early
if (std::fabs(a[i][i]-a[i-1][i-1]) >= safe*m_err) break;
}
return (ans);
}
int
main(int argc, char *argv[])
{
FunctionDerivative<double> dfsin(std::sin);
double argument = 0;
std::cout << "sin (" << argument << ") == " << std::sin(argument)
<< " - dfsin(" << argument << ") == " << dfsin(argument)
<< " +/- " << dfsin.get_err() << std::endl;
return 0;
}
]]></plaintext></span>