RombergIntegrator.java

  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.math3.analysis.integration;

  18. import org.apache.commons.math3.exception.MaxCountExceededException;
  19. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  20. import org.apache.commons.math3.exception.NumberIsTooLargeException;
  21. import org.apache.commons.math3.exception.NumberIsTooSmallException;
  22. import org.apache.commons.math3.exception.TooManyEvaluationsException;
  23. import org.apache.commons.math3.util.FastMath;

  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.  * @since 1.2
  35.  */
  36. public class RombergIntegrator extends BaseAbstractUnivariateIntegrator {

  37.     /** Maximal number of iterations for Romberg. */
  38.     public static final int ROMBERG_MAX_ITERATIONS_COUNT = 32;

  39.     /**
  40.      * Build a Romberg integrator with given accuracies and iterations counts.
  41.      * @param relativeAccuracy relative accuracy of the result
  42.      * @param absoluteAccuracy absolute accuracy of the result
  43.      * @param minimalIterationCount minimum number of iterations
  44.      * @param maximalIterationCount maximum number of iterations
  45.      * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
  46.      * @exception NotStrictlyPositiveException if minimal number of iterations
  47.      * is not strictly positive
  48.      * @exception NumberIsTooSmallException if maximal number of iterations
  49.      * is lesser than or equal to the minimal number of iterations
  50.      * @exception NumberIsTooLargeException if maximal number of iterations
  51.      * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT}
  52.      */
  53.     public RombergIntegrator(final double relativeAccuracy,
  54.                              final double absoluteAccuracy,
  55.                              final int minimalIterationCount,
  56.                              final int maximalIterationCount)
  57.         throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {
  58.         super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
  59.         if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) {
  60.             throw new NumberIsTooLargeException(maximalIterationCount,
  61.                                                 ROMBERG_MAX_ITERATIONS_COUNT, false);
  62.         }
  63.     }

  64.     /**
  65.      * Build a Romberg integrator with given iteration counts.
  66.      * @param minimalIterationCount minimum number of iterations
  67.      * @param maximalIterationCount maximum number of iterations
  68.      * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
  69.      * @exception NotStrictlyPositiveException if minimal number of iterations
  70.      * is not strictly positive
  71.      * @exception NumberIsTooSmallException if maximal number of iterations
  72.      * is lesser than or equal to the minimal number of iterations
  73.      * @exception NumberIsTooLargeException if maximal number of iterations
  74.      * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT}
  75.      */
  76.     public RombergIntegrator(final int minimalIterationCount,
  77.                              final int maximalIterationCount)
  78.         throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {
  79.         super(minimalIterationCount, maximalIterationCount);
  80.         if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) {
  81.             throw new NumberIsTooLargeException(maximalIterationCount,
  82.                                                 ROMBERG_MAX_ITERATIONS_COUNT, false);
  83.         }
  84.     }

  85.     /**
  86.      * Construct a Romberg integrator with default settings
  87.      * (max iteration count set to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
  88.      */
  89.     public RombergIntegrator() {
  90.         super(DEFAULT_MIN_ITERATIONS_COUNT, ROMBERG_MAX_ITERATIONS_COUNT);
  91.     }

  92.     /** {@inheritDoc} */
  93.     @Override
  94.     protected double doIntegrate()
  95.         throws TooManyEvaluationsException, MaxCountExceededException {

  96.         final int m = iterations.getMaximalCount() + 1;
  97.         double previousRow[] = new double[m];
  98.         double currentRow[]  = new double[m];

  99.         TrapezoidIntegrator qtrap = new TrapezoidIntegrator();
  100.         currentRow[0] = qtrap.stage(this, 0);
  101.         iterations.incrementCount();
  102.         double olds = currentRow[0];
  103.         while (true) {

  104.             final int i = iterations.getCount();

  105.             // switch rows
  106.             final double[] tmpRow = previousRow;
  107.             previousRow = currentRow;
  108.             currentRow = tmpRow;

  109.             currentRow[0] = qtrap.stage(this, i);
  110.             iterations.incrementCount();
  111.             for (int j = 1; j <= i; j++) {
  112.                 // Richardson extrapolation coefficient
  113.                 final double r = (1L << (2 * j)) - 1;
  114.                 final double tIJm1 = currentRow[j - 1];
  115.                 currentRow[j] = tIJm1 + (tIJm1 - previousRow[j - 1]) / r;
  116.             }
  117.             final double s = currentRow[i];
  118.             if (i >= getMinimalIterationCount()) {
  119.                 final double delta  = FastMath.abs(s - olds);
  120.                 final double rLimit = getRelativeAccuracy() * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5;
  121.                 if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) {
  122.                     return s;
  123.                 }
  124.             }
  125.             olds = s;
  126.         }

  127.     }

  128. }