Trapèzes et intégrales

Voila je dois faire un programme calculant une intégrale par la méthode des trapèzes, plus précisément celle de Romberg. Et il faudrait qu'au lancement du programme, l'utilisateur puisse choisir les fonctions qui pourront être utiliées dans le programme.
exemple : l'utilisateur lance le programme, il a le choix entre : exponentielle, cosinus, ln, sinus, tangeante, ch,sh et th.
C'est assez pressé, merci d'avance.

-------------------------------
Réponse au message :
-------------------------------

> Voila je dois faire un programme calculant une intégrale par la méthode des trapèzes, plus précisément celle de Romberg. Et il faudrait qu'au lancement du programme, l'utilisateur puisse choisir les fonctions qui pourront être utiliées dans le programme.
> exemple : l'utilisateur lance le programme, il a le choix entre : exponentielle, cosinus, ln, sinus, tangeante, ch,sh et th.
> C'est assez pressé, merci d'avance.
>
>
>

Voici le code que mon collègue Michel Bernard a écrit,
Cette méthode converge très rapidememt.

En espérant cette réponse pas  trop ardive (je viens de m'inscrire).
Si vous avez entretemps une autre méthode qui convergerait encore plus ramidement,
n'hesitez pas à nous  la transmettre !

Bien à vous,

Yves

#include <iostream>

#include <string>

#include <sstream>

#include <vector>

#include <math.h>

using namespace std;

double RombergIntegrator( double func(double), ///< Function to integrate

double a, ///< Lower bound of the integral

double b, ///< Upper bound of the integral

double eps = 1.E-9 ///< Accuracy

)

{

double h = (b-a)/2.; // initial step

int n = 0; // iteration counter

int pow2 = 1; // powers of 2

double i00, iprev;

vector< double > I; // store extrapolated values from the Neville's algorithm

int k;

// compute I(0,0)

double i = h*( (*func)(a) + (*func)(b) );

double i0 = i;

while ( n<20 )

{

I.push_back(i);

// save latest estimation

iprev = i;

// save I(n,0)

i00 = i0;

// compute I(n+1,0) with recursive trapezoidal rule

i0 = 0;

for (k=1;k<=pow2;k++)

{

i0 += (*func)( a + h*(2.*k - 1) );

}

i = i0 = i00/2 + h * i0;

// Richardson extrapolation using Neville's algorithm

int pow4 = 1;

for (k=0;k<=n;k++)

{

pow4 *= 4;

double itemp = I[k];

cout << "I(" << n << "," << k << ") = " << itemp << '\t';

I[k] = i;

i = (pow4 * i - itemp) / (pow4 - 1);

}

cout << endl;

// exit if wanted accuracy has been raached

if ( fabs( i - iprev ) <= eps*fabs(i) ) return i;

n++;

pow2 *= 2;

h /= 2;

}

cout << "Too many steps in RombergIntegrator" << endl;

return 0.;

}

double f(double x)

{

return ( x*x + x + 1) * cos(x);

}

double g(double x)

{

return ( x*x*x*x*log( x + sqrt( x*x + 1 ) ) );

}

main()

{

double pi=3.1415926535;

cout.precision(15);

cout << "I(f)=" << RombergIntegrator( f, 0., pi/2. ) << endl;

cout << "I(g)=" << RombergIntegrator( g, 0., 2. ) << endl;

return 0;

}

Voici le code que mon collègue Michel Bernard a écrit,
Cette méthode converge très rapidememt.

En espérant cette réponse pas  trop ardive (je viens de m'inscrire).
Si vous avez entretemps une autre méthode qui convergerait encore plus ramidement,
n'hesitez pas à nous  la transmettre !

Bien à vous,

Yves

#include <iostream>

#include <string>

#include <sstream>

#include <vector>

#include <math.h>

using namespace std;

double RombergIntegrator( double func(double), ///< Function to integrate

                                     double a, ///< Lower bound of the integral

                                     double b, ///< Upper bound of the integral

                                    double eps = 1.E-9 ///< Accuracy

                                    )

{

double h = (b-a)/2.; // initial step

int n = 0; // iteration counter

int pow2 = 1; // powers of 2

double i00, iprev;

vector< double > I; // store extrapolated values from the Neville's algorithm

int k;

// compute I(0,0)

double i = h*( (*func)(a) + (*func)(b) );

double i0 = i;

while ( n<20 )

{

I.push_back(i);

// save latest estimation

iprev = i;

// save I(n,0)

i00 = i0;

// compute I(n+1,0) with recursive trapezoidal rule

i0 = 0;

for (k=1;k<=pow2;k++)

{

i0 += (*func)( a + h*(2.*k - 1) );

}

i = i0 = i00/2 + h * i0;

// Richardson extrapolation using Neville's algorithm

int pow4 = 1;

for (k=0;k<=n;k++)

{

pow4 *= 4;

double itemp = I[k];

cout << "I(" << n << "," << k << ") = " << itemp << '\t';

I[k] = i;

i = (pow4 * i - itemp) / (pow4 - 1);

}

cout << endl;

// exit if wanted accuracy has been raached

if ( fabs( i - iprev ) <= eps*fabs(i) ) return i;

n++;

pow2 *= 2;

h /= 2;

}

cout << "Too many steps in RombergIntegrator" << endl;

return 0.;

}

double f(double x)

{

return ( x*x + x + 1) * cos(x);

}

double g(double x)

{

return ( x*x*x*x*log( x + sqrt( x*x + 1 ) ) );

}

main()

{

double pi=3.1415926535;

cout.precision(15);

cout << "I(f)=" << RombergIntegrator( f, 0., pi/2. ) << endl;

cout << "I(g)=" << RombergIntegrator( g, 0., 2. ) << endl;

return 0;

}