inte_qawf_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_QAWF_H
24 #define O2SCL_GSL_INTE_QAWF_H
25 
26 /** \file inte_qawf_gsl.h
27  \brief File defining \ref o2scl::inte_qawf_gsl_sin and
28  \ref o2scl::inte_qawf_gsl_cos
29 */
30 
31 #include <o2scl/inte.h>
32 #include <o2scl/inte_qawo_gsl.h>
33 #include <o2scl/inte_qagiu_gsl.h>
34 
35 #ifndef DOXYGEN_NO_O2NS
36 namespace o2scl {
37 #endif
38 
39  /** \brief Adaptive integration for oscillatory integrals (GSL)
40 
41  The Fourier integral
42  \f[
43  \int_a^{\infty} f(x) \sin(\omega x)~dx
44  \f]
45  is computed for some frequency parameter \f$ \omega \f$,
46  stored in \ref inte_qawo_gsl_sin::omega .
47 
48  The integral is computed using the same method as \ref
49  inte_qawo_gsl_sin and \ref inte_qawo_gsl_cos over each of the
50  subintervals,
51  \f{eqnarray*}{
52  C_1 &=& [a, a+c] \\
53  C_2 &=& [a+c, a+2c] \\
54  &\vdots & \\
55  C_k &=& [a +(k-1)c,\, a+kc],
56  \f}
57  where \f$ c = (2\mathrm{floor}(|\omega|)+1)\pi/|\omega|\f$. This
58  width is chosen to cover an odd number of periods so that the
59  contributions from the intervals alternate in sign and are
60  monotonically decreasing when \f$ f \f$ is positive and
61  monotonically decreasing. The sum of this sequence of
62  contributions is accelerated using the \f$ \varepsilon \f$
63  algorithm.
64 
65  The algorithm uses zero for the relative tolerance \ref
66  inte::tol_rel and attempts to compute the integral to an
67  overall absolute tolerance set by \ref inte::tol_abs. The
68  following strategy is used: on each interval \f$ C_k\f$, the
69  algorithm tries to achieve the tolerance
70  \f[
71  \mathrm{TOL}_k = u_k\cdot \epsilon_{\mathrm{abs}}
72  \f]
73  where \f$ u_k = (1-p)p^{k-1} \f$ and \f$ p = 0.9\f$. The sum of
74  the geometric series of contributions from each interval gives
75  an overall tolerance of \f$ \epsilon_{\mathrm{abs}}\f$. If the
76  integration of a subinterval leads to difficulties then the accu
77  racy requirement for subsequent intervals is relaxed,
78  \f[
79  \mathrm{TOL}_k =
80  u_k\cdot \max\{\epsilon_{\mathrm{abs}}, E_1, \ldots, E_{k-1} \}
81  \f]
82  where \f$ E_k\f$ is the estimated error on the interval \f$ C_k\f$.
83 
84  See \ref gslinte_subsect in the User's guide for general
85  information about the GSL integration classes.
86 
87  When verbose output is enabled, this class outputs information
88  from both the subintegrations performed by \ref
89  inte_qawo_gsl_sin and the overall integration progress in this
90  class.
91 
92  \todo More documentation and examples for the
93  qawf, qawo and qawc integrators.
94  */
95  template<class func_t> class inte_qawf_gsl_sin :
96  public inte_qawo_gsl_sin<func_t> {
97 
98  public:
99 
101  }
102 
103  virtual ~inte_qawf_gsl_sin() {}
104 
105  /** \brief Integrate function \c func from \c a to \c b and place
106  the result in \c res and the error in \c err
107  */
108  virtual int integ_err(func_t &func, double a, double b,
109  double &res, double &err) {
110 
111  this->otable=gsl_integration_qawo_table_alloc
112  (this->omega,1.0,GSL_INTEG_SINE,this->n_levels);
113  this->cyclew=new inte_workspace_gsl;
114  this->cyclew->allocate(this->w->limit);
115 
116  int status=qawf(func,a,this->tol_abs,&res,&err);
117 
118  gsl_integration_qawo_table_free(this->otable);
119  this->cyclew->free();
120  delete this->cyclew;
121 
122  return status;
123  }
124 
125 #ifndef DOXYGEN_INTERNAL
126 
127  protected:
128 
129  /// The integration workspace
131 
132  /** \brief The full GSL integration routine called by integ_err()
133  */
134  int qawf(func_t &func, const double a,
135  const double epsabs, double *result, double *abserr) {
136 
137  double area, errsum;
138  double res_ext, err_ext;
139  double correc, total_error = 0.0, truncation_error;
140 
141  size_t ktmin = 0;
142  size_t iteration = 0;
143 
145 
146  double cycle;
147  //double omega = this->otable->omega;
148 
149  const double p = 0.9;
150  double factor = 1;
151  double initial_eps, eps;
152  int error_type = 0;
153 
154  /* Initialize results */
155 
156  this->w->initialise(a,a);
157 
158  *result = 0;
159  *abserr = 0;
160 
161  size_t limit=this->w->limit;
162  /*
163  if (limit > this->w->limit) {
164  std::string estr="Iteration limit exceeds workspace ";
165  estr+="in inte_qawf_gsl::qawf().";
166  O2SCL_ERR(estr.c_str(),exc_einval);
167  }
168  */
169 
170  /* Test on accuracy */
171 
172  if (epsabs <= 0) {
173  std::string estr="The absolute tolerance must be positive ";
174  estr+="in inte_qawf_gsl::qawf().";
175  O2SCL_ERR(estr.c_str(),exc_ebadtol);
176  }
177 
178  if (this->omega == 0.0) {
179  if (this->otable->sine == GSL_INTEG_SINE) {
180  /* The function sin(w x) f(x) is always zero for w = 0 */
181 
182  *result = 0;
183  *abserr = 0;
184 
185  return success;
186  } else {
187  /* The function cos(w x) f(x) is always f(x) for w = 0 */
188 
190 
191  int status=iu.integ_err(func,a,0.0,*result,*abserr);
192 
193  return status;
194  }
195  }
196 
197  if (epsabs > GSL_DBL_MIN / (1 - p)) {
198  eps = epsabs * (1 - p);
199  } else {
200  eps = epsabs;
201  }
202 
203  initial_eps = eps;
204 
205  area = 0;
206  errsum = 0;
207 
208  res_ext = 0;
209  err_ext = GSL_DBL_MAX;
210  correc = 0;
211 
212  cycle = (2 * floor (fabs (this->omega)) + 1) *
213  M_PI / fabs (this->omega);
214 
215  gsl_integration_qawo_table_set_length (this->otable, cycle);
216 
217  this->initialise_table (&table);
218 
219  for (iteration = 0; iteration < limit; iteration++) {
220  double area1, error1, reseps, erreps;
221 
222  double a1 = a + iteration * cycle;
223  double b1 = a1 + cycle;
224 
225  double epsabs1 = eps * factor;
226 
227  int status=this->qawo(func,a1,epsabs1,0.0,cyclew,this->otable,
228  &area1,&error1);
229 
230  this->w->append_interval (a1, b1, area1, error1);
231 
232  factor *= p;
233 
234  area = area + area1;
235  errsum = errsum + error1;
236 
237  /* estimate the truncation error as 50 times the final term */
238 
239  truncation_error = 50 * fabs (area1);
240 
241  total_error = errsum + truncation_error;
242 
243  if (total_error < epsabs && iteration > 4) {
244  goto compute_result;
245  }
246 
247  if (error1 > correc) {
248  correc = error1;
249  }
250 
251  if (status) {
252  eps = GSL_MAX_DBL (initial_eps, correc * (1.0 - p));
253  }
254 
255  if (status && total_error < 10 * correc && iteration > 3) {
256  goto compute_result;
257  }
258 
259  this->append_table (&table, area);
260 
261  if (table.n < 2) {
262  continue;
263  }
264 
265  this->qelg (&table, &reseps, &erreps);
266 
267  ktmin++;
268 
269  if (ktmin >= 15 && err_ext < 0.001 * total_error) {
270  error_type = 4;
271  }
272 
273  if (erreps < err_ext) {
274  ktmin = 0;
275  err_ext = erreps;
276  res_ext = reseps;
277 
278  if (err_ext + 10 * correc <= epsabs)
279  break;
280  if (err_ext <= epsabs && 10 * correc >= epsabs)
281  break;
282  }
283 
284  if (this->verbose>0) {
285  std::cout << "inte_qawf_gsl Iter: " << iteration;
286  std::cout.setf(std::ios::showpos);
287  std::cout << " Res: " << area;
288  std::cout.unsetf(std::ios::showpos);
289  std::cout << " Err: " << total_error
290  << " Tol: " << epsabs << std::endl;
291  if (this->verbose>1) {
292  char ch;
293  std::cout << "Press a key and type enter to continue. " ;
294  std::cin >> ch;
295  }
296  }
297 
298  }
299 
300  if (iteration == limit) error_type = 1;
301 
302  if (err_ext == GSL_DBL_MAX)
303  goto compute_result;
304 
305  err_ext = err_ext + 10 * correc;
306 
307  *result = res_ext;
308  *abserr = err_ext;
309 
310  if (error_type == 0) {
311  return success ;
312  }
313 
314  if (res_ext != 0.0 && area != 0.0) {
315  if (err_ext / fabs (res_ext) > errsum / fabs (area))
316  goto compute_result;
317  } else if (err_ext > errsum) {
318  goto compute_result;
319  } else if (area == 0.0) {
320  goto return_error;
321  }
322 
323  if (error_type == 4) {
324  err_ext = err_ext + truncation_error;
325  }
326 
327  goto return_error;
328 
329  compute_result:
330 
331  *result = area;
332  *abserr = total_error;
333 
334  return_error:
335 
336  if (error_type > 2)
337  error_type--;
338 
339  if (error_type == 0) {
340  return success;
341  } else if (error_type == 1) {
342  std::string estr="Number of iterations was insufficient ";
343  estr+=" in inte_qawf_gsl::qawf().";
344  O2SCL_ERR(estr.c_str(),exc_emaxiter);
345  } else if (error_type == 2) {
346  std::string estr="Roundoff error prevents tolerance ";
347  estr+="from being achieved in inte_qawf_gsl::qawf().";
348  O2SCL_ERR(estr.c_str(),exc_eround);
349  } else if (error_type == 3) {
350  std::string estr="Bad integrand behavior ";
351  estr+=" in inte_qawf_gsl::qawf().";
352  O2SCL_ERR(estr.c_str(),exc_esing);
353  } else if (error_type == 4) {
354  std::string estr="Roundoff error detected in extrapolation table ";
355  estr+="in inte_qawf_gsl::qawf().";
356  O2SCL_ERR(estr.c_str(),exc_eround);
357  } else if (error_type == 5) {
358  std::string estr="Integral is divergent or slowly convergent ";
359  estr+="in inte_qawf_gsl::qawf().";
360  O2SCL_ERR(estr.c_str(),exc_ediverge);
361  } else {
362  std::string estr="Could not integrate function in inte_qawf_gsl";
363  estr+="::qawf() (it may have returned a non-finite result).";
364  O2SCL_ERR(estr.c_str(),exc_efailed);
365  }
366  // No return statement needed since the above if statement
367  // always forces a return, but some compilers like having one
368  // anyway.
369  return o2scl::success;
370  }
371 
372  /// Add the oscillating part to the integrand
373  virtual double transform(double t, func_t &func) {
374  return func(t)*sin(this->omega*t);
375  }
376 
377 #endif
378 
379  /// Return string denoting type ("inte_qawf_gsl_sin")
380  const char *type() { return "inte_qawf_gsl_sin"; }
381 
382  };
383 
384  /** \brief Adaptive integration a function with finite limits of
385  integration (GSL)
386 
387  The Fourier integral
388  \f[
389  \int_a^{\infty} f(x) \cos(\omega x)~dx
390  \f]
391  is computed for some frequency parameter \f$ \omega \f$ .
392 
393  This class is exactly analogous to \ref inte_qawf_gsl_sin .
394  See that class documentation for more details.
395  */
396  template<class func_t> class inte_qawf_gsl_cos :
397  public inte_qawf_gsl_sin<func_t> {
398 
399  public:
400 
402  }
403 
404  virtual ~inte_qawf_gsl_cos() {}
405 
406  /** \brief Integrate function \c func from \c a to \c b and place
407  the result in \c res and the error in \c err
408  */
409  virtual int integ_err(func_t &func, double a, double b,
410  double &res, double &err) {
411 
412  this->otable=gsl_integration_qawo_table_alloc
413  (this->omega,b-a,GSL_INTEG_COSINE,this->n_levels);
414  this->cyclew=new inte_workspace_gsl;
415  this->cyclew->allocate(this->w->limit);
416 
417  int status=this->qawf(func,a,this->tol_abs,&res,&err);
418 
419  gsl_integration_qawo_table_free(this->otable);
420  this->cyclew->free();
421  delete this->cyclew;
422 
423  return status;
424 
425  }
426 
427 #ifndef DOXYGEN_INTERNAL
428 
429  protected:
430 
431  /// Add the oscillating part to the integrand
432  virtual double transform(double t, func_t &func) {
433  return func(t)*cos(this->omega*t);
434  }
435 
436 #endif
437 
438  /// Return string denoting type ("inte_qawf_gsl_cos")
439  const char *type() { return "inte_qawf_gsl_cos"; }
440 
441  };
442 
443 #ifndef DOXYGEN_NO_O2NS
444 }
445 #endif
446 
447 #endif
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
Data table table class.
Definition: table.h:49
double omega
The user-specified frequency (default 1.0)
Definition: inte_qawo_gsl.h:71
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 integration a function with finite limits of integration (GSL)
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
const char * type()
Return string denoting type ("inte_qawf_gsl_sin")
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.
size_t n_levels
The number of bisection levels (default 10)
Definition: inte_qawo_gsl.h:74
generic failure
Definition: err_hnd.h:61
inte_workspace_gsl * cyclew
The integration workspace.
double tol_abs
The maximum absolute uncertainty in the value of the integral (default )
Definition: inte.h:73
Integrate a function over the interval (GSL)
int initialise(double a, double b)
Initialize the workspace for an integration with limits a and b.
virtual int integ_err(func_t &func, double a, double b, double &res, double &err)
Integrate a function over the interval giving result res and error err.
int allocate(size_t sz)
Allocate a workspace.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
void append_interval(double a1, double b1, double area1, double error1)
Push a new interval to the workspace stack.
void qelg(struct extrapolation_table *table, double *result, double *abserr)
Determines the limit of a given sequence of approximations.
Adaptive integration for oscillatory integrals (GSL)
Definition: inte_qawf_gsl.h:95
inte_workspace_gsl * w
The integration workspace.
A structure for extrapolation for inte_qags_gsl.
const char * type()
Return string denoting type ("inte_qawf_gsl_cos")
size_t n
Index of new element in the first column.
size_t limit
Maximum number of subintervals allocated.
user specified an invalid tolerance
Definition: err_hnd.h:77
Success.
Definition: err_hnd.h:47
int free()
Free allocated workspace memory.
virtual double transform(double t, func_t &func)
Add the oscillating part to the integrand.
Adaptive integration for oscillatory integrals (GSL)
Definition: inte_qawo_gsl.h:58
gsl_integration_qawo_table * otable
The integration workspace.
Definition: inte_qawo_gsl.h:98
int verbose
Verbosity.
Definition: inte.h:60
virtual double transform(double t, func_t &func)
Add the oscillating part to the integrand.
failed because of roundoff error
Definition: err_hnd.h:87
integral or series is divergent
Definition: err_hnd.h:95
int qawf(func_t &func, const double a, const double epsabs, double *result, double *abserr)
The full GSL integration routine called by integ_err()

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