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;
18  
19  import org.apache.commons.math.FunctionEvaluationException;
20  import org.apache.commons.math.MaxIterationsExceededException;
21  
22  /**
23   * Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html">
24   * Romberg Algorithm</a> for integration of real univariate functions. For
25   * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
26   * chapter 3.
27   * <p>
28   * Romberg integration employs k successvie refinements of the trapezoid
29   * rule to remove error terms less than order O(N^(-2k)). Simpson's rule
30   * is a special case of k = 2.</p>
31   *  
32   * @version $Revision: 620312 $ $Date: 2008-02-10 12:28:59 -0700 (Sun, 10 Feb 2008) $
33   * @since 1.2
34   */
35  public class RombergIntegrator extends UnivariateRealIntegratorImpl {
36  
37      /** serializable version identifier */
38      private static final long serialVersionUID = -1058849527738180243L;
39  
40      /**
41       * Construct an integrator for the given function.
42       * 
43       * @param f function to integrate
44       */
45      public RombergIntegrator(UnivariateRealFunction f) {
46          super(f, 32);
47      }
48  
49      /**
50       * Integrate the function in the given interval.
51       * 
52       * @param min the lower bound for the interval
53       * @param max the upper bound for the interval
54       * @return the value of integral
55       * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
56       * or the integrator detects convergence problems otherwise
57       * @throws FunctionEvaluationException if an error occurs evaluating the
58       * function
59       * @throws IllegalArgumentException if any parameters are invalid
60       */
61      public double integrate(double min, double max) throws MaxIterationsExceededException,
62          FunctionEvaluationException, IllegalArgumentException {
63          
64          int i = 1, j, m = maximalIterationCount + 1;
65          // Array strcture here can be improved for better space
66          // efficiency because only the lower triangle is used.
67          double r, t[][] = new double[m][m], s, olds;
68  
69          clearResult();
70          verifyInterval(min, max);
71          verifyIterationCount();
72  
73          TrapezoidIntegrator qtrap = new TrapezoidIntegrator(this.f);
74          t[0][0] = qtrap.stage(min, max, 0);
75          olds = t[0][0];
76          while (i <= maximalIterationCount) {
77              t[i][0] = qtrap.stage(min, max, i);
78              for (j = 1; j <= i; j++) {
79                  // Richardson extrapolation coefficient
80                  r = (1L << (2 * j)) -1;
81                  t[i][j] = t[i][j-1] + (t[i][j-1] - t[i-1][j-1]) / r;
82              }
83              s = t[i][i];
84              if (i >= minimalIterationCount) {
85                  if (Math.abs(s - olds) <= Math.abs(relativeAccuracy * olds)) {
86                      setResult(s, i);
87                      return result;
88                  }
89              }
90              olds = s;
91              i++;
92          }
93          throw new MaxIterationsExceededException(maximalIterationCount);
94      }
95  
96      /**
97       * Verifies that the iteration limits are valid and within the range.
98       * 
99       * @throws IllegalArgumentException if not
100      */
101     protected void verifyIterationCount() throws IllegalArgumentException {
102         super.verifyIterationCount();
103         // at most 32 bisection refinements due to higher order divider
104         if (maximalIterationCount > 32) {
105             throw new IllegalArgumentException
106                 ("Iteration upper limit out of [0, 32] range: " +
107                 maximalIterationCount);
108         }
109     }
110 }