//
//
//
//
// Lab. Calcolo II - Esempio di codice
//
//
//
//
//
//
// <![CDATA[
// Example code: Integral (Romberg algorithm) functor.
// (francesco.prelz@mi.infn.it 20060911)
#include <iostream>
#include <iomanip>
#include <functional>
#include <cmath>
template <typename DT>
class FunctionIntegral
{
public:
typedef DT value_type;
typedef DT (*function_ptr_type)(DT); // Normal - nonmember function pointer
typedef std::function<DT (DT arg)> functor_type;
// Creator
FunctionIntegral(functor_type &f);
FunctionIntegral(function_ptr_type f);
// Destructor
~FunctionIntegral();
DT operator() (DT argument);
DT integral(DT a, DT b);
DT get_err() { return m_err; };
private:
void polint(DT xa[], DT ya[], int n, DT x, DT &y, DT &dy);
DT trapzd(DT a, DT b, int n);
functor_type m_f;
DT m_err;
};
template <typename DT>
FunctionIntegral<DT>::FunctionIntegral(functor_type &f) : m_f(f), m_err(0) {}
template <typename DT>
FunctionIntegral<DT>::FunctionIntegral(function_ptr_type f) : m_f(f), m_err(0) {}
template <typename DT>
FunctionIntegral<DT>::~FunctionIntegral() {}
template <typename DT>
DT FunctionIntegral<DT>::operator() (DT argument)
{
// Assume integral function starting from 0.
return integral(0,argument);
}
template <typename DT>
DT FunctionIntegral<DT>::integral(DT a, DT b)
{
const DT EPS=1e-6;
const int JMAX=20;
const int JMAXP=(JMAX+1);
const int K=5;
// Here EPS is the fractional accuracy desired, as determined by the
// extrapolation error estimate; JMAX limits the total number of steps;
// K is the number of points used in the extrapolation.
// This method returns the integral of the function func from a to b.
// Integration is performed by Romberg's method of order 2K, where, e.g.,
// K=2 is Simpson's rule.
DT ss;
DT s[JMAXP],h[JMAXP+1]; // These store the successive trapezoidal
// approximations and their relative stepsizes.
int j;
h[1]=1.0;
for (j=1;j<=JMAX;j++)
{
s[j]=trapzd(a,b,j);
if (j >= K)
{
polint(&h[j-K],&s[j-K],K,0.0,ss,m_err);
if (std::fabs(m_err) <= EPS*std::fabs(ss)) return ss;
}
h[j+1] = 0.25*h[j];
// This is a key step: The factor is 0.25 even though the stepsize is
// decreased by only 0.5. This makes the extrapolation a polynomial in
// h^2, not just a polynomial in h.
}
std::cerr << "Too many steps in method FunctionIntegral<>::integral."
<< std::endl;
return 0.0; // Never get here.
}
template <typename DT>
DT FunctionIntegral<DT>::trapzd(DT a, DT b, int n)
{
// This routine computes the nth stage of refinement of an extended
// trapezoidal rule. func is input as a pointer to the function to be
// integrated between limits a and b, also input. When called
// with n=1, the routine returns the crudest estimate of the integral.
// Subsequent calls with n=2,3,... (in that sequential order) will
// improve the accuracy by adding 2^(n-2) additional interior points */
DT x, tnm, sum, del;
static DT s;
int it,j;
if (n == 1)
{
return(s=0.5*(b-a)*(m_f(a)+m_f(b)));
}
else
{
for(it=1,j=1;j<n-1;j++) it <<= 1;
tnm = it;
del=(b-a)/tnm; // This is the spacing of the points to be added.
x=a+0.5*del;
for(sum=0.0,j=1;j<=it;j++,x+=del) sum+=m_f(x);
s=0.5*(s+(b-a)*sum/tnm); // This replaces s by its refined value
return(s);
}
}
template <typename DT>
void FunctionIntegral<DT>::polint(DT xa[], DT ya[], int n, DT x, DT &y, DT &dy)
{
// Given arrays xa[1..n] and ya[1..n], and given a value x, this routine
// returns a value y, and an error estimate dy. If P(x) is the polynomial
// of degree N - 1 such that P(xai) = yai, i=1, .... n, then the returned
// value y = P(x). */
int i,m,ns=1;
double den,dif,dift,ho,hp,w;
double *c,*d;
dif=std::fabs(x-xa[1]);
c = new double[n+1];
if (!c)
{
std::cerr << "Out of memory." << std::endl;
return;
}
d = new double[n+1];
if (!d)
{
std::cerr << "Out of memory." << std::endl;
delete[] c;
return;
}
for (i=1;i<=n;i++) // Here we find the index ns of the closest tabie entry,
{
if ( (dift=std::fabs(x-xa[i])) < dif)
{
ns=i;
dif=dift;
}
c[i] = ya[i]; // and initialize the tableau of c's and d's.
d[i] = ya[i];
}
y=ya[ns--]; // This is the initial approximation to y.
for (m=1;m<n;m++) // For each column of the tableau,
{
for (i=1;i<=n-m;i++) // we loop over the current c's and d's
{ // and update them.
ho=xa[i]-x;
hp=xa[i+m]-x;
w=c[i+1]-d[i];
if ( (den=ho-hp) == 0.0)
{
// This error can occur only if two input xa's are
// (to within roundoff) identical.
std::cerr << "Error in method FunctionIntegral<>::polint"
<< std::endl;
delete[] c;
delete[] d;
return;
}
den=w/den;
d[i]=hp*den; // Here the c's and d's are updated.
c[i]=ho*den;
}
y += (dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));
// After each column in the tableau is completed, we decide which
// correction, c or d, we want to add to our accumulating value of
// y, i.e., which path to take through the tableau-forking up or down.
// We do this in such a way as to take the most "straight line" route
// through the tableau to its apex, updating ns accordingly to keep track
// of where we are. This route keeps the partial approximations centered
// (insofar as possible) on the target x. The last dy added is thus
// the error indication.
}
delete[] d;
delete[] c;
}
int
main(int argc, char *argv[])
{
FunctionIntegral<double> intsin(std::sin);
double argument = M_PI;
std::cout << "sin (" << argument << ") == " << std::sin(argument)
<< " - intsin(" << argument << ") == " << intsin(argument)
<< " +/- " << intsin.get_err() << std::endl;
return 0;
}
]]></plaintext></span>