|
Packit |
67cb25 |
/* integration/qags.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_integration.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "initialise.c"
|
|
Packit |
67cb25 |
#include "set_initial.c"
|
|
Packit |
67cb25 |
#include "qpsrt.c"
|
|
Packit |
67cb25 |
#include "util.c"
|
|
Packit |
67cb25 |
#include "reset.c"
|
|
Packit |
67cb25 |
#include "qpsrt2.c"
|
|
Packit |
67cb25 |
#include "qelg.c"
|
|
Packit |
67cb25 |
#include "positivity.c"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int qags (const gsl_function * f, const double a, const double
|
|
Packit |
67cb25 |
b, const double epsabs, const double epsrel, const size_t limit,
|
|
Packit |
67cb25 |
gsl_integration_workspace * workspace, double *result, double *abserr,
|
|
Packit |
67cb25 |
gsl_integration_rule * q);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_integration_qags (const gsl_function *f,
|
|
Packit |
67cb25 |
double a, double b,
|
|
Packit |
67cb25 |
double epsabs, double epsrel, size_t limit,
|
|
Packit |
67cb25 |
gsl_integration_workspace * workspace,
|
|
Packit |
67cb25 |
double * result, double * abserr)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status = qags (f, a, b, epsabs, epsrel, limit,
|
|
Packit |
67cb25 |
workspace,
|
|
Packit |
67cb25 |
result, abserr,
|
|
Packit |
67cb25 |
&gsl_integration_qk21) ;
|
|
Packit |
67cb25 |
return status ;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* QAGI: evaluate an integral over an infinite range using the
|
|
Packit |
67cb25 |
transformation
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
integrate(f(x),-Inf,Inf) = integrate((f((1-t)/t) + f(-(1-t)/t))/t^2,0,1)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double i_transform (double t, void *params);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_integration_qagi (gsl_function * f,
|
|
Packit |
67cb25 |
double epsabs, double epsrel, size_t limit,
|
|
Packit |
67cb25 |
gsl_integration_workspace * workspace,
|
|
Packit |
67cb25 |
double *result, double *abserr)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_function f_transform;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
f_transform.function = &i_transform;
|
|
Packit |
67cb25 |
f_transform.params = f;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = qags (&f_transform, 0.0, 1.0,
|
|
Packit |
67cb25 |
epsabs, epsrel, limit,
|
|
Packit |
67cb25 |
workspace,
|
|
Packit |
67cb25 |
result, abserr,
|
|
Packit |
67cb25 |
&gsl_integration_qk15);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double
|
|
Packit |
67cb25 |
i_transform (double t, void *params)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_function *f = (gsl_function *) params;
|
|
Packit |
67cb25 |
double x = (1 - t) / t;
|
|
Packit |
67cb25 |
double y = GSL_FN_EVAL (f, x) + GSL_FN_EVAL (f, -x);
|
|
Packit |
67cb25 |
return (y / t) / t;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* QAGIL: Evaluate an integral over an infinite range using the
|
|
Packit |
67cb25 |
transformation,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
integrate(f(x),-Inf,b) = integrate(f(b-(1-t)/t)/t^2,0,1)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
struct il_params { double b ; gsl_function * f ; } ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double il_transform (double t, void *params);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_integration_qagil (gsl_function * f,
|
|
Packit |
67cb25 |
double b,
|
|
Packit |
67cb25 |
double epsabs, double epsrel, size_t limit,
|
|
Packit |
67cb25 |
gsl_integration_workspace * workspace,
|
|
Packit |
67cb25 |
double *result, double *abserr)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_function f_transform;
|
|
Packit |
67cb25 |
struct il_params transform_params ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
transform_params.b = b ;
|
|
Packit |
67cb25 |
transform_params.f = f ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
f_transform.function = &il_transform;
|
|
Packit |
67cb25 |
f_transform.params = &transform_params;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = qags (&f_transform, 0.0, 1.0,
|
|
Packit |
67cb25 |
epsabs, epsrel, limit,
|
|
Packit |
67cb25 |
workspace,
|
|
Packit |
67cb25 |
result, abserr,
|
|
Packit |
67cb25 |
&gsl_integration_qk15);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double
|
|
Packit |
67cb25 |
il_transform (double t, void *params)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
struct il_params *p = (struct il_params *) params;
|
|
Packit |
67cb25 |
double b = p->b;
|
|
Packit |
67cb25 |
gsl_function * f = p->f;
|
|
Packit |
67cb25 |
double x = b - (1 - t) / t;
|
|
Packit |
67cb25 |
double y = GSL_FN_EVAL (f, x);
|
|
Packit |
67cb25 |
return (y / t) / t;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* QAGIU: Evaluate an integral over an infinite range using the
|
|
Packit |
67cb25 |
transformation
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
integrate(f(x),a,Inf) = integrate(f(a+(1-t)/t)/t^2,0,1)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
struct iu_params { double a ; gsl_function * f ; } ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double iu_transform (double t, void *params);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_integration_qagiu (gsl_function * f,
|
|
Packit |
67cb25 |
double a,
|
|
Packit |
67cb25 |
double epsabs, double epsrel, size_t limit,
|
|
Packit |
67cb25 |
gsl_integration_workspace * workspace,
|
|
Packit |
67cb25 |
double *result, double *abserr)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int status;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_function f_transform;
|
|
Packit |
67cb25 |
struct iu_params transform_params ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
transform_params.a = a ;
|
|
Packit |
67cb25 |
transform_params.f = f ;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
f_transform.function = &iu_transform;
|
|
Packit |
67cb25 |
f_transform.params = &transform_params;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
status = qags (&f_transform, 0.0, 1.0,
|
|
Packit |
67cb25 |
epsabs, epsrel, limit,
|
|
Packit |
67cb25 |
workspace,
|
|
Packit |
67cb25 |
result, abserr,
|
|
Packit |
67cb25 |
&gsl_integration_qk15);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double
|
|
Packit |
67cb25 |
iu_transform (double t, void *params)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
struct iu_params *p = (struct iu_params *) params;
|
|
Packit |
67cb25 |
double a = p->a;
|
|
Packit |
67cb25 |
gsl_function * f = p->f;
|
|
Packit |
67cb25 |
double x = a + (1 - t) / t;
|
|
Packit |
67cb25 |
double y = GSL_FN_EVAL (f, x);
|
|
Packit |
67cb25 |
return (y / t) / t;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Main integration function */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
qags (const gsl_function * f,
|
|
Packit |
67cb25 |
const double a, const double b,
|
|
Packit |
67cb25 |
const double epsabs, const double epsrel,
|
|
Packit |
67cb25 |
const size_t limit,
|
|
Packit |
67cb25 |
gsl_integration_workspace * workspace,
|
|
Packit |
67cb25 |
double *result, double *abserr,
|
|
Packit |
67cb25 |
gsl_integration_rule * q)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double area, errsum;
|
|
Packit |
67cb25 |
double res_ext, err_ext;
|
|
Packit |
67cb25 |
double result0, abserr0, resabs0, resasc0;
|
|
Packit |
67cb25 |
double tolerance;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double ertest = 0;
|
|
Packit |
67cb25 |
double error_over_large_intervals = 0;
|
|
Packit |
67cb25 |
double reseps = 0, abseps = 0, correc = 0;
|
|
Packit |
67cb25 |
size_t ktmin = 0;
|
|
Packit |
67cb25 |
int roundoff_type1 = 0, roundoff_type2 = 0, roundoff_type3 = 0;
|
|
Packit |
67cb25 |
int error_type = 0, error_type2 = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
size_t iteration = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int positive_integrand = 0;
|
|
Packit |
67cb25 |
int extrapolate = 0;
|
|
Packit |
67cb25 |
int disallow_extrapolation = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
struct extrapolation_table table;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Initialize results */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
initialise (workspace, a, b);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*result = 0;
|
|
Packit |
67cb25 |
*abserr = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (limit > workspace->limit)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Test on accuracy */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("tolerance cannot be achieved with given epsabs and epsrel",
|
|
Packit |
67cb25 |
GSL_EBADTOL);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Perform the first integration */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
q (f, a, b, &result0, &abserr0, &resabs0, &resasc0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
set_initial_result (workspace, result0, abserr0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (abserr0 <= 100 * GSL_DBL_EPSILON * resabs0 && abserr0 > tolerance)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*result = result0;
|
|
Packit |
67cb25 |
*abserr = abserr0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_ERROR ("cannot reach tolerance because of roundoff error"
|
|
Packit |
67cb25 |
"on first attempt", GSL_EROUND);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if ((abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*result = result0;
|
|
Packit |
67cb25 |
*abserr = abserr0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (limit == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*result = result0;
|
|
Packit |
67cb25 |
*abserr = abserr0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_ERROR ("a maximum of one iteration was insufficient", GSL_EMAXITER);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Initialization */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
initialise_table (&table);
|
|
Packit |
67cb25 |
append_table (&table, result0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
area = result0;
|
|
Packit |
67cb25 |
errsum = abserr0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
res_ext = result0;
|
|
Packit |
67cb25 |
err_ext = GSL_DBL_MAX;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
positive_integrand = test_positivity (result0, resabs0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
iteration = 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
do
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t current_level;
|
|
Packit |
67cb25 |
double a1, b1, a2, b2;
|
|
Packit |
67cb25 |
double a_i, b_i, r_i, e_i;
|
|
Packit |
67cb25 |
double area1 = 0, area2 = 0, area12 = 0;
|
|
Packit |
67cb25 |
double error1 = 0, error2 = 0, error12 = 0;
|
|
Packit |
67cb25 |
double resasc1, resasc2;
|
|
Packit |
67cb25 |
double resabs1, resabs2;
|
|
Packit |
67cb25 |
double last_e_i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Bisect the subinterval with the largest error estimate */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
current_level = workspace->level[workspace->i] + 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
a1 = a_i;
|
|
Packit |
67cb25 |
b1 = 0.5 * (a_i + b_i);
|
|
Packit |
67cb25 |
a2 = b1;
|
|
Packit |
67cb25 |
b2 = b_i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
iteration++;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
q (f, a1, b1, &area1, &error1, &resabs1, &resasc1);
|
|
Packit |
67cb25 |
q (f, a2, b2, &area2, &error2, &resabs2, &resasc2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
area12 = area1 + area2;
|
|
Packit |
67cb25 |
error12 = error1 + error2;
|
|
Packit |
67cb25 |
last_e_i = e_i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Improve previous approximations to the integral and test for
|
|
Packit |
67cb25 |
accuracy.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
We write these expressions in the same way as the original
|
|
Packit |
67cb25 |
QUADPACK code so that the rounding errors are the same, which
|
|
Packit |
67cb25 |
makes testing easier. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
errsum = errsum + error12 - e_i;
|
|
Packit |
67cb25 |
area = area + area12 - r_i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (resasc1 != error1 && resasc2 != error2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double delta = r_i - area12;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (fabs (delta) <= 1.0e-5 * fabs (area12) && error12 >= 0.99 * e_i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (!extrapolate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
roundoff_type1++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
roundoff_type2++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
if (iteration > 10 && error12 > e_i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
roundoff_type3++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Test for roundoff and eventually set error flag */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (roundoff_type1 + roundoff_type2 >= 10 || roundoff_type3 >= 20)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
error_type = 2; /* round off error */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (roundoff_type2 >= 5)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
error_type2 = 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* set error flag in the case of bad integrand behaviour at
|
|
Packit |
67cb25 |
a point of the integration range */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (subinterval_too_small (a1, a2, b2))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
error_type = 4;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* append the newly-created intervals to the list */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (errsum <= tolerance)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
goto compute_result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (error_type)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (iteration >= limit - 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
error_type = 1;
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (iteration == 2) /* set up variables on first iteration */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
error_over_large_intervals = errsum;
|
|
Packit |
67cb25 |
ertest = tolerance;
|
|
Packit |
67cb25 |
append_table (&table, area);
|
|
Packit |
67cb25 |
continue;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (disallow_extrapolation)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
continue;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
error_over_large_intervals += -last_e_i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (current_level < workspace->maximum_level)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
error_over_large_intervals += error12;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (!extrapolate)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* test whether the interval to be bisected next is the
|
|
Packit |
67cb25 |
smallest interval. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (large_interval (workspace))
|
|
Packit |
67cb25 |
continue;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
extrapolate = 1;
|
|
Packit |
67cb25 |
workspace->nrmax = 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (!error_type2 && error_over_large_intervals > ertest)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (increase_nrmax (workspace))
|
|
Packit |
67cb25 |
continue;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Perform extrapolation */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
append_table (&table, area);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
qelg (&table, &reseps, &abseps);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
ktmin++;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (ktmin > 5 && err_ext < 0.001 * errsum)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
error_type = 5;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (abseps < err_ext)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ktmin = 0;
|
|
Packit |
67cb25 |
err_ext = abseps;
|
|
Packit |
67cb25 |
res_ext = reseps;
|
|
Packit |
67cb25 |
correc = error_over_large_intervals;
|
|
Packit |
67cb25 |
ertest = GSL_MAX_DBL (epsabs, epsrel * fabs (reseps));
|
|
Packit |
67cb25 |
if (err_ext <= ertest)
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Prepare bisection of the smallest interval. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (table.n == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
disallow_extrapolation = 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (error_type == 5)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* work on interval with largest error */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
reset_nrmax (workspace);
|
|
Packit |
67cb25 |
extrapolate = 0;
|
|
Packit |
67cb25 |
error_over_large_intervals = errsum;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
while (iteration < limit);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*result = res_ext;
|
|
Packit |
67cb25 |
*abserr = err_ext;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (err_ext == GSL_DBL_MAX)
|
|
Packit |
67cb25 |
goto compute_result;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (error_type || error_type2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (error_type2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
err_ext += correc;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (error_type == 0)
|
|
Packit |
67cb25 |
error_type = 3;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (res_ext != 0.0 && area != 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (err_ext / fabs (res_ext) > errsum / fabs (area))
|
|
Packit |
67cb25 |
goto compute_result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (err_ext > errsum)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
goto compute_result;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (area == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
goto return_error;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Test on divergence. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double max_area = GSL_MAX_DBL (fabs (res_ext), fabs (area));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (!positive_integrand && max_area < 0.01 * resabs0)
|
|
Packit |
67cb25 |
goto return_error;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ratio = res_ext / area;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (ratio < 0.01 || ratio > 100.0 || errsum > fabs (area))
|
|
Packit |
67cb25 |
error_type = 6;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
goto return_error;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
compute_result:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*result = sum_results (workspace);
|
|
Packit |
67cb25 |
*abserr = errsum;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return_error:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (error_type > 2)
|
|
Packit |
67cb25 |
error_type--;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (error_type == 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (error_type == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("number of iterations was insufficient", GSL_EMAXITER);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (error_type == 2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("cannot reach tolerance because of roundoff error",
|
|
Packit |
67cb25 |
GSL_EROUND);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (error_type == 3)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("bad integrand behavior found in the integration interval",
|
|
Packit |
67cb25 |
GSL_ESING);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (error_type == 4)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("roundoff error detected in the extrapolation table",
|
|
Packit |
67cb25 |
GSL_EROUND);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (error_type == 5)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("integral is divergent, or slowly convergent",
|
|
Packit |
67cb25 |
GSL_EDIVERGE);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("could not integrate function", GSL_EFAILED);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|