inte_qng_gsl.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner and Jerry Gagelman
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 #ifndef O2SCL_INTE_QNG_GSL_H
24 #define O2SCL_INTE_QNG_GSL_H
25 
26 /** \file inte_qng_gsl.h
27  \brief File defining \ref o2scl::inte_qng_gsl
28 */
29 
30 #include <o2scl/string_conv.h>
31 #include <o2scl/inte.h>
32 #include <o2scl/inte_gsl.h>
33 #include <o2scl/funct.h>
34 
35 /** \brief A namespace for the quadrature coefficients for
36  non-adaptive integration
37 
38  <b>Documentation from GSL</b>: \n
39  Gauss-Kronrod-Patterson quadrature coefficients for use in
40  quadpack routine qng. These coefficients were calculated with
41  101 decimal digit arithmetic by L. W. Fullerton, Bell Labs, Nov
42  1981.
43 
44 */
46 
47  /** Abcissae common to the 10-, 21-, 43- and 87-point rule */
48  static const double x1[5] = {
49  0.973906528517171720077964012084452,
50  0.865063366688984510732096688423493,
51  0.679409568299024406234327365114874,
52  0.433395394129247190799265943165784,
53  0.148874338981631210884826001129720
54  };
55 
56  /** Weights of the 10-point formula */
57  static const double w10[5] = {
58  0.066671344308688137593568809893332,
59  0.149451349150580593145776339657697,
60  0.219086362515982043995534934228163,
61  0.269266719309996355091226921569469,
62  0.295524224714752870173892994651338
63  };
64 
65  /** Abcissae common to the 21-, 43- and 87-point rule */
66  static const double x2[5] = {
67  0.995657163025808080735527280689003,
68  0.930157491355708226001207180059508,
69  0.780817726586416897063717578345042,
70  0.562757134668604683339000099272694,
71  0.294392862701460198131126603103866
72  };
73 
74  /** Weights of the 21-point formula for abcissae x1 */
75  static const double w21a[5] = {
76  0.032558162307964727478818972459390,
77  0.075039674810919952767043140916190,
78  0.109387158802297641899210590325805,
79  0.134709217311473325928054001771707,
80  0.147739104901338491374841515972068
81  };
82 
83  /** Weights of the 21-point formula for abcissae x2 */
84  static const double w21b[6] = {
85  0.011694638867371874278064396062192,
86  0.054755896574351996031381300244580,
87  0.093125454583697605535065465083366,
88  0.123491976262065851077958109831074,
89  0.142775938577060080797094273138717,
90  0.149445554002916905664936468389821
91  };
92 
93  /** Abscissae common to the 43- and 87-point rule */
94  static const double x3[11] = {
95  0.999333360901932081394099323919911,
96  0.987433402908088869795961478381209,
97  0.954807934814266299257919200290473,
98  0.900148695748328293625099494069092,
99  0.825198314983114150847066732588520,
100  0.732148388989304982612354848755461,
101  0.622847970537725238641159120344323,
102  0.499479574071056499952214885499755,
103  0.364901661346580768043989548502644,
104  0.222254919776601296498260928066212,
105  0.074650617461383322043914435796506
106  };
107 
108  /** Weights of the 43-point formula for abscissae x1, x3 */
109  static const double w43a[10] = {
110  0.016296734289666564924281974617663,
111  0.037522876120869501461613795898115,
112  0.054694902058255442147212685465005,
113  0.067355414609478086075553166302174,
114  0.073870199632393953432140695251367,
115  0.005768556059769796184184327908655,
116  0.027371890593248842081276069289151,
117  0.046560826910428830743339154433824,
118  0.061744995201442564496240336030883,
119  0.071387267268693397768559114425516
120  };
121 
122  /** Weights of the 43-point formula for abscissae x3 */
123  static const double w43b[12] = {
124  0.001844477640212414100389106552965,
125  0.010798689585891651740465406741293,
126  0.021895363867795428102523123075149,
127  0.032597463975345689443882222526137,
128  0.042163137935191811847627924327955,
129  0.050741939600184577780189020092084,
130  0.058379395542619248375475369330206,
131  0.064746404951445885544689259517511,
132  0.069566197912356484528633315038405,
133  0.072824441471833208150939535192842,
134  0.074507751014175118273571813842889,
135  0.074722147517403005594425168280423
136  };
137 
138  /** Abscissae of the 87-point rule */
139  static const double x4[22] = {
140  0.999902977262729234490529830591582,
141  0.997989895986678745427496322365960,
142  0.992175497860687222808523352251425,
143  0.981358163572712773571916941623894,
144  0.965057623858384619128284110607926,
145  0.943167613133670596816416634507426,
146  0.915806414685507209591826430720050,
147  0.883221657771316501372117548744163,
148  0.845710748462415666605902011504855,
149  0.803557658035230982788739474980964,
150  0.757005730685495558328942793432020,
151  0.706273209787321819824094274740840,
152  0.651589466501177922534422205016736,
153  0.593223374057961088875273770349144,
154  0.531493605970831932285268948562671,
155  0.466763623042022844871966781659270,
156  0.399424847859218804732101665817923,
157  0.329874877106188288265053371824597,
158  0.258503559202161551802280975429025,
159  0.185695396568346652015917141167606,
160  0.111842213179907468172398359241362,
161  0.037352123394619870814998165437704
162  };
163 
164  /** Weights of the 87-point formula for abscissae x1, x2, x3 */
165  static const double w87a[21] = {
166  0.008148377384149172900002878448190,
167  0.018761438201562822243935059003794,
168  0.027347451050052286161582829741283,
169  0.033677707311637930046581056957588,
170  0.036935099820427907614589586742499,
171  0.002884872430211530501334156248695,
172  0.013685946022712701888950035273128,
173  0.023280413502888311123409291030404,
174  0.030872497611713358675466394126442,
175  0.035693633639418770719351355457044,
176  0.000915283345202241360843392549948,
177  0.005399280219300471367738743391053,
178  0.010947679601118931134327826856808,
179  0.016298731696787335262665703223280,
180  0.021081568889203835112433060188190,
181  0.025370969769253827243467999831710,
182  0.029189697756475752501446154084920,
183  0.032373202467202789685788194889595,
184  0.034783098950365142750781997949596,
185  0.036412220731351787562801163687577,
186  0.037253875503047708539592001191226
187  };
188 
189  /** Weights of the 87-point formula for abscissae x4 */
190  static const double w87b[23] = {
191  0.000274145563762072350016527092881,
192  0.001807124155057942948341311753254,
193  0.004096869282759164864458070683480,
194  0.006758290051847378699816577897424,
195  0.009549957672201646536053581325377,
196  0.012329447652244853694626639963780,
197  0.015010447346388952376697286041943,
198  0.017548967986243191099665352925900,
199  0.019938037786440888202278192730714,
200  0.022194935961012286796332102959499,
201  0.024339147126000805470360647041454,
202  0.026374505414839207241503786552615,
203  0.028286910788771200659968002987960,
204  0.030052581128092695322521110347341,
205  0.031646751371439929404586051078883,
206  0.033050413419978503290785944862689,
207  0.034255099704226061787082821046821,
208  0.035262412660156681033782717998428,
209  0.036076989622888701185500318003895,
210  0.036698604498456094498018047441094,
211  0.037120549269832576114119958413599,
212  0.037334228751935040321235449094698,
213  0.037361073762679023410321241766599
214  };
215 }
216 
217 #ifndef DOXYGEN_NO_O2NS
218 namespace o2scl {
219 #endif
220 
221  /** \brief Non-adaptive integration from a to b (GSL)
222 
223  The function \ref integ() uses 10-point, 21-point, 43-point, and
224  87-point Gauss-Kronrod integration successively until the
225  integral is returned within the accuracy specified by \ref
226  inte::tol_abs and \ref inte::tol_rel. The 10-point rule is only
227  used to estimate the error for the 21-point rule. The result of
228  the 21-point calculation is used to estimate the error for the
229  43-point rule, and so forth.
230 
231  The error handler is called if the 87-point integration fails
232  to produce the desired accuracy. If \ref inte::err_nonconv is
233  false (the default is true), then the error handler is never
234  called and when the desired accuracy is not obtained the
235  result of the 87-point integration is returned along with
236  the associated error.
237 
238  The return value of the function to be integrated is ignored.
239 
240  The integration fails and the error handler is called if the
241  tolerances are too small, i.e. if either \f$ \mathrm{tol_{abs}}
242  \leq 0 \f$ or if \f$ \mathrm{tol_{rel}} < 50 \cdot
243  \epsilon_\mathrm{double} \f$ where \f$ \epsilon_\mathrm{double}
244  \approx 1.11 \times 10^{-14}) \f$.
245 
246  The integration coefficients are stored in the \ref
247  o2scl_inte_qng_coeffs namespace.
248 
249  \comment
250  \todo Shouldn't feval be last_iter ?
251  7/27/09 - Actually I think no, because the concepts are somewhat
252  different.
253  \comment
254  */
255  template<class func_t=funct11> class inte_qng_gsl :
256  public inte<func_t>, public inte_gsl {
257 
258  public:
259 
260  inte_qng_gsl() {
261  min_rel_tol=0.5e-28;
262  feval=0;
263  }
264 
265  /** \brief The number of function evalutions for the last integration
266 
267  Set to either 0, 21, 43, or 87, depending on the number of
268  function evaluations that were used. This variable is zero if
269  an error occurs before any function evaluations were performed
270  and is never equal 10, since in the 10-point method, the
271  21-point result is used to estimate the error. If the function
272  fails to achieve the desired precision, feval is 87.
273  */
274  size_t feval;
275 
276  /** \brief Minimum relative tolerance for integration
277  (default \f$ 5 \times 10^{-29} \f$ )
278  */
279  double min_rel_tol;
280 
281  /** \brief Integrate function \c func from \c a to \c b
282  giving result \c res and error \c err
283  */
284  virtual int integ_err(func_t &func, double a, double b,
285  double &res, double &err2) {
286 
287  // Array of function values which have been computed
288  double savfun[21];
289 
290  // 10, 21, 43 and 87 point results
291  double res10, res21, res43, res87;
292 
293  double result_kronrod, err;
294 
295  // Approximation to the integral of abs(f)
296  double resabs;
297 
298  // Approximation to the integral of abs(f-i/(b-a))
299  double resasc;
300 
301  // Function values for the computation of resabs and resasc
302  double fv1[5], fv2[5], fv3[5], fv4[5];
303 
304  const double half_length = 0.5 * (b - a);
305  const double abs_half_length = fabs(half_length);
306  const double center = 0.5 * (b + a);
307 
308  double f_center;
309  f_center=func(center);
310 
311  double dbl_eps=std::numeric_limits<double>::epsilon();
312 
313  if (this->tol_abs <= 0 && (this->tol_rel < 50 * dbl_eps ||
314  this->tol_rel < min_rel_tol)) {
315  res = 0;
316  err2 = 0;
317  feval = 0;
318 
319  std::string estr=((std::string)"Tolerance cannot be achieved ")+
320  "with given value of tol_abs, "+o2scl::dtos(this->tol_abs)+
321  ", and tol_rel, "+o2scl::dtos(this->tol_rel)+
322  ", in inte_qng_gsl::integ_err().";
323  O2SCL_ERR(estr.c_str(),exc_ebadtol);
324  };
325 
326  // Compute the integral using the 10- and 21-point formula.
327  // Also, compute resabs and resasc.
328 
329  res10 = 0;
330  res21 = o2scl_inte_qng_coeffs::w21b[5] * f_center;
331  resabs = o2scl_inte_qng_coeffs::w21b[5] * fabs(f_center);
332 
333  for(int k=0; k < 5; k++) {
334  const double abscissa = half_length * o2scl_inte_qng_coeffs::x1[k];
335  double fval1, fval2;
336  fval1=func(center+abscissa);
337  fval2=func(center-abscissa);
338  const double fval = fval1 + fval2;
339  res10 += o2scl_inte_qng_coeffs::w10[k] * fval;
340  res21 += o2scl_inte_qng_coeffs::w21a[k] * fval;
341  resabs += o2scl_inte_qng_coeffs::w21a[k] *
342  (fabs(fval1) + fabs(fval2));
343  savfun[k] = fval;
344  fv1[k] = fval1;
345  fv2[k] = fval2;
346  }
347 
348  for(int k=0; k < 5; k++) {
349  const double abscissa = half_length * o2scl_inte_qng_coeffs::x2[k];
350  double fval1, fval2;
351  fval1=func(center+abscissa);
352  fval2=func(center-abscissa);
353  const double fval = fval1 + fval2;
354  res21 += o2scl_inte_qng_coeffs::w21b[k] * fval;
355  resabs += o2scl_inte_qng_coeffs::w21b[k] *
356  (fabs(fval1) + fabs(fval2));
357  savfun[k + 5] = fval;
358  fv3[k] = fval1;
359  fv4[k] = fval2;
360  }
361 
362  resabs *= abs_half_length;
363 
364  {
365  const double mean = 0.5 * res21;
366 
367  resasc = o2scl_inte_qng_coeffs::w21b[5] * fabs(f_center - mean);
368 
369  for(int k=0; k < 5; k++) {
370  resasc+=(fabs(fv1[k]-mean)+fabs(fv2[k]-mean))*
372  +(fabs(fv3[k]-mean)+fabs(fv4[k]-mean))*
374  }
375  resasc*=abs_half_length;
376  }
377 
378  // Test for convergence.
379 
380  result_kronrod = res21 * half_length;
381  err = rescale_error ((res21 - res10) * half_length, resabs, resasc);
382 
383  if (this->verbose>0) {
384  std::cout << "inte_qng_gsl Iter: " << 1;
385  std::cout.setf(std::ios::showpos);
386  std::cout << " Res: " << result_kronrod;
387  std::cout.unsetf(std::ios::showpos);
388 
389  double ttol=this->tol_rel * fabs(result_kronrod);
390  if (ttol<this->tol_abs) ttol=this->tol_abs;
391 
392  std::cout << " Err: " << err
393  << " Tol: " << ttol << std::endl;
394  if (this->verbose>1) {
395  char ch;
396  std::cout << "Press a key and type enter to continue. " ;
397  std::cin >> ch;
398  }
399  }
400 
401  if (err < this->tol_abs || err < this->tol_rel * fabs(result_kronrod)) {
402  res = result_kronrod;
403  err2 = err;
404  feval = 21;
405  return success;
406  }
407 
408  // Compute the integral using the 43-point formula.
409 
410  res43 = o2scl_inte_qng_coeffs::w43b[11] * f_center;
411 
412  for(int k=0; k < 10; k++) {
413  res43 += savfun[k] * o2scl_inte_qng_coeffs::w43a[k];
414  }
415 
416  for(int k=0; k < 11; k++) {
417  const double abscissa = half_length * o2scl_inte_qng_coeffs::x3[k];
418  double fval1, fval2;
419  fval1=func(center+abscissa);
420  fval2=func(center-abscissa);
421  const double fval = fval1+fval2;
422  res43 += fval * o2scl_inte_qng_coeffs::w43b[k];
423  savfun[k + 10] = fval;
424  }
425 
426  // Test for convergence.
427 
428  result_kronrod = res43 * half_length;
429  err = rescale_error ((res43 - res21) * half_length, resabs, resasc);
430 
431  if (this->verbose>0) {
432  std::cout << "inte_qng_gsl Iter: " << 2;
433  std::cout.setf(std::ios::showpos);
434  std::cout << " Res: " << result_kronrod;
435  std::cout.unsetf(std::ios::showpos);
436 
437  double ttol=this->tol_rel * fabs(result_kronrod);
438  if (ttol<this->tol_abs) ttol=this->tol_abs;
439 
440  std::cout << " Err: " << err
441  << " Tol: " << ttol << std::endl;
442  if (this->verbose>1) {
443  char ch;
444  std::cout << "Press a key and type enter to continue. " ;
445  std::cin >> ch;
446  }
447  }
448 
449  if (err < this->tol_abs || err < this->tol_rel * fabs(result_kronrod)) {
450  res = result_kronrod;
451  err2 = err;
452  feval = 43;
453  return success;
454  }
455 
456  // Compute the integral using the 87-point formula.
457 
458  res87 = o2scl_inte_qng_coeffs::w87b[22] * f_center;
459 
460  for(int k=0; k < 21; k++) {
461  res87 += savfun[k] * o2scl_inte_qng_coeffs::w87a[k];
462  }
463 
464  for(int k=0; k < 22; k++) {
465  const double abscissa = half_length * o2scl_inte_qng_coeffs::x4[k];
466  double fval1, fval2;
467  fval1=func(center+abscissa);
468  fval2=func(center-abscissa);
469  res87 += o2scl_inte_qng_coeffs::w87b[k] * (fval1+fval2);
470  }
471 
472  // Test for convergence.
473 
474  result_kronrod = res87 * half_length;
475  err = rescale_error ((res87 - res43) * half_length, resabs, resasc);
476 
477  if (this->verbose>0) {
478  std::cout << "inte_qng_gsl Iter: " << 3;
479  std::cout.setf(std::ios::showpos);
480  std::cout << " Res: " << result_kronrod;
481  std::cout.unsetf(std::ios::showpos);
482 
483  double ttol=this->tol_rel * fabs(result_kronrod);
484  if (ttol<this->tol_abs) ttol=this->tol_abs;
485 
486  std::cout << " Err: " << err
487  << " Tol: " << ttol << std::endl;
488  if (this->verbose>1) {
489  char ch;
490  std::cout << "Press a key and type enter to continue. " ;
491  std::cin >> ch;
492  }
493  }
494 
495  if (err < this->tol_abs || err < this->tol_rel * fabs(result_kronrod)) {
496  res = result_kronrod;
497  err2 = err;
498  feval = 87;
499  return success;
500  }
501 
502  // Failed to converge
503 
504  res = result_kronrod;
505  err2 = err;
506  feval = 87;
507 
508  std::string estr="Failed to reach tolerance ";
509  estr+="in inte_qng_gsl::integ_err().";
510  O2SCL_CONV_RET(estr.c_str(),exc_etol,this->err_nonconv);
511  }
512 
513  /// Return string denoting type ("inte_qng_gsl")
514  const char *type() { return "inte_qng_gsl"; }
515 
516  };
517 
518 #ifndef DOXYGEN_NO_O2NS
519 }
520 #endif
521 
522 #endif
static const double w21b[6]
Definition: inte_qng_gsl.h:84
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
static const double w43a[10]
Definition: inte_qng_gsl.h:109
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
Definition: err_hnd.h:292
Non-adaptive integration from a to b (GSL)
Definition: inte_qng_gsl.h:255
static const double w10[5]
Definition: inte_qng_gsl.h:57
const char * type()
Return string denoting type ("inte_qng_gsl")
Definition: inte_qng_gsl.h:514
failed to reach the specified tolerance
Definition: err_hnd.h:79
A namespace for the quadrature coefficients for non-adaptive integration.
Definition: inte_qng_gsl.h:45
virtual int integ_err(func_t &func, double a, double b, double &res, double &err2)
Integrate function func from a to b giving result res and error err.
Definition: inte_qng_gsl.h:284
Base integration class [abstract base].
Definition: inte.h:44
static const double x3[11]
Definition: inte_qng_gsl.h:94
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
static const double w87b[23]
Definition: inte_qng_gsl.h:190
GSL integration base.
Definition: inte_gsl.h:62
size_t feval
The number of function evalutions for the last integration.
Definition: inte_qng_gsl.h:274
static const double w87a[21]
Definition: inte_qng_gsl.h:165
static const double x4[22]
Definition: inte_qng_gsl.h:139
static const double x2[5]
Definition: inte_qng_gsl.h:66
static const double w21a[5]
Definition: inte_qng_gsl.h:75
static const double x1[5]
Definition: inte_qng_gsl.h:48
double min_rel_tol
Minimum relative tolerance for integration (default )
Definition: inte_qng_gsl.h:279
user specified an invalid tolerance
Definition: err_hnd.h:77
Success.
Definition: err_hnd.h:47
static const double w43b[12]
Definition: inte_qng_gsl.h:123

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