|
Packit |
67cb25 |
/* integration/qaws.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 <math.h>
|
|
Packit |
67cb25 |
#include <float.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.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 "append.c"
|
|
Packit |
67cb25 |
#include "qpsrt.c"
|
|
Packit |
67cb25 |
#include "util.c"
|
|
Packit |
67cb25 |
#include "qc25s.c"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_integration_qaws (gsl_function * f,
|
|
Packit |
67cb25 |
const double a, const double b,
|
|
Packit |
67cb25 |
gsl_integration_qaws_table * t,
|
|
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 |
{
|
|
Packit |
67cb25 |
double area, errsum;
|
|
Packit |
67cb25 |
double result0, abserr0;
|
|
Packit |
67cb25 |
double tolerance;
|
|
Packit |
67cb25 |
size_t iteration = 0;
|
|
Packit |
67cb25 |
int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;
|
|
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 |
if (b <= a)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("limits must form an ascending sequence, a < b", GSL_EINVAL) ;
|
|
Packit |
67cb25 |
}
|
|
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 |
{
|
|
Packit |
67cb25 |
double area1, area2;
|
|
Packit |
67cb25 |
double error1, error2;
|
|
Packit |
67cb25 |
int err_reliable1, err_reliable2;
|
|
Packit |
67cb25 |
double a1 = a;
|
|
Packit |
67cb25 |
double b1 = 0.5 * (a + b);
|
|
Packit |
67cb25 |
double a2 = b1;
|
|
Packit |
67cb25 |
double b2 = b;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
qc25s (f, a, b, a1, b1, t, &area1, &error1, &err_reliable1);
|
|
Packit |
67cb25 |
qc25s (f, a, b, a2, b2, t, &area2, &error2, &err_reliable2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (error1 > error2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
append_interval (workspace, a1, b1, area1, error1);
|
|
Packit |
67cb25 |
append_interval (workspace, a2, b2, area2, error2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
append_interval (workspace, a2, b2, area2, error2);
|
|
Packit |
67cb25 |
append_interval (workspace, a1, b1, area1, error1);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result0 = area1 + area2;
|
|
Packit |
67cb25 |
abserr0 = error1 + error2;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Test on accuracy */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Test on accuracy, use 0.01 relative error as an extra safety
|
|
Packit |
67cb25 |
margin on the first iteration (ignored for subsequent iterations) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (abserr0 < tolerance && abserr0 < 0.01 * fabs(result0))
|
|
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 |
area = result0;
|
|
Packit |
67cb25 |
errsum = abserr0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
iteration = 2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
do
|
|
Packit |
67cb25 |
{
|
|
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 |
int err_reliable1, err_reliable2;
|
|
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 |
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 |
qc25s (f, a, b, a1, b1, t, &area1, &error1, &err_reliable1);
|
|
Packit |
67cb25 |
qc25s (f, a, b, a2, b2, t, &area2, &error2, &err_reliable2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
area12 = area1 + area2;
|
|
Packit |
67cb25 |
error12 = error1 + error2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
errsum += (error12 - e_i);
|
|
Packit |
67cb25 |
area += area12 - r_i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (err_reliable1 && err_reliable2)
|
|
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 |
roundoff_type1++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
if (iteration >= 10 && error12 > e_i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
roundoff_type2++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (errsum > tolerance)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (roundoff_type1 >= 6 || roundoff_type2 >= 20)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
error_type = 2; /* round off error */
|
|
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 = 3;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
iteration++;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
while (iteration < limit && !error_type && errsum > tolerance);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*result = sum_results (workspace);
|
|
Packit |
67cb25 |
*abserr = errsum;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (errsum <= tolerance)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (error_type == 2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("roundoff error prevents tolerance from being achieved",
|
|
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 (iteration == limit)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("maximum number of subdivisions reached", GSL_EMAXITER);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
GSL_ERROR ("could not integrate function", GSL_EFAILED);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|