Go to the documentation of this file.
41 #ifndef O2SCL_GSL_INTE_QAG_H
42 #define O2SCL_GSL_INTE_QAG_H
47 #include <o2scl/inte.h>
48 #include <o2scl/inte_kronrod_gsl.h>
49 #include <o2scl/funct.h>
50 #include <o2scl/string_conv.h>
52 #ifndef DOXYGEN_NO_O2NS
95 virtual int integ_err(func_t &func,
double a,
double b,
96 double &res,
double &err) {
100 #ifndef DOXYGEN_INTERNAL
109 int qag(func_t &func,
const double a,
const double b,
110 const double l_epsabs,
const double l_epsrel,
111 double *result,
double *abserr) {
114 double result0, abserr0, resabs0, resasc0;
116 size_t iteration = 0;
117 int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;
128 double dbl_eps=std::numeric_limits<double>::epsilon();
131 (l_epsrel < 50 * dbl_eps || l_epsrel < 0.5e-28)) {
133 std::string estr=
"Tolerance cannot be achieved with given ";
134 estr+=
"value of tol_abs, "+
dtos(l_epsabs)+
", and tol_rel, "+
135 dtos(l_epsrel)+
", in inte_qag_gsl::qag().";
141 this->
gauss_kronrod(func,a,b,&result0,&abserr0,&resabs0,&resasc0);
147 tolerance = GSL_MAX_DBL(l_epsabs, l_epsrel * fabs (result0));
151 round_off=gsl_coerce_double(50 * dbl_eps * resabs0);
153 if (abserr0 <= round_off && abserr0 > tolerance) {
162 std::string estr=
"Cannot reach tolerance because of roundoff ";
163 estr+=
"error on first attempt in inte_qag_gsl::qag().";
166 }
else if ((abserr0 <= tolerance &&
167 abserr0 != resasc0) || abserr0 == 0.0) {
177 }
else if (this->
w->
limit == 1) {
187 "in inte_qag_gsl::qag().",
196 double a1, b1, a2, b2;
197 double a_i, b_i, r_i, e_i;
198 double area1 = 0, area2 = 0, area12 = 0;
199 double error1 = 0, error2 = 0, error12 = 0;
200 double resasc1, resasc2;
201 double resabs1, resabs2;
205 this->
w->
retrieve (&a_i, &b_i, &r_i, &e_i);
208 b1 = 0.5 * (a_i + b_i);
212 this->
gauss_kronrod(func,a1,b1,&area1,&error1,&resabs1,&resasc1);
213 this->
gauss_kronrod(func,a2,b2,&area2,&error2,&resabs2,&resasc2);
215 area12 = area1 + area2;
216 error12 = error1 + error2;
218 errsum += (error12 - e_i);
219 area += area12 - r_i;
221 if (resasc1 != error1 && resasc2 != error2) {
222 double delta = r_i - area12;
224 if (fabs (delta) <= 1.0e-5 * fabs (area12) &&
225 error12 >= 0.99 * e_i) {
228 if (iteration >= 10 && error12 > e_i) {
233 tolerance = GSL_MAX_DBL (l_epsabs, l_epsrel * fabs (area));
235 if (errsum > tolerance) {
236 if (roundoff_type1 >= 6 || roundoff_type2 >= 20) {
248 this->
w->
update (a1, b1, area1, error1, a2, b2, area2, error2);
250 this->
w->
retrieve (&a_i, &b_i, &r_i, &e_i);
253 std::cout <<
"inte_qag_gsl Iter: " << iteration;
254 std::cout.setf(std::ios::showpos);
255 std::cout <<
" Res: " << area;
256 std::cout.unsetf(std::ios::showpos);
257 std::cout <<
" Err: " << errsum
258 <<
" Tol: " << tolerance << std::endl;
261 std::cout <<
"Press a key and type enter to continue. " ;
268 }
while (iteration < this->
w->
limit && !error_type &&
276 if (errsum <= tolerance) {
278 }
else if (error_type == 2) {
279 std::string estr=
"Roundoff error prevents tolerance ";
280 estr+=
"from being achieved in inte_qag_gsl::qag().";
282 }
else if (error_type == 3) {
283 std::string estr=
"Bad integrand behavior ";
284 estr+=
" in inte_qag_gsl::qag().";
286 }
else if (iteration == this->
w->
limit) {
287 std::string estr=
"Maximum number of subdivisions ("+
itos(iteration);
288 estr+=
") reached in inte_qag_gsl::qag().";
291 std::string estr=
"Could not integrate function in inte_qag_gsl::";
292 estr+=
"qag() (it may have returned a non-finite result).";
307 const char *
type() {
return "inte_qag_gsl"; }
311 #ifndef DOXYGEN_NO_O2NS
@ exc_esing
apparent singularity detected
virtual int integ_err(func_t &func, double a, double b, double &res, double &err)
Integrate function func from a to b and place the result in res and the error in err.
inte_qag_gsl()
Create an integrator with the specified rule.
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
std::string dtos(const fp_t &x, int prec=6, bool auto_prec=false)
Convert a double to a string.
inte_workspace_gsl * w
The integration workspace.
Adaptive numerical integration of a function (without singularities) on a bounded interval (GSL)
int retrieve(double *a, double *b, double *r, double *e) const
Retrieve the ith result from the workspace stack.
int set_initial_result(double result, double error)
Update the workspace with the result and error from the first integration.
@ exc_efailed
generic failure
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
@ exc_emaxiter
exceeded max number of iterations
int initialise(double a, double b)
Initialize the workspace for an integration with limits a and b.
int qag(func_t &func, const double a, const double b, const double l_epsabs, const double l_epsrel, double *result, double *abserr)
Perform an adaptive integration given the coefficients, and returning result.
double sum_results()
Add up all of the contributions to construct the final result.
size_t last_iter
The most recent number of iterations taken.
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
int update(double a1, double b1, double area1, double error1, double a2, double b2, double area2, double error2)
Determine which new subinterval to add to the workspace stack and perform update.
virtual void gauss_kronrod(funct &func, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
Integration wrapper for user-specified function type.
double tol_abs
The maximum absolute uncertainty in the value of the integral (default )
std::string itos(int x)
Convert an integer to a string.
double tol_rel
The maximum relative uncertainty in the value of the integral (default )
@ exc_eround
failed because of roundoff error
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
size_t limit
Maximum number of subintervals allocated.
Basic Gauss-Kronrod integration class (GSL)
bool err_nonconv
If true, call the error handler if the routine does not converge or reach the desired tolerance (defa...
const char * type()
Return string denoting type ("inte_qag_gsl")
@ exc_ebadtol
user specified an invalid tolerance
int subinterval_too_small(double a1, double a2, double b2)
Test whether the proposed subdivision falls before floating-point precision.
Documentation generated with Doxygen. Provided under the
GNU Free Documentation License (see License Information).