23 #ifndef O2SCL_GSL_INTE_QAWO_H 24 #define O2SCL_GSL_INTE_QAWO_H 31 #include <o2scl/inte.h> 32 #include <o2scl/inte_qawc_gsl.h> 34 #ifndef DOXYGEN_NO_O2NS 79 virtual int integ_err(func_t &func,
double a,
double b,
80 double &res,
double &err) {
82 otable=gsl_integration_qawo_table_alloc
83 (omega,b-a,GSL_INTEG_SINE,n_levels);
86 this->
w,this->
otable,&res,&err);
88 gsl_integration_qawo_table_free(
otable);
93 #ifndef DOXYGEN_INTERNAL 104 int qawo(func_t &func,
const double a,
const double epsabs,
106 gsl_integration_qawo_table *wf,
double *result,
double *abserr) {
109 double res_ext, err_ext;
110 double result0=0.0, abserr0=0.0, resabs0=0.0, resasc0=0.0;
114 double error_over_large_intervals = 0;
115 double reseps = 0, abseps = 0, correc = 0;
117 int roundoff_type1 = 0, roundoff_type2 = 0, roundoff_type3 = 0;
118 int error_type = 0, error_type2 = 0;
120 size_t iteration = 0;
122 int positive_integrand = 0;
125 int disallow_extrapolation = 0;
129 double b = a + wf->L ;
130 double abs_omega = fabs (wf->omega) ;
139 size_t limit=this->
w->
limit;
143 double dbl_eps=std::numeric_limits<double>::epsilon();
145 if (epsabs <= 0 && (epsrel < 50 * dbl_eps || epsrel < 0.5e-28)) {
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().";
155 this->
qc25f(func, a, b, wf, 0,
156 &result0, &abserr0, &resabs0, &resasc0);
160 tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
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;
171 std::cout <<
"Press a key and type enter to continue. " ;
176 if (abserr0 <= 100 * GSL_DBL_EPSILON * resabs0 &&
177 abserr0 > tolerance) {
182 std::string estr=
"Cannot reach tolerance because of roundoff error ";
183 estr+=
"on first attempt in inte_qawo_gsl_sin::qawo().";
185 }
else if ((abserr0 <= tolerance && abserr0 != resasc0) ||
192 }
else if (limit == 1) {
197 std::string estr=
"A maximum of 1 iteration was insufficient ";
198 estr+=
"in inte_qawo_gsl_sin::qawo().";
206 if (0.5 * abs_omega * fabs(b - a) <= 2) {
215 err_ext = GSL_DBL_MAX;
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;
234 loc_w->
retrieve (&a_i, &b_i, &r_i, &e_i);
236 current_level = loc_w->
level[loc_w->
i] + 1;
238 if (current_level >= wf->n) {
244 b1 = 0.5 * (a_i + b_i);
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);
255 area12 = area1 + area2;
256 error12 = error1 + error2;
266 errsum = errsum + error12 - e_i;
267 area = area + area12 - r_i;
269 tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
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;
280 std::cout <<
"Press a key and type enter to continue. " ;
285 if (resasc1 != error1 && resasc2 != error2) {
287 double delta = r_i - area12;
289 if (fabs (delta) <= 1.0e-5 * fabs (area12) &&
290 error12 >= 0.99 * e_i) {
299 if (iteration > 10 && error12 > e_i) {
306 if (roundoff_type1 + roundoff_type2 >= 10 || roundoff_type3 >= 20) {
310 if (roundoff_type2 >= 5) {
323 loc_w->
update (a1, b1, area1, error1, a2, b2, area2, error2);
325 if (errsum <= tolerance) {
333 if (iteration >= limit - 1) {
340 if (iteration == 2 && extall) {
341 error_over_large_intervals = errsum;
347 if (disallow_extrapolation) {
352 error_over_large_intervals += -last_e_i;
354 if (current_level < loc_w->maximum_level)
356 error_over_large_intervals += error12;
375 double width = loc_w->
blist[i] - loc_w->
alist[i];
377 if (0.25 * fabs(width) * abs_omega > 2) {
382 error_over_large_intervals = errsum;
389 if (!error_type2 && error_over_large_intervals > ertest) {
402 error_over_large_intervals = errsum;
406 this->
qelg (&table, &reseps, &abseps);
410 if (ktmin > 5 && err_ext < 0.001 * errsum) {
414 if (abseps < err_ext) {
418 correc = error_over_large_intervals;
419 ertest = GSL_MAX_DBL (epsabs, epsrel * fabs (reseps));
420 if (err_ext <= ertest) {
428 disallow_extrapolation = 1;
431 if (error_type == 5) {
439 error_over_large_intervals = errsum;
441 }
while (iteration < limit);
446 if (err_ext == GSL_DBL_MAX)
449 if (error_type || error_type2) {
458 if (result != 0 && area != 0) {
459 if (err_ext / fabs (res_ext) > errsum / fabs (area)) {
462 }
else if (err_ext > errsum) {
464 }
else if (area == 0.0) {
472 double max_area = GSL_MAX_DBL (fabs (res_ext), fabs (area));
474 if (!positive_integrand && max_area < 0.01 * resabs0) {
480 double ratio = res_ext / area;
482 if (ratio < 0.01 || ratio > 100 || errsum > fabs (area)) {
496 if (error_type > 2) {
502 if (error_type == 0) {
504 }
else if (error_type == 1) {
505 std::string estr=
"Number of iterations was insufficient ";
506 estr+=
" in inte_qawo_gsl_sin::qawo().";
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().";
512 }
else if (error_type == 3) {
513 std::string estr=
"Bad integrand behavior ";
514 estr+=
" in inte_qawo_gsl_sin::qawo().";
516 }
else if (error_type == 4) {
517 std::string estr=
"Roundoff error detected in extrapolation table ";
518 estr+=
"in inte_qawo_gsl_sin::qawo().";
520 }
else if (error_type == 5) {
521 std::string estr=
"Integral is divergent or slowly convergent ";
522 estr+=
"in inte_qawo_gsl_sin::qawo().";
524 }
else if (error_type == -1) {
525 std::string estr=
"Exceeded limit of trigonometric table ";
526 estr+=
"inte_qawo_gsl_sin::qawo()";
529 std::string estr=
"Could not integrate function in inte_qawo_gsl";
530 estr+=
"::qawo() (it may have returned a non-finite result).";
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,
546 const double center = 0.5 * (a + b);
547 const double half_length = 0.5 * (b - a);
549 const double par = omega * half_length;
551 if (fabs (par) < 2) {
560 double cheb12[13], cheb24[25];
561 double result_abs, res12_cos, res12_sin, res24_cos, res24_sin;
562 double est_cos, est_sin;
568 if (level >= wf->n) {
570 O2SCL_ERR(
"Table overflow in inte_qawo_gsl::qc25f().",
577 moment = wf->chebmo + 25 * level;
579 res12_cos = cheb12[12] * moment[12];
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];
588 res24_cos = cheb24[24] * moment[24];
591 result_abs = fabs(cheb24[24]) ;
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]);
600 est_cos = fabs(res24_cos - res12_cos);
601 est_sin = fabs(res24_sin - res12_sin);
603 c = half_length * cos(center * omega);
604 s = half_length * sin(center * omega);
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);
610 *result = c * res24_cos - s * res24_sin;
611 *abserr = fabs(c * est_cos) + fabs(s * est_sin);
614 *resabs = result_abs * half_length;
615 *resasc = GSL_DBL_MAX;
623 return func(t)*sin(this->omega*t);
629 const char *
type() {
return "inte_qawo_gsl_sin"; }
659 double &res,
double &err) {
661 this->
otable=gsl_integration_qawo_table_alloc
665 this->
w,this->
otable,&res,&err);
667 gsl_integration_qawo_table_free(this->
otable);
672 #ifndef DOXYGEN_INTERNAL 678 return func(t)*cos(this->
omega*t);
684 const char *
type() {
return "inte_qawo_gsl_cos"; }
688 #ifndef DOXYGEN_NO_O2NS 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.
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='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.
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.
double omega
The user-specified frequency (default 1.0)
sanity check failed - shouldn't happen
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
exceeded max number of iterations
int large_interval(inte_workspace_gsl *workspace)
Determine if an interval is large.
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...
table table limit exceeded
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)
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.
Chebyshev integration base class (GSL)
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 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.
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 )
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.
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
Adaptive integration for oscillatory integrals (GSL)
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.
failed because of roundoff error
integral or series is divergent