file "integration.h"
#ifndef INTEGRATION_H_
#define INTEGRATION_H_
class Integration
{
private:
double (*nfunc)(double,double);
double xsav;
double f1(double);
double f2(double);
double x1,y1,x2,y2;
int accuracy;
public:
Integration();
Integration(int naccuracy,double xx1,double xx2,double yy1,double yy2);
virtual ~Integration() { }
double qgaus(double (*f)(double),double min,double max) const;
double quad2d(double (*func)(double,double));
void SetIntLim(const double xx1,const double xx2,const double yy1,const
double yy2);
void SetAccuracy(const int naccuracy);
};
#endif
file "integration.cpp"
#include "stdafx.h"
#include "Integration.h"
Integration::Integration()
{
x1 = 0;
x2 = 1;
y1 = 0;
y2 = 1;
accuracy = 20;
}
Integration::Integration(int naccuracy,double xx1,double xx2,double
yy1,double yy2)
{
x1 = xx1 ;
x2 = xx2 ;
y1 = yy1;
y2 = yy2;
accuracy = naccuracy;
}
void Integration::SetAccuracy(const int naccuracy)
{
accuracy = naccuracy;
}
void Integration::SetIntLim(const double xx1,const double xx2,const double
yy1,const double yy2)
{
x1 = xx1 ;
x2 = xx2 ;
y1 = yy1;
y2 = yy2;
}
double Integration::qgaus(double (*f)(double),double min,double max) const
{
double sqrt3 = 1.7320508075688772935;
double step = (max - min)/accuracy;
double a,b = min;
double integral = 0;
if ((max - min) <= 0)
{
return 0 ;
}
for (int i=0;i<accuracy;i++)
{
a = b;
b +=step;
integral += ( (b-a)/2) * ( (*f)(-(b-a)/(2*sqrt3)+ (a+b)/2 ) + (*f)(
(b-a)/(sqrt3*2) + (a+b)/2 ));
}
return integral;
}
double Integration::quad2d(double (*func)(double,double))
{
nfunc = func ;
return qgaus(f1,x1,x2);
}
double Integration::f1(double x)
{
xsav = x ;
return qgaus(f2,y1,y2);
}
double Integration::f2(double y)
{
return nfunc(xsav,y);
}