Home | History | Annotate | Download | only in integration 
      1 /*
      2  * Licensed to the Apache Software Foundation (ASF) under one or more
      3  * contributor license agreements.  See the NOTICE file distributed with
      4  * this work for additional information regarding copyright ownership.
      5  * The ASF licenses this file to You under the Apache License, Version 2.0
      6  * (the "License"); you may not use this file except in compliance with
      7  * the License.  You may obtain a copy of the License at
      8  *
      9  *      http://www.apache.org/licenses/LICENSE-2.0
     10  *
     11  * Unless required by applicable law or agreed to in writing, software
     12  * distributed under the License is distributed on an "AS IS" BASIS,
     13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
     14  * See the License for the specific language governing permissions and
     15  * limitations under the License.
     16  */
     17 package org.apache.commons.math.analysis.integration;
     18 
     19 import org.apache.commons.math.FunctionEvaluationException;
     20 import org.apache.commons.math.MathRuntimeException;
     21 import org.apache.commons.math.MaxIterationsExceededException;
     22 import org.apache.commons.math.analysis.UnivariateRealFunction;
     23 import org.apache.commons.math.exception.util.LocalizedFormats;
     24 import org.apache.commons.math.util.FastMath;
     25 
     26 /**
     27  * Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html">
     28  * Romberg Algorithm</a> for integration of real univariate functions. For
     29  * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
     30  * chapter 3.
     31  * <p>
     32  * Romberg integration employs k successive refinements of the trapezoid
     33  * rule to remove error terms less than order O(N^(-2k)). Simpson's rule
     34  * is a special case of k = 2.</p>
     35  *
     36  * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 fvr. 2011) $
     37  * @since 1.2
     38  */
     39 public class RombergIntegrator extends UnivariateRealIntegratorImpl {
     40 
     41     /**
     42      * Construct an integrator for the given function.
     43      *
     44      * @param f function to integrate
     45      * @deprecated as of 2.0 the integrand function is passed as an argument
     46      * to the {@link #integrate(UnivariateRealFunction, double, double)}method.
     47      */
     48     @Deprecated
     49     public RombergIntegrator(UnivariateRealFunction f) {
     50         super(f, 32);
     51     }
     52 
     53     /**
     54      * Construct an integrator.
     55      */
     56     public RombergIntegrator() {
     57         super(32);
     58     }
     59 
     60     /** {@inheritDoc} */
     61     @Deprecated
     62     public double integrate(final double min, final double max)
     63         throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException {
     64         return integrate(f, min, max);
     65     }
     66 
     67     /** {@inheritDoc} */
     68     public double integrate(final UnivariateRealFunction f, final double min, final double max)
     69         throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException {
     70 
     71         final int m = maximalIterationCount + 1;
     72         double previousRow[] = new double[m];
     73         double currentRow[]  = new double[m];
     74 
     75         clearResult();
     76         verifyInterval(min, max);
     77         verifyIterationCount();
     78 
     79         TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
     80         currentRow[0] = qtrap.stage(f, min, max, 0);
     81         double olds = currentRow[0];
     82         for (int i = 1; i <= maximalIterationCount; ++i) {
     83 
     84             // switch rows
     85             final double[] tmpRow = previousRow;
     86             previousRow = currentRow;
     87             currentRow = tmpRow;
     88 
     89             currentRow[0] = qtrap.stage(f, min, max, i);
     90             for (int j = 1; j <= i; j++) {
     91                 // Richardson extrapolation coefficient
     92                 final double r = (1L << (2 * j)) - 1;
     93                 final double tIJm1 = currentRow[j - 1];
     94                 currentRow[j] = tIJm1 + (tIJm1 - previousRow[j - 1]) / r;
     95             }
     96             final double s = currentRow[i];
     97             if (i >= minimalIterationCount) {
     98                 final double delta  = FastMath.abs(s - olds);
     99                 final double rLimit = relativeAccuracy * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5;
    100                 if ((delta <= rLimit) || (delta <= absoluteAccuracy)) {
    101                     setResult(s, i);
    102                     return result;
    103                 }
    104             }
    105             olds = s;
    106         }
    107         throw new MaxIterationsExceededException(maximalIterationCount);
    108     }
    109 
    110     /** {@inheritDoc} */
    111     @Override
    112     protected void verifyIterationCount() throws IllegalArgumentException {
    113         super.verifyIterationCount();
    114         // at most 32 bisection refinements due to higher order divider
    115         if (maximalIterationCount > 32) {
    116             throw MathRuntimeException.createIllegalArgumentException(
    117                     LocalizedFormats.INVALID_ITERATIONS_LIMITS,
    118                     0, 32);
    119         }
    120     }
    121 }
    122