| /* |
| * 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.math4.analysis.integration; |
| |
| import org.apache.commons.math4.exception.MaxCountExceededException; |
| import org.apache.commons.math4.exception.NotStrictlyPositiveException; |
| import org.apache.commons.math4.exception.NumberIsTooLargeException; |
| import org.apache.commons.math4.exception.NumberIsTooSmallException; |
| import org.apache.commons.math4.exception.TooManyEvaluationsException; |
| import org.apache.commons.math4.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) { |
| 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) { |
| 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() { |
| final int m = iterations.getMaximalCount() + 1; |
| double previousRow[] = new double[m]; |
| double currentRow[] = new double[m]; |
| |
| TrapezoidIntegrator qtrap = new TrapezoidIntegrator(); |
| currentRow[0] = qtrap.stage(this, 0); |
| iterations.increment(); |
| double olds = currentRow[0]; |
| while (true) { |
| |
| final int i = iterations.getCount(); |
| |
| // switch rows |
| final double[] tmpRow = previousRow; |
| previousRow = currentRow; |
| currentRow = tmpRow; |
| |
| currentRow[0] = qtrap.stage(this, i); |
| iterations.increment(); |
| 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; |
| } |
| |
| } |
| |
| } |