/************************************************************* * Program to demonstrate integration of a real function F(x) * * by Romberg's method * * ---------------------------------------------------------- * * REFERENCE: "Mathematiques en Turbo-Pascal (Part 1) by * * Marc Ducamp and Alain Reverchon, Eyrolles, * * Paris, 1987" [BIBLI 03]. * * C++ version by J-P Moreau * * (www.jpmoreau.fr) * * ---------------------------------------------------------- * * SAMPLE RUN: * * (Integrate sin(x) from x=0 to x=1) * * * * Integral of a function F(X) by Romberg's method * * * * Input begin and end values for x variable: * * * * X0 = 0 * * X1 = 1 * * * * Input desired precision: 1e-10 * * * * * * Value of integral : 4.596977e-001 * * * * Obtained precision : 9.599083e-011 * * * * Number of iterations: 4 * * * *************************************************************/ #include #include #define MAXITER 15 double x0, // begin x value x1, // end x value prec, // desired precision integral, // result of integral obtprec; // obtained precision int niter; // number of actual iterations // Given function to integrate double FUNC(double x) { return sin(x); } /****************************************************** * Integral of a function FUNC(X) by Romberg's method * * --------------------------------------------------- * * INPUTS: * * a begin value of x variable * * b end value of x variable * * prec desired precision * * itermin minimum number of ietrations * * itermax maximum number of iterations * * * * OUTPUTS: * * obtprec obtained precision for integral * * n number of iterations done * * * * RETURNED VALUE the integral of FUNC(X) from a to b * * * ******************************************************/ double RombergIntegral(double a,double b,double prec, double *obtprec, int *n, int itermin, int itermax) { int i,j; double pas,r,s,ta; double t[MAXITER][MAXITER]; if (itermax>MAXITER) itermax=MAXITER; r = FUNC(a); ta = (r + FUNC(b)) / 2; *n=0; pas=b-a; t[0][0]=ta*pas; do { *n = *n + 1; pas = pas/2; s=ta; for (i=1; i=itermax) return t[*n][0]; } while (*obtprec >prec || *n