inte_qawo_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 #ifndef O2SCL_GSL_INTE_QAWO_H
24 #define O2SCL_GSL_INTE_QAWO_H
25 
26 /** \file inte_qawo_gsl.h
27  \brief File defining \ref o2scl::inte_qawo_gsl_sin and
28  \ref o2scl::inte_qawo_gsl_cos
29 */
30 
31 #include <o2scl/inte.h>
32 #include <o2scl/inte_qawc_gsl.h>
33 
34 #ifndef DOXYGEN_NO_O2NS
35 namespace o2scl {
36 #endif
37 
38  /** \brief Adaptive integration for oscillatory integrals (GSL)
39 
40  The integral
41  \f[
42  \int_a^b f(x) \sin (\omega x)~dx
43  \f]
44  is computed for some frequency parameter \f$ \omega \f$,
45  stored in \ref inte_qawo_gsl_sin::omega .
46 
47  An adaptive algorithm together with an series-acceleration
48  method like that of \ref inte_qags_gsl is used. Those
49  subintervals with "large" widths \f$ d \equiv b-a \f$ where \f$
50  d\omega > 4 \f$ are computed using a 25-point Clenshaw-Curtis
51  integration rule to handle the oscillatory behavior. In order to
52  work efficiently, the Chebyshev moments for the particular
53  weight function \f$ W \f$ are computed in advance.
54 
55  See \ref gslinte_subsect in the User's guide for general
56  information about the GSL integration classes.
57  */
58  template<class func_t> class inte_qawo_gsl_sin :
59  public inte_cheb_gsl<func_t> {
60 
61  public:
62 
64  n_levels=10;
65  omega=1.0;
66  }
67 
68  virtual ~inte_qawo_gsl_sin() {}
69 
70  /// The user-specified frequency (default 1.0)
71  double omega;
72 
73  /// The number of bisection levels (default 10)
74  size_t n_levels;
75 
76  /** \brief Integrate function \c func from \c a to \c b and place
77  the result in \c res and the error in \c err
78  */
79  virtual int integ_err(func_t &func, double a, double b,
80  double &res, double &err) {
81 
82  otable=gsl_integration_qawo_table_alloc
83  (omega,b-a,GSL_INTEG_SINE,n_levels);
84 
85  int status=qawo(func,a,this->tol_abs,this->tol_rel,
86  this->w,this->otable,&res,&err);
87 
88  gsl_integration_qawo_table_free(otable);
89 
90  return status;
91  }
92 
93 #ifndef DOXYGEN_INTERNAL
94 
95  protected:
96 
97  /// The integration workspace
98  gsl_integration_qawo_table *otable;
99 
100  /** \brief The full GSL integration routine called by integ_err()
101 
102  \future Remove goto statements.
103  */
104  int qawo(func_t &func, const double a, const double epsabs,
105  const double epsrel, inte_workspace_gsl *loc_w,
106  gsl_integration_qawo_table *wf, double *result, double *abserr) {
107 
108  double area, errsum;
109  double res_ext, err_ext;
110  double result0=0.0, abserr0=0.0, resabs0=0.0, resasc0=0.0;
111  double tolerance;
112 
113  double ertest = 0;
114  double error_over_large_intervals = 0;
115  double reseps = 0, abseps = 0, correc = 0;
116  size_t ktmin = 0;
117  int roundoff_type1 = 0, roundoff_type2 = 0, roundoff_type3 = 0;
118  int error_type = 0, error_type2 = 0;
119 
120  size_t iteration = 0;
121 
122  int positive_integrand = 0;
123  int extrapolate = 0;
124  int extall = 0;
125  int disallow_extrapolation = 0;
126 
128 
129  double b = a + wf->L ;
130  double abs_omega = fabs (wf->omega) ;
131 
132  /* Initialize results */
133 
134  loc_w->initialise(a,b);
135 
136  *result = 0;
137  *abserr = 0;
138 
139  size_t limit=this->w->limit;
140 
141  /* Test on accuracy */
142 
143  double dbl_eps=std::numeric_limits<double>::epsilon();
144 
145  if (epsabs <= 0 && (epsrel < 50 * dbl_eps || epsrel < 0.5e-28)) {
146  this->last_iter=0;
147  std::string estr="Tolerance cannot be achieved with given ";
148  estr+="value of tol_abs, "+dtos(epsabs)+", and tol_rel, "+
149  dtos(epsrel)+", in inte_qawo_gsl_sin::qawo().";
150  O2SCL_ERR(estr.c_str(),exc_ebadtol);
151  }
152 
153  /* Perform the first integration */
154 
155  this->qc25f(func, a, b, wf, 0,
156  &result0, &abserr0, &resabs0, &resasc0);
157 
158  loc_w->set_initial_result (result0, abserr0);
159 
160  tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
161 
162  if (this->verbose>0) {
163  std::cout << "inte_qawo_gsl Iter: " << 1;
164  std::cout.setf(std::ios::showpos);
165  std::cout << " Res: " << result0;
166  std::cout.unsetf(std::ios::showpos);
167  std::cout << " Err: " << abserr0
168  << " Tol: " << tolerance << std::endl;
169  if (this->verbose>1) {
170  char ch;
171  std::cout << "Press a key and type enter to continue. " ;
172  std::cin >> ch;
173  }
174  }
175 
176  if (abserr0 <= 100 * GSL_DBL_EPSILON * resabs0 &&
177  abserr0 > tolerance) {
178  *result = result0;
179  *abserr = abserr0;
180 
181  this->last_iter=1;
182  std::string estr="Cannot reach tolerance because of roundoff error ";
183  estr+="on first attempt in inte_qawo_gsl_sin::qawo().";
184  O2SCL_CONV_RET(estr.c_str(),exc_eround,this->err_nonconv);
185  } else if ((abserr0 <= tolerance && abserr0 != resasc0) ||
186  abserr0 == 0.0) {
187  *result = result0;
188  *abserr = abserr0;
189 
190  this->last_iter=1;
191  return success;
192  } else if (limit == 1) {
193  *result = result0;
194  *abserr = abserr0;
195 
196  this->last_iter=1;
197  std::string estr="A maximum of 1 iteration was insufficient ";
198  estr+="in inte_qawo_gsl_sin::qawo().";
199  O2SCL_CONV_RET(estr.c_str(),exc_emaxiter,this->err_nonconv);
200  }
201 
202  /* Initialization */
203 
204  this->initialise_table(&table);
205 
206  if (0.5 * abs_omega * fabs(b - a) <= 2) {
207  this->append_table (&table, result0);
208  extall = 1;
209  }
210 
211  area = result0;
212  errsum = abserr0;
213 
214  res_ext = result0;
215  err_ext = GSL_DBL_MAX;
216 
217  positive_integrand = this->test_positivity (result0, resabs0);
218 
219  iteration = 1;
220 
221  do {
222 
223  size_t current_level;
224  double a1, b1, a2, b2;
225  double a_i, b_i, r_i, e_i;
226  double area1 = 0, area2 = 0, area12 = 0;
227  double error1 = 0, error2 = 0, error12 = 0;
228  double resasc1=0.0, resasc2=0.0;
229  double resabs1=0.0, resabs2=0.0;
230  double last_e_i;
231 
232  /* Bisect the subinterval with the largest error estimate */
233 
234  loc_w->retrieve (&a_i, &b_i, &r_i, &e_i);
235 
236  current_level = loc_w->level[loc_w->i] + 1;
237 
238  if (current_level >= wf->n) {
239  error_type = -1 ; /* exceeded limit of table */
240  break ;
241  }
242 
243  a1 = a_i;
244  b1 = 0.5 * (a_i + b_i);
245  a2 = b1;
246  b2 = b_i;
247 
248  iteration++;
249 
250  qc25f(func, a1, b1, wf, current_level,
251  &area1, &error1, &resabs1, &resasc1);
252  qc25f(func, a2, b2, wf, current_level,
253  &area2, &error2, &resabs2, &resasc2);
254 
255  area12 = area1 + area2;
256  error12 = error1 + error2;
257  last_e_i = e_i;
258 
259  /* Improve previous approximations to the integral and test for
260  accuracy.
261 
262  We write these expressions in the same way as the original
263  QUADPACK code so that the rounding errors are the same, which
264  makes testing easier. */
265 
266  errsum = errsum + error12 - e_i;
267  area = area + area12 - r_i;
268 
269  tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
270 
271  if (this->verbose>0) {
272  std::cout << "inte_qawo_gsl Iter: " << iteration;
273  std::cout.setf(std::ios::showpos);
274  std::cout << " Res: " << area;
275  std::cout.unsetf(std::ios::showpos);
276  std::cout << " Err: " << errsum
277  << " Tol: " << tolerance << std::endl;
278  if (this->verbose>1) {
279  char ch;
280  std::cout << "Press a key and type enter to continue. " ;
281  std::cin >> ch;
282  }
283  }
284 
285  if (resasc1 != error1 && resasc2 != error2) {
286 
287  double delta = r_i - area12;
288 
289  if (fabs (delta) <= 1.0e-5 * fabs (area12) &&
290  error12 >= 0.99 * e_i) {
291 
292  if (!extrapolate) {
293  roundoff_type1++;
294  } else {
295  roundoff_type2++;
296  }
297  }
298 
299  if (iteration > 10 && error12 > e_i) {
300  roundoff_type3++;
301  }
302  }
303 
304  /* Test for roundoff and eventually set error flag */
305 
306  if (roundoff_type1 + roundoff_type2 >= 10 || roundoff_type3 >= 20) {
307  error_type = 2; /* round off error */
308  }
309 
310  if (roundoff_type2 >= 5) {
311  error_type2 = 1;
312  }
313 
314  /* set error flag in the case of bad integrand behaviour at
315  a point of the integration range */
316 
317  if (loc_w->subinterval_too_small (a1, a2, b2)) {
318  error_type = 4;
319  }
320 
321  /* append the newly-created intervals to the list */
322 
323  loc_w->update (a1, b1, area1, error1, a2, b2, area2, error2);
324 
325  if (errsum <= tolerance) {
326  goto compute_result;
327  }
328 
329  if (error_type) {
330  break;
331  }
332 
333  if (iteration >= limit - 1) {
334  error_type = 1;
335  break;
336  }
337 
338  /* set up variables on first iteration */
339 
340  if (iteration == 2 && extall) {
341  error_over_large_intervals = errsum;
342  ertest = tolerance;
343  this->append_table (&table, area);
344  continue;
345  }
346 
347  if (disallow_extrapolation) {
348  continue;
349  }
350 
351  if (extall) {
352  error_over_large_intervals += -last_e_i;
353 
354  if (current_level < loc_w->maximum_level)
355  {
356  error_over_large_intervals += error12;
357  }
358 
359  if (extrapolate) {
360  goto label70;
361  }
362  }
363 
364  if (this->large_interval(loc_w)) {
365  continue;
366  }
367 
368  if (extall) {
369  extrapolate = 1;
370  loc_w->nrmax = 1;
371  } else {
372  /* test whether the interval to be bisected next is the
373  smallest interval. */
374  size_t i = loc_w->i;
375  double width = loc_w->blist[i] - loc_w->alist[i];
376 
377  if (0.25 * fabs(width) * abs_omega > 2) {
378  continue;
379  }
380 
381  extall = 1;
382  error_over_large_intervals = errsum;
383  ertest = tolerance;
384  continue;
385  }
386 
387  label70:
388 
389  if (!error_type2 && error_over_large_intervals > ertest) {
390  if (this->increase_nrmax (loc_w)) {
391  continue;
392  }
393  }
394 
395  /* Perform extrapolation */
396 
397  this->append_table (&table, area);
398 
399  if (table.n < 3) {
400  this->reset_nrmax(loc_w);
401  extrapolate = 0;
402  error_over_large_intervals = errsum;
403  continue;
404  }
405 
406  this->qelg (&table, &reseps, &abseps);
407 
408  ktmin++;
409 
410  if (ktmin > 5 && err_ext < 0.001 * errsum) {
411  error_type = 5;
412  }
413 
414  if (abseps < err_ext) {
415  ktmin = 0;
416  err_ext = abseps;
417  res_ext = reseps;
418  correc = error_over_large_intervals;
419  ertest = GSL_MAX_DBL (epsabs, epsrel * fabs (reseps));
420  if (err_ext <= ertest) {
421  break;
422  }
423  }
424 
425  /* Prepare bisection of the smallest interval. */
426 
427  if (table.n == 1) {
428  disallow_extrapolation = 1;
429  }
430 
431  if (error_type == 5) {
432  break;
433  }
434 
435  /* work on interval with largest error */
436 
437  this->reset_nrmax (loc_w);
438  extrapolate = 0;
439  error_over_large_intervals = errsum;
440 
441  } while (iteration < limit);
442 
443  *result = res_ext;
444  *abserr = err_ext;
445 
446  if (err_ext == GSL_DBL_MAX)
447  goto compute_result;
448 
449  if (error_type || error_type2) {
450 
451  if (error_type2) {
452  err_ext += correc;
453  }
454 
455  if (error_type == 0)
456  error_type = 3;
457 
458  if (result != 0 && area != 0) {
459  if (err_ext / fabs (res_ext) > errsum / fabs (area)) {
460  goto compute_result;
461  }
462  } else if (err_ext > errsum) {
463  goto compute_result;
464  } else if (area == 0.0) {
465  goto return_error;
466  }
467  }
468 
469  /* Test on divergence. */
470 
471  {
472  double max_area = GSL_MAX_DBL (fabs (res_ext), fabs (area));
473 
474  if (!positive_integrand && max_area < 0.01 * resabs0) {
475  goto return_error;
476  }
477  }
478 
479  {
480  double ratio = res_ext / area;
481 
482  if (ratio < 0.01 || ratio > 100 || errsum > fabs (area)) {
483  error_type = 6;
484  }
485  }
486 
487  goto return_error;
488 
489  compute_result:
490 
491  *result = loc_w->sum_results();
492  *abserr = errsum;
493 
494  return_error:
495 
496  if (error_type > 2) {
497  error_type--;
498  }
499 
500  this->last_iter=iteration;
501 
502  if (error_type == 0) {
503  return success;
504  } else if (error_type == 1) {
505  std::string estr="Number of iterations was insufficient ";
506  estr+=" in inte_qawo_gsl_sin::qawo().";
507  O2SCL_CONV_RET(estr.c_str(),exc_emaxiter,this->err_nonconv);
508  } else if (error_type == 2) {
509  std::string estr="Roundoff error prevents tolerance ";
510  estr+="from being achieved in inte_qawo_gsl_sin::qawo().";
511  O2SCL_CONV_RET(estr.c_str(),exc_eround,this->err_nonconv);
512  } else if (error_type == 3) {
513  std::string estr="Bad integrand behavior ";
514  estr+=" in inte_qawo_gsl_sin::qawo().";
515  O2SCL_CONV_RET(estr.c_str(),exc_esing,this->err_nonconv);
516  } else if (error_type == 4) {
517  std::string estr="Roundoff error detected in extrapolation table ";
518  estr+="in inte_qawo_gsl_sin::qawo().";
519  O2SCL_CONV_RET(estr.c_str(),exc_eround,this->err_nonconv);
520  } else if (error_type == 5) {
521  std::string estr="Integral is divergent or slowly convergent ";
522  estr+="in inte_qawo_gsl_sin::qawo().";
523  O2SCL_CONV_RET(estr.c_str(),exc_ediverge,this->err_nonconv);
524  } else if (error_type == -1) {
525  std::string estr="Exceeded limit of trigonometric table ";
526  estr+="inte_qawo_gsl_sin::qawo()";
527  O2SCL_ERR(estr.c_str(),exc_etable);
528  } else {
529  std::string estr="Could not integrate function in inte_qawo_gsl";
530  estr+="::qawo() (it may have returned a non-finite result).";
531  O2SCL_ERR(estr.c_str(),exc_efailed);
532  }
533 
534  // No return statement needed since the above if statement
535  // always forces a return, but some compilers like having one
536  // anyway.
537  return o2scl::success;
538  }
539 
540  /// 25-point quadrature for oscillating functions
541  void qc25f(func_t &func, double a, double b,
542  gsl_integration_qawo_table *wf, size_t level,
543  double *result, double *abserr, double *resabs,
544  double *resasc) {
545 
546  const double center = 0.5 * (a + b);
547  const double half_length = 0.5 * (b - a);
548 
549  const double par = omega * half_length;
550 
551  if (fabs (par) < 2) {
552 
553  this->gauss_kronrod(func,a,b,result,abserr,resabs,resasc);
554 
555  return;
556 
557  } else {
558 
559  double *moment;
560  double cheb12[13], cheb24[25];
561  double result_abs, res12_cos, res12_sin, res24_cos, res24_sin;
562  double est_cos, est_sin;
563  double c, s;
564  size_t i;
565 
566  this->inte_cheb_series(func, a, b, cheb12, cheb24);
567 
568  if (level >= wf->n) {
569  /* table overflow should not happen, check before calling */
570  O2SCL_ERR("Table overflow in inte_qawo_gsl::qc25f().",
571  exc_esanity);
572  return;
573  }
574 
575  /* obtain moments from the table */
576 
577  moment = wf->chebmo + 25 * level;
578 
579  res12_cos = cheb12[12] * moment[12];
580  res12_sin = 0 ;
581 
582  for (i = 0; i < 6; i++) {
583  size_t k = 10 - 2 * i;
584  res12_cos += cheb12[k] * moment[k];
585  res12_sin += cheb12[k + 1] * moment[k + 1];
586  }
587 
588  res24_cos = cheb24[24] * moment[24];
589  res24_sin = 0 ;
590 
591  result_abs = fabs(cheb24[24]) ;
592 
593  for (i = 0; i < 12; i++) {
594  size_t k = 22 - 2 * i;
595  res24_cos += cheb24[k] * moment[k];
596  res24_sin += cheb24[k + 1] * moment[k + 1];
597  result_abs += fabs(cheb24[k]) + fabs(cheb24[k+1]);
598  }
599 
600  est_cos = fabs(res24_cos - res12_cos);
601  est_sin = fabs(res24_sin - res12_sin);
602 
603  c = half_length * cos(center * omega);
604  s = half_length * sin(center * omega);
605 
606  if (wf->sine == GSL_INTEG_SINE) {
607  *result = c * res24_sin + s * res24_cos;
608  *abserr = fabs(c * est_sin) + fabs(s * est_cos);
609  } else {
610  *result = c * res24_cos - s * res24_sin;
611  *abserr = fabs(c * est_cos) + fabs(s * est_sin);
612  }
613 
614  *resabs = result_abs * half_length;
615  *resasc = GSL_DBL_MAX;
616 
617  return;
618  }
619  }
620 
621  /// Add the oscillating part to the integrand
622  virtual double transform(double t, func_t &func) {
623  return func(t)*sin(this->omega*t);
624  }
625 
626 #endif
627 
628  /// Return string denoting type ("inte_qawo_gsl_sin")
629  const char *type() { return "inte_qawo_gsl_sin"; }
630 
631  };
632 
633  /** \brief Adaptive integration a function with finite limits of
634  integration (GSL)
635 
636  The integral
637  \f[
638  \int_a^b f(x) \cos (\omega x)~dx
639  \f]
640  is computed for some frequency parameter \f$ \omega \f$ .
641 
642  This class is exactly analogous to \ref inte_qawo_gsl_sin .
643  See that class documentation for more details.
644  */
645  template<class func_t> class inte_qawo_gsl_cos :
646  public inte_qawo_gsl_sin<func_t> {
647 
648  public:
649 
651  }
652 
653  virtual ~inte_qawo_gsl_cos() {}
654 
655  /** \brief Integrate function \c func from \c a to \c b and place
656  the result in \c res and the error in \c err
657  */
658  virtual int integ_err(func_t &func, double a, double b,
659  double &res, double &err) {
660 
661  this->otable=gsl_integration_qawo_table_alloc
662  (this->omega,b-a,GSL_INTEG_COSINE,this->n_levels);
663 
664  int status=this->qawo(func,a,this->tol_abs,this->tol_rel,
665  this->w,this->otable,&res,&err);
666 
667  gsl_integration_qawo_table_free(this->otable);
668 
669  return status;
670  }
671 
672 #ifndef DOXYGEN_INTERNAL
673 
674  protected:
675 
676  /// Add the oscillating part to the integrand
677  virtual double transform(double t, func_t &func) {
678  return func(t)*cos(this->omega*t);
679  }
680 
681 #endif
682 
683  /// Return string denoting type ("inte_qawo_gsl_cos")
684  const char *type() { return "inte_qawo_gsl_cos"; }
685 
686  };
687 
688 #ifndef DOXYGEN_NO_O2NS
689 }
690 #endif
691 
692 #endif
const char * type()
Return string denoting type ("inte_qawo_gsl_cos")
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_qawo_gsl.h:79
Integration workspace for the GSL integrators.
int qawo(func_t &func, const double a, const double epsabs, const double epsrel, inte_workspace_gsl *loc_w, gsl_integration_qawo_table *wf, double *result, double *abserr)
The full GSL integration routine called by integ_err()
void initialise_table(struct extrapolation_table *table)
Initialize the table.
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
double * alist
Left endpoints of subintervals.
void inte_cheb_series(func2_t &f, double a, double b, double *cheb12, double *cheb24)
Compute Chebyshev series expansion using a FFT method.
double * blist
Right endpoints of subintervals.
Data table table class.
Definition: table.h:49
double omega
The user-specified frequency (default 1.0)
Definition: inte_qawo_gsl.h:71
sanity check failed - shouldn&#39;t happen
Definition: err_hnd.h:65
size_t * level
Numbers of subdivisions made.
void append_table(struct extrapolation_table *table, double y)
Append a result to the table.
apparent singularity detected
Definition: err_hnd.h:93
exceeded max number of iterations
Definition: err_hnd.h:73
int large_interval(inte_workspace_gsl *workspace)
Determine if an interval is large.
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
table table limit exceeded
Definition: err_hnd.h:103
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. ...
size_t n_levels
The number of bisection levels (default 10)
Definition: inte_qawo_gsl.h:74
size_t i
Index of current subinterval.
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.
generic failure
Definition: err_hnd.h:61
Chebyshev integration base class (GSL)
Definition: inte_qawc_gsl.h:46
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
virtual void gauss_kronrod(func_t &func, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
Integration wrapper for internal transformed function type.
int test_positivity(double result, double resabs)
Test if the integrand satisfies .
int initialise(double a, double b)
Initialize the workspace for an integration with limits a and b.
int set_initial_result(double result, double error)
Update the workspace with the result and error from the first integration.
const char * type()
Return string denoting type ("inte_qawo_gsl_sin")
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
virtual double transform(double t, func_t &func)
Add the oscillating part to the integrand.
size_t nrmax
Counter for extrapolation routine.
void qelg(struct extrapolation_table *table, double *result, double *abserr)
Determines the limit of a given sequence of approximations.
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 double transform(double t, func_t &func)
Add the oscillating part to the integrand.
void reset_nrmax(inte_workspace_gsl *workspace)
Reset workspace to work on the interval with the largest error.
int increase_nrmax(inte_workspace_gsl *workspace)
Increase workspace.
A structure for extrapolation for inte_qags_gsl.
size_t n
Index of new element in the first column.
Adaptive integration a function with finite limits of integration (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.
user specified an invalid tolerance
Definition: err_hnd.h:77
Success.
Definition: err_hnd.h:47
Adaptive integration for oscillatory integrals (GSL)
Definition: inte_qawo_gsl.h:58
void qc25f(func_t &func, double a, double b, gsl_integration_qawo_table *wf, size_t level, double *result, double *abserr, double *resabs, double *resasc)
25-point quadrature for oscillating functions
gsl_integration_qawo_table * otable
The integration workspace.
Definition: inte_qawo_gsl.h:98
int verbose
Verbosity.
Definition: inte.h:60
failed because of roundoff error
Definition: err_hnd.h:87
integral or series is divergent
Definition: err_hnd.h:95

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