diff options
Diffstat (limited to 'gsl-1.9/integration/qawf.c')
| -rw-r--r-- | gsl-1.9/integration/qawf.c | 281 |
1 files changed, 281 insertions, 0 deletions
diff --git a/gsl-1.9/integration/qawf.c b/gsl-1.9/integration/qawf.c new file mode 100644 index 0000000..ba72f53 --- /dev/null +++ b/gsl-1.9/integration/qawf.c @@ -0,0 +1,281 @@ +/* integration/qawf.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 "initialise.c" +#include "append.c" +#include "qelg.c" + +int +gsl_integration_qawf (gsl_function * f, + const double a, + const double epsabs, + const size_t limit, + gsl_integration_workspace * workspace, + gsl_integration_workspace * cycle_workspace, + gsl_integration_qawo_table * wf, + double *result, double *abserr) +{ + double area, errsum; + double res_ext, err_ext; + double correc, total_error = 0.0, truncation_error; + + size_t ktmin = 0; + size_t iteration = 0; + + struct extrapolation_table table; + + double cycle; + double omega = wf->omega; + + const double p = 0.9; + double factor = 1; + double initial_eps, eps; + int error_type = 0; + + /* Initialize results */ + + initialise (workspace, a, a); + + *result = 0; + *abserr = 0; + + if (limit > workspace->limit) + { + GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ; + } + + /* Test on accuracy */ + + if (epsabs <= 0) + { + GSL_ERROR ("absolute tolerance epsabs must be positive", GSL_EBADTOL) ; + } + + if (omega == 0.0) + { + if (wf->sine == GSL_INTEG_SINE) + { + /* The function sin(w x) f(x) is always zero for w = 0 */ + + *result = 0; + *abserr = 0; + + return GSL_SUCCESS; + } + else + { + /* The function cos(w x) f(x) is always f(x) for w = 0 */ + + int status = gsl_integration_qagiu (f, a, epsabs, 0.0, + cycle_workspace->limit, + cycle_workspace, + result, abserr); + return status; + } + } + + if (epsabs > GSL_DBL_MIN / (1 - p)) + { + eps = epsabs * (1 - p); + } + else + { + eps = epsabs; + } + + initial_eps = eps; + + area = 0; + errsum = 0; + + res_ext = 0; + err_ext = GSL_DBL_MAX; + correc = 0; + + cycle = (2 * floor (fabs (omega)) + 1) * M_PI / fabs (omega); + + gsl_integration_qawo_table_set_length (wf, cycle); + + initialise_table (&table); + + for (iteration = 0; iteration < limit; iteration++) + { + double area1, error1, reseps, erreps; + + double a1 = a + iteration * cycle; + double b1 = a1 + cycle; + + double epsabs1 = eps * factor; + + int status = gsl_integration_qawo (f, a1, epsabs1, 0.0, limit, + cycle_workspace, wf, + &area1, &error1); + + append_interval (workspace, a1, b1, area1, error1); + + factor *= p; + + area = area + area1; + errsum = errsum + error1; + + /* estimate the truncation error as 50 times the final term */ + + truncation_error = 50 * fabs (area1); + + total_error = errsum + truncation_error; + + if (total_error < epsabs && iteration > 4) + { + goto compute_result; + } + + if (error1 > correc) + { + correc = error1; + } + + if (status) + { + eps = GSL_MAX_DBL (initial_eps, correc * (1.0 - p)); + } + + if (status && total_error < 10 * correc && iteration > 3) + { + goto compute_result; + } + + append_table (&table, area); + + if (table.n < 2) + { + continue; + } + + qelg (&table, &reseps, &erreps); + + ktmin++; + + if (ktmin >= 15 && err_ext < 0.001 * total_error) + { + error_type = 4; + } + + if (erreps < err_ext) + { + ktmin = 0; + err_ext = erreps; + res_ext = reseps; + + if (err_ext + 10 * correc <= epsabs) + break; + if (err_ext <= epsabs && 10 * correc >= epsabs) + break; + } + + } + + if (iteration == limit) + error_type = 1; + + if (err_ext == GSL_DBL_MAX) + goto compute_result; + + err_ext = err_ext + 10 * correc; + + *result = res_ext; + *abserr = err_ext; + + if (error_type == 0) + { + return GSL_SUCCESS ; + } + + if (res_ext != 0.0 && area != 0.0) + { + if (err_ext / fabs (res_ext) > errsum / fabs (area)) + goto compute_result; + } + else if (err_ext > errsum) + { + goto compute_result; + } + else if (area == 0.0) + { + goto return_error; + } + + if (error_type == 4) + { + err_ext = err_ext + truncation_error; + } + + goto return_error; + +compute_result: + + *result = area; + *abserr = total_error; + +return_error: + + if (error_type > 2) + error_type--; + + if (error_type == 0) + { + return GSL_SUCCESS; + } + else if (error_type == 1) + { + GSL_ERROR ("number of iterations was insufficient", GSL_EMAXITER); + } + else if (error_type == 2) + { + GSL_ERROR ("cannot reach tolerance because of roundoff error", + GSL_EROUND); + } + else if (error_type == 3) + { + GSL_ERROR ("bad integrand behavior found in the integration interval", + GSL_ESING); + } + else if (error_type == 4) + { + GSL_ERROR ("roundoff error detected in the extrapolation table", + GSL_EROUND); + } + else if (error_type == 5) + { + GSL_ERROR ("integral is divergent, or slowly convergent", + GSL_EDIVERGE); + } + else + { + GSL_ERROR ("could not integrate function", GSL_EFAILED); + } + +} + |