|
|
Boost : |
From: Dan Mcleran (dan.mcleran_at_[hidden])
Date: 2002-11-04 10:50:15
Martin,
I've modified my code by incorporating some of yours and making some
changes. There are obviously plenty of details yet to be done. What I've got
now allows the user to:
* specify argument and return type.
* specify which Functor to integrate over
* specify which method of integration to use
All of this is done when the _integrate function is called. Here'e the code,
let me know what you think. I compiled and ran the code with gcc 2.95.2.
Borland Builder 5 could not hack it. I would like to make some changes in
the future to allow Borland to compile this code, maybe with some checks for
which compiler the programmer is using.
I've included an implementation of Richardsons extrapolation and the Romberg
rule. I've used this web page as a reference:
http://www2.egr.uh.edu/~rbarton/Classes/ECE2331/Files.pdf/Notes/Notes_2-11-0
2.pdf
The algortihm builds a Romberg table until the relative tolerance is met.
I've placed all of my templates in a file called NumericalIntegration.hpp
//NumericalIntegration.hpp
#ifndef _NUMERICAL_INTEGRATION_HPP_INCLUDED_
#define _NUMERICAL_INTEGRATION_HPP_INCLUDED_
#include <cmath>
#include <algorithm>
#include <vector>
#include <boost/utility.hpp>
#include <iostream>
template<typename argument_type,typename return_type = argument_type>struct
OneOver
{
public:
typedef argument_type ArgumentType;
typedef return_type ReturnType;
public:
ReturnType operator()(const ArgumentType& arg)
{
return ((ReturnType)1/(ReturnType)arg);
}
};
template<typename argument_type,typename return_type = argument_type>struct
XSquared
{
public:
typedef argument_type ArgumentType;
typedef return_type ReturnType;
public:
ReturnType operator()(const ArgumentType& arg)
{
return (arg * arg);
}
};
template<typename argument_type,typename return_type = argument_type>struct
Poly
{
public:
typedef argument_type ArgumentType;
typedef return_type ReturnType;
public:
ReturnType operator()(const ArgumentType& arg)
{
return (0.2 + (25 * arg) - (200 * pow(arg,2.0)) + (675 *
pow(arg,3.0)) - (900 * pow(arg,4.0)) + (400 * pow(arg,5.0)));
}
};
template<class Functor> class Integrator
{
public:
Integrator(){};
Integrator(typename Functor::ArgumentType a,typename
Functor::ArgumentType b,typename Functor::ArgumentType tol) :
_a(a),_b(b),_tol(tol){};
~Integrator(){};
public:
typedef typename Functor::ArgumentType ArgumentType;
typedef typename Functor::ReturnType ReturnType;
protected:
Functor _f;
typename Functor::ArgumentType _a;
typename Functor::ArgumentType _b;
typename Functor::ArgumentType _tol;
};
template<class Functor> class TrapezoidalIntegrator : public
Integrator<Functor>
{
public:
TrapezoidalIntegrator(){};
explicit TrapezoidalIntegrator(typename
Integrator<Functor>::ArgumentType a,typename
Integrator<Functor>::ArgumentType b,typename
Integrator<Functor>::ArgumentType tol):Integrator<Functor>(a,b,tol){};
~TrapezoidalIntegrator(){};
public:
typedef typename Integrator<Functor>::ArgumentType ArgumentType;
typedef typename Integrator<Functor>::ReturnType ReturnType;
public:
ReturnType operator ()()
{
ReturnType area(0);
int const n =
static_cast<int>(std::ceil(std::sqrt(std::fabs(_b-_a)/_tol)));
ArgumentType xl = _a;
ReturnType fl = _f(xl);
for(int i = 0;i < n;++i)
{
ArgumentType xr = _a + (_b-_a)*(i+1)/n;
ReturnType fr = _f(xr);
area += (_b-_a)/n * (fl+fr)/2;
xl=xr;
fl=fr;
}
return area;
}
};
template<class Functor> class MidpointIntegrator : public
Integrator<Functor>
{
public:
MidpointIntegrator(){};
explicit MidpointIntegrator(typename
Integrator<Functor>::ArgumentType a,typename
Integrator<Functor>::ArgumentType b,typename
Integrator<Functor>::ArgumentType tol):Integrator<Functor>(a,b,tol){};
~MidpointIntegrator(){};
public:
typedef typename Integrator<Functor>::ArgumentType ArgumentType;
typedef typename Integrator<Functor>::ReturnType ReturnType;
public:
ReturnType operator()()
{
ReturnType area(0);
int const n =
static_cast<int>(std::ceil(std::sqrt(std::fabs(_b-_a)/_tol)));
ArgumentType xl = _a;
for (int i=0; i<n; ++i) {
ArgumentType xr = _a + (_b-_a)*(i+1)/n;
ReturnType f = _f((xl+xr)/2);
area += (_b-_a)/n * f;
xl = xr;
}
return area;
}
};
template<class Functor> class RombergIntegrator : public Integrator<Functor>
{
public:
RombergIntegrator(){};
explicit RombergIntegrator(typename
Integrator<Functor>::ArgumentType a,typename
Integrator<Functor>::ArgumentType b,typename
Integrator<Functor>::ArgumentType tol):Integrator<Functor>(a,b,tol){};
~RombergIntegrator(){};
public:
typedef typename Integrator<Functor>::ArgumentType ArgumentType;
typedef typename Integrator<Functor>::ReturnType ReturnType;
typedef std::vector<ReturnType> ResultVector;
typedef ResultVector::iterator itResultVector;
typedef std::vector<ResultVector> RombergTable;
typedef RombergTable::iterator itRombergTable;
public:
ReturnType operator()()
{
RombergTable vRombergTable;
itRombergTable itk;
ResultVector vResults;
itResultVector itj;
ReturnType area(0);
int step_size(0);
int j(0);
int k(0);
int const MAX_N =
static_cast<int>(std::ceil(std::sqrt(std::fabs((_b - _a)/_tol))));
double power(0);
double factor1(0);
double factor2(0);
ArgumentType tolerance(0);
bool bToleranceMet(false);
vResults.clear();
//calculate area via trapezoidal rule for n = 1,2,4,etc.
for(int n = 1;n <= MAX_N;n *= 2)
{
tolerance = std::fabs(_b-_a)/(n*n);
TrapezoidalIntegrator<Functor> func(_a,_b,tolerance);
area = func();
vResults.push_back(area);
}
//build Romberg table
vRombergTable.push_back(vResults);
for(itk = vRombergTable.begin();itk !=
vRombergTable.end();++itk)
{
vResults.clear();
k = vRombergTable.size() + 1;
power = pow(4,(k-1));
factor1 = (power/(power - 1));
factor2 = (1/(power - 1));
for(itj = (*itk).begin();itj != (*itk).end() - 1;++itj)
{
j = static_cast<int>(std::distance((*itk).begin(),itj) +
1);
area = (factor1 * (*(boost::next(itj))) - (factor2 *
(*itj)));
vResults.push_back(area);
//check for meeting error tolerance
if(vResults.size() > 1)
{
if(std::fabs(vResults[vResults.size() - 1] -
vResults[vResults.size() - 2]) <= tolerance)
{
bToleranceMet = true;
break;//we have met the error tol
}
}
}
if(bToleranceMet == true)
{
break;
}
else
{
vRombergTable.push_back(vResults);
if(vResults.size() == 1)
{
break;
}
}
}
return area;
}
};
template<typename argument_type,typename
return_type,template<typename,typename>class Functor,template<class>class
IntegrationType> IntegrationType<Functor<argument_type,return_type> >
_integrate(typename Functor<argument_type,return_type>::ArgumentType
a,typename Functor<argument_type,return_type>::ArgumentType b,typename
Functor<argument_type,return_type>::ArgumentType tol)
{
IntegrationType<Functor<argument_type,return_type> > mine(a,b,tol);
return mine;
}
#endif //_NUMERICAL_INTEGRATION_HPP_INCLUDED_
Here's the main.cpp:
//--------------------------------------------------------------------------
-
#include "NumericalIntegration.hpp"
#include <iostream>
#include <iomanip>
//--------------------------------------------------------------------------
-
int main(int argc, char* argv[])
{
using namespace std;
unsigned char key;
double area(0.0);
double a(1.0);
double b(2.0);
double c(0.0);
double d(0.8);
cout<<setprecision(12)<<endl;
cout<<"****** Trapezoidal Rule *******"<<endl
<<"Area (1e-3) =
"<<_integrate<double,double,OneOver,TrapezoidalIntegrator>(a,b,1e-3)()<<endl
<<"Area (1e-4) =
"<<_integrate<double,double,OneOver,TrapezoidalIntegrator>(a,b,1e-4)()<<endl
<<"Area (1e-5) =
"<<_integrate<double,double,OneOver,TrapezoidalIntegrator>(a,b,1e-5)()<<endl
<<"Area (1e-6) =
"<<_integrate<double,double,OneOver,TrapezoidalIntegrator>(a,b,1e-6)()<<endl
<<"Area (1e-7) =
"<<_integrate<double,double,OneOver,TrapezoidalIntegrator>(a,b,1e-7)()<<endl
;
cout<<"****** Midpoint Rule *******"<<endl
<<"Area (1e-3) =
"<<_integrate<double,double,OneOver,MidpointIntegrator>(a,b,1e-3)()<<endl
<<"Area (1e-4) =
"<<_integrate<double,double,OneOver,MidpointIntegrator>(a,b,1e-4)()<<endl
<<"Area (1e-5) =
"<<_integrate<double,double,OneOver,MidpointIntegrator>(a,b,1e-5)()<<endl
<<"Area (1e-6) =
"<<_integrate<double,double,OneOver,MidpointIntegrator>(a,b,1e-6)()<<endl
<<"Area (1e-7) =
"<<_integrate<double,double,OneOver,MidpointIntegrator>(a,b,1e-7)()<<endl;
cout<<"****** Romberg Rule *******"<<endl
<<"Area (1e-3) =
"<<_integrate<double,double,Poly,RombergIntegrator>(c,d,1e-3)()<<endl
<<"Area (1e-4) =
"<<_integrate<double,double,Poly,RombergIntegrator>(c,d,1e-4)()<<endl
<<"Area (1e-5) =
"<<_integrate<double,double,Poly,RombergIntegrator>(c,d,1e-5)()<<endl
<<"Area (1e-6) =
"<<_integrate<double,double,Poly,RombergIntegrator>(c,d,1e-6)()<<endl
<<"Area (1e-7) =
"<<_integrate<double,double,Poly,RombergIntegrator>(c,d,1e-7)()<<endl;
cout<<"Press any key";
cin>>key;
return 0;
}
//--------------------------------------------------------------------------
-
Boost list run by bdawes at acm.org, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk