// file: qags_test.cpp // // C++ Program to test the qags automatic integrator from // the gsl numerical library. // // Programmer: Dick Furnstahl furnstahl.1@osu.edu // // Revision history: // 12/26/03 original C++ version, modified from C version // // Notes: // * Example taken from the GNU Scientific Library Reference Manual // Edition 1.1, for GSL Version 1.1 9 January 2002 // URL: gsl/ref/gsl-ref_23.html#SEC364 // * Compile and link with: // g++ -Wall -o qags_test qags_test.cpp -lgsl -lgslcblas // * gsl routines have built-in // extern "C" { //
// } // so they can be called from C++ programs without modification // //*********************************************************************// // The following details are taken from the GSL documentation // // Each algorithm computes an approximation to a definite integral of // the form, // // I = \int_a^b f(x) w(x) dx // // where w(x) is a weight function (for general integrands w(x)=1). The // user provides absolute and relative error bounds (epsabs, epsrel) // which specify the following accuracy requirement, // // |RESULT - I| <= max(epsabs, epsrel |I|) // // where RESULT is the numerical approximation obtained by the // algorithm. The algorithms attempt to estimate the absolute error // ABSERR = |RESULT - I| in such a way that the following inequality // holds, // // |RESULT - I| <= ABSERR <= max(epsabs, epsrel |I|) // // The routines will fail to converge if the error bounds are too // stringent, but always return the best approximation obtained up to // that stage. // // // QAGS adaptive integration with singularities // // Function: int gsl_integration_qags (const gsl_function * f, // double a, double b, // double epsabs, double epsrel, // size_t limit, // gsl_integration_workspace * workspace, // double *result, // double *abserr) // // This function applies the Gauss-Kronrod 21-point integration rule // adaptively until an estimate of the integral of f over (a,b) is // achieved within the desired absolute and relative error limits, // epsabs and epsrel. The results are extrapolated using the // epsilon-algorithm, which accelerates the convergence of the integral // in the presence of discontinuities and integrable singularities. The // function returns the final approximation from the extrapolation, // result, and an estimate of the absolute error, abserr. The // subintervals and their results are stored in the memory provided by // workspace. The maximum number of subintervals is given by limit, // which may not exceed the allocated size of the workspace. // // // The integrator QAGS will handle a large class of definite integrals. // For example, consider the following integral, which has a // algebraic-logarithmic singularity at the origin, // // \int_0^1 x^{-1/2} log(x) dx = -4 // // The program below computes this integral to a relative accuracy bound // of 1e-8. // // // // The results below show that the desired accuracy is achieved after 8 // subdivisions. // // result = -3.999999999999973799 // exact result = -4.000000000000000000 // estimated error = 0.000000000000499600 // actual error = 0.000000000000026201 // intervals = 8 // // In fact, the extrapolation procedure used by QAGS produces an // accuracy of many more digits. The error estimate returned // by the extrapolation procedure is larger than the actual error, // giving a margin of safety of one order of magnitude. // //*********************************************************************// // include files #include #include #include #include using namespace std; #include // function prototypes double my_integrand (double x, void *params); //*********************************************************************// int main (void) { gsl_integration_workspace *work_ptr = gsl_integration_workspace_alloc (1000); double lower_limit = 0; /* lower limit a */ double upper_limit = 1; /* upper limit b */ double abs_error = 1.0e-8; /* to avoid round-off problems */ double rel_error = 1.0e-8; /* the result will usually be much better */ double result; /* the result from the integration */ double error; /* the estimated error from the integration */ double alpha = 1.0; // parameter in integrand double expected = -4.0; // exact answer gsl_function My_function; void *params_ptr = α My_function.function = &my_integrand; My_function.params = params_ptr; gsl_integration_qags (&My_function, lower_limit, upper_limit, abs_error, rel_error, 1000, work_ptr, &result, &error); cout.setf (ios::fixed, ios::floatfield); // output in fixed format cout.precision (18); // 18 digits in doubles int width = 20; // setw width for output cout << "result = " << setw(width) << result << endl; cout << "exact result = " << setw(width) << expected << endl; cout << "estimated error = " << setw(width) << error << endl; cout << "actual error = " << setw(width) << result - expected << endl; cout << "intervals = " << work_ptr->size << endl; return 0; } //*********************************************************************// double my_integrand (double x, void *params) { // Mathematica form: Log[alpha*x]/Sqrt[x] // The next line recovers alpha from the passed params pointer double alpha = *(double *) params; return (log (alpha * x) / sqrt (x)); }