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;

  19. import org.apache.commons.math3.exception.MathIllegalArgumentException;
  20. import org.apache.commons.math3.exception.MaxCountExceededException;
  21. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  22. import org.apache.commons.math3.exception.NumberIsTooSmallException;
  23. import org.apache.commons.math3.exception.TooManyEvaluationsException;
  24. import org.apache.commons.math3.exception.util.LocalizedFormats;
  25. import org.apache.commons.math3.util.FastMath;

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

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

  63.     /** Weights for the 2 points method. */
  64.     private static final double[] WEIGHTS_2 = {
  65.         1.0,
  66.         1.0
  67.     };

  69.     /** Abscissas for the 3 points method. */
  70.     private static final double[] ABSCISSAS_3 = {
  71.         -FastMath.sqrt(0.6),
  72.          0.0,
  73.          FastMath.sqrt(0.6)
  74.     };

  76.     /** Weights for the 3 points method. */
  77.     private static final double[] WEIGHTS_3 = {
  78.         5.0 / 9.0,
  79.         8.0 / 9.0,
  80.         5.0 / 9.0
  81.     };

  83.     /** Abscissas for the 4 points method. */
  84.     private static final double[] ABSCISSAS_4 = {
  85.         -FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0),
  86.         -FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0),
  87.          FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0),
  88.          FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0)
  89.     };

  91.     /** Weights for the 4 points method. */
  92.     private static final double[] WEIGHTS_4 = {
  93.         (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0,
  94.         (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0,
  95.         (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0,
  96.         (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0
  97.     };

  99.     /** Abscissas for the 5 points method. */
  100.     private static final double[] ABSCISSAS_5 = {
  101.         -FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0),
  102.         -FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0),
  103.          0.0,
  104.          FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0),
  105.          FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0)
  106.     };

  108.     /** Weights for the 5 points method. */
  109.     private static final double[] WEIGHTS_5 = {
  110.         (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0,
  111.         (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0,
  112.         128.0 / 225.0,
  113.         (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0,
  114.         (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0
  115.     };

  117.     /** Abscissas for the current method. */
  118.     private final double[] abscissas;

  120.     /** Weights for the current method. */
  121.     private final double[] weights;

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

  166.     }

  168.     /**
  169.      * Build a Legendre-Gauss integrator with given accuracies.
  170.      * @param n number of points desired (must be between 2 and 5 inclusive)
  171.      * @param relativeAccuracy relative accuracy of the result
  172.      * @param absoluteAccuracy absolute accuracy of the result
  173.      * @exception MathIllegalArgumentException if number of points is out of [2; 5]
  174.      */
  175.     public LegendreGaussIntegrator(final int n,
  176.                                    final double relativeAccuracy,
  177.                                    final double absoluteAccuracy)
  178.         throws MathIllegalArgumentException {
  179.         this(n, relativeAccuracy, absoluteAccuracy,
  180.              DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT);
  181.     }

  183.     /**
  184.      * Build a Legendre-Gauss integrator with given iteration counts.
  185.      * @param n number of points desired (must be between 2 and 5 inclusive)
  186.      * @param minimalIterationCount minimum number of iterations
  187.      * @param maximalIterationCount maximum number of iterations
  188.      * @exception MathIllegalArgumentException if number of points is out of [2; 5]
  189.      * @exception NotStrictlyPositiveException if minimal number of iterations
  190.      * is not strictly positive
  191.      * @exception NumberIsTooSmallException if maximal number of iterations
  192.      * is lesser than or equal to the minimal number of iterations
  193.      */
  194.     public LegendreGaussIntegrator(final int n,
  195.                                    final int minimalIterationCount,
  196.                                    final int maximalIterationCount)
  197.         throws MathIllegalArgumentException {
  198.         this(n, DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY,
  199.              minimalIterationCount, maximalIterationCount);
  200.     }

  202.     /** {@inheritDoc} */
  203.     @Override
  204.     protected double doIntegrate()
  205.         throws MathIllegalArgumentException, TooManyEvaluationsException, MaxCountExceededException {

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

  210.         int n = 2;
  211.         while (true) {

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

  216.             // estimate error
  217.             final double delta = FastMath.abs(t - oldt);
  218.             final double limit =
  219.                 FastMath.max(getAbsoluteAccuracy(),
  220.                              getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5);

  222.             // check convergence
  223.             if ((iterations.getCount() + 1 >= getMinimalIterationCount()) && (delta <= limit)) {
  224.                 return t;
  225.             }

  227.             // prepare next iteration
  228.             double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / abscissas.length));
  229.             n = FastMath.max((int) (ratio * n), n + 1);
  230.             oldt = t;
  231.             iterations.incrementCount();

  233.         }

  235.     }

  237.     /**
  238.      * Compute the n-th stage integral.
  239.      * @param n number of steps
  240.      * @return the value of n-th stage integral
  241.      * @throws TooManyEvaluationsException if the maximum number of evaluations
  242.      * is exceeded.
  243.      */
  244.     private double stage(final int n)
  245.         throws TooManyEvaluationsException {

  247.         // set up the step for the current stage
  248.         final double step     = (getMax() - getMin()) / n;
  249.         final double halfStep = step / 2.0;

  251.         // integrate over all elementary steps
  252.         double midPoint = getMin() + halfStep;
  253.         double sum = 0.0;
  254.         for (int i = 0; i < n; ++i) {
  255.             for (int j = 0; j < abscissas.length; ++j) {
  256.                 sum += weights[j] * computeObjectiveValue(midPoint + halfStep * abscissas[j]);
  257.             }
  258.             midPoint += step;
  259.         }

  261.         return halfStep * sum;

  263.     }

  265. }
Created with JaCoCo 0.7.2.201409121644