123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266 |
- /* integration/qag.c
- *
- * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 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 3 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 "gsl__config.h"
- #include <stdlib.h>
- #include "gsl_math.h"
- #include "gsl_errno.h"
- #include "gsl_integration.h"
- #include "gsl_integration__initialise.c"
- #include "gsl_integration__set_initial.c"
- #include "gsl_integration__qpsrt.c"
- #include "gsl_integration__util.c"
- static int
- qag (const gsl_function *f,
- const double a, const double b,
- const double epsabs, const double epsrel,
- const size_t limit,
- gsl_integration_workspace * workspace,
- double * result, double * abserr,
- gsl_integration_rule * q) ;
- int
- gsl_integration_qag (const gsl_function *f,
- double a, double b,
- double epsabs, double epsrel, size_t limit,
- int key,
- gsl_integration_workspace * workspace,
- double * result, double * abserr)
- {
- int status ;
- gsl_integration_rule * integration_rule = gsl_integration_qk15 ;
- if (key < GSL_INTEG_GAUSS15)
- {
- key = GSL_INTEG_GAUSS15 ;
- }
- else if (key > GSL_INTEG_GAUSS61)
- {
- key = GSL_INTEG_GAUSS61 ;
- }
- switch (key)
- {
- case GSL_INTEG_GAUSS15:
- integration_rule = gsl_integration_qk15 ;
- break ;
- case GSL_INTEG_GAUSS21:
- integration_rule = gsl_integration_qk21 ;
- break ;
- case GSL_INTEG_GAUSS31:
- integration_rule = gsl_integration_qk31 ;
- break ;
- case GSL_INTEG_GAUSS41:
- integration_rule = gsl_integration_qk41 ;
- break ;
- case GSL_INTEG_GAUSS51:
- integration_rule = gsl_integration_qk51 ;
- break ;
- case GSL_INTEG_GAUSS61:
- integration_rule = gsl_integration_qk61 ;
- break ;
- default:
- GSL_ERROR("value of key does specify a known integration rule",
- GSL_EINVAL) ;
- }
- status = qag (f, a, b, epsabs, epsrel, limit,
- workspace,
- result, abserr,
- integration_rule) ;
-
- return status ;
- }
- static int
- qag (const gsl_function * f,
- const double a, const double b,
- const double epsabs, const double epsrel,
- const size_t limit,
- gsl_integration_workspace * workspace,
- double *result, double *abserr,
- gsl_integration_rule * q)
- {
- double area, errsum;
- double result0, abserr0, resabs0, resasc0;
- double tolerance;
- size_t iteration = 0;
- int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;
- double round_off;
- /* Initialize results */
- initialise (workspace, a, b);
- *result = 0;
- *abserr = 0;
- if (limit > workspace->limit)
- {
- GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
- }
- if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
- {
- GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
- GSL_EBADTOL);
- }
- /* perform the first integration */
- q (f, a, b, &result0, &abserr0, &resabs0, &resasc0);
- set_initial_result (workspace, result0, abserr0);
- /* Test on accuracy */
- tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
- /* need IEEE rounding here to match original quadpack behavior */
- round_off = GSL_COERCE_DBL (50 * GSL_DBL_EPSILON * resabs0);
- if (abserr0 <= round_off && abserr0 > tolerance)
- {
- *result = result0;
- *abserr = abserr0;
- GSL_ERROR ("cannot reach tolerance because of roundoff error "
- "on first attempt", GSL_EROUND);
- }
- else if ((abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0.0)
- {
- *result = result0;
- *abserr = abserr0;
- return GSL_SUCCESS;
- }
- else if (limit == 1)
- {
- *result = result0;
- *abserr = abserr0;
- GSL_ERROR ("a maximum of one iteration was insufficient", GSL_EMAXITER);
- }
- area = result0;
- errsum = abserr0;
- iteration = 1;
- do
- {
- double a1, b1, a2, b2;
- double a_i, b_i, r_i, e_i;
- double area1 = 0, area2 = 0, area12 = 0;
- double error1 = 0, error2 = 0, error12 = 0;
- double resasc1, resasc2;
- double resabs1, resabs2;
- /* Bisect the subinterval with the largest error estimate */
- retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
- a1 = a_i;
- b1 = 0.5 * (a_i + b_i);
- a2 = b1;
- b2 = b_i;
- q (f, a1, b1, &area1, &error1, &resabs1, &resasc1);
- q (f, a2, b2, &area2, &error2, &resabs2, &resasc2);
- area12 = area1 + area2;
- error12 = error1 + error2;
- errsum += (error12 - e_i);
- area += area12 - r_i;
- if (resasc1 != error1 && resasc2 != error2)
- {
- double delta = r_i - area12;
- if (fabs (delta) <= 1.0e-5 * fabs (area12) && error12 >= 0.99 * e_i)
- {
- roundoff_type1++;
- }
- if (iteration >= 10 && error12 > e_i)
- {
- roundoff_type2++;
- }
- }
- tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
- if (errsum > tolerance)
- {
- if (roundoff_type1 >= 6 || roundoff_type2 >= 20)
- {
- error_type = 2; /* round off error */
- }
- /* set error flag in the case of bad integrand behaviour at
- a point of the integration range */
- if (subinterval_too_small (a1, a2, b2))
- {
- error_type = 3;
- }
- }
- update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);
- retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
- iteration++;
- }
- while (iteration < limit && !error_type && errsum > tolerance);
- *result = sum_results (workspace);
- *abserr = errsum;
- if (errsum <= tolerance)
- {
- return GSL_SUCCESS;
- }
- else if (error_type == 2)
- {
- GSL_ERROR ("roundoff error prevents tolerance from being achieved",
- GSL_EROUND);
- }
- else if (error_type == 3)
- {
- GSL_ERROR ("bad integrand behavior found in the integration interval",
- GSL_ESING);
- }
- else if (iteration == limit)
- {
- GSL_ERROR ("maximum number of subdivisions reached", GSL_EMAXITER);
- }
- else
- {
- GSL_ERROR ("could not integrate function", GSL_EFAILED);
- }
- }
|