23 #ifndef O2SCL_GSL_INTE_QAWF_H 24 #define O2SCL_GSL_INTE_QAWF_H 31 #include <o2scl/inte.h> 32 #include <o2scl/inte_qawo_gsl.h> 33 #include <o2scl/inte_qagiu_gsl.h> 35 #ifndef DOXYGEN_NO_O2NS 109 double &res,
double &err) {
111 this->
otable=gsl_integration_qawo_table_alloc
118 gsl_integration_qawo_table_free(this->
otable);
125 #ifndef DOXYGEN_INTERNAL 134 int qawf(func_t &func,
const double a,
135 const double epsabs,
double *result,
double *abserr) {
138 double res_ext, err_ext;
139 double correc, total_error = 0.0, truncation_error;
142 size_t iteration = 0;
149 const double p = 0.9;
151 double initial_eps, eps;
161 size_t limit=this->
w->
limit;
173 std::string estr=
"The absolute tolerance must be positive ";
174 estr+=
"in inte_qawf_gsl::qawf().";
178 if (this->
omega == 0.0) {
179 if (this->
otable->sine == GSL_INTEG_SINE) {
191 int status=iu.
integ_err(func,a,0.0,*result,*abserr);
197 if (epsabs > GSL_DBL_MIN / (1 - p)) {
198 eps = epsabs * (1 - p);
209 err_ext = GSL_DBL_MAX;
212 cycle = (2 * floor (fabs (this->
omega)) + 1) *
213 M_PI / fabs (this->
omega);
215 gsl_integration_qawo_table_set_length (this->
otable, cycle);
219 for (iteration = 0; iteration < limit; iteration++) {
220 double area1, error1, reseps, erreps;
222 double a1 = a + iteration * cycle;
223 double b1 = a1 + cycle;
225 double epsabs1 = eps * factor;
227 int status=this->
qawo(func,a1,epsabs1,0.0,cyclew,this->
otable,
235 errsum = errsum + error1;
239 truncation_error = 50 * fabs (area1);
241 total_error = errsum + truncation_error;
243 if (total_error < epsabs && iteration > 4) {
247 if (error1 > correc) {
252 eps = GSL_MAX_DBL (initial_eps, correc * (1.0 - p));
255 if (status && total_error < 10 * correc && iteration > 3) {
265 this->
qelg (&table, &reseps, &erreps);
269 if (ktmin >= 15 && err_ext < 0.001 * total_error) {
273 if (erreps < err_ext) {
278 if (err_ext + 10 * correc <= epsabs)
280 if (err_ext <= epsabs && 10 * correc >= epsabs)
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;
293 std::cout <<
"Press a key and type enter to continue. " ;
300 if (iteration == limit) error_type = 1;
302 if (err_ext == GSL_DBL_MAX)
305 err_ext = err_ext + 10 * correc;
310 if (error_type == 0) {
314 if (res_ext != 0.0 && area != 0.0) {
315 if (err_ext / fabs (res_ext) > errsum / fabs (area))
317 }
else if (err_ext > errsum) {
319 }
else if (area == 0.0) {
323 if (error_type == 4) {
324 err_ext = err_ext + truncation_error;
332 *abserr = total_error;
339 if (error_type == 0) {
341 }
else if (error_type == 1) {
342 std::string estr=
"Number of iterations was insufficient ";
343 estr+=
" in inte_qawf_gsl::qawf().";
345 }
else if (error_type == 2) {
346 std::string estr=
"Roundoff error prevents tolerance ";
347 estr+=
"from being achieved in inte_qawf_gsl::qawf().";
349 }
else if (error_type == 3) {
350 std::string estr=
"Bad integrand behavior ";
351 estr+=
" in inte_qawf_gsl::qawf().";
353 }
else if (error_type == 4) {
354 std::string estr=
"Roundoff error detected in extrapolation table ";
355 estr+=
"in inte_qawf_gsl::qawf().";
357 }
else if (error_type == 5) {
358 std::string estr=
"Integral is divergent or slowly convergent ";
359 estr+=
"in inte_qawf_gsl::qawf().";
362 std::string estr=
"Could not integrate function in inte_qawf_gsl";
363 estr+=
"::qawf() (it may have returned a non-finite result).";
374 return func(t)*sin(this->
omega*t);
380 const char *
type() {
return "inte_qawf_gsl_sin"; }
410 double &res,
double &err) {
412 this->
otable=gsl_integration_qawo_table_alloc
417 int status=this->
qawf(func,a,this->
tol_abs,&res,&err);
419 gsl_integration_qawo_table_free(this->
otable);
427 #ifndef DOXYGEN_INTERNAL 433 return func(t)*cos(this->
omega*t);
439 const char *
type() {
return "inte_qawf_gsl_cos"; }
443 #ifndef DOXYGEN_NO_O2NS double tol_abs
The maximum absolute uncertainty in the value of the integral (default )
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...
double omega
The user-specified frequency (default 1.0)
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
exceeded max number of iterations
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)
inte_workspace_gsl * cyclew
The integration workspace.
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.
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)
inte_workspace_gsl * w
The integration workspace.
const char * type()
Return string denoting type ("inte_qawf_gsl_cos")
size_t limit
Maximum number of subintervals allocated.
user specified an invalid tolerance
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)
gsl_integration_qawo_table * otable
The integration workspace.
virtual double transform(double t, func_t &func)
Add the oscillating part to the integrand.
failed because of roundoff error
integral or series is divergent
int qawf(func_t &func, const double a, const double epsabs, double *result, double *abserr)
The full GSL integration routine called by integ_err()