Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

improved array structure in Romberg integrator 

since the integrator uses only two rows at a time, we allocate only two arrays and reuse them

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@824822 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2009-10-13 15:56:51 +00:00
parent e975825996
commit dd9b7b477a
1 changed files with 20 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
src/main/java/org/apache/commons/math/analysis/integration/RombergIntegrator.java
View File

@ -67,37 +67,41 @@ public class RombergIntegrator extends UnivariateRealIntegratorImpl {
final double min, final double max)
throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException {
int i = 1, j, m = maximalIterationCount + 1;
// Array structure here can be improved for better space
// efficiency because only the lower triangle is used.
double r, t[][] = new double[m][m], s, olds;
final int m = maximalIterationCount + 1;
double previousRow[] = new double[m];
double currentRow[] = new double[m];
clearResult();
verifyInterval(min, max);
verifyIterationCount();
TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
t[0][0] = qtrap.stage(f, min, max, 0);
olds = t[0][0];
while (i <= maximalIterationCount) {
t[i][0] = qtrap.stage(f, min, max, i);
for (j = 1; j <= i; j++) {
currentRow[0] = qtrap.stage(f, min, max, 0);
double olds = currentRow[0];
for (int i = 1; i <= maximalIterationCount; ++i) {
// switch rows
final double[] tmpRow = previousRow;
previousRow = currentRow;
currentRow = tmpRow;
currentRow[0] = qtrap.stage(f, min, max, i);
for (int j = 1; j <= i; j++) {
// Richardson extrapolation coefficient
r = (1L << (2 * j)) -1;
t[i][j] = t[i][j-1] + (t[i][j-1] - t[i-1][j-1]) / r;
final double r = (1L << (2 * j)) - 1;
final double tIJm1 = currentRow[j - 1];
currentRow[j] = tIJm1 + (tIJm1 - previousRow[j - 1]) / r;
}
s = t[i][i];
final double s = currentRow[i];
if (i >= minimalIterationCount) {
final double delta = Math.abs(s - olds);
final double rLimit =
relativeAccuracy * (Math.abs(olds) + Math.abs(s)) * 0.5;
final double delta = Math.abs(s - olds);
final double rLimit = relativeAccuracy * (Math.abs(olds) + Math.abs(s)) * 0.5;
if ((delta <= rLimit) || (delta <= absoluteAccuracy)) {
setResult(s, i);
return result;
}
}
olds = s;
i++;
}
throw new MaxIterationsExceededException(maximalIterationCount);
}