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

[poincare] Fix recursion bug in adaptiveQuadrature() integration

Browse Source 
The recursive adaptiveQuadrature() function is initially called with the
maximum number of recursion levels to perform, so adaptiveQuadrature()
should decrement numberOfRecursions at each recursion level.

When adaptiveQuadrature() needs to recurse, it splits the interval in
two, calculates the integral for each half by calling itself
recursively, and adds the result from each half.  However, it then
ignores this result and returns NAN because it is checking a stale
quadKG.absoluteError value.  Rearrange the code to avoid that.
This commit is contained in:
Ian Abbott
2017-09-15 12:07:41 +01:00
committed by EmilieNumworks
parent 6392735e70
commit 28df37813c
1 changed files with 5 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
poincare/src/integral.cpp
View File

@@ -171,14 +171,14 @@ template<typename T>
T Integral::adaptiveQuadrature(T a, T b, T eps, int numberOfIterations, VariableContext<T> xContext, AngleUnit angleUnit) const {
DetailedResult<T> quadKG = kronrodGaussQuadrature(a, b, xContext, angleUnit);
T result = quadKG.integral;
if (numberOfIterations < k_maxNumberOfIterations && quadKG.absoluteError > eps) {
if (quadKG.absoluteError <= eps) {
return result;
} else if (--numberOfIterations > 0) {
T m = (a+b)/2;
result = adaptiveQuadrature<T>(a, m, eps/2, numberOfIterations+1, xContext, angleUnit) + adaptiveQuadrature<T>(m, b, eps/2, numberOfIterations+1, xContext, angleUnit);
}
if (quadKG.absoluteError > eps) {
return adaptiveQuadrature<T>(a, m, eps/2, numberOfIterations, xContext, angleUnit) + adaptiveQuadrature<T>(m, b, eps/2, numberOfIterations, xContext, angleUnit);
} else {
return NAN;
}
return result;
}
#endif