FieldRombergIntegrator.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.  *      https://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. /*
  18.  * This is not the original file distributed by the Apache Software Foundation
  19.  * It has been modified by the Hipparchus project
  20.  */
  21. package org.hipparchus.analysis.integration;

  22. import org.hipparchus.CalculusFieldElement;
  23. import org.hipparchus.Field;
  24. import org.hipparchus.exception.LocalizedCoreFormats;
  25. import org.hipparchus.exception.MathIllegalArgumentException;
  26. import org.hipparchus.exception.MathIllegalStateException;
  27. import org.hipparchus.util.FastMath;
  28. import org.hipparchus.util.MathArrays;

  29. /**
  30.  * Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html">
  31.  * Romberg Algorithm</a> for integration of real univariate functions. For
  32.  * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X,
  33.  * chapter 3.
  34.  * <p>
  35.  * Romberg integration employs k successive refinements of the trapezoid
  36.  * rule to remove error terms less than order O(N^(-2k)). Simpson's rule
  37.  * is a special case of k = 2.</p>
  38.  * @param <T> Type of the field elements.
  39.  * @since 2.0
  40.  */
  41. public class FieldRombergIntegrator<T extends CalculusFieldElement<T>> extends BaseAbstractFieldUnivariateIntegrator<T> {

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

  44.     /**
  45.      * Build a Romberg integrator with given accuracies and iterations counts.
  46.      * @param field field to which function argument and value belong
  47.      * @param relativeAccuracy relative accuracy of the result
  48.      * @param absoluteAccuracy absolute accuracy of the result
  49.      * @param minimalIterationCount minimum number of iterations
  50.      * @param maximalIterationCount maximum number of iterations
  51.      * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
  52.      * @exception MathIllegalArgumentException if minimal number of iterations
  53.      * is not strictly positive
  54.      * @exception MathIllegalArgumentException if maximal number of iterations
  55.      * is lesser than or equal to the minimal number of iterations
  56.      * @exception MathIllegalArgumentException if maximal number of iterations
  57.      * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT}
  58.      */
  59.     public FieldRombergIntegrator(final Field<T> field,
  60.                                   final double relativeAccuracy,
  61.                                   final double absoluteAccuracy,
  62.                                   final int minimalIterationCount,
  63.                                   final int maximalIterationCount)
  64.         throws MathIllegalArgumentException {
  65.         super(field, relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
  66.         if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) {
  67.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE_BOUND_EXCLUDED,
  68.                                                    maximalIterationCount, ROMBERG_MAX_ITERATIONS_COUNT);
  69.         }
  70.     }

  71.     /**
  72.      * Build a Romberg integrator with given iteration counts.
  73.      * @param field field to which function argument and value belong
  74.      * @param minimalIterationCount minimum number of iterations
  75.      * @param maximalIterationCount maximum number of iterations
  76.      * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
  77.      * @exception MathIllegalArgumentException if minimal number of iterations
  78.      * is not strictly positive
  79.      * @exception MathIllegalArgumentException if maximal number of iterations
  80.      * is lesser than or equal to the minimal number of iterations
  81.      * @exception MathIllegalArgumentException if maximal number of iterations
  82.      * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT}
  83.      */
  84.     public FieldRombergIntegrator(final Field<T> field,
  85.                                   final int minimalIterationCount,
  86.                                   final int maximalIterationCount)
  87.         throws MathIllegalArgumentException {
  88.         super(field, minimalIterationCount, maximalIterationCount);
  89.         if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) {
  90.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE_BOUND_EXCLUDED,
  91.                                                    maximalIterationCount, ROMBERG_MAX_ITERATIONS_COUNT);
  92.         }
  93.     }

  94.     /**
  95.      * Construct a Romberg integrator with default settings
  96.      * @param field field to which function argument and value belong
  97.      * (max iteration count set to {@link #ROMBERG_MAX_ITERATIONS_COUNT})
  98.      */
  99.     public FieldRombergIntegrator(final Field<T> field) {
  100.         super(field, DEFAULT_MIN_ITERATIONS_COUNT, ROMBERG_MAX_ITERATIONS_COUNT);
  101.     }

  102.     /** {@inheritDoc} */
  103.     @Override
  104.     protected T doIntegrate()
  105.         throws MathIllegalStateException {

  106.         final int m = iterations.getMaximalCount() + 1;
  107.         T[] previousRow = MathArrays.buildArray(getField(), m);
  108.         T[] currentRow = MathArrays.buildArray(getField(), m);

  109.         FieldTrapezoidIntegrator<T> qtrap = new FieldTrapezoidIntegrator<>(getField());
  110.         currentRow[0] = qtrap.stage(this, 0);
  111.         iterations.increment();
  112.         T olds = currentRow[0];
  113.         while (true) {

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

  115.             // switch rows
  116.             final T[] tmpRow = previousRow;
  117.             previousRow = currentRow;
  118.             currentRow = tmpRow;

  119.             currentRow[0] = qtrap.stage(this, i);
  120.             iterations.increment();
  121.             for (int j = 1; j <= i; j++) {
  122.                 // Richardson extrapolation coefficient
  123.                 final double r = (1L << (2 * j)) - 1;
  124.                 final T tIJm1 = currentRow[j - 1];
  125.                 currentRow[j] = tIJm1.add(tIJm1.subtract(previousRow[j - 1]).divide(r));
  126.             }
  127.             final T s = currentRow[i];
  128.             if (i >= getMinimalIterationCount()) {
  129.                 final double delta  = FastMath.abs(s.subtract(olds)).getReal();
  130.                 final double rLimit = FastMath.abs(olds).add(FastMath.abs(s)).multiply(0.5 * getRelativeAccuracy()).getReal();
  131.                  if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) {
  132.                     return s;
  133.                 }
  134.             }
  135.             olds = s;
  136.         }

  137.     }

  138. }