/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.analysis.integration;
import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.exception.NotStrictlyPositiveException;
import org.apache.commons.math3.exception.NumberIsTooLargeException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.exception.TooManyEvaluationsException;
import org.apache.commons.math3.util.FastMath;
/**
* Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html">
* Romberg Algorithm</a> for integration of real univariate functions. For
* reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
* chapter 3.
* <p>
* Romberg integration employs k successive refinements of the trapezoid
* rule to remove error terms less than order O(N^(-2k)). Simpson's rule
* is a special case of k = 2.</p>
*
* @since 1.2
*/
public class RombergIntegrator extends BaseAbstractUnivariateIntegrator {
/** Maximal number of iterations for Romberg. */
public static final int ROMBERG_MAX_ITERATIONS_COUNT = 32;
/**
* Build a Romberg integrator with given accuracies and iterations counts.
* @param relativeAccuracy relative accuracy of the result
* @param absoluteAccuracy absolute accuracy of the result
* @param minimalIterationCount minimum number of iterations
* @param maximalIterationCount maximum number of iterations
* (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
* @exception NotStrictlyPositiveException if minimal number of iterations
* is not strictly positive
* @exception NumberIsTooSmallException if maximal number of iterations
* is lesser than or equal to the minimal number of iterations
* @exception NumberIsTooLargeException if maximal number of iterations
* is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT}
*/
public RombergIntegrator(final double relativeAccuracy,
final double absoluteAccuracy,
final int minimalIterationCount,
final int maximalIterationCount)
throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {
super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) {
throw new NumberIsTooLargeException(maximalIterationCount,
ROMBERG_MAX_ITERATIONS_COUNT, false);
}
}
/**
* Build a Romberg integrator with given iteration counts.
* @param minimalIterationCount minimum number of iterations
* @param maximalIterationCount maximum number of iterations
* (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
* @exception NotStrictlyPositiveException if minimal number of iterations
* is not strictly positive
* @exception NumberIsTooSmallException if maximal number of iterations
* is lesser than or equal to the minimal number of iterations
* @exception NumberIsTooLargeException if maximal number of iterations
* is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT}
*/
public RombergIntegrator(final int minimalIterationCount,
final int maximalIterationCount)
throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {
super(minimalIterationCount, maximalIterationCount);
if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) {
throw new NumberIsTooLargeException(maximalIterationCount,
ROMBERG_MAX_ITERATIONS_COUNT, false);
}
}
/**
* Construct a Romberg integrator with default settings
* (max iteration count set to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
*/
public RombergIntegrator() {
super(DEFAULT_MIN_ITERATIONS_COUNT, ROMBERG_MAX_ITERATIONS_COUNT);
}
/** {@inheritDoc} */
@Override
protected double doIntegrate()
throws TooManyEvaluationsException, MaxCountExceededException {
final int m = getMaximalIterationCount() + 1;
double previousRow[] = new double[m];
double currentRow[] = new double[m];
TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
currentRow[0] = qtrap.stage(this, 0);
incrementCount();
double olds = currentRow[0];
while (true) {
final int i = getIterations();
// switch rows
final double[] tmpRow = previousRow;
previousRow = currentRow;
currentRow = tmpRow;
currentRow[0] = qtrap.stage(this, i);
incrementCount();
for (int j = 1; j <= i; j++) {
// Richardson extrapolation coefficient
final double r = (1L << (2 * j)) - 1;
final double tIJm1 = currentRow[j - 1];
currentRow[j] = tIJm1 + (tIJm1 - previousRow[j - 1]) / r;
}
final double s = currentRow[i];
if (i >= getMinimalIterationCount()) {
final double delta = FastMath.abs(s - olds);
final double rLimit = getRelativeAccuracy() * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5;
if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) {
return s;
}
}
olds = s;
}
}
}