ROMBERG INTEGRATION

#include "iostream.h"
#include "math.h"
#include "Romberg.h"

int main()
{
        // tableau contenant des paramêtres qui n'entre pas dans l'intégration
	double * Tab;
	Tab = new double[2];
	*Tab=2.0;
	*(Tab+1)=1.0;
	
        //constructeur 
	Romberg R(1,2,Tab);	
	
	cout << endl
	<<R.RI();
	cout << endl;
	delete [] Tab;
	
return 0;
}
//voila pour le .h qui contient la définition de la classe
#ifndef Romberg_H
#define Romberg_H

#include  "math.h"
#define Tol 0.000000001

class Romberg
{
public: double a; //borne inf
		double b; //borne sup
		double * para; //pointeur de tableau contenant les paramétres additionnelles de la fonction

		Romberg(double A, double B, double *Para)
		{
			a=A;
			b=B;
			para=Para; //ici il n'y a qu'1 paramétre
		}

		//fonction à intégrer et à editer (ici f(x)=1/(x^k) + g
		//on intégre par rapport à x, k et g sont des paramêtres additionnelles 
		
		double ff(double i)
		{return pow(i,- *para)+ *(para+1);}
		//fin de la fonction


		double TR(double n)
		{
			double i, temp, H;
			temp=0;
			H = (b - a) / n;

			for(i=a; i<=b; i+=H)
				{
				if(i==a || i==b )
				{temp +=ff(i) / 2;}
				else
				{temp +=ff(i);}
				}

			return temp * H;
		}
		
		double Ink(double R1,double R2,int k)
		{
			return (pow(4,k) * R1 - R2) / (pow(4,k) - 1);
		}

		double RI()
		{
			const int MaxIter = 11;
			double temp;
			double Mresult[MaxIter][MaxIter];
			int n, k, i, m;

			m = 0;

			for(k=0;k<MaxIter;k++)
				{
				n = static_cast<int>(pow(2,k));
				if(n==1)	{Mresult[m][0] = TR(n);}
				else    {Mresult[m][0] = TR(n);}

				for(i=1;i<=k;i++)
				{Mresult[m][i] = Ink(Mresult[m][i-1], Mresult[m - 1][i-1], i);}

				temp= fabs(Mresult[m-1][m-1] - Mresult[m][m]);
	
				if(temp<Tol && m>=2){k = 10;}
				temp= Mresult[m][m];
				m++;
				}
		return temp;
		}

};


#endif