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 The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
apparent singularity detected
exceeded max number of iterations
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
size_t last_iter
The most recent number of iterations taken.
bool err_nonconv
If true, call the error handler if the routine does not converge or reach the desired tolerance (defa...
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(func_t &func, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
Integration wrapper for user-specified function type.
inte_qag_gsl()
Create an integrator with the specified rule.
int subinterval_too_small(double a1, double a2, double b2)
Test whether the proposed subdivision falls before floating-point precision.
double tol_abs
The maximum absolute uncertainty in the value of the integral (default )
int initialise(double a, double b)
Initialize the workspace for an integration with limits a and b.
const char * type()
Return string denoting type ("inte_qag_gsl")
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.
int set_initial_result(double result, double error)
Update the workspace with the result and error from the first integration.
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Basic Gauss-Kronrod integration class (GSL)
inte_workspace_gsl * w
The integration workspace.
double sum_results()
Add up all of the contributions to construct the final result.
double tol_rel
The maximum relative uncertainty in the value of the integral (default )
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.
Adaptive numerical integration of a function (without singularities) on a bounded interval (GSL) ...
size_t limit
Maximum number of subintervals allocated.
int retrieve(double *a, double *b, double *r, double *e) const
Retrieve the ith result from the workspace stack.
std::string itos(int x)
Convert an integer to a string.
user specified an invalid tolerance
failed because of roundoff error