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.math3.analysis.integration;
018    
019    import org.apache.commons.math3.exception.MaxCountExceededException;
020    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
021    import org.apache.commons.math3.exception.NumberIsTooLargeException;
022    import org.apache.commons.math3.exception.NumberIsTooSmallException;
023    import org.apache.commons.math3.exception.TooManyEvaluationsException;
024    import org.apache.commons.math3.util.FastMath;
025    
026    /**
027     * Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html">
028     * Romberg Algorithm</a> for integration of real univariate functions. For
029     * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
030     * chapter 3.
031     * <p>
032     * Romberg integration employs k successive refinements of the trapezoid
033     * rule to remove error terms less than order O(N^(-2k)). Simpson's rule
034     * is a special case of k = 2.</p>
035     *
036     * @version $Id: RombergIntegrator.java 1364387 2012-07-22 18:14:11Z tn $
037     * @since 1.2
038     */
039    public class RombergIntegrator extends BaseAbstractUnivariateIntegrator {
040    
041        /** Maximal number of iterations for Romberg. */
042        public static final int ROMBERG_MAX_ITERATIONS_COUNT = 32;
043    
044        /**
045         * Build a Romberg integrator with given accuracies and iterations counts.
046         * @param relativeAccuracy relative accuracy of the result
047         * @param absoluteAccuracy absolute accuracy of the result
048         * @param minimalIterationCount minimum number of iterations
049         * @param maximalIterationCount maximum number of iterations
050         * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
051         * @exception NotStrictlyPositiveException if minimal number of iterations
052         * is not strictly positive
053         * @exception NumberIsTooSmallException if maximal number of iterations
054         * is lesser than or equal to the minimal number of iterations
055         * @exception NumberIsTooLargeException if maximal number of iterations
056         * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT}
057         */
058        public RombergIntegrator(final double relativeAccuracy,
059                                 final double absoluteAccuracy,
060                                 final int minimalIterationCount,
061                                 final int maximalIterationCount)
062            throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {
063            super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
064            if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) {
065                throw new NumberIsTooLargeException(maximalIterationCount,
066                                                    ROMBERG_MAX_ITERATIONS_COUNT, false);
067            }
068        }
069    
070        /**
071         * Build a Romberg integrator with given iteration counts.
072         * @param minimalIterationCount minimum number of iterations
073         * @param maximalIterationCount maximum number of iterations
074         * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
075         * @exception NotStrictlyPositiveException if minimal number of iterations
076         * is not strictly positive
077         * @exception NumberIsTooSmallException if maximal number of iterations
078         * is lesser than or equal to the minimal number of iterations
079         * @exception NumberIsTooLargeException if maximal number of iterations
080         * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT}
081         */
082        public RombergIntegrator(final int minimalIterationCount,
083                                 final int maximalIterationCount)
084            throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {
085            super(minimalIterationCount, maximalIterationCount);
086            if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) {
087                throw new NumberIsTooLargeException(maximalIterationCount,
088                                                    ROMBERG_MAX_ITERATIONS_COUNT, false);
089            }
090        }
091    
092        /**
093         * Construct a Romberg integrator with default settings
094         * (max iteration count set to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
095         */
096        public RombergIntegrator() {
097            super(DEFAULT_MIN_ITERATIONS_COUNT, ROMBERG_MAX_ITERATIONS_COUNT);
098        }
099    
100        /** {@inheritDoc} */
101        @Override
102        protected double doIntegrate()
103            throws TooManyEvaluationsException, MaxCountExceededException {
104    
105            final int m = iterations.getMaximalCount() + 1;
106            double previousRow[] = new double[m];
107            double currentRow[]  = new double[m];
108    
109            TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
110            currentRow[0] = qtrap.stage(this, 0);
111            iterations.incrementCount();
112            double olds = currentRow[0];
113            while (true) {
114    
115                final int i = iterations.getCount();
116    
117                // switch rows
118                final double[] tmpRow = previousRow;
119                previousRow = currentRow;
120                currentRow = tmpRow;
121    
122                currentRow[0] = qtrap.stage(this, i);
123                iterations.incrementCount();
124                for (int j = 1; j <= i; j++) {
125                    // Richardson extrapolation coefficient
126                    final double r = (1L << (2 * j)) - 1;
127                    final double tIJm1 = currentRow[j - 1];
128                    currentRow[j] = tIJm1 + (tIJm1 - previousRow[j - 1]) / r;
129                }
130                final double s = currentRow[i];
131                if (i >= getMinimalIterationCount()) {
132                    final double delta  = FastMath.abs(s - olds);
133                    final double rLimit = getRelativeAccuracy() * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5;
134                    if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) {
135                        return s;
136                    }
137                }
138                olds = s;
139            }
140    
141        }
142    
143    }