inte_qag_gsl.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Jerry Gagelman and Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /*
24  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
25  *
26  * This program is free software; you can redistribute it and/or modify
27  * it under the terms of the GNU General Public License as published by
28  * the Free Software Foundation; either version 3 of the License, or (at
29  * your option) any later version.
30  *
31  * This program is distributed in the hope that it will be useful, but
32  * WITHOUT ANY WARRANTY; without even the implied warranty of
33  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
34  * General Public License for more details.
35  *
36  * You should have received a copy of the GNU General Public License
37  * along with this program; if not, write to the Free Software
38  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
39  * 02110-1301, USA.
40  */
41 #ifndef O2SCL_GSL_INTE_QAG_H
42 #define O2SCL_GSL_INTE_QAG_H
43 
44 /** \file inte_qag_gsl.h
45  \brief File defining \ref o2scl::inte_qag_gsl
46 */
47 #include <o2scl/inte.h>
48 #include <o2scl/inte_kronrod_gsl.h>
49 #include <o2scl/funct.h>
50 #include <o2scl/string_conv.h>
51 
52 #ifndef DOXYGEN_NO_O2NS
53 namespace o2scl {
54 #endif
55 
56  /** \brief Adaptive numerical integration of a function (without
57  singularities) on a bounded interval (GSL)
58 
59  Adaptive integration of a univariate function requires two main
60  procedures: approximating the integral on a bounded interval,
61  and estimating the approximation error. The algorithm
62  recursively refines the interval, computing the integral and its
63  error estimate on each subinterval, until the total error
64  estimate over \b all subintervals falls within the
65  user-specified tolerance. The value returned is the sum of the
66  (approximated) integrals over all subintervals.
67 
68  See \ref gslinte_subsect in the User's guide for general
69  information about the GSL integration classes.
70 
71  \future There are a few fine-tuned parameters which should
72  be re-expressed as data members in the convergence tests.
73  \future Should QUADPACK parameters in round-off tests be subject
74  to adjustments by the end user?
75  \future Add functions to examine the contents of the workspace
76  to detect regions where the integrand may be problematic;
77  possibly call these functions automatically depending on
78  verbosity settings.
79  */
80  template<class func_t=funct11> class inte_qag_gsl :
81  public inte_kronrod_gsl<func_t> {
82 
83  public:
84 
85  /// Create an integrator with the specified rule
87  }
88 
89  virtual ~inte_qag_gsl() {}
90 
91  /** \brief Integrate function \c func from \c a to \c b and place
92  the result in \c res and the error in \c err
93 
94  */
95  virtual int integ_err(func_t &func, double a, double b,
96  double &res, double &err) {
97  return qag(func,a,b,this->tol_abs,this->tol_rel,&res,&err);
98  }
99 
100 #ifndef DOXYGEN_INTERNAL
101 
102  protected:
103 
104  /** \brief Perform an adaptive integration given the coefficients,
105  and returning \c result
106 
107  \future Just move this function to integ_err().
108  */
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) {
112 
113  double area, errsum;
114  double result0, abserr0, resabs0, resasc0;
115  double tolerance;
116  size_t iteration = 0;
117  int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;
118 
119  double round_off;
120 
121  /* Initialize results */
122 
123  this->w->initialise(a,b);
124 
125  *result = 0;
126  *abserr = 0;
127 
128  double dbl_eps=std::numeric_limits<double>::epsilon();
129 
130  if (l_epsabs <= 0 &&
131  (l_epsrel < 50 * dbl_eps || l_epsrel < 0.5e-28)) {
132  this->last_iter=0;
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().";
136  O2SCL_ERR(estr.c_str(),exc_ebadtol);
137  }
138 
139  /* perform the first integration */
140 
141  this->gauss_kronrod(func,a,b,&result0,&abserr0,&resabs0,&resasc0);
142 
143  this->w->set_initial_result(result0,abserr0);
144 
145  /* Test on accuracy */
146 
147  tolerance = GSL_MAX_DBL(l_epsabs, l_epsrel * fabs (result0));
148 
149  /* need IEEE rounding here to match original quadpack behavior */
150 
151  round_off=gsl_coerce_double(50 * dbl_eps * resabs0);
152 
153  if (abserr0 <= round_off && abserr0 > tolerance) {
154 
155  *result = result0;
156  *abserr = abserr0;
157 
158  // We start with 1 here, because an integration
159  // was already performed above
160  this->last_iter=1;
161 
162  std::string estr="Cannot reach tolerance because of roundoff ";
163  estr+="error on first attempt in inte_qag_gsl::qag().";
164  O2SCL_CONV_RET(estr.c_str(),exc_eround,this->err_nonconv);
165 
166  } else if ((abserr0 <= tolerance &&
167  abserr0 != resasc0) || abserr0 == 0.0) {
168  *result = result0;
169  *abserr = abserr0;
170 
171  // We start with 1 here, because an integration
172  // was already performed above
173  this->last_iter=1;
174 
175  return success;
176 
177  } else if (this->w->limit == 1) {
178 
179  *result = result0;
180  *abserr = abserr0;
181 
182  // We start with 1 here, because an integration
183  // was already performed above
184  this->last_iter=1;
185 
186  O2SCL_CONV2_RET("A maximum of 1 iteration was insufficient ",
187  "in inte_qag_gsl::qag().",
188  exc_emaxiter,this->err_nonconv);
189  }
190 
191  area = result0;
192  errsum = abserr0;
193 
194  iteration = 1;
195  do {
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;
202 
203  /* Bisect the subinterval with the largest error estimate */
204 
205  this->w->retrieve (&a_i, &b_i, &r_i, &e_i);
206 
207  a1 = a_i;
208  b1 = 0.5 * (a_i + b_i);
209  a2 = b1;
210  b2 = b_i;
211 
212  this->gauss_kronrod(func,a1,b1,&area1,&error1,&resabs1,&resasc1);
213  this->gauss_kronrod(func,a2,b2,&area2,&error2,&resabs2,&resasc2);
214 
215  area12 = area1 + area2;
216  error12 = error1 + error2;
217 
218  errsum += (error12 - e_i);
219  area += area12 - r_i;
220 
221  if (resasc1 != error1 && resasc2 != error2) {
222  double delta = r_i - area12;
223 
224  if (fabs (delta) <= 1.0e-5 * fabs (area12) &&
225  error12 >= 0.99 * e_i) {
226  roundoff_type1++;
227  }
228  if (iteration >= 10 && error12 > e_i) {
229  roundoff_type2++;
230  }
231  }
232 
233  tolerance = GSL_MAX_DBL (l_epsabs, l_epsrel * fabs (area));
234 
235  if (errsum > tolerance) {
236  if (roundoff_type1 >= 6 || roundoff_type2 >= 20) {
237  // round off error
238  error_type = 2;
239  }
240  /* set error flag in the case of bad integrand behaviour at
241  a point of the integration range */
242 
243  if (this->w->subinterval_too_small (a1, a2, b2)) {
244  error_type = 3;
245  }
246  }
247 
248  this->w->update (a1, b1, area1, error1, a2, b2, area2, error2);
249 
250  this->w->retrieve (&a_i, &b_i, &r_i, &e_i);
251 
252  if (this->verbose>0) {
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;
259  if (this->verbose>1) {
260  char ch;
261  std::cout << "Press a key and type enter to continue. " ;
262  std::cin >> ch;
263  }
264  }
265 
266  iteration++;
267 
268  } while (iteration < this->w->limit && !error_type &&
269  errsum > tolerance);
270 
271  *result = this->w->sum_results();
272  *abserr = errsum;
273 
274  this->last_iter=iteration;
275 
276  if (errsum <= tolerance) {
277  return success;
278  } else if (error_type == 2) {
279  std::string estr="Roundoff error prevents tolerance ";
280  estr+="from being achieved in inte_qag_gsl::qag().";
281  O2SCL_CONV_RET(estr.c_str(),exc_eround,this->err_nonconv);
282  } else if (error_type == 3) {
283  std::string estr="Bad integrand behavior ";
284  estr+=" in inte_qag_gsl::qag().";
285  O2SCL_CONV_RET(estr.c_str(),exc_esing,this->err_nonconv);
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().";
289  O2SCL_CONV_RET(estr.c_str(),exc_emaxiter,this->err_nonconv);
290  } else {
291  std::string estr="Could not integrate function in inte_qag_gsl::";
292  estr+="qag() (it may have returned a non-finite result).";
293  O2SCL_ERR(estr.c_str(),exc_efailed);
294  }
295 
296  // No return statement needed since the above if statement
297  // always forces a return, but some compilers like having one
298  // anyway.
299  return o2scl::success;
300  }
301 
302 #endif
303 
304  public:
305 
306  /// Return string denoting type ("inte_qag_gsl")
307  const char *type() { return "inte_qag_gsl"; }
308 
309  };
310 
311 #ifndef DOXYGEN_NO_O2NS
312 }
313 #endif
314 
315 #endif
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
Definition: err_hnd.h:292
apparent singularity detected
Definition: err_hnd.h:93
exceeded max number of iterations
Definition: err_hnd.h:73
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
Definition: err_hnd.h:298
size_t last_iter
The most recent number of iterations taken.
Definition: inte.h:63
bool err_nonconv
If true, call the error handler if the routine does not converge or reach the desired tolerance (defa...
Definition: inte.h:81
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.
generic failure
Definition: err_hnd.h:61
inte_qag_gsl()
Create an integrator with the specified rule.
Definition: inte_qag_gsl.h:86
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 )
Definition: inte.h:73
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")
Definition: inte_qag_gsl.h:307
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.
Definition: inte_qag_gsl.h:109
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.
Definition: err_hnd.h:273
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 )
Definition: inte.h:68
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.
Definition: inte_qag_gsl.h:95
Adaptive numerical integration of a function (without singularities) on a bounded interval (GSL) ...
Definition: inte_qag_gsl.h:80
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
Definition: err_hnd.h:77
Success.
Definition: err_hnd.h:47
int verbose
Verbosity.
Definition: inte.h:60
failed because of roundoff error
Definition: err_hnd.h:87

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).