View Javadoc

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  
24  /**
25   * Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html">
26   * Romberg Algorithm</a> for integration of real univariate functions. For
27   * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
28   * chapter 3.
29   * <p>
30   * Romberg integration employs k successive refinements of the trapezoid
31   * rule to remove error terms less than order O(N^(-2k)). Simpson's rule
32   * is a special case of k = 2.</p>
33   *  
34   * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
35   * @since 1.2
36   */
37  public class RombergIntegrator extends UnivariateRealIntegratorImpl {
38  
39      /**
40       * Construct an integrator for the given function.
41       * 
42       * @param f function to integrate
43       * @deprecated as of 2.0 the integrand function is passed as an argument
44       * to the {@link #integrate(UnivariateRealFunction, double, double)}method.
45       */
46      @Deprecated
47      public RombergIntegrator(UnivariateRealFunction f) {
48          super(f, 32);
49      }
50  
51      /**
52       * Construct an integrator.
53       */
54      public RombergIntegrator() {
55          super(32);
56      }
57  
58      /** {@inheritDoc} */
59      @Deprecated
60      public double integrate(final double min, final double max)
61          throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException {
62          return integrate(f, min, max);
63      }
64  
65      /** {@inheritDoc} */
66      public double integrate(final UnivariateRealFunction f,
67                              final double min, final double max)
68          throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException {
69          
70          int i = 1, j, m = maximalIterationCount + 1;
71          // Array structure here can be improved for better space
72          // efficiency because only the lower triangle is used.
73          double r, t[][] = new double[m][m], s, olds;
74  
75          clearResult();
76          verifyInterval(min, max);
77          verifyIterationCount();
78  
79          TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
80          t[0][0] = qtrap.stage(f, min, max, 0);
81          olds = t[0][0];
82          while (i <= maximalIterationCount) {
83              t[i][0] = qtrap.stage(f, min, max, i);
84              for (j = 1; j <= i; j++) {
85                  // Richardson extrapolation coefficient
86                  r = (1L << (2 * j)) -1;
87                  t[i][j] = t[i][j-1] + (t[i][j-1] - t[i-1][j-1]) / r;
88              }
89              s = t[i][i];
90              if (i >= minimalIterationCount) {
91                  final double delta = Math.abs(s - olds);
92                  final double rLimit =
93                      relativeAccuracy * (Math.abs(olds) + Math.abs(s)) * 0.5; 
94                  if ((delta <= rLimit) || (delta <= absoluteAccuracy)) {
95                      setResult(s, i);
96                      return result;
97                  }
98              }
99              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 }