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     */
017    package org.apache.commons.math.analysis.integration;
018    
019    import org.apache.commons.math.FunctionEvaluationException;
020    import org.apache.commons.math.MathRuntimeException;
021    import org.apache.commons.math.MaxIterationsExceededException;
022    import org.apache.commons.math.analysis.UnivariateRealFunction;
023    
024    /**
025     * Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html">
026     * Romberg Algorithm</a> for integration of real univariate functions. For
027     * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
028     * chapter 3.
029     * <p>
030     * Romberg integration employs k successive refinements of the trapezoid
031     * rule to remove error terms less than order O(N^(-2k)). Simpson's rule
032     * is a special case of k = 2.</p>
033     *  
034     * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
035     * @since 1.2
036     */
037    public class RombergIntegrator extends UnivariateRealIntegratorImpl {
038    
039        /**
040         * Construct an integrator for the given function.
041         * 
042         * @param f function to integrate
043         * @deprecated as of 2.0 the integrand function is passed as an argument
044         * to the {@link #integrate(UnivariateRealFunction, double, double)}method.
045         */
046        @Deprecated
047        public RombergIntegrator(UnivariateRealFunction f) {
048            super(f, 32);
049        }
050    
051        /**
052         * Construct an integrator.
053         */
054        public RombergIntegrator() {
055            super(32);
056        }
057    
058        /** {@inheritDoc} */
059        @Deprecated
060        public double integrate(final double min, final double max)
061            throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException {
062            return integrate(f, min, max);
063        }
064    
065        /** {@inheritDoc} */
066        public double integrate(final UnivariateRealFunction f,
067                                final double min, final double max)
068            throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException {
069            
070            int i = 1, j, m = maximalIterationCount + 1;
071            // Array structure here can be improved for better space
072            // efficiency because only the lower triangle is used.
073            double r, t[][] = new double[m][m], s, olds;
074    
075            clearResult();
076            verifyInterval(min, max);
077            verifyIterationCount();
078    
079            TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
080            t[0][0] = qtrap.stage(f, min, max, 0);
081            olds = t[0][0];
082            while (i <= maximalIterationCount) {
083                t[i][0] = qtrap.stage(f, min, max, i);
084                for (j = 1; j <= i; j++) {
085                    // Richardson extrapolation coefficient
086                    r = (1L << (2 * j)) -1;
087                    t[i][j] = t[i][j-1] + (t[i][j-1] - t[i-1][j-1]) / r;
088                }
089                s = t[i][i];
090                if (i >= minimalIterationCount) {
091                    final double delta = Math.abs(s - olds);
092                    final double rLimit =
093                        relativeAccuracy * (Math.abs(olds) + Math.abs(s)) * 0.5; 
094                    if ((delta <= rLimit) || (delta <= absoluteAccuracy)) {
095                        setResult(s, i);
096                        return result;
097                    }
098                }
099                olds = s;
100                i++;
101            }
102            throw new MaxIterationsExceededException(maximalIterationCount);
103        }
104    
105        /** {@inheritDoc} */
106        @Override
107        protected void verifyIterationCount() throws IllegalArgumentException {
108            super.verifyIterationCount();
109            // at most 32 bisection refinements due to higher order divider
110            if (maximalIterationCount > 32) {
111                throw MathRuntimeException.createIllegalArgumentException(
112                        "invalid iteration limits: min={0}, max={1}",
113                        0, 32);
114            }
115        }
116    }