24 #ifndef O2SCL_GSL_INTE_QAWC_H 25 #define O2SCL_GSL_INTE_QAWC_H 31 #include <o2scl/err_hnd.h> 32 #include <o2scl/inte.h> 33 #include <o2scl/inte_singular_gsl.h> 35 #ifndef DOXYGEN_NO_O2NS 55 double a0 = log (fabs ((1.0 - cc) / (1.0 + cc)));
56 double a1 = 2 + a0 * cc;
61 for (k = 2; k < 25; k++) {
65 a2 = 2.0 * cc * a1 - a0;
67 const double km1 = k - 1.0;
68 a2 = 2.0 * cc * a1 - a0 - 4.0 / (km1 * km1 - 1.0);
114 template<
class func2_t>
116 double *cheb12,
double *cheb24) {
118 double fval[25], v[12];
123 const double x[11] = { 0.9914448613738104,
133 0.1305261922200516 };
135 const double center = 0.5 * (b + a);
136 const double half_length = 0.5 * (b - a);
146 for (i = 1; i < 12; i++) {
147 const size_t j = 24 - i;
148 const double u = half_length * x[i-1];
156 for (i = 0; i < 12; i++) {
157 const size_t j = 24 - i;
158 v[i] = fval[i] - fval[j];
159 fval[i] = fval[i] + fval[j];
163 const double alam1 = v[0] - v[8];
164 const double alam2 = x[5] * (v[2] - v[6] - v[10]);
166 cheb12[3] = alam1 + alam2;
167 cheb12[9] = alam1 - alam2;
171 const double alam1 = v[1] - v[7] - v[9];
172 const double alam2 = v[3] - v[5] - v[11];
174 const double alam = x[2] * alam1 + x[8] * alam2;
176 cheb24[3] = cheb12[3] + alam;
177 cheb24[21] = cheb12[3] - alam;
181 const double alam = x[8] * alam1 - x[2] * alam2;
182 cheb24[9] = cheb12[9] + alam;
183 cheb24[15] = cheb12[9] - alam;
188 const double part1 = x[3] * v[4];
189 const double part2 = x[7] * v[8];
190 const double part3 = x[5] * v[6];
193 const double alam1 = v[0] + part1 + part2;
194 const double alam2 = x[1] * v[2] + part3 + x[9] * v[10];
196 cheb12[1] = alam1 + alam2;
197 cheb12[11] = alam1 - alam2;
201 const double alam1 = v[0] - part1 + part2;
202 const double alam2 = x[9] * v[2] - part3 + x[1] * v[10];
203 cheb12[5] = alam1 + alam2;
204 cheb12[7] = alam1 - alam2;
209 const double alam = (x[0] * v[1] + x[2] * v[3] + x[4] * v[5]
210 + x[6] * v[7] + x[8] * v[9] + x[10] * v[11]);
211 cheb24[1] = cheb12[1] + alam;
212 cheb24[23] = cheb12[1] - alam;
216 const double alam = (x[10] * v[1] - x[8] * v[3] + x[6] * v[5]
217 - x[4] * v[7] + x[2] * v[9] - x[0] * v[11]);
218 cheb24[11] = cheb12[11] + alam;
219 cheb24[13] = cheb12[11] - alam;
223 const double alam = (x[4] * v[1] - x[8] * v[3] - x[0] * v[5]
224 - x[10] * v[7] + x[2] * v[9] + x[6] * v[11]);
225 cheb24[5] = cheb12[5] + alam;
226 cheb24[19] = cheb12[5] - alam;
230 const double alam = (x[6] * v[1] - x[2] * v[3] - x[10] * v[5]
231 + x[0] * v[7] - x[8] * v[9] - x[4] * v[11]);
232 cheb24[7] = cheb12[7] + alam;
233 cheb24[17] = cheb12[7] - alam;
236 for (i = 0; i < 6; i++) {
237 const size_t j = 12 - i;
238 v[i] = fval[i] - fval[j];
239 fval[i] = fval[i] + fval[j];
243 const double alam1 = v[0] + x[7] * v[4];
244 const double alam2 = x[3] * v[2];
246 cheb12[2] = alam1 + alam2;
247 cheb12[10] = alam1 - alam2;
250 cheb12[6] = v[0] - v[4];
253 const double alam = x[1] * v[1] + x[5] * v[3] + x[9] * v[5];
254 cheb24[2] = cheb12[2] + alam;
255 cheb24[22] = cheb12[2] - alam;
259 const double alam = x[5] * (v[1] - v[3] - v[5]);
260 cheb24[6] = cheb12[6] + alam;
261 cheb24[18] = cheb12[6] - alam;
265 const double alam = x[9] * v[1] - x[5] * v[3] + x[1] * v[5];
266 cheb24[10] = cheb12[10] + alam;
267 cheb24[14] = cheb12[10] - alam;
270 for (i = 0; i < 3; i++) {
271 const size_t j = 6 - i;
272 v[i] = fval[i] - fval[j];
273 fval[i] = fval[i] + fval[j];
276 cheb12[4] = v[0] + x[7] * v[2];
277 cheb12[8] = fval[0] - x[7] * fval[2];
280 const double alam = x[3] * v[1];
281 cheb24[4] = cheb12[4] + alam;
282 cheb24[20] = cheb12[4] - alam;
286 const double alam = x[7] * fval[1] - fval[3];
287 cheb24[8] = cheb12[8] + alam;
288 cheb24[16] = cheb12[8] - alam;
291 cheb12[0] = fval[0] + fval[2];
294 const double alam = fval[1] + fval[3];
295 cheb24[0] = cheb12[0] + alam;
296 cheb24[24] = cheb12[0] - alam;
299 cheb12[12] = v[0] - v[2];
300 cheb24[12] = cheb12[12];
302 for (i = 1; i < 12; i++) {
303 cheb12[i] *= 1.0 / 6.0;
306 cheb12[0] *= 1.0 / 12.0;
307 cheb12[12] *= 1.0 / 12.0;
309 for (i = 1; i < 24; i++) {
310 cheb24[i] *= 1.0 / 12.0;
313 cheb24[0] *= 1.0 / 24.0;
314 cheb24[24] *= 1.0 / 24.0;
368 double &res,
double &err) {
370 return this->qawc(func,a,b,s,this->
tol_abs,this->
tol_rel,&res,&err);
373 #ifndef DOXYGEN_INTERNAL 379 int qawc(func_t &func,
const double a,
const double b,
const double c,
380 const double epsabs,
const double epsrel,
381 double *result,
double *abserr) {
384 double result0, abserr0;
386 size_t iteration = 0;
387 int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;
390 double lower, higher;
397 size_t limit=this->
w->
limit;
410 double dbl_eps=std::numeric_limits<double>::epsilon();
412 if (epsabs <= 0 && (epsrel < 50 * dbl_eps || epsrel < 0.5e-28)) {
414 std::string estr=
"Tolerance cannot be achieved with given ";
415 estr+=
"value of tol_abs, "+
dtos(epsabs)+
", and tol_rel, "+
416 dtos(epsrel)+
", in inte_qawc_gsl::qawc().";
420 if (c == a || c == b) {
422 std::string estr=
"Cannot integrate with singularity on endpoint ";
423 estr+=
"in inte_qawc_gsl::qawc().";
429 this->qc25c(func,lower,higher,c,&result0,&abserr0, &err_reliable);
436 tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
438 if (abserr0 < tolerance && abserr0 < 0.01 * fabs(result0)) {
441 *result = sign * result0;
445 }
else if (limit == 1) {
447 *result = sign * result0;
451 std::string estr=
"A maximum of 1 iteration was insufficient ";
452 estr+=
"in inte_qawc_gsl::qawc().";
463 double a1, b1, a2, b2;
464 double a_i, b_i, r_i, e_i;
465 double area1 = 0, area2 = 0, area12 = 0;
466 double error1 = 0, error2 = 0, error12 = 0;
467 int err_reliable1, err_reliable2;
471 this->
w->
retrieve (&a_i, &b_i, &r_i, &e_i);
474 b1 = 0.5 * (a_i + b_i);
478 if (c > a1 && c <= b1) {
479 b1 = 0.5 * (c + b2) ;
481 }
else if (c > b1 && c < b2) {
482 b1 = 0.5 * (a1 + c) ;
486 qc25c (func, a1, b1, c, &area1, &error1, &err_reliable1);
487 qc25c (func, a2, b2, c, &area2, &error2, &err_reliable2);
489 area12 = area1 + area2;
490 error12 = error1 + error2;
492 errsum += (error12 - e_i);
493 area += area12 - r_i;
495 if (err_reliable1 && err_reliable2) {
496 double delta = r_i - area12;
498 if (fabs (delta) <= 1.0e-5 * fabs (area12) &&
499 error12 >= 0.99 * e_i) {
502 if (iteration >= 10 && error12 > e_i) {
507 tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
509 if (errsum > tolerance) {
510 if (roundoff_type1 >= 6 || roundoff_type2 >= 20) {
522 this->
w->
update (a1, b1, area1, error1, a2, b2, area2, error2);
524 this->
w->
retrieve (&a_i, &b_i, &r_i, &e_i);
527 std::cout <<
"inte_qawc_gsl Iter: " << iteration;
528 std::cout.setf(std::ios::showpos);
529 std::cout <<
" Res: " << area;
530 std::cout.unsetf(std::ios::showpos);
531 std::cout <<
" Err: " << errsum
532 <<
" Tol: " << tolerance << std::endl;
535 std::cout <<
"Press a key and type enter to continue. " ;
542 }
while (iteration < limit && !error_type && errsum > tolerance);
548 if (errsum <= tolerance) {
550 }
else if (error_type == 2) {
551 std::string estr=
"Roundoff error prevents tolerance ";
552 estr+=
"from being achieved in inte_qawc_gsl::qawc().";
554 }
else if (error_type == 3) {
555 std::string estr=
"Bad integrand behavior ";
556 estr+=
" in inte_qawc_gsl::qawc().";
558 }
else if (iteration == limit) {
559 std::string estr=
"Maximum number of subdivisions ("+
itos(iteration);
560 estr+=
") reached in inte_qawc_gsl::qawc().";
563 std::string estr=
"Could not integrate function in inte_qawc_gsl::";
564 estr+=
"qawc() (it may have returned a non-finite result).";
575 void qc25c(func_t &func,
double a,
double b,
double c,
576 double *result,
double *abserr,
int *err_reliable) {
578 double cc = (2 * c - b - a) / (b - a);
580 if (fabs (cc) > 1.1) {
581 double resabs, resasc;
585 if (*abserr == resasc) {
595 double cheb12[13], cheb24[25], moment[25];
596 double res12 = 0, res24 = 0;
601 for (i = 0; i < 13; i++) {
602 res12 += cheb12[i] * moment[i];
605 for (i = 0; i < 25; i++) {
606 res24 += cheb24[i] * moment[i];
610 *abserr = fabs(res24 - res12) ;
627 const char *
type() {
return "inte_qawc_gsl"; }
631 #ifndef DOXYGEN_NO_O2NS Adaptive Cauchy principal value integration (GSL)
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.
virtual double transform(double t, func_t &func)
Add the singularity to the function.
void inte_cheb_series(func2_t &f, double a, double b, double *cheb12, double *cheb24)
Compute Chebyshev series expansion using a FFT method.
invalid argument supplied by user
apparent singularity detected
exceeded max number of iterations
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...
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. ...
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 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.
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
void compute_moments(double cc, double *moment)
Compute the Chebyshev moments.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
int qawc(func_t &func, const double a, const double b, const double c, const double epsabs, const double epsrel, double *result, double *abserr)
The full GSL integration routine called by integ_err()
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 )
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.
std::string itos(int x)
Convert an integer to a string.
void qc25c(func_t &func, double a, double b, double c, double *result, double *abserr, int *err_reliable)
25-point quadrature for Cauchy principal values
user specified an invalid tolerance
const char * type()
Return string denoting type ("inte_qawc_gsl")
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.
failed because of roundoff error