#include "integ.h" #include "generalfit.h" // A faire : // 1 dans GLInteg, optimiser recalcul en utilisant des poids // calcules sur [0,1] et chgt de variable // 2 dans GLInteg, distinction nStep = subdivisions de l'intervalle, // et on applique GL d'ordre "ordre" sur chaque subdivision, en // utilisant le (1) //++ // Class Integrator // Lib Outils++ // include integ.h // // Classe abstraite d'intégration numérique 1D. // On fournit une fonction double f(double) au constructeur, ou // une GeneralFunction avec des paramètres définis. // L'objet Integrator est convertible en valeur double qui est la valeur // de l'intégrale. Diverses méthodes permettent de choisir des options // de calcul, et ces méthodes retournent une référence sur l'objet, pour // permettre une notation chaînée. //-- //++ // Links Implementations // TrpzInteg // GLInteg //-- //++ // Titre Constructeurs //-- //++ Integrator::Integrator() // // Constructeur par défaut. L'objet n'est pas utilisable en l'état. //-- : mFunc(NULL), mGFF(NULL), mGFFParm(NULL), mNStep(50), mDX(-1), mReqPrec(-1), mXMin(0), mXMax(1) {} //++ Integrator::Integrator(FUNC f, double xmin, double xmax) // // Constructeur à partir de la fonction double->double, et des // bornes d'intégration. //-- : mFunc(f), mGFF(NULL), mGFFParm(NULL), mNStep(50), mDX(-1), mReqPrec(-1), mXMin(xmin), mXMax(xmax) {} //++ Integrator::Integrator(FUNC f) // // Constructeur sans spécifier les bornes. Elles sont positionnées // à [0,1], et on pourra les modifier plus tard. //-- : mFunc(f), mGFF(NULL), mGFFParm(NULL), mNStep(50), mDX(-1), mReqPrec(-1), mXMin(0), mXMax(1) {} //++ Integrator::Integrator(GeneralFunction* gff, double* par, double xmin, double xmax) // // Constructeur à partir d'une GeneralFunction. La fonction doit être une // fonction de une variable, et on fournit les valeurs des paramètres ainsi // que les bornes d'intégration. //-- : mFunc(NULL), mGFF(gff), mGFFParm(par), mNStep(50), mDX(-1), mReqPrec(-1), mXMin(xmin), mXMax(xmax) { DBASSERT(gff->NVar() == 1); } //++ Integrator::Integrator(GeneralFunction* gff, double* par) // // Constructeur à partir d'une GeneralFunction. La fonction doit être une // fonction de une variable, et on fournit les valeurs des paramètres. // On ne spécifie pas les bornes. Elles sont positionnées // à [0,1], et on pourra les modifier plus tard. //-- : mFunc(NULL), mGFF(gff), mGFFParm(par), mNStep(50), mDX(-1), mReqPrec(-1), mXMin(0), mXMax(1) { DBASSERT(gff->NVar() == 1); } Integrator::~Integrator() {} //++ // Titre Méthodes //-- //++ Integrator& Integrator::NStep(int n) // // Spécifie le nombre de pas pour l'intégration numérique. // La signification peut dépendre de la méthode d'intégration. //-- { mNStep = n; mDX = mReqPrec = -1; StepsChanged(); return *this; } //++ Integrator& Integrator::DX(double d) // // Spécifie le pas d'intégration. // La signification peut dépendre de la méthode d'intégration. //-- { mDX = d; mNStep = mReqPrec = -1; StepsChanged(); return *this; } //++ Integrator& Integrator::ReqPrec(double p) // // Spécifie la précision souhaitée. // La signification peut dépendre de la méthode d'intégration. // Non disponible dans toutes les méthodes d'intégration. //-- { DBASSERT( !"Pas encore implemente !"); mReqPrec = p; mDX = mNStep = -1; StepsChanged(); return *this; } //++ Integrator& Integrator::Limits(double xmin, double xmax) // // Spécifie les bornes de l'intégrale. //-- { mXMin = xmin; mXMax = xmax; LimitsChanged(); return *this; } //++ Integrator& Integrator::Func(FUNC f) // // Spécifie la fonction à intégrer, sous forme double f(double). //-- { mFunc = f; mGFFParm = NULL; mFunc = NULL; FuncChanged(); return *this; } //++ Integrator& Integrator::Func(GeneralFunction* gff, double* par) // // Spécifie la fonction à intégrer, sous forme de GeneralFunction // à une variable, et les paramètres sont fournis. //-- { mGFF = gff; mGFFParm = par; mFunc = NULL; DBASSERT(gff->NVar() == 1); FuncChanged(); return *this; } double Integrator::FVal(double x) const { DBASSERT( mFunc || (mGFF && mGFFParm) ); return mFunc ? (*mFunc)(x) : mGFF->Value(&x, mGFFParm); } double Integrator::ValueBetween(double xmin, double xmax) { Limits(xmin,xmax); return Value(); } //++ // Class TrpzInteg // Lib Outils++ // include integ.h // // Classe d'intégration par la méthode des trapèzes. // Voir Integrator pour les méthodes. Le nombre de pas // est le nombre de trapèze, le pas d'intégration est // la largeur des trapèzez. Impossible de demander une // précision. // //-- //++ // Links Parents // Integrator //-- //++ // Titre Constructeurs // Voir Integrator pour les détails. //-- //++ TrpzInteg::TrpzInteg() // //-- {} //++ TrpzInteg::TrpzInteg(FUNC f, double xmin, double xmax) // //-- : Integrator(f, xmin, xmax) {} //++ TrpzInteg::TrpzInteg(FUNC f) // //-- : Integrator(f) {} //++ TrpzInteg::TrpzInteg(GeneralFunction* gff, double* par, double xmin, double xmax) // //-- : Integrator(gff, par, xmin, xmax) {} //++ TrpzInteg::TrpzInteg(GeneralFunction* gff, double* par) // //-- : Integrator(gff, par) {} TrpzInteg::~TrpzInteg() {} double TrpzInteg::Value() { double dx = mDX; double nstep = mNStep; if (dx > 0) nstep = (mXMax - mXMin)/dx; if (nstep <= 0) nstep = 10; dx = (mXMax - mXMin) / nstep; double s = (FVal(mXMin) + FVal(mXMax))/2; double x = dx + mXMin; for (int i=1; i=mXMax) THROW(rangeCheckErr); // On introduira les classes d'exections apres reflexion et de maniere systematique (Rz+cmv) // if (x<=mXMin || x>=mXMax) throw range_error("GLInteg::AddBound bound outside interval"); mBounds.insert(x); return *this; } GLInteg& GLInteg::SetOrder(int order) { mOrder = order; InvalWeights(); return *this; } double GLInteg::Value() { if (!mXPos) ComputeWeights(); double s = 0; set::iterator i=mBounds.begin(); set::iterator j=i; j++; while (j != mBounds.end()) { double s1 = 0; double x1 = *i; double x2 = *j; for (int k=0; k EPS_gauleg); mXPos[i-1] = xm-xl*z; mXPos[mOrder-i] = xm+xl*z; mWeights[i-1] = 2.0*xl/((1.0-z*z)*pp*pp); mWeights[mOrder-i] = mWeights[i-1]; } }