    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.analysis.UnivariateFunction;
    020    import org.apache.commons.math3.analysis.integration.gauss.GaussIntegratorFactory;
    021    import org.apache.commons.math3.analysis.integration.gauss.GaussIntegrator;
    022    import org.apache.commons.math3.exception.MaxCountExceededException;
    023    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    024    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    025    import org.apache.commons.math3.exception.TooManyEvaluationsException;
    026    import org.apache.commons.math3.util.FastMath;
    027
    028    /**
    029     * This algorithm divides the integration interval into equally-sized
    030     * sub-interval and on each of them performs a
    031     * <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html">
    032     * Legendre-Gauss</a> quadrature.
    033     *
    034     * @version $Id: IterativeLegendreGaussIntegrator.java 1416643 2012-12-03 19:37:14Z tn $
    035     * @since 3.1
    036     */
    037
    038    public class IterativeLegendreGaussIntegrator
    039        extends BaseAbstractUnivariateIntegrator {
    040        /** Factory that computes the points and weights. */
    041        private static final GaussIntegratorFactory FACTORY
    042            = new GaussIntegratorFactory();
    043        /** Number of integration points (per interval). */
    044        private final int numberOfPoints;
    045
    046        /**
    047         * Builds an integrator with given accuracies and iterations counts.
    048         *
    049         * @param n Number of integration points.
    050         * @param relativeAccuracy Relative accuracy of the result.
    051         * @param absoluteAccuracy Absolute accuracy of the result.
    052         * @param minimalIterationCount Minimum number of iterations.
    053         * @param maximalIterationCount Maximum number of iterations.
    054         * @throws NotStrictlyPositiveException if minimal number of iterations
    055         * is not strictly positive.
    056         * @throws NumberIsTooSmallException if maximal number of iterations
    057         * is smaller than or equal to the minimal number of iterations.
    058         */
    059        public IterativeLegendreGaussIntegrator(final int n,
    060                                                final double relativeAccuracy,
    061                                                final double absoluteAccuracy,
    062                                                final int minimalIterationCount,
    063                                                final int maximalIterationCount)
    064            throws NotStrictlyPositiveException, NumberIsTooSmallException {
    065            super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount);
    066            numberOfPoints = n;
    067        }
    068
    069        /**
    070         * Builds an integrator with given accuracies.
    071         *
    072         * @param n Number of integration points.
    073         * @param relativeAccuracy Relative accuracy of the result.
    074         * @param absoluteAccuracy Absolute accuracy of the result.
    075         */
    076        public IterativeLegendreGaussIntegrator(final int n,
    077                                                final double relativeAccuracy,
    078                                                final double absoluteAccuracy) {
    079            this(n, relativeAccuracy, absoluteAccuracy,
    080                 DEFAULT_MIN_ITERATIONS_COUNT, DEFAULT_MAX_ITERATIONS_COUNT);
    081        }
    082
    083        /**
    084         * Builds an integrator with given iteration counts.
    085         *
    086         * @param n Number of integration points.
    087         * @param minimalIterationCount Minimum number of iterations.
    088         * @param maximalIterationCount Maximum number of iterations.
    089         * @throws NotStrictlyPositiveException if minimal number of iterations
    090         * is not strictly positive.
    091         * @throws NumberIsTooSmallException if maximal number of iterations
    092         * is smaller than or equal to the minimal number of iterations.
    093         */
    094        public IterativeLegendreGaussIntegrator(final int n,
    095                                                final int minimalIterationCount,
    096                                                final int maximalIterationCount) {
    097            this(n, DEFAULT_RELATIVE_ACCURACY, DEFAULT_ABSOLUTE_ACCURACY,
    098                 minimalIterationCount, maximalIterationCount);
    099        }
    100
    101        /** {@inheritDoc} */
    102        @Override
    103        protected double doIntegrate()
    104            throws TooManyEvaluationsException, MaxCountExceededException {
    105            // Compute first estimate with a single step.
    106            double oldt = stage(1);
    107
    108            int n = 2;
    109            while (true) {
    110                // Improve integral with a larger number of steps.
    111                final double t = stage(n);
    112
    113                // Estimate the error.
    114                final double delta = FastMath.abs(t - oldt);
    115                final double limit =
    116                    FastMath.max(getAbsoluteAccuracy(),
    117                                 getRelativeAccuracy() * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5);
    118
    119                // check convergence
    120                if (iterations.getCount() + 1 >= getMinimalIterationCount() &&
    121                    delta <= limit) {
    122                    return t;
    123                }
    124
    125                // Prepare next iteration.
    126                final double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / numberOfPoints));
    127                n = FastMath.max((int) (ratio * n), n + 1);
    128                oldt = t;
    129                iterations.incrementCount();
    130            }
    131        }
    132
    133        /**
    134         * Compute the n-th stage integral.
    135         *
    136         * @param n Number of steps.
    137         * @return the value of n-th stage integral.
    138         * @throws TooManyEvaluationsException if the maximum number of evaluations
    139         * is exceeded.
    140         */
    141        private double stage(final int n)
    142            throws TooManyEvaluationsException {
    143            // Function to be integrated is stored in the base class.
    144            final UnivariateFunction f = new UnivariateFunction() {
    145                    public double value(double x) {
    146                        return computeObjectiveValue(x);
    147                    }
    148                };
    149
    150            final double min = getMin();
    151            final double max = getMax();
    152            final double step = (max - min) / n;
    153
    154            double sum = 0;
    155            for (int i = 0; i < n; i++) {
    156                // Integrate over each sub-interval [a, b].
    157                final double a = min + i * step;
    158                final double b = a + step;
    159                final GaussIntegrator g = FACTORY.legendreHighPrecision(numberOfPoints, a, b);
    160                sum += g.integrate(f);
    161            }
    162
    163            return sum;
    164        }
    165    }