001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math4.legacy.analysis.integration;
018
019import org.apache.commons.math4.legacy.exception.NumberIsTooLargeException;
020import org.apache.commons.math4.core.jdkmath.JdkMath;
021
022/**
023 * Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html">
024 * Romberg Algorithm</a> for integration of real univariate functions. For
025 * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
026 * chapter 3.
027 * <p>
028 * Romberg integration employs k successive refinements of the trapezoid
029 * rule to remove error terms less than order O(N^(-2k)). Simpson's rule
030 * is a special case of k = 2.</p>
031 *
032 * @since 1.2
033 */
034public class RombergIntegrator extends BaseAbstractUnivariateIntegrator {
035
036    /** Maximal number of iterations for Romberg. */
037    public static final int ROMBERG_MAX_ITERATIONS_COUNT = 32;
038
039    /**
040     * Build a Romberg integrator with given accuracies and iterations counts.
041     * @param relativeAccuracy relative accuracy of the result
042     * @param absoluteAccuracy absolute accuracy of the result
043     * @param minimalIterationCount minimum number of iterations
044     * @param maximalIterationCount maximum number of iterations
045     * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
046     * @exception org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException if minimal number of iterations
047     * is not strictly positive
048     * @exception org.apache.commons.math4.legacy.exception.NumberIsTooSmallException if maximal number of iterations
049     * is lesser than or equal to the minimal number of iterations
050     * @exception NumberIsTooLargeException if maximal number of iterations
051     * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT}
052     */
053    public RombergIntegrator(final double relativeAccuracy,
054                             final double absoluteAccuracy,
055                             final int minimalIterationCount,
056                             final int maximalIterationCount) {
057        super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
058        if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) {
059            throw new NumberIsTooLargeException(maximalIterationCount,
060                                                ROMBERG_MAX_ITERATIONS_COUNT, false);
061        }
062    }
063
064    /**
065     * Build a Romberg integrator with given iteration counts.
066     * @param minimalIterationCount minimum number of iterations
067     * @param maximalIterationCount maximum number of iterations
068     * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
069     * @exception org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException if minimal number of iterations
070     * is not strictly positive
071     * @exception org.apache.commons.math4.legacy.exception.NumberIsTooSmallException if maximal number of iterations
072     * is lesser than or equal to the minimal number of iterations
073     * @exception NumberIsTooLargeException if maximal number of iterations
074     * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT}
075     */
076    public RombergIntegrator(final int minimalIterationCount,
077                             final int maximalIterationCount) {
078        super(minimalIterationCount, maximalIterationCount);
079        if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) {
080            throw new NumberIsTooLargeException(maximalIterationCount,
081                                                ROMBERG_MAX_ITERATIONS_COUNT, false);
082        }
083    }
084
085    /**
086     * Construct a Romberg integrator with default settings
087     * (max iteration count set to {@link #ROMBERG_MAX_ITERATIONS_COUNT}).
088     */
089    public RombergIntegrator() {
090        super(DEFAULT_MIN_ITERATIONS_COUNT, ROMBERG_MAX_ITERATIONS_COUNT);
091    }
092
093    /** {@inheritDoc} */
094    @Override
095    protected double doIntegrate() {
096        final int m = iterations.getMaximalCount() + 1;
097        double[] previousRow = new double[m];
098        double[] currentRow = new double[m];
099
100        TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
101        currentRow[0] = qtrap.stage(this, 0);
102        iterations.increment();
103        double olds = currentRow[0];
104        while (true) {
105
106            final int i = iterations.getCount();
107
108            // switch rows
109            final double[] tmpRow = previousRow;
110            previousRow = currentRow;
111            currentRow = tmpRow;
112
113            currentRow[0] = qtrap.stage(this, i);
114            iterations.increment();
115            for (int j = 1; j <= i; j++) {
116                // Richardson extrapolation coefficient
117                final double r = (1L << (2 * j)) - 1;
118                final double tIJm1 = currentRow[j - 1];
119                currentRow[j] = tIJm1 + (tIJm1 - previousRow[j - 1]) / r;
120            }
121            final double s = currentRow[i];
122            if (i >= getMinimalIterationCount()) {
123                final double delta  = JdkMath.abs(s - olds);
124                final double rLimit = getRelativeAccuracy() * (JdkMath.abs(olds) + JdkMath.abs(s)) * 0.5;
125                if (delta <= rLimit || delta <= getAbsoluteAccuracy()) {
126                    return s;
127                }
128            }
129            olds = s;
130        }
131    }
132}