    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math3.analysis.integration;
    018
    019    import org.apache.commons.math3.exception.MathIllegalArgumentException;
    020    import org.apache.commons.math3.exception.MaxCountExceededException;
    021    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    022    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    023    import org.apache.commons.math3.exception.TooManyEvaluationsException;
    024    import org.apache.commons.math3.exception.util.LocalizedFormats;
    025    import org.apache.commons.math3.util.FastMath;
    026
    027    /**
    028     * Implements the <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html">
    029     * Legendre-Gauss</a> quadrature formula.
    030     * <p>
    031     * Legendre-Gauss integrators are efficient integrators that can
    032     * accurately integrate functions with few function evaluations. A
    033     * Legendre-Gauss integrator using an n-points quadrature formula can
    034     * integrate 2n-1 degree polynomials exactly.
    035     * </p>
    036     * <p>
    037     * These integrators evaluate the function on n carefully chosen
    038     * abscissas in each step interval (mapped to the canonical [-1,1] interval).
    039     * The evaluation abscissas are not evenly spaced and none of them are
    040     * at the interval endpoints. This implies the function integrated can be
    041     * undefined at integration interval endpoints.
    042     * </p>
    043     * <p>
    044     * The evaluation abscissas x<sub>i</sub> are the roots of the degree n
    045     * Legendre polynomial. The weights a<sub>i</sub> of the quadrature formula
    046     * integrals from -1 to +1 &int; Li<sup>2</sup> where Li (x) =
    047     * &prod; (x-x<sub>k</sub>)/(x<sub>i</sub>-x<sub>k</sub>) for k != i.
    048     * </p>
    049     * <p>
    050     * @version $Id: LegendreGaussIntegrator.java 1364452 2012-07-22 22:30:01Z erans $
    051     * @since 1.2
    052     * @deprecated As of 3.1 (to be removed in 4.0). Please use
    053     * {@link IterativeLegendreGaussIntegrator} instead.
    054     */
    055    @Deprecated
    056    public class LegendreGaussIntegrator extends BaseAbstractUnivariateIntegrator {
    057
    058        /** Abscissas for the 2 points method. */
    059        private static final double[] ABSCISSAS_2 = {
    060            -1.0 / FastMath.sqrt(3.0),
    061             1.0 / FastMath.sqrt(3.0)
    062        };
    063
    064        /** Weights for the 2 points method. */
    065        private static final double[] WEIGHTS_2 = {
    066            1.0,
    067            1.0
    068        };
    069
    070        /** Abscissas for the 3 points method. */
    071        private static final double[] ABSCISSAS_3 = {
    072            -FastMath.sqrt(0.6),
    073             0.0,
    074             FastMath.sqrt(0.6)
    075        };
    076
    077        /** Weights for the 3 points method. */
    078        private static final double[] WEIGHTS_3 = {
    079            5.0 / 9.0,
    080            8.0 / 9.0,
    081            5.0 / 9.0
    082        };
    083
    084        /** Abscissas for the 4 points method. */
    085        private static final double[] ABSCISSAS_4 = {
    086            -FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0),
    087            -FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0),
    088             FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0),
    089             FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0)
    090        };
    091
    092        /** Weights for the 4 points method. */
    093        private static final double[] WEIGHTS_4 = {
    094            (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0,
    095            (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0,
    096            (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0,
    097            (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0
    098        };
    099
    100        /** Abscissas for the 5 points method. */
    101        private static final double[] ABSCISSAS_5 = {
    102            -FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0),
    103            -FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0),
    104             0.0,
    105             FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0),
    106             FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0)
    107        };
    108
    109        /** Weights for the 5 points method. */
    110        private static final double[] WEIGHTS_5 = {
    111            (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0,
    112            (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0,
    113            128.0 / 225.0,
    114            (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0,
    115            (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0
    116        };
    117
    118        /** Abscissas for the current method. */
    119        private final double[] abscissas;
    120
    121        /** Weights for the current method. */
    122        private final double[] weights;
    123
    124        /**
    125         * Build a Legendre-Gauss integrator with given accuracies and iterations counts.
    126         * @param n number of points desired (must be between 2 and 5 inclusive)
    127         * @param relativeAccuracy relative accuracy of the result
    128         * @param absoluteAccuracy absolute accuracy of the result
    129         * @param minimalIterationCount minimum number of iterations
    130         * @param maximalIterationCount maximum number of iterations
    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 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            }
    165
    166        }
    167
    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         */
    174        public LegendreGaussIntegrator(final int n,
    175                                       final double relativeAccuracy,
    176                                       final double absoluteAccuracy) {
    177            this(n, relativeAccuracy, absoluteAccuracy,
    178                 DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT);
    179        }
    180
    181        /**
    182         * Build a Legendre-Gauss integrator with given iteration counts.
    183         * @param n number of points desired (must be between 2 and 5 inclusive)
    184         * @param minimalIterationCount minimum number of iterations
    185         * @param maximalIterationCount maximum number of iterations
    186         * @exception NotStrictlyPositiveException if minimal number of iterations
    187         * is not strictly positive
    188         * @exception NumberIsTooSmallException if maximal number of iterations
    189         * is lesser than or equal to the minimal number of iterations
    190         */
    191        public LegendreGaussIntegrator(final int n,
    192                                       final int minimalIterationCount,
    193                                       final int maximalIterationCount) {
    194            this(n, DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY,
    195                 minimalIterationCount, maximalIterationCount);
    196        }
    197
    198        /** {@inheritDoc} */
    199        @Override
    200        protected double doIntegrate()
    201            throws TooManyEvaluationsException, MaxCountExceededException {
    202
    203            // compute first estimate with a single step
    204            double oldt = stage(1);
    205
    206            int n = 2;
    207            while (true) {
    208
    209                // improve integral with a larger number of steps
    210                final double t = stage(n);
    211
    212                // estimate error
    213                final double delta = FastMath.abs(t - oldt);
    214                final double limit =
    215                    FastMath.max(getAbsoluteAccuracy(),
    216                                 getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5);
    217
    218                // check convergence
    219                if ((iterations.getCount() + 1 >= getMinimalIterationCount()) && (delta <= limit)) {
    220                    return t;
    221                }
    222
    223                // prepare next iteration
    224                double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / abscissas.length));
    225                n = FastMath.max((int) (ratio * n), n + 1);
    226                oldt = t;
    227                iterations.incrementCount();
    228
    229            }
    230
    231        }
    232
    233        /**
    234         * Compute the n-th stage integral.
    235         * @param n number of steps
    236         * @return the value of n-th stage integral
    237         * @throws TooManyEvaluationsException if the maximum number of evaluations
    238         * is exceeded.
    239         */
    240        private double stage(final int n)
    241            throws TooManyEvaluationsException {
    242
    243            // set up the step for the current stage
    244            final double step     = (getMax() - getMin()) / n;
    245            final double halfStep = step / 2.0;
    246
    247            // integrate over all elementary steps
    248            double midPoint = getMin() + halfStep;
    249            double sum = 0.0;
    250            for (int i = 0; i < n; ++i) {
    251                for (int j = 0; j < abscissas.length; ++j) {
    252                    sum += weights[j] * computeObjectiveValue(midPoint + halfStep * abscissas[j]);
    253                }
    254                midPoint += step;
    255            }
    256
    257            return halfStep * sum;
    258
    259        }
    260
    261    }