diff options
Diffstat (limited to 'gsl-1.9/integration/qng.c')
| -rw-r--r-- | gsl-1.9/integration/qng.c | 190 |
1 files changed, 190 insertions, 0 deletions
diff --git a/gsl-1.9/integration/qng.c b/gsl-1.9/integration/qng.c new file mode 100644 index 0000000..ecf8e62 --- /dev/null +++ b/gsl-1.9/integration/qng.c @@ -0,0 +1,190 @@ +/* integration/qng.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <math.h> +#include <float.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_integration.h> + +#include "err.c" +#include "qng.h" + +int +gsl_integration_qng (const gsl_function *f, + double a, double b, + double epsabs, double epsrel, + double * result, double * abserr, size_t * neval) +{ + double fv1[5], fv2[5], fv3[5], fv4[5]; + double savfun[21]; /* array of function values which have been computed */ + double res10, res21, res43, res87; /* 10, 21, 43 and 87 point results */ + double result_kronrod, err ; + double resabs; /* approximation to the integral of abs(f) */ + double resasc; /* approximation to the integral of abs(f-i/(b-a)) */ + + const double half_length = 0.5 * (b - a); + const double abs_half_length = fabs (half_length); + const double center = 0.5 * (b + a); + const double f_center = GSL_FN_EVAL(f, center); + + int k ; + + if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28)) + { + * result = 0; + * abserr = 0; + * neval = 0; + GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel", + GSL_EBADTOL); + }; + + /* Compute the integral using the 10- and 21-point formula. */ + + res10 = 0; + res21 = w21b[5] * f_center; + resabs = w21b[5] * fabs (f_center); + + for (k = 0; k < 5; k++) + { + const double abscissa = half_length * x1[k]; + const double fval1 = GSL_FN_EVAL(f, center + abscissa); + const double fval2 = GSL_FN_EVAL(f, center - abscissa); + const double fval = fval1 + fval2; + res10 += w10[k] * fval; + res21 += w21a[k] * fval; + resabs += w21a[k] * (fabs (fval1) + fabs (fval2)); + savfun[k] = fval; + fv1[k] = fval1; + fv2[k] = fval2; + } + + for (k = 0; k < 5; k++) + { + const double abscissa = half_length * x2[k]; + const double fval1 = GSL_FN_EVAL(f, center + abscissa); + const double fval2 = GSL_FN_EVAL(f, center - abscissa); + const double fval = fval1 + fval2; + res21 += w21b[k] * fval; + resabs += w21b[k] * (fabs (fval1) + fabs (fval2)); + savfun[k + 5] = fval; + fv3[k] = fval1; + fv4[k] = fval2; + } + + resabs *= abs_half_length ; + + { + const double mean = 0.5 * res21; + + resasc = w21b[5] * fabs (f_center - mean); + + for (k = 0; k < 5; k++) + { + resasc += + (w21a[k] * (fabs (fv1[k] - mean) + fabs (fv2[k] - mean)) + + w21b[k] * (fabs (fv3[k] - mean) + fabs (fv4[k] - mean))); + } + resasc *= abs_half_length ; + } + + result_kronrod = res21 * half_length; + + err = rescale_error ((res21 - res10) * half_length, resabs, resasc) ; + + /* test for convergence. */ + + if (err < epsabs || err < epsrel * fabs (result_kronrod)) + { + * result = result_kronrod ; + * abserr = err ; + * neval = 21; + return GSL_SUCCESS; + } + + /* compute the integral using the 43-point formula. */ + + res43 = w43b[11] * f_center; + + for (k = 0; k < 10; k++) + { + res43 += savfun[k] * w43a[k]; + } + + for (k = 0; k < 11; k++) + { + const double abscissa = half_length * x3[k]; + const double fval = (GSL_FN_EVAL(f, center + abscissa) + + GSL_FN_EVAL(f, center - abscissa)); + res43 += fval * w43b[k]; + savfun[k + 10] = fval; + } + + /* test for convergence */ + + result_kronrod = res43 * half_length; + err = rescale_error ((res43 - res21) * half_length, resabs, resasc); + + if (err < epsabs || err < epsrel * fabs (result_kronrod)) + { + * result = result_kronrod ; + * abserr = err ; + * neval = 43; + return GSL_SUCCESS; + } + + /* compute the integral using the 87-point formula. */ + + res87 = w87b[22] * f_center; + + for (k = 0; k < 21; k++) + { + res87 += savfun[k] * w87a[k]; + } + + for (k = 0; k < 22; k++) + { + const double abscissa = half_length * x4[k]; + res87 += w87b[k] * (GSL_FN_EVAL(f, center + abscissa) + + GSL_FN_EVAL(f, center - abscissa)); + } + + /* test for convergence */ + + result_kronrod = res87 * half_length ; + + err = rescale_error ((res87 - res43) * half_length, resabs, resasc); + + if (err < epsabs || err < epsrel * fabs (result_kronrod)) + { + * result = result_kronrod ; + * abserr = err ; + * neval = 87; + return GSL_SUCCESS; + } + + /* failed to converge */ + + * result = result_kronrod ; + * abserr = err ; + * neval = 87; + + GSL_ERROR("failed to reach tolerance with highest-order rule", GSL_ETOL) ; +} |