LegendreGaussIntegrator.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.MathIllegalArgumentException;
  19. import org.apache.commons.math3.exception.MaxCountExceededException;
  20. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  21. import org.apache.commons.math3.exception.NumberIsTooSmallException;
  22. import org.apache.commons.math3.exception.TooManyEvaluationsException;
  23. import org.apache.commons.math3.exception.util.LocalizedFormats;
  24. import org.apache.commons.math3.util.FastMath;

  25. /**
  26.  * Implements the <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html">
  27.  * Legendre-Gauss</a> quadrature formula.
  28.  * <p>
  29.  * Legendre-Gauss integrators are efficient integrators that can
  30.  * accurately integrate functions with few function evaluations. A
  31.  * Legendre-Gauss integrator using an n-points quadrature formula can
  32.  * integrate 2n-1 degree polynomials exactly.
  33.  * </p>
  34.  * <p>
  35.  * These integrators evaluate the function on n carefully chosen
  36.  * abscissas in each step interval (mapped to the canonical [-1,1] interval).
  37.  * The evaluation abscissas are not evenly spaced and none of them are
  38.  * at the interval endpoints. This implies the function integrated can be
  39.  * undefined at integration interval endpoints.
  40.  * </p>
  41.  * <p>
  42.  * The evaluation abscissas x<sub>i</sub> are the roots of the degree n
  43.  * Legendre polynomial. The weights a<sub>i</sub> of the quadrature formula
  44.  * integrals from -1 to +1 &int; Li<sup>2</sup> where Li (x) =
  45.  * &prod; (x-x<sub>k</sub>)/(x<sub>i</sub>-x<sub>k</sub>) for k != i.
  46.  * </p>
  47.  * <p>
  48.  * @since 1.2
  49.  * @deprecated As of 3.1 (to be removed in 4.0). Please use
  50.  * {@link IterativeLegendreGaussIntegrator} instead.
  51.  */
  52. @Deprecated
  53. public class LegendreGaussIntegrator extends BaseAbstractUnivariateIntegrator {

  54.     /** Abscissas for the 2 points method. */
  55.     private static final double[] ABSCISSAS_2 = {
  56.         -1.0 / FastMath.sqrt(3.0),
  57.          1.0 / FastMath.sqrt(3.0)
  58.     };

  59.     /** Weights for the 2 points method. */
  60.     private static final double[] WEIGHTS_2 = {
  61.         1.0,
  62.         1.0
  63.     };

  64.     /** Abscissas for the 3 points method. */
  65.     private static final double[] ABSCISSAS_3 = {
  66.         -FastMath.sqrt(0.6),
  67.          0.0,
  68.          FastMath.sqrt(0.6)
  69.     };

  70.     /** Weights for the 3 points method. */
  71.     private static final double[] WEIGHTS_3 = {
  72.         5.0 / 9.0,
  73.         8.0 / 9.0,
  74.         5.0 / 9.0
  75.     };

  76.     /** Abscissas for the 4 points method. */
  77.     private static final double[] ABSCISSAS_4 = {
  78.         -FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0),
  79.         -FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0),
  80.          FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0),
  81.          FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0)
  82.     };

  83.     /** Weights for the 4 points method. */
  84.     private static final double[] WEIGHTS_4 = {
  85.         (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0,
  86.         (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0,
  87.         (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0,
  88.         (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0
  89.     };

  90.     /** Abscissas for the 5 points method. */
  91.     private static final double[] ABSCISSAS_5 = {
  92.         -FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0),
  93.         -FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0),
  94.          0.0,
  95.          FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0),
  96.          FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0)
  97.     };

  98.     /** Weights for the 5 points method. */
  99.     private static final double[] WEIGHTS_5 = {
  100.         (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0,
  101.         (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0,
  102.         128.0 / 225.0,
  103.         (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0,
  104.         (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0
  105.     };

  106.     /** Abscissas for the current method. */
  107.     private final double[] abscissas;

  108.     /** Weights for the current method. */
  109.     private final double[] weights;

  110.     /**
  111.      * Build a Legendre-Gauss integrator with given accuracies and iterations counts.
  112.      * @param n number of points desired (must be between 2 and 5 inclusive)
  113.      * @param relativeAccuracy relative accuracy of the result
  114.      * @param absoluteAccuracy absolute accuracy of the result
  115.      * @param minimalIterationCount minimum number of iterations
  116.      * @param maximalIterationCount maximum number of iterations
  117.      * @exception MathIllegalArgumentException if number of points is out of [2; 5]
  118.      * @exception NotStrictlyPositiveException if minimal number of iterations
  119.      * is not strictly positive
  120.      * @exception NumberIsTooSmallException if maximal number of iterations
  121.      * is lesser than or equal to the minimal number of iterations
  122.      */
  123.     public LegendreGaussIntegrator(final int n,
  124.                                    final double relativeAccuracy,
  125.                                    final double absoluteAccuracy,
  126.                                    final int minimalIterationCount,
  127.                                    final int maximalIterationCount)
  128.         throws MathIllegalArgumentException, NotStrictlyPositiveException, NumberIsTooSmallException {
  129.         super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
  130.         switch(n) {
  131.         case 2 :
  132.             abscissas = ABSCISSAS_2;
  133.             weights   = WEIGHTS_2;
  134.             break;
  135.         case 3 :
  136.             abscissas = ABSCISSAS_3;
  137.             weights   = WEIGHTS_3;
  138.             break;
  139.         case 4 :
  140.             abscissas = ABSCISSAS_4;
  141.             weights   = WEIGHTS_4;
  142.             break;
  143.         case 5 :
  144.             abscissas = ABSCISSAS_5;
  145.             weights   = WEIGHTS_5;
  146.             break;
  147.         default :
  148.             throw new MathIllegalArgumentException(
  149.                     LocalizedFormats.N_POINTS_GAUSS_LEGENDRE_INTEGRATOR_NOT_SUPPORTED,
  150.                     n, 2, 5);
  151.         }

  152.     }

  153.     /**
  154.      * Build a Legendre-Gauss integrator with given accuracies.
  155.      * @param n number of points desired (must be between 2 and 5 inclusive)
  156.      * @param relativeAccuracy relative accuracy of the result
  157.      * @param absoluteAccuracy absolute accuracy of the result
  158.      * @exception MathIllegalArgumentException if number of points is out of [2; 5]
  159.      */
  160.     public LegendreGaussIntegrator(final int n,
  161.                                    final double relativeAccuracy,
  162.                                    final double absoluteAccuracy)
  163.         throws MathIllegalArgumentException {
  164.         this(n, relativeAccuracy, absoluteAccuracy,
  165.              DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT);
  166.     }

  167.     /**
  168.      * Build a Legendre-Gauss integrator with given iteration counts.
  169.      * @param n number of points desired (must be between 2 and 5 inclusive)
  170.      * @param minimalIterationCount minimum number of iterations
  171.      * @param maximalIterationCount maximum number of iterations
  172.      * @exception MathIllegalArgumentException if number of points is out of [2; 5]
  173.      * @exception NotStrictlyPositiveException if minimal number of iterations
  174.      * is not strictly positive
  175.      * @exception NumberIsTooSmallException if maximal number of iterations
  176.      * is lesser than or equal to the minimal number of iterations
  177.      */
  178.     public LegendreGaussIntegrator(final int n,
  179.                                    final int minimalIterationCount,
  180.                                    final int maximalIterationCount)
  181.         throws MathIllegalArgumentException {
  182.         this(n, DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY,
  183.              minimalIterationCount, maximalIterationCount);
  184.     }

  185.     /** {@inheritDoc} */
  186.     @Override
  187.     protected double doIntegrate()
  188.         throws MathIllegalArgumentException, TooManyEvaluationsException, MaxCountExceededException {

  189.         // compute first estimate with a single step
  190.         double oldt = stage(1);

  191.         int n = 2;
  192.         while (true) {

  193.             // improve integral with a larger number of steps
  194.             final double t = stage(n);

  195.             // estimate error
  196.             final double delta = FastMath.abs(t - oldt);
  197.             final double limit =
  198.                 FastMath.max(getAbsoluteAccuracy(),
  199.                              getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5);

  200.             // check convergence
  201.             if ((iterations.getCount() + 1 >= getMinimalIterationCount()) && (delta <= limit)) {
  202.                 return t;
  203.             }

  204.             // prepare next iteration
  205.             double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / abscissas.length));
  206.             n = FastMath.max((int) (ratio * n), n + 1);
  207.             oldt = t;
  208.             iterations.incrementCount();

  209.         }

  210.     }

  211.     /**
  212.      * Compute the n-th stage integral.
  213.      * @param n number of steps
  214.      * @return the value of n-th stage integral
  215.      * @throws TooManyEvaluationsException if the maximum number of evaluations
  216.      * is exceeded.
  217.      */
  218.     private double stage(final int n)
  219.         throws TooManyEvaluationsException {

  220.         // set up the step for the current stage
  221.         final double step     = (getMax() - getMin()) / n;
  222.         final double halfStep = step / 2.0;

  223.         // integrate over all elementary steps
  224.         double midPoint = getMin() + halfStep;
  225.         double sum = 0.0;
  226.         for (int i = 0; i < n; ++i) {
  227.             for (int j = 0; j < abscissas.length; ++j) {
  228.                 sum += weights[j] * computeObjectiveValue(midPoint + halfStep * abscissas[j]);
  229.             }
  230.             midPoint += step;
  231.         }

  232.         return halfStep * sum;

  233.     }

  234. }