UpsilonNumworks/Upsilon
1
0
Fork 0
mirror of https://github.com/UpsilonNumworks/Upsilon.git synced 2026-03-24 00:00:44 +01:00
Code Issues Packages Projects Releases Wiki Activity

[poincare] Clearer variable names in Gauss-Kronrod algorithm

Browse Source
This commit is contained in:
Léa Saviot
2019-01-21 16:18:51 +01:00
committed by EmilieNumworks
parent 0a5290d5a6
commit b97873473b
1 changed files with 41 additions and 32 deletions
Download Patch File Download Diff File Expand all files Collapse all files

73
poincare/src/integral.cpp
View File

@@ -113,72 +113,81 @@ IntegralNode::DetailedResult<T> IntegralNode::kronrodGaussQuadrature(T a, T b, C
static T max = sizeof(T) == sizeof(double) ? DBL_MAX : FLT_MAX;
/* We here use Kronrod-Legendre quadrature with n = 21
* The abscissa and weights are taken from QUADPACK library. */
const static T wg[5]= {0.066671344308688137593568809893332, 0.149451349150580593145776339657697,
0.219086362515982043995534934228163, 0.269266719309996355091226921569469, 0.295524224714752870173892994651338};
const static T xgk[11]= {0.995657163025808080735527280689003, 0.973906528517171720077964012084452,
// Abscissae for the gauss (odd weights) and kronrod rules (all weights)
const static T x[11]= {0.995657163025808080735527280689003, 0.973906528517171720077964012084452,
0.930157491355708226001207180059508, 0.865063366688984510732096688423493, 0.780817726586416897063717578345042,
0.679409568299024406234327365114874, 0.562757134668604683339000099272694, 0.433395394129247190799265943165784,
0.294392862701460198131126603103866, 0.148874338981631210884826001129720, 0.000000000000000000000000000000000};
const static T wgk[11]= {0.011694638867371874278064396062192, 0.032558162307964727478818972459390,
// Weights for the gauss integral
const static T wGauss[5]= {0.066671344308688137593568809893332, 0.149451349150580593145776339657697,
0.219086362515982043995534934228163, 0.269266719309996355091226921569469, 0.295524224714752870173892994651338};
// Weights for the kronrod rule
const static T wKronrod[11]= {0.011694638867371874278064396062192, 0.032558162307964727478818972459390,
0.054755896574351996031381300244580, 0.075039674810919952767043140916190, 0.093125454583697605535065465083366,
0.109387158802297641899210590325805, 0.123491976262065851077958109831074, 0.134709217311473325928054001771707,
0.142775938577060080797094273138717, 0.147739104901338491374841515972068, 0.149445554002916905664936468389821};
T fv1[10];
T fv2[10];
T centr = 0.5*(a+b);
T hlgth = 0.5*(b-a);
T dhlgth = std::fabs(hlgth);
T center = 0.5 * (a+b);
T halfLength = 0.5 * (b-a);
T absHalfLength = std::fabs(halfLength);
DetailedResult<T> errorResult;
errorResult.integral = NAN;
errorResult.absoluteError = 0;
T resg = 0;
T fc = functionValueAtAbscissa(centr, context, complexFormat, angleUnit);
if (std::isnan(fc)) {
T gaussIntegral = 0;
T fCenter = functionValueAtAbscissa(center, context, complexFormat, angleUnit);
if (std::isnan(fCenter)) {
return errorResult;
}
T resk = wgk[10]*fc;
T resabs = std::fabs(resk);
T kronrodIntegral = wKronrod[10] * fCenter;
T absKronrodIntegral = std::fabs(kronrodIntegral);
for (int j = 0; j < 10; j++) {
T absc = hlgth*xgk[j];
T fval1 = functionValueAtAbscissa(centr-absc, context, complexFormat, angleUnit);
T xDelta = halfLength * x[j];
T fval1 = functionValueAtAbscissa(center - xDelta, context, complexFormat, angleUnit);
if (std::isnan(fval1)) {
return errorResult;
}
T fval2 = functionValueAtAbscissa(centr+absc, context, complexFormat, angleUnit);
T fval2 = functionValueAtAbscissa(center + xDelta, context, complexFormat, angleUnit);
if (std::isnan(fval2)) {
return errorResult;
}
fv1[j] = fval1;
fv2[j] = fval2;
T fsum = fval1+fval2;
if (j%2 == 1) {
resg += wg[j/2]*fsum;
T fsum = fval1 + fval2;
if (j % 2 == 1) {
gaussIntegral += wGauss[j/2] * fsum;
}
resk += wgk[j]*fsum;
resabs += wgk[j]*(std::fabs(fval1)+std::fabs(fval2));
kronrodIntegral += wKronrod[j] * fsum;
absKronrodIntegral += wKronrod[j] * (std::fabs(fval1) + std::fabs(fval2));
}
T reskh = resk*0.5;
T resasc = wgk[10]*std::fabs(fc-reskh);
T halfKronrodIntegral = 0.5 * kronrodIntegral;
T kronrodIntegralDifference = wKronrod[10] * std::fabs(fCenter - halfKronrodIntegral);
for (int j = 0; j < 10; j++) {
resasc += wgk[j]*(std::fabs(fv1[j]-reskh)+std::fabs(fv2[j]-reskh));
kronrodIntegralDifference += wKronrod[j] * (std::fabs(fv1[j] - halfKronrodIntegral) + std::fabs(fv2[j] - halfKronrodIntegral));
}
T integral = resk*hlgth;
resabs = resabs*dhlgth;
resasc = resasc*dhlgth;
T abserr = std::fabs((resk-resg)*hlgth);
if (resasc != 0 && abserr != 0) {
abserr = 1 > std::pow((T)(200*abserr/resasc), (T)1.5)? resasc*std::pow((T)(200*abserr/resasc), (T)1.5) : resasc;
T integral = kronrodIntegral * halfLength;
absKronrodIntegral = absKronrodIntegral * absHalfLength;
kronrodIntegralDifference = kronrodIntegralDifference * absHalfLength;
T absError = std::fabs((kronrodIntegral - gaussIntegral) * halfLength);
if (kronrodIntegralDifference != 0 && absError != 0) {
T errorCoefficient = std::pow((T)(200*absError/kronrodIntegralDifference), (T)1.5);
absError = 1 > errorCoefficient ? kronrodIntegralDifference * errorCoefficient : kronrodIntegralDifference;
}
if (resabs > max/(50.0*epsilon)) {
abserr = abserr > epsilon*50*resabs ? abserr : epsilon*50*resabs;
if (absKronrodIntegral > max/(50.0 * epsilon)) {
T minError = epsilon * 50 * absKronrodIntegral;
absError = absError > minError ? absError : minError;
}
DetailedResult<T> result;
result.integral = integral;
result.absoluteError = abserr;
result.absoluteError = absError;
return result;
}